Решение систем дифференциальных уравнений методом Рунге-Кутты 4 порядка

2. Суть метода

Разбор и рассмотрение методов, применяемых на практике для решения дифференциальных уравнений, мы начнем с их широкой категории, известной под общим названием методов Рунге-Кутта.

Методы Рунге-Кутта обладают следующими свойствами:

1. Эти методы являются одноступенчатыми: чтобы найти уm+1, нужна

информация о предыдущей точке xm, ym.

2. Они согласуются с рядом Тейлора вплоть до членов порядка hp, где степень р

различна для различных методов и называется порядковым номером или

порядком метода.

3. Они не требуют вычисления производных от f (x, y), а требуют вычисления

самой функции.

Рассмотрим сначала геометрическое построение и выведем некоторые формулы на основе геометрических аналогий. После этого мы подтвердим полученные результаты аналитически.

Предположим, нам известна точка xm, ym на искомой кривой. Тогда мы можем провести прямую линию с тангенсом угла наклона уў m=f (xm, ym), которая пройдет через точку xm, ym. Это построение показано на рис. 1, где кривая представляет собой точное, но конечно неизвестное решение уравнения, а прямая линия L1 построена так, как это только что описано. {4}

Тогда следующей точкой решения можно считать ту, где прямая L1 пересечет ординату, проведенную через точку x=xm+1=xm+h.

Уравнение прямой L1 выглядит так: y=ym+yў m(x-xm) так как yў =f (xm, ym) и кроме того, xm+1=xm+h тогда уравнение примет вид

ym+1=ym+h*f (xm, ym) 1. 1

Ошибка при x=xm+1 показана в виде отрезка е. Очевидно, найденное таким образом приближенное значение согласуется с разложением в ряд Тейлора вплоть до членов порядка h, так что ошибка ограничения равна et=Кh2

Заметим, что хотя точка на графике 1 была показана на кривой, в действительности ym является приближенным значением и не лежит точно на кривой.

Формула 1. 1 описывает метод Эйлера, один из самых старых и широко известных методов численного интегрирования дифференциальных уравнений. Отметим, что метод Эйлера является одним из методов Рунге-Кутта первого порядка.

Рассмотрим исправленный метод Эйлера и модификационный метод Эйлера. В исправленном методе Эйлера мы находим средний тангенс угла наклона касательной для двух точек: xm, ym и xm+h, ym+hyў m. Последняя точка есть та самая, которая в методе Эйлера обозначалась xm+1, ym+1. Геометрический процесс нахождения точки xm+1, ym+1 можно проследить по рис. 2. С помощью метода Эйлера находится точка xm+h, ym+hyў m, лежащая на прямой L1. В этой точке снова вычисляется тангенс, дает прямую L. Наконец, через точку xm, ym мы проводим прямую L, параллельную L. Точка, в которой прямая L пересечется с ординатой, восстановленной из x=xm+1=xm+h, и будет искомой точкой xm+1, ym+1.

Тангенс угла наклона прямой L и прямой L равен

Ф (xm, ym, h)=Ѕ [f (xm, ym)+f (xm+h, ym+yў mh)] 1. 2

где yў m=f (xm, ym) 1. 3

Уравнение линии L при этом записывается в виде

y=ym+(x-xm)Ф (xm, ym, h),

так что

ym+1=ym+hФ (xm, ym, h). 1. 4

Соотношения 1. 2, 1. 3, 1. 4 описывают исправленный метод Эйлера. {5}

Чтобы выяснить, насколько хорошо этот метод согласуется с разложением в ряд Тейлора, вспомним, что разложение в ряд функции f (x, y) можно записать следующим образом:

f (x, y)=f (xm, ym)+(x-xm)¶ f/¶ x+(y-ym)¶ f/¶ x+ј 1. 5

где частные производные вычисляются при x=xm и y=ym.

Подставляя в формулу 1. 5 x=xm+h и y=ym+hyў m и используя выражение 1. 3 для yў m, получаем

f (xm+h, ym+hyў m)=f+hfx+hffy+O (h2),

где снова функция f и ее производные вычисляются в точке xm, ym. Подставляя результат в 1. 2 и производя необходимые преобразования, получаем

Ф (xm, ym, h)=f+h/2(fx+ffy)+O (h2).

Подставим полученное выражение в 1. 4 и сравним с рядом Тейлора

ym+1=ym+hf+h2/2(fx+ffy)+O (h3).

Как видим, исправленный метод Эйлера согласуется с разложением в ряд Тейлора вплоть до членов степени h2, являясь, таким образом, методом Рунге-Кутты второго порядка.

Рассмотрим модификационный метод Эйлера. Рассмотрим рис. 3 где первоначальное построение сделано так же, как и на рис. 2. Но на этот раз мы берем точку, лежащую на пересечении этой прямой и ординатой x=x+h/2. На рисунке эта точка образована через Р, а ее ордината равна y=ym+(h/2)yў m. Вычислим тангенс угла наклона касательной в этой точке

Ф (xm, ym, h)=f+(xm+h/2, ym+h/2*yў m), 1. 6

где yў m=f (xm, ym) 1. 7

Прямая с таким наклоном, проходящая через Р, обозначена через L *. Вслед за тем, мы проводим через точку xm, ym прямую параллельную L *, и обозначаем ее через L0. Пересечение этой прямой с ординатой x=xm+h и даст искомую точку xm+1, ym+1. Уравнение прямой можно записать в виде y=ym+(x-xm)Ф (xm, ym, h),

где Ф задается формулой 1. 6. Поэтому

ym+1=ym+hФ (xm, ym, h) 1. 8

Соотношения 1. 6, 1. 7, 1. 8 описывают так называемый модификационный метод Эйлера и является еще одним методом Рунге-Кутта второго порядка. Обобщим оба метода. Заметим, что оба метода описываются формулами вида

ym+1=ym+hФ (xm, ym, h) 1. 9

и в обоих случаях Ф имеет вид

Ф (xm, ym, h)=a1f (xm, ym)+a2f (xm+b1h, ym+b2hyў m), 1. 10

где yў m=f (xm, ym) 1. 11

В частности, для исправленного метода Эйлера

a1=a2=½;

b1=b2=1. {6}

В то время как для модификационного метода Эйлера

a1=0, a2=1,

b1=b2=½.

Формулы 1. 9, 1. 10, 1. 11 описывают некоторый метод типа Рунге-Кутты. Посмотрим, какого порядка метод можно рассчитывать получить в лучшем случае и каковы допустимые значения параметров a1, a2, b1 и b2.

Чтобы получить соответствие ряду Тейлора вплоть до членов степени h, в общем случае достаточно одного параметра. Чтобы получить согласование вплоть до членов степени h2, потребуется еще два параметра, так как необходимо учитывать члены h2fx и h2ffy. Так как у нас имеется всего четыре параметра, три из которых потребуются для создания согласования с рядом Тейлора вплоть до членов порядка h2, то самое лучшее, на что здесь можно рассчитывать — это метод второго порядка.

В разложении f (x, y) в ряд 1. 5 в окрестности точки xm, ym положим x=xm+b1h,

y=ym+b2hf.

Тогда f (xm+b1h, ym+b2hf)=f+b1hfx+b2hffy+O (h2), где функция и производные в правой части равенства вычислены в точке xm, ym.

Тогда 1. 9 можно переписать в виде ym+1=ym+h[a1f+a2f+h (a2b1fx+a2b2ffy)]+O (h3).

Сравнив эту формулу с разложением в ряд Тейлора, можно переписать в виде

ym+1=ym+h[a1f+a2f+h (a2b1fx+a2b2ffy)]+O (h3).

Если потребовать совпадения членов hf, то a1+a2=1.

Сравнивая члены, содержащие h2fx, получаем a2b1=½.

Сравнивая члены, содержащие h2ffy, получаем a2b2=½.

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

Положим, например, a2=w № 0. тогда a1=1-w, b1=b2=½w и соотношения 1. 9, 1. 10, 1. 11 сведутся к

ym+1=ym+h[(1-w)f (xm, ym)+w f (xm+h/2w, ym+h/2w f (xm, ym))]+O (h3) 1. 12

Это наиболее общая форма записи метода Рунге-Кутта второго порядка. При w =½ мы получаем исправленный метод Эйлера, при w =1 получаем модификационный метод Эйлера. Для всех w, отличных от нуля, ошибка ограничения равна