ODU2

Доказательство. Схема рассуждение здесь та же, что и в одномерном случае. Сначала

единственность неподвижной точки интегрального оператора A , откуда следует существование и единственность решения задачи Коши (3’), а также − его гладкость.

Теорема 2. (Следствие теоремы 1).Предположим, что в условиях задачи Коши (1)

Выберем положительное число M такое, что на множестве Q будут выполняться нера-

ди функций класса Cn .

Замечание. Рассмотрим задачу Коши (1) для линейного уравнения вида y n p1 x y n 1 ... pn 1 x y ' pn x y F x ,

30

Действительно, сейчас f x; y1, y2,..., yn F x p1 x yn ... pn 1 x y2 pn x y1 и

31

§10. Понятие о сеточных методах решения дифференциальных уравнений.

В предыдущих параграфах для получения приближенного решения ДУ мы уже использовали сначала изоклины затем − метод последовательных приближений.

Суть сеточных методов состоит в том, что пространство гладких функций на от-

ностных уравнений, а u x − точное решение исходной задачи Коши, следует доби-

Рассмотрим более подробно простейшую и самую старую из разностных схем − схему Эйлера (метод Эйлера) приближенного решения задачи Коши y ' f x, y , y x0 y0 :

С геометрической точки зрения на i-м шаге вычислений отрезок интегральной кривой заменяется отрезком касательной к той интегральной кривой, которая проходит через

точку xi , yi , а вся искомая кривая − ломанной. Данный метод обычно называют методом ломанных Эйлера.

32

Применим сначала эту схему к функции u x , являющейся точным решением задачи Коши, считая, что правая часть уравнения f x, y непрерывно дифференцируема в некоторой замкнутой области G , заведомо содержащей графики точного и при-

Доказательство. Отметим, прежде всего, что функция x ex 1 x имеет при x 0

единственную точку минимума, причём 0 0 и, следовательно, 1 x ex при x 0 .

Применим для доказательства леммы метод математической индукции. При k=0 доказываемое неравенство является очевидным. Предположим, что при некотором целом

Предположим еще, что ошибка округления при вычислениях удовлетворяет условию

Таким образом, скорость сходимости для схемы Эйлера равна O h (в предполо-

жении, что ошибка округления имеет порядок O h2 ).

33

Существует множество разностных схем с более высокой скоростью сходимости.

Схемы Рунге – Кутта. )

Построим схему вида yi 1 yi h f xi , yi f xi h, yi f xi , yi h , имеющую на точном решении задачи Коши u x порядок аппроксимации, равный O h3 . Если

Сравнивая разложения (1) и (2), видим, нужный порядок аппроксимации мы получим,

Приведем схему Рунге – Кутта, имеющую четвёртый порядок точности:

yi 1 yi h6 k1 2k2 2k3 k4 , где

k1 f xi , yi , k2 f xi h2 , yi h2 k1 , k3 f xi h2 , yi h2 k2 , k4 f xi h, yi hk3 .

Априорные оценки точности в последнем случае слишком громоздки, поэтому в прак-

тических расчётах эти оценки не используются. Удобнее сравнивать между собой ре-

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

_____________________________________________________

) Рунге (Runge - Tolmé, Carl David) − немецкий математик: 1856 – 1927. Кутта (Kutta, Martin Wilhelm) – немецкий инженер: 1867 – 1944.

34

Методы прогноза и коррекции.

Простейший из методов прогноза и коррекции состоит в том, что сначала для

относительно yi 1 . Если требуется более двух-трёх итераций, вместо этого приходится уменьшать шаг расчёта.

Замечание 1. Для того чтобы начать вычисления, требуется знать, кроме y0 , ещё и y1 .

Вычисление y1 с приемлемой точностью можно выполнить, например, методом Рунге – Кутта 4-го порядка.

Этот метод удобен ещё и тем, что при достаточной гладкости правой части уравнения, т.е. функции f x, y , оценка погрешности может быть оценена с помощью следующей приближенной формулы:

Приведем более точную схему прогноза и коррекции − схему Адамса: )

_________________________________________________________

) Адамс (Adams, John Couch) − английский астроном: 1819 – 1892.

35

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

Замечание 3. Для решения уравнения, разрешенного относительно старшей производной, можно заменить уравнение нормальной системой уравнений первого порядка.

Пример. На отрезке 0; 3 с шагом 0.1 найти приближённое решение задачи Коши

y '' e 2 y ; y 0 0, y ' 0 0

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

Приведём простую программу ode, где для решения данной задачи на первом шаге используется метод Рунге – Кутта, а на остальных шагах − схема прогноз - коррекция (в виде обращения к подпрограммам rk и prcor). Вычисление вектора

помещен в приложении к параграфу.

program ode

dimension y0(2),y1(2),y(2),f1(2),f(2) external p

open(unit=1, file='ode.res', mode='write',status='unknown') eps=0.1e-03

y1(1)=0.

y1(2)=0.

y0(1)=0.

y0(2)=0.

x=0.

h=0.1

call rk(2,x,y0,h,p,f) write(1,2) x,(y0(j) , j=1,2) do 1 k=1,29

call prcor(2,x,y1,y0,h,p,f,eps,f1,y) write(1,2) x,(y0(j) ,j=1,2)

1continue

2format(2x,f5.2,2x,f9.6,2x,f9.6)

close(1)

stop

end program ode

Таблица, образовавшаяся в файле ‘ode.res’ и содержащая три столбца результатов

ещё, что точным решением задачи является функция ln ch x . Её график в рассматриваемом диапазоне значений аргумента полностью сливается с графиком построенного приближенного решения.

36

Приложение. Подрограммы rk, prcor и p.

subroutine rk(n,x,y,h,p,f) dimension f(1),y(1),y1(10),fk(4,10) call p(x,y,f)

do 1 i=1,n fk(1,i)=h*f(i)

1y1(i)=y(i)+0.5*fk(1,i)

x=x+0.5*h call p(x,y1,f) do 2 i=1,n fk(2,i)=h*f(i)

2y1(i)=y(i)+0.5*fk(2,i) call p(x,y1,f)

do 3 i=1,n fk(3,i)=h*f(i)

3y1(i)=y(i)+fk(3,i)

x=x+0.5*h call p(x,y1,f) do 4 i=1,n fk(4,i)=h*f(i)

d=fk(1,i)+2.*fk(2,i)+2.*fk(3,i)+fk(4,i) d=d/6.

4y(i)=y(i)+d

return end

37

subroutine prcor(n,x,y1,y0,h,p,f,eps,f1,y) dimension y0(1),y1(1),y(1),f1(1),f(1) call p(x,y0,f)

x=x+h

h2=0.5*h do 1 j=1,n

y(j)=y1(j)+2.*h*f(j)

1y1(j)=y0(j)

3err=0.

call p(x,y,f1) do 2 j=1,n

y0(j)=y1(j)+h2*(f(j)+f1(j)) ex=abs(y0(j)-y(j)) if(ex.gt.err) err=ex

2y(j)=y0(j)

if(err.gt.eps) go to 3 return

end

subroutine p(x,y,f) dimension y(2),f(2) f(1)=y(2) f(2)=exp(-2*y(1)) return

end

38