Численные методы решения дифференциальных уравнений

Численное дифференцирование  применяется тогда, когда функция  f(x), для которой нужно найти производную, задана таблично или же определяется сложной формулой. В первом случае методы дифференциального исчисления неприемлемы, а во втором случае их использование вызывают значительные трудности. В обоих случаях вместо функции f(x) рассматривают интерполирующую функцию ϕ(х) и считают производную f(xИзм.

Лист

№ докум.

Подпись

Дата

Лист

9

 

КП.230100. 08 .ПЗ

) приближенно равной производной  ϕ(х).

Вычисление первой производной

По определению, первая производная  гладкой функции f(x) в точке x вычисляется  как

 

При вычислении первой производной  функции f(x) на компьютере мы заменяем бесконечно малое h ® ¥ на малое, но конечное значение h:

 

где O(h) — ошибка вычисления производной, естественно зависящая  от h. Формула  называется разностной схемой для вычисления первой производной (более точно — правой разностной схемой или просто правой разностью). Аналогично, может быть записана левая разностная схема.

Разложим функцию f(x) в  ряд Тейлора в точке x + h:

 

откуда следует, что в  первом порядке разложения

 

При выборе очень малого h ошибки округления при вычислении на компьютере могут быть сравнимы или больше h. Следовательно мы заинтересованы в алгоритме, дающим меньшую величину ошибки при той же величине h.

Такой улучшенный алгоритм легко получить разлагая функцию f(x) в ряд Тейлора в точках x + h и x - h, вычитая затем один результат из другого, что дает

 

где погрешность вычисления первой производной

 

Это — центральная разностная схема (центральная разность).

Можно пойти по пути улучшения точности метода вычисления первой производной и дальше. Например, рассматривая разложение функции f(x) в ряд Тейлора в точках x + h, x + 2h, x - h и x - 2h, можно получить четырехточечную сИзм.

Лист

№ докум.

Подпись

Дата

Лист

10

 

КП.230100. 08 .ПЗ

хему и т.д.

Вычисление второй производной

По аналогии с получением разностных схем для вычисления первой производной могут быть получены разностные схемы для второй производной функции f(x). Наиболее часто используемая схема:

 

с погрешностью

 

2.1   Интерполяционная формула Ньютона

Пусть в равноотстоящих точках , где h – шаг интерполяции, заданы значения для функции . Требуется подобрать полином степени не выше n , удовлетворяющий условиям (1).

Введем конечные разности для последовательности значений  :

          (2)

Условия (1) эквивалентны равенствам:

при

Получим первую интерполяционную формулу Ньютона:

        (3)

где    – число шагов интерполяции от начальной точки  до точки х. 
Формулу (3) целесообразно использовать для интерполяции функции  в окрестности начальной точки , где  q по абсолютной величине мало.

 

3. Численные методы решения дифференциальных уравнений

Обыкновенным  дифференциальным уравнением n-го порядка называется следующее уравнение, которое содеИзм.

Лист

№ докум.

Подпись

Дата

Лист

11

 

КП.230100. 08 .ПЗ

ржит одну или несколько производных от искомой функции y(x):

, здесь y⁽ⁿ⁾ обозначает производную порядка n некоторой функции y(x), x – это независимая переменная.

 

3.1   Метод Эйлера

Наиболее простым способом численного решения задачи Коши для ОДУ первого порядка является метод Эйлера. В его основе лежит аппроксимация производной отношением конечных приращений зависимой (y) и независимой (x) переменных между узлами равномерной сетки:

где y_i+1  это искомое значение функции в точке x_i+1.

Если теперь преобразовать это  уравнение, и учесть равномерность  сетки интегрирования, то получится  итерационная формула, по которой можно  вычислить  y_i+1 , если известно y_i  в точке х_i:

                                                                                           

Сравнивая формулу Эйлера с общим выражением, полученным ранее, видно, что для  приближенного вычисления интеграла  в в методе Эйлера используется простейшая формула интегрирования - формула прямоугольников по левому краю отрезка.

Графическая интерпретация метода Эйлера, исходя из вида решаемого уравнения следует, что значение  есть значение производной функции y(x) в точке x=x_i -    , и, таким образом, равно тангенсу угла наклона касательной, проведенной к графику функции y(x) в точке x=x_i.

Рис.1 Графическая интерпретация метода Эйлера.  

Из прямоугольного треугольника на рисунке можно найти 

,

откуда и получается фИзм.

Лист

№ докум.

Подпись

Дата

Лист

12

 

КП.230100. 08 .ПЗ

ормула Эйлера. Таким образом, суть метода Эйлера заключается в замене функции y(x) на отрезке интегрирования прямой линией, касательной к графику в точке x=x_i. Если искомая функция сильно отличается от линейной на отрезке интегрирования, то погрешность вычисления будет значительной. Ошибка метода Эйлера прямо пропорциональна шагу интегрирования:

При больших h метод Эйлера весьма неточен. Он дает все более точное приближение при уменьшении шага интегрирования. Если отрезок [x_i, x_i+1] слишком велик, то каждый участок [x_i, x_i+1] разбивается на N отрезков интегрирования и к каждому их них применяется формула Эйлера с шагом , то есть шаг интегрирования h берется меньше шага сетки, на которой определяется решение.

3.2   Методы Рунге-Кутты

Улучшение точности решения ОДУ  первого порядка возможно за счет увеличения точности приближенного  вычисления интеграла в выражении .

Воспользовавшись формулой Симпсона, можно получить еще более точную формулу для решения задачи Коши для ОДУ первого порядка -  широко используемого в вычислительной практике метода Рунге-Кутты.

В формуле Симпсона для приближенного  вычисления определенного интеграла  используются значения под интегрального  выражения в  трех точках. В интеграле их всего две, поэтому введем дополнительную точку в середине отрезка [x_i+1 x_i]

тогда можно переписать так:

Полученное выражение является неявным, так как в правой части  содержатся  еще не определенные значения функции y_i+h/2 и y_i+1. Чтобы воспользоваться этоИзм.

Лист

№ докум.

Подпись

Дата

Лист

13

 

КП.230100. 08 .ПЗ

й формулой, надо использовать некоторое приближение для вычисления этих значений

 

При использовании различных методов  приближенного вычисления этих величин, получаются выражения для методов  Рунге-Кутты различного порядка  точности.

Алгоритм Рунге-Кутты третьего порядка - РК3 (погрешность порядка h³):

                                                                       

где                                                     

Алгоритм Рунге-Кутты четвертого порядка- РК4 (погрешность порядка h⁴):

                                                

где                                                     

Алгоритмы третьего и четвертого порядков требуют на каждом шаге трех и четырех  вычислений функции соответственно, но являются весьма точными.

Практическая часть

Задание 1. Необходимо вычислить значение интеграла от заданной функции f(x) на отрезке [a,b] при делении отрезка на 10 равных частей по формуле левых, правых и средних прямоугольников; по формуле трапеций; по формуле Симпсона. Произвести оценку погрешности методов интегрирования и сравниИзм.

Лист

№ докум.

Подпись

Дата

Лист

14

 

КП.230100. 08 .ПЗ

ть точность полученных результатов. Написать программу вычисления интеграла этой же функции по формулам левых, правых и средних прямоугольников, трапеций, Симпсона.

Рис.2. Задание 1 а).

 

Рис.3. Задание 1 б).

Сравнивая полученные результаты, сделанныеИзм.

Лист

№ докум.

Подпись

Дата

Лист

15

 

КП.230100. 08 .ПЗ

 в Exel и результаты разработанной программы можно прийти к выводу, что программа написана, верно. Анализируя погрешность методов можно заключить, что метод Симпсона самый точный из представленных методов.

Блок схема задания 1:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Задание 2. ВычисИзм.

Лист

№ докум.

Подпись

Дата

Лист

16

 

КП.230100. 08 .ПЗ

лить значение производной функции, заданной таблично, используя любые инструментальные средства. Написать программу вычисления значения производной этой же функции, используя интерполяционную формулу Ньютона.

Рис.4. Задание 2 а).

Рис.5. Задание 2 б).

Проанализировав данные, приходим к выводу, что интерполяционная формула  Ньютона достаточно точная.

 

Изм.

Лист

№ докум.

Подпись

Дата

Лист

17

 

КП.230100. 08 .ПЗ

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Задание 3. РешиИзм.

Лист

№ докум.

Подпись

Дата

Лист

18

 

КП.230100. 08 .ПЗ

ть задачу Коши для дифференциального уравнения y=f(x,y) на отрезке [a,b] при заданном начальном условии и шаге интегрирования h, используя любые инструментальные средства, методами Эйлера и Рунге-Кутта 4-го порядка. Написать программу решения обыкновенных дифференциальных уравнений методами Эйлера и Рунге-Кутта 4-го порядка.

Рис.6. Задание 3 а).

Рис.7. Задание 3 б).

Результаты вычисленные в Exel совпадают с результатами программы. Погрешность вычисления методом Эйлера достаточно велика, для уменьшения погрешносИзм.

Лист

№ докум.

Подпись

Дата

Лист