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

ПРИБЛИЖЕННОЕ РЕШЕНИЕ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ И СИСТЕМ. Дифференциальные уравнения являются основным математическим инструментом моделирования и анализа разнообразных явлений и процессов в науке и технике. Методы их решения подразделяются на два класса: 1) аналитические методы, в которых решение получается в виде аналитических функций; 2) численные (приближенные) методы, где искомые интегральные кривые получают в виде таблиц их числовых значений. Применение аналитических методов позволяет исследовать полученные решения методами математического анализа и сделать соответствующие выводы о свойствах моделируемого явления или процесса. К сожалению, с помощью таких методов можно решать достаточно ограниченный круг реальных задач. Численные методы позволяют получить с определенной точностью приближенное решение практически любой задачи. Решить дифференциальное уравнение y  f  x,y  (8.1) численным методом означает, что для заданной последовательности аргументов x0 , x1 , x2 , . . . , xn и числа y0  y  x0  , не определяя аналитического вида функции y  F  x  , найти значения y1 , y2 , . . . ,yn , удовлетворяющие условиям: F  x0   y0 , yk  F  xk  , k  1 , n . Рассмотрим три наиболее распространенных при решении практических задач численных метода интегрирования: метод Эйлера, метод РунгеКутта и метод Адамса. МЕТОД ЭЙЛЕРА Этот метод обладает малой точностью и применяется в основном для ориентировочных расчетов. Однако идеи, положенные в основу метода Эйлера, являются исходными для ряда других численных методов. Пусть дано дифференциальное уравнение с начальным условием (задача Коши) y  f  x,y  , y0  y  x0  (8.2) и выполняются условия существования и единственности решения. Требуется найти решение y  x  задачи Коши (8.2). Выбрав шаг h - достаточно малый, равный h   b  a  n , строим систему равноотстоящих точек x0 , x1 , . . . , xn , xi  x0  ih, i  0,n . 1 Искомую интегральную кривую y  y  x  , проходящую через точку M 0  x0 , y0  , приближенно заменим ломаной Эйлера с вершинами M i  xi yi  , i  0,1,2... (Рис.8.1). y y2 M 1 y1 M0  y0 h x0 tg  y0  y1  y0 h h x1 x2 x Рис.8.1 Звено ломаной M i M i 1 , заключенное между xi и xi 1 , наклонено к оси OX под углом  . Тангенс этого угла вычисляется по формуле: y y tg i  i 1 i  y  xi   f  xi , yi  . h Сделав преобразование, получим формулу Эйлера: yi 1  yi  h f  xi , yi  , i =0, n . (8.3) Вычисление значений y1, y2 , . . . , yn осуществляется с использованием формулы (8.3) следующим образом. По заданным начальным условиям x0  a и y0 полагая i  0 в выражении (8.3) вычисляется значение y1  y0  h f  x0 , y0  . (8.4) Далее определяя значение аргумента x по формуле x1  x0  h , используя найденное значение y1 и полагая в формуле (8.3) i  1 вычисляем следующее приближенное значение интегральной кривой y  F  x  : y2  y1  h f ( x1, y1) . (8.5) Поступая аналогичным образом при i  2, n  1 определяем все yi , в том числе последнее значение остальные значения yn  yn 1  h f  xn 1, yn 1  , которое соответствует значению аргумента xn  b . Таким образом, соединяя на координатной плоскости точки  x0 , y0  ,  x1, y1  , . . . ,  xn ,yn  отрезками прямых, получаем ломаную линию с вершинами в точках M 0  x0 , y0  , M1  x1, y1  , . . . ,M n  xn , yn  . Запишем разложение yi 1 в ряд Тейлора: Как правило, шаг h выбирают таким образом, чтобы h2   , где  2 заданная точность. Метод Эйлера может быть применен к решению систем дифференциальных уравнений. Пусть задана система двух дифференциальных уравнений первого порядка  y  f1 ( x, y, z ), (8.6)   z  f 2 ( x, y, z ), с начальными условиями y  x0   y0 , z  x0   z0 . Необходимо найти решение этой задачи Коши. Проводя аналогичные рассуждения, получаем расчетные формулы вида: yi 1  yi  hf1  xi , yi , zi  , zi 1  zi  hf 2  xi , yi , zi  , (8.7) xi 1  xi  h , i  0,1,2,..., где h  шаг интегрирования. В результате применения расчетной схемы (8.7) получается приближенное представление интегральных кривых y  F1  x  и z  F2  x  в форме двух ломаных Эйлера, построенных по полученным точкам  xi , yi  ,  xi , zi  , i  0,1,2,... . Достоинством метода Эйлера является его простота и высокая скорость поиска решения. Недостатками метода Эйлера являются малая точность и систематическое накопление ошибок, так как при вычислении значений на каждом последующем шаге исходные данные не являются точными и содержат погрешности, зависящие от неточности предшествующих вычислений. МЕТОД РУНГЕ-КУТТА(Ы) Данный метод является одним из наиболее распространенных численных методов интегрирования обыкновенных дифференциальных уравнений. По сравнению с описанным выше методом Эйлера метод РунгеКутта имеет более высокую точность, но невысокую скорость поиска решения, так как метод относится к классу многошаговых методов. Пусть дано дифференциальное уравнение первого порядка y  f  x,y  , (8.8) с начальным условием y0  y  x0  . (8.9) Выберем шаг h и для краткости введем обозначения xi  x0  ih , yi  y  xi  , где i  0,1,... . Рассмотрим числа: 3 k1i  hf  xi , yi  , k  h  k2i  hf  xi  , yi  1i  , 2 2   k  h  k3i  hf  xi  , yi  2i  , 2 2   (8.10) k4i  hf  xi  h, yi  k3i  . По методу Рунге-Кутта последовательные значения yi искомой функции y определяются по формуле: 1 yi 1  yi   k1i  2k2i  2k3i  k4i  . (8.11) 6 Погрешность метода Рунге-Кутта, заданного формулой (8.11), на каждом шаге есть величина порядка h 5 (в предположении, что f  x, y   C   ). Формулу (8.11) еще называют формулой Рунге-Кутта четвертого порядка точности. Помимо формулы (8.11) существуют еще другие формулы типа Рунге-Кутта(ы) с иными порядками точности. В частности формула yi 1  yi  k2i  формула Рунге-Кутта второго порядка точности. Эта фор5 мула на каждом шаге дает погрешность порядка h 3 . Для определения правильности выбора шага h на практике обычно на каждом этапе из двух шагов применяют двойной пересчет. Исходя из текущего верного значения y  xi  , вычисляют y  xi  2h  двумя способами: вначале с шагом h , а затем с шагом 2h . Если расхождение полученных результатов не превышает допустимой погрешности, то шаг h для данного этапа выбран правильно и полученное с его помощью значение можно принять за y  xi  2h  . В противном случае шаг уменьшается в два раза. Эту вычислительную схему легко запрограммировать на ЭВМ. Метод Рунге-Кутта(ы) может быть использован и при решении систем дифференциальных уравнений. Рассмотрим задачу Коши для системы двух дифференциальных уравнений:  y  f1 ( x, y, z ),   z  f 2 ( x, y, z ), с начальными условиями y  x0   y0 , z  x0   z0 . Формулы метода Рунге-Кутта для данной системы примут вид: 4 1  k1i  2k2i  2k3i  k4i  , 6 1 zi 1  zi   m1i  2m2i  2m3i  m4i  , i  0, n. 6 yi 1  yi  где k1i  hf1  xi , yi , zi  , m1i  hf 2  xi , yi , zi  , k2i  hf1  xi  0.5h ,yi  0.5k1i ,zi  0.5m1i  , m2i  hf 2  xi  0.5h ,yi  0.5k1i ,zi  0.5m1i  , k3i  hf1  xi  0.5h ,yi  0.5k2i ,zi  0.5m2i  , m3i  hf 2  xi  0.5h ,yi  0.5k2i ,zi  0.5m2i  , k4i  hf1  xi  h ,yi  k3i ,zi  m3i  , m4i  hf 2  xi  h ,yi  k3i ,zi  m3i  . Метод Рунге-Кутта(ы) обладает значительной точностью и, несмотря на свою трудоемкость, широко используется при численном решении дифференциальных уравнений и систем. Важным преимуществом этого метода является возможность применения переменного шага, что позволяет учитывать локальные особенности искомой функции. МЕТОД АДАМСА Этот метод численного интегрирования разработан Адамсом в 1855 году по просьбе известного английского артиллериста Башфорта, который занимался баллистикой. В последствии этот метод был забыт и вновь открыт в начале XX века норвежским математиком Штермером. Популяризация метода Адамса и дальнейшее его усовершенствование связано с именем А.Н. Крылова. Пусть дано дифференциальное уравнение первого порядка y  f  x,y  , (8.12) с начальным условием y0  y  x0  . (8.13) Пусть xi , i  0,1,...  система равноотстоящих значений с шагом h и yi  y  xi  . Очевидно, что xi 1 yi   y  x  dx . (8.14) xi Запишем вторую интерполяционную формулу Ньютона с точностью 5 до разностей четвертого порядка: q  q  1 2 q  q  1 q  2  3 y  x   yi  qyi 1   yi  2   yi 3 , (8.15) 2! 3! x  xn где q  . h В формуле (8.15) функцию y заменим на производную y , получим: q  q  1 2 q  q  1 q  2  3 y  x   yi  qyi1   yi 2   yi3 . (8.16) 2! 3! Так как dx  h  dq , то подставив (8.16) в (8.14), получим: 1  q2  q 2 q3  3q 2  2q 3 yi  h   yi  qyi1   yi 2   yi3  dq . 2 6 0  После преобразований будем иметь: 1 5 3 yi 1  yi  hyi    hyi1    2  hyi 2    3  hyi3  . (8.17) 2 12 8 Формула (8.17) называется экстраполяционной формулой Адамса. Для начала итерационного процесса нужно знать начальные значения y0 , y1, y2 , y3 , так называемый начальный отрезок, который определяют, исходя из начального условия (8.13), каким-либо численным методом (например, методом Рунге-Кутта). Зная значения y0 , y1, y2 , y3 из (8.12) находят y0 , y1 , y2 , y3 и составляют таблицу разностей:   hy0  ,   hy1  ,   hy2  ,  2  hy0  ,  2  hy1  , 3  hy0  . (8.18) Дальнейшие значения yi , i  4,5,... искомого решения можно шаг за шагом вычислять по формуле Адамса (8.17), пополняя по мере необходимости таблицу разностей (8.18). Для работы на ЭВМ формулу Адамса применяют в раскрытом виде. Так как yi1  yi  yi1,  2 yi 2  yi  2 yi1  yi 2 , 3 yi3  yi  3 yi1  3 yi 2  yi3 , то после приведения подобных членов имеем: h yi 1  yi   55 yi  59 yi1  37 yi 2  9 yi3  , (8.19) 24 xi 1  xi  h. На практике шаг h выбирают так, чтобы можно было пренебречь величиной 1 3   hyi 2  . 24 6 Метод Адамса легко распространяется на системы дифференциальных уравнений. Погрешность метода Адамса имеет тот же порядок, что и метод Рунге-Кутта. 7