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

Тема 1. Численные методы решения задачи Коши для обыкновенных дифференциальных уравнений Лекция 1.1. Постановка задачи численного решения обыкновенных дифференциальных уравнений Рассмотрим основные понятия, связанные с дифференциальными уравнениями вообще и с процессом поиска их решения в частности. Дифференциальные уравнения, как известно, являются функциональными уравнениями, т.е. их решением является функция, которая задает уравнение своими производными. Порядок дифференциального уравнения определяется максимальным порядком таких производных функций, входящих в это уравнений. Дифференциальные уравнения бывают двух видов: обыкновенные и в частных производных. Обыкновенным называется уравнение, у которого искомая функция – это «обыкновенная» функция, функция от одного аргумента y  yx  . Математическая запись ее производных имеет вид: y, y, ..., y ( k ) или dy d 2 y dk y . , , ... , dx dx 2 dx k В дифференциальном уравнении, заданном в частных производных, искомая функция зависит от нескольких (двух и более) аргументов. Например z  ux, y  – функция двух аргументов x и y . В этом случае ее дифференциалы обозначаются:  , uxy  , uyy  ,... или ux , uy , uxx u u  2u  2u  2u , , , , , ... . x y x 2 xy y 2 Рассмотрим обыкновенное дифференциальное уравнение первого порядка в общем виде: dyx   f  x, y ( x )  dx или в более кратко: dy  f  x, y  , dx где правая часть уравнения f x, y  – выражение, зависящее в общем случае и от аргумента x , и от самой функции y (x) . Общим решением данного уравнения является не одна функция, а целое семейство функций, которые отличаются друг от друга константой С . В самом деле, если от обеих частей уравнения взять неопределенный интеграл, то искомую функцию можно выразить в явном виде:  dy   f x, y dx  y( x)  F ( x)  C , где F (x) – первообразная для f x, y( x) . Интерпретировать графически общее решение можно как множество интегральных кривых на плоскости в Декартовой системе координат (см. рис. 1.1) Рис. 1.1. Геометрическая интерпретация общего решения ДУ Для того чтобы из общего множества решений выбрать одно единственное частное решение (соответствующее конкретному значению С ), необходимо задать начальное условие для искомой функции y (x) , т.е. определить значение y (x) в начальной точке x  x0 : yx0   y0 В этом случае возникает начальная задача (задача Коши) для дифференциального уравнения:  dy   f x, y ;  dx  y  x0   y0 . Такая задача предполагает нахождение одного частного решения, удовлетворяющего заданному начальному условию. Графически это означает, что из множества интегральных кривых нужно выбрать ту, которая проходит через точку x0 , y0  на плоскости (см. рис.1.2). Рис. 1.2. Геометрическая интерпретация частного решения ДУ При использовании численных методов решение задачи Коши – функцию yx  находят не в явном виде, представленном как алгебраическое выражение или формула, а в дискретном виде – как вектор значений y   y0 , y1 , ..., yn  , полученных в предварительно построенных узлах x0 , x1 , ..., xn  из отрезка x0 , xn  , левый конец которого задан значением x0 из начального условия. Таким образом, численно можно вычислять только одно частное решение, т.е. решать начальную задачу для дифференциального уравнения. На практике чаще всего встречаются задачи для систем дифференциальных уравнений (ДУ). Если рассматривается система обыкновенных ДУ, то имеем дело с векторфункциями одного аргумента. В векторном виде начальная задача для системы ДУ имеет вид: y x   f x, y x , y x0   y0 ,  y10   y1  x    f1         y 20   y2  x    f2  где y  x      , f   ...  , y 0   ...  ...       y  x   f   0  n   n  yn  Для решения систем ДУ с начальными условиями численными методами используются алгоритмы, в основе которых лежат те же принципы, что и для вычисления начальной задачи для одного ДУ. С другой стороны известно, что путем простой замены одно обыкновенное ДУ nого порядка можно свести к системе из n уравнений первого порядка. Пусть имеется ДУ порядка n:   dny  f x, y, y, ... , y ( n1) , n dx С начальными условиями (кроме самой функции в начальной точке нужно знать значения всех ее производных до n-1–го порядка включительно): yx0   y00 , yx0   y10 , …, y ( n1) x0   yn1,0 Путем (n-1) -ой замены: y  z1 , y  z 2 , … , y ( n1)  z n1 Получим систему из n обыкновенных ДУ первого порядка:  y   z1 z  z 2  1  z 2  z3  ...  z n 2  z n1   z n 1  f ( x, y, z1 , z 2 , ... , z n1 , ) с начальными условиями: yx0   y00 , z1 x0   y10 , …, z n1 x0   yn1,0 Таким образом, достаточно рассмотреть алгоритмы численных методов вычисления задачи Коши для обыкновенного ДУ первого порядка. А при необходимости выполнить их модификацию, для решения системы ДУ первого порядка или ДУ порядка n. Будем рассматривать одномерные случаи, то есть задачу (1)-(2) Задача (1-2) имеет единственное решение (выясним когда?) Теорема Пикара :   Пусть дана область G  x, y  x  x0  a, y  y0  b Если: Пусть в области G функция f x, y  удовлетворяет по переменной y условию Липшица: Существует константа Липшица L , что f x, y1   f x, y2   L  y1  y2 Липшица Если функция f по переменной y дифференцируема , то L  max G  f y Тогда: в области G решение задачи (1) - (2) существует и единственно. Методы решения задачи Коши. Делятся на следующие типы: - аналитические - численные Начнем рассмотрения с - неравенство Метод решения основанный на разложении функции в ряд yx   f x, y , yx0   y0 Разложим функцию yx  yx0  yx0  y k  x0  2 x  x0   x  x0   .....  x  x0 k  ..... yx   yx0   1! 2! k! y  x0   y 0 y x0   f x0 , y x0   f x0 , y0         y    y    f  x, y x   f x  f y y x  f x  f y f                   y   y    f x  f y f   f xx  f xy f  f y  f x  f y f    f xy  f yy        Выбираем определённое количество слагаемых ряда. Заменим в них производные соответствующими выражениями и получаем расчетную формулу. Пример:  y  x 2  y 2   y 0  2 Решение: y  x0   y 0  2 y  x0   f  x0 , y  x0   f 0,2  0 2  2 2  4    y  x    y  f x  f y f  2 x  2 y x 2  y 2     y  x0   2  0  2  2 0 2  2 2  4  4  16 yx   2  4 x  0  16 x  02  2  4 x  8 x 2 1 2 Метод Пикара (итерационный метод) y k 1 x   A  y k  x  Идея метода: Функция правой части f x, y   f1 x, y   f 2 x, y  так чтобы представление задачи Коши имело решение в произвольной достаточно малой функции g x   y x   f1 x, y   g x    y  x0   y 0 Итерационный процесс метода строится следующим образом y 0   x   y 0 Решается задача Коши, где в функцию f 2 подставляется начальное приближение:     1 1 0    y  x   f 1 x, y  f 2 x , y  x   y 1 x  - решение, полученное на первой итерации  1   y x0   y0 Снова решаем задачу Коши   2  2    y x   f1 x, y  f 2 x, yx   y 2  x  - решение, полученное на второй итерации  2    y x0   y0     k 1 k 1   f 2 x, y  k   x   y  x   f 1 x, y В общем случае будем иметь:  k 1   y  x0   y 0 В результате получим функциональную последовательность: y0  y 0  x , y 1 x , y 2  x ,........., y k  x ,.....  y * x  , которая сходится к решению y * x  Если функцию f x, y  не расщеплять на 2 слагаемых , (на практике часто это так), то получим:    y k 1 x   f x, y k  x    f  x0   y 0 Проинтегрируем полученное равенство: x x x0 x0   x   k 1 k  k 1 k 1 k   y t dt   f t, y t  dt  y x   y x0    f t, y t  dt x0 Получим формулу Пикара: x  k 1 t  k t t  y x   y0   f x , y x  dx  x0  0   y x   y0  Пример:   y  x 2  y 2   y 0  2 y 0   x   2 x  t3  x3 y 1 x   2   t 2  4 dt  2    4t  0x  2  4 x  3 3  2 x  t3   2  2 y x   2   t   2  4t    dt  ............ 3   0     Выясним сходимость метода: x   y k 1 x   y0   f x, y k  x  dx функциональная - последовательность: x0 y1 x , y2 x ,........., yk x  Оценим:  x  x x0 x     y0   f x, y k  x  dx  y0   f x, y k 1 x  dx = y k 1 x   y k  x     f x, y k   f x, y k 1 x0  dx      f x, y   f x, y  x  k 1 k x0 dx  (для f применим x0 x неравенство Липшица)  L  y k  x   y k 1 x  dx x0 То есть разность на k – ом шаге может быть определена через разность на (k-1) Оценим разность при k=1 y 1 x   y x   0  y 1  x   y0 x   f x, y dx x0 x   f  x, y  dx  x0 A  max  x  f x, y 0  x   Ax  x0  При k=2 y 2  x x   y x   L  y 1 x0 При k=3 2  x  x0  dx  A  L   x  x0 dx  A  L  x 1 y x0 2  A L2 x  x0  L 2! x x x0 x0 y 3 x   y 2  x   L  y 2  y1 dx  A  L2   y k 1 A Lk 1 x  x0  x   y x   k  1! L k  k 1 x  x0 2 dx  2! A A L3 x  x0  3  L2  x  x0   23 L 3! - (k+1) – й член ряда Тейлора для функции e L  x x0   x   x 2 x3  e  1  x    ........ k  0 2! 3!   k  Таким образом y k 1 x   y k  x    0 , то есть метод Пикара сходится Замечания: 1) При малом шаге метод Пикара достаточно быстро сходится, решение можно получить за 5 – 10 итераций с точностью до 15 – ти знаков после запятой. 2) Интегралы в формуле вносят погрешность необходимо согласовать сетку вычислений с сеткой узлов интегрирования. Метод Эйлера (ломаных) Решаем задачу Коши:  y x   f x, y 1   y x0   y0 2 Метод заключается в построении ломаной, которая приближалась бы к интегральной кривой. Пытаемся приблизить нашу нелинейную функцию к линейной, получим ломаную составленную из кусочков прямых. Существует несколько подходов к построению расчетной формулы метода Эйлера: 1 –й подход yx   yx0   yx0 x  x0  - уравнение касательной yx   yx0   f x0 , yx0 x  x0  yx   y0  f x0 , y0 x  x0  - уравнение касательной в точке x0 , y0  y1  y0  f x0 , y0 x1  x0  ----------в точке x1 , y1  Точка x1 , y1  берем за приближенную точку, так как через неё проходит S – ая кривая, близкая к исходной кривой. В ней появится касательная: y2  y1  f x1 , y1 h2 и так далее yk  yk 1  f xk 1 , yk 1 hk - общая формула Эйлера 2 – й подход (разностный) y  x   y0   yx0  x  x0   y x0  x  x0 2  .......... 1! 2! Возьмем 2 сл – х: yx   y0  f x, y x  x0  - формула Говорят, что формула Эйлера имеет 2 – ой порядок точности в точке x0 (то есть локально). Глобальный – (на единицу меньше) 1 – й порядок точности 3 – й подход (конечно - разностный) yx   yx   yx0   yx   yx0   yx  x0  x  x0 4 – й подход (интегральный) Возьмем формулу Пикара x yx   y0   f x, yx dx и на отрезке от x0 до x при? x0 а) Формулу левых прямоугольников yx   y0  f x0 , y0 x1  x0  - явная схема б) Формулу правых прямоугольников yx1   y0  f x1 , y1 x1  x0  - неявная формула в) Формула прямоугольников на двойном отрезке y2  y0  f x1 , y1 x2  x0  - формула будет работать на 2 – ном интервале Рассмотрим неявную формулу Эйлера в сочетании с явной y1  y0  f x0 , y0 h1 - подставим - y1  y0  f x, y0  f x0 , y0 h1  Можно снова подставить в неявную формулу и так далее Получаем следующую процедуру (метод Эйлера – Коши) y10   y0  h1  f x0 , y0  , y1k   y0  h1  f x1 , y1  По ней производить итерации по неявной формуле (несколько уточнений). Затем снова по явной и так далее. Рассмотрим формулу трапеций для интегрирования: y  x   y0  y1  y0  h1  f x0 , y0   f x1 , y1  2 h1  f x0 , y0   f x1 , y0  h1 f x0 , y0  2   - формула Хьюна Формула Хьюна имеет глобально 2 – порядок точности. Геометрическая интерпретация формулы Хьюна: Вычислительные действия формулы Хьюна: 1) k1  f x0 , y0  - тангенс угла наклона касательной в точке x0 , y0  2) k2  f x1 , y0  h1k1  1 1 3) k  k1  k 2 - средневзвешенный коэффициент 2 2 4) y1  y0  h1 k Методы Рунге – Кутта Формулы Эйлера и Хьюна относятся к группе методов Рунге – Кутта, группа работает по аналогичной схеме. Определение: Количество S  точек, в которых получают вычисления углов наклона k1 , k2 ,....., k s  называется стадийностью (или порядком), а метод S – стадийным. Так метод Хьюна является 2 – х стадийным методом Рунге – Кутта. Схема: 1) Вычисление k1 , k 2 ,....., k s  2) Вычисление средне – ого коэффициента по формуле: k  b1k1  b2 k 2  .......  bs k s , b1  b2  ......  bs  1 3) y1  y0  h1k Выбор y1 зависит от коэффициентов bi и в каких точках брать решение. Семейство 2 – х стадийных методов Рунге – Кутта В общем виде 2 – х стадийная схема имеет вид: yx   yx0   hb1 f x0 , y0   b2 f x0  ch, y  ahf x0 , y0  Схематически ее удобно представлять в виде: нулевая строка – 0 , то есть вычисляется в точке x0 , y0  Выпишем формулу для 4 – х параметрового семейства 2 – х стадийных методов (для увеличения точности расчетов). Будем считать, что величина шага h достаточно мала, чтобы считать ch - приращением x Разложим функцию  f x0  ch, y  ahf x0 , y0  в ряд Тейлора , ограничиваясь линейными   членами f x0 , y0   f xx0 , y ch  f y x, y ahf x, y   o h 2 подставим в (*):     hb  b  f x , y   hb cf   b af  f  yx   y0  h b1 f  b2  f  f xch  f yahf   O h3 y1  y0 1 2 2 x 2 y Но из ряда Тейлора: yx1   y0  fh   f x  f y f    h2  O h3 2 Соберем при одинаковых множителях коэффициенты:   b1  b2  1 b1  b2  1  2  h 1 2  h b2 c    b2 c   2 2 2 h 1  ab2  h 2 ab2   2  2 Получили 3 уравнения относительно 4 – х неизвестных парметров Пусть b2   - свободный параметр Тогда b1 1   , c  1 2 , a 1 2 Подставим эти коэффициенты в формулу парам-ого семейства (*):  h h   y1  y0  h 1    f x0 , y0   f  x0  , y0  f x0 , y0   2 2    Однопараметричное семейство h h   При   1 y1  y0  hf  x0  , y0  f x0 , y0  - формула 2 – го пр. 2 2    h 1 y1  y0   f x0 , y0   f x0  h, y0  hf x0 , y0  - формула Хьюна 2 2 Аналогичным образом построены S – стадийные методы Рунге – Кутта по следующему алгоритму: 1) Поиск k i : k1  f xk , yk  k 2  f xk  c1h, yk  ha11k1  k3  f xk  c2 h, yk  ha21k1  a22k2  …………………….. k s  f xk  cs1h, yk  has1,1k1  as1, 2 k 2  ......  as1,s1k s1  2) Вычисление средневзвешенного коэффициента: k  s  k  b1k1  b2 k2  .......  bs ks   b j  1;  akj  ck  j 1  j 1  3) yk 1  yk  hk Построение методов более высокого порядка проводится аналогично. Методы малого порядка на практике используются редко. Часто используются методы 4 – го порядка и более Геометрическая интерпретация: k1  f  xk , y k  : A h h   k 2  f  xk  , y k  k1  : B 2 2   : h h   k 3  f  xk  , y k  k 2  : C 2 2   k 4  f  xk  h, y k  hk3  : D 1 2 2 1 k  k1  k 2  k3  k 4 6 6 6 6 yk 1  yk  hk Методы Адамса Идея метода xk 1  y   f  x, y 1   y  xk   y0 2  y dx  xk xk 1  f  x, y  xk y k 1  y k  xk 1  f x, yx dx3 xk Заменить значения подинтегральной функции интерполяционным множеством, построенным по известным значениям функции f . Осуществить экстранополяцию многочлена (построить продолжение) и выполнить интегрирование. То есть методы Адамса основаны на интегрировании многочлена. Они делятся на 2 группы: 1) методы Адамса – Башфорта (экстраполяция для точки ) 2) методы Адамса – Моутона (экстраполяция для точки xk 1 ) (Явные методы) Для замены подинтегральной функции в (3) используем Pn x  - 2 – ую формулу Вспомним многочлен Ньютона I формула: Ищем многочлен Ньютона для f x, yx  f x  : 1 f x0  2 f x0  k f x0       x  x0 x  x1 .......x  xk 1  x  x  x  x x  x  .....  1 1!h1 2!h 2 k!h k Рассмотрим слагаемые: Pn x   f x0   0: f x0   f 0  0 f 0 - кон. раз 0 – го порядка 0 f 0  f 0 - раздел. к.р. 0 – го порядка (обозначается f x0  ) h0 1: 1 f 0  f1  f 0 - кон. раз 1 – го порядка f x0 , x  1 f 0 f1  f 0  - раздел. к.р. 1 – го порядка h h1 f f    f x0 , x, f x1 , x2 , f x2 , x3 ,..... f xn1 , xn   n1 n  h   2: 2 f 0  1 f1  1 f 2   f 2  f1    f1  f 0  кон. раз 2 – го порядка f1 f 2  2 f 0 f1  f 0 h h  f x, x2   f x0 , x f x0 , x1 , x2   2   2 h h h h …………………………………………………………………….. n: n f 0  n1 f1  n1 f 0 кон. раз n – го порядка n1 f1 n1 f 0  n1 n 1 n f 0 f x,..., xn   f x0 ,.., xn  h h f x0 , x1 ,..., xn   n   h h h Тогда f x0 , x1  x  x0   f x0 , x1 , x2  x  x0 x  x1   .....  f x0 , x1 ,..xn  x  x0 x  x1 ..x  xn1  1! 2! n! Pn x   f 0  II формула:  xn  : f xn1  2 f xn2  n f x0   x  xn   x  xn x  xn1   .....  x  xn x  xn1 .......x  x1  Pn x   f xn   1!h1 2!h 2 n!h n Через разд. к.р.: f xn1 , xn  x  xn   1 f xn2 , xn1 , xn x  xn x  xn1   ......  1 f x0 ,.., xn x  xn ....x  x1   1 2! n! f n  f n1 f n1  f n2  f n  f n1 1 h  x  xn    h x  xn x  xn1    fn  h 2! h f n  f n1 f n1  f n2 f n1  f n2 f n2  f n3   h h h h  1 h h    ................... 3! h Вспомнили. Pn  x   f n  Используем для интерполяции подинтегральной функции в (3) 2 – ую формулу многочлена Ньютона относительно точки xk yk 1  y k  xk 1  P x dx k xk Точность формулы Адамса – Башфорта (ее порядок) зависит от точности Pk x  (то есть от выбранного числа слагаемых) Если возьмем 1 слагаемое в (4): yk 1  yk  xk 1  f k dx  yxk   xk 1 xk  xk xk 1 f xk , yxk dx  yxk   f xk , yxk   dx  yk  f xk , yk xk 1  xk  то есть yk 1  yk  hf xk , yk  - формула Эйлера xk Таким образом, метод Эйлера – есть метод Адамса – Башфорта 1 – го порядка Возьмем 2 слагаемых в (4): y k 1  y k  xk 1  xk f  f k 1  x  xk dx  yk   fk  k h    y k  h1 f xk , y k    yk  f k f k 1  x  x k   h  2 2 xk 1 xk xk 1 xk 1 xk xk  P1 x dx   f k  f k 1  x  xk   h  f  xk , y k   f xk 1 , y k 1     y k  h f xk , y k     2    h 3 f k  f k 1  2 yk 1  yk  h 3 f xk , yk   f xk 1 , yk 1  - формула Адамса – Башфорта 2 – го порядка 2 3 слагаемых: y k 1  y k  xk 1  P x dx  y 3 xk k  xk 1 xk 1  P x dx   2 xk f  f k 1  f k 2 h  y k  3 f k  f k 1   k 2 2h 2  yk  xk  f k  f k 1 f k 1  f k 2   h h  2 h    x  x  xk 1 2 k    x  xk  x  xk 1 dx      h x  x k   xk 3 2   h 3 f k  f k 1   f k  2 f k 12  f k 2  x  xk   hx  xk   2 3 2 2h   xk 1 xk  3 3  3  h 3 f k  f k 1   f k  2 f k 12  f k 2  h  h   yk  h 3 f k  f k 1    f k  2 f k 1  2f k 2 5h  2 2 2 2h 2  6h  3 h  y k  23 f k  16 f k 1  5 f k 2  12  yk  Формула Адамса – Башфорта 3 – го порядка.