Модифицированный метод Эйлера с итерационным уточнением и переменным шагом Текст научной статьи по специальности «Математика»

УДК 512.1

МОДИФИЦИРОВАННЫЙ МЕТОД ЭЙЛЕРА С ИТЕРАЦИОННЫМ УТОЧНЕНИЕМ

И ПЕРЕМЕННЫМ ШАГОМ

С.В. Трубников

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

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

Введение

Значительную часть прикладных математических задач составляют краевые задачи для дифференциальных уравнений. Поэтому одной из центральных проблем современной прикладной математики является разработка и исследование численных методов решения подобных краевых задач. Этой проблеме посвящена обширная литература (см., например, библиографию в [1], [2] и [3]). В статье [5] описан новый подход к построению численных методов решения задач Коши для обыкновенных дифференциальных уравнений, основанный на аппроксимации приближенного решения с помощью кусочно-полиномиальной интерполяции многочленами Эрмита, а также на принципе минимизации невязки. С помощью этого подхода был построен новый численный метод для решения задачи Коши, названный исправленным методом Эйлера с итерационным уточнением, который можно отнести к итерационно-разностным методам. В нем используется сетка с постоянным шагом. В данной статье описан новый численный метод с переменным шагом сетки и описана процедура автоматического выбора узлов сетки, обеспечивающих заданную точность. Этот метод также относится к итерационно-разностным методам.

1. Составная функциональная кинематическая кривая и ее свойства

Кинематические кривые введены и описаны в статье [4]. Функции, задающие составные кинематические кривые, представляют собой результат кусочно-многочленной интерполяции многочленами Эрмита 5 порядка с двумя трехкратными узлами. Уравнение составной кинематической кривой на плоскости хОу можно записать в виде:

2 5

г = г(г) = ^ Т5 (г - г +1)-Оу,_1 + ^ Т5 (г - г +1)-Q5-ji при г е [/-1,/], г = 1,2,* , т . (1)

}=0 j=3

Здесь т - заданное натуральное число, О j г = Qxj г -1 + Qyji' '] ( j = 0,1,2, г = 0,1, * , т ) -

заданные векторы, Т5 (г) - интерполяционные многочлены Эрмита 5 порядка [2], которые можно записать в виде:

Т05 () = 1 -10*3 +15*4 - бг5, Т5(г)=г -бг3 + 8*4 -3*5, Т25(г) = 1 г2 -2г3 + 2г4 -2г5,

Т35 (г ) = 2 г3 - г4 + 2 г5, Т45 (г ) = -4г3 + 7г4 - 3г5, Т55 (г ) = 10г3 - 15г4 + бг5. (2)

Если ввести обозначения компонент векторной функции

г(г) = Гх (г)- 1 + Гу (г)-] , (3)

то равенство (1) можно записать в виде:

2 5

ГХ(і) = XТ! (і -і +1)‘ °х/г-1 + XГ/(і -і + !)• °х5-іі при іє[і -1,і], і = ІД* ,т • (4)

І=0 і=з

гу (і) = XТ/ (і-*' +^ 6_УІ і-1 + XТ/ (і-/ +^ Оу5-і і при і є[/-1,/], і = 1,2,‘ , т • (5)

j=0 j=3

Одно из главных свойств кинематической кривой [4] состоит в том, что

2

Гх (і) = 6x0 і , _ 6x1 і, ё [Х2 ) _ 6х2 і, і = 0,1,‘ ,т ,

йі йГ

Гу (і ) = Оу 0 і ,

йгу (і ) _^ &\ (і )

йі

= 6у1 і ,

йі

уу> _ 2

Оу 2 і, і = 0,1,‘ ,т •

Введем новую кривую, задаваемую уравнением (3) и уравнениями

Гх (і ) = Хі-1 +(хі - Хі-1 У(і - і + 1) при і є[і -1> і] , і = 1,2

(6)

(7)

(8)

гу (г) = X Т](г - г + 1)-^М + ХТ7(г - г + 1)-Qy/ 5-ji при г е[ - ^ г], г = 1,2,к , т . (9)

j=0 j=3

Здесь х{ (i = 0,1,* , т ) - заданные числа, такие, что

х0 < Х1 < ...Хт , (10)

а Qyrji-1 и Qylji - заданные постоянные.

Для функции, заданной формулой (9), будут справедливы условия, аналогичные условиям (7). Если считать, что i е [/' -1, /'], то

/\ йгу (і - 0) й2г (і -1)

(і)=Оуі0 і, і;---= <2уі 1 і, ----Т2-_ °уі2 і , і = 0,1,‘ ,т •

Гу\і ] =

йі ” йі А если считать, что і є [і, і +1], то

(11)

Гу (і ) = 6

уг 0 і

йгу(і + 0) _ й 2Гу (і +1)

'_ 6уг 1 і ,

= буг2 і, і = 0,1,К ,т •

(12)

йі йґ

Отсюда видно, что функция у = Гу (і) и её производные определяются неоднозначно в точках і = 1,2,_, т -1

Из выражения (8) и условия (10) следует, что функция х = гх (і) является возрастающей и, следовательно, обратимой Обратная функция определяется элементарно

і = гх^(х) = (і-1)+——при х є[хі-1 ,хі], і = 1,2,* ,т • (13)

хі - хі-1

Поэтому уравнения (8), (9) неявно задают функцию

у(х) = гу (гх( 1)(х))_ X ТІ

І=0

х-х

і-1

хі - хі-1

V ' ' 1 У

• 6уг/і-1 + X ТІ І=3

х-х

і-1

хі - хі-1 V ' ' 1

• °уІ 5-

Іі

при х е [хг-1 , х{ ], г = 1,2, • , т . (14)

Эта функция является однозначной во всех точках [х0 , хт ], кроме, вообще говоря, точек х{

(г = 1,2,* , т -1). Значения этой функции и её производных выражаются через значения гх (г)

и гу (г ) .

ёу(х)

ёх

= ^ /ёЛ), ёг / ёг

(15)

(16)

ё 2 у(х)

ёх

' ё2 Гу (г) ёгх () ёГу (г) ё2 гх () / ( ёгх ( )3

х = 1 а2 ёг ёг ёг2 V / 1 л )

(17)

х

х=Г„ х

х

Из (8) следует, что

Гх (г) = х , г = 0,1,К ,т . (18)

Если считать, что г е [г -1, г ], то

ёгх (г - 0)

\ = х - хг-1, г = 1,2,‘ , т . (19)

ёг

Если считать, что г е [г, г +1], то

т'х |/ + ^ = хг+1 - хг , г = 0,1,К , т - 1. (20)

ёг

Вторая производная

ё2гх (г - 0) . ё2гх (г + 0) .

----^—- =0, г = 1,2,* ,т; -----^—- =0, г = 0,1,* ,т -1. (21)

сИ2 ёг2

Из условий (11), (12), (15) - (21) следует, что

У(х ) = б ёУ(хг - 0) = ^уИ г ё2У(хг - 0) = @у12 г

ПХ)=Уу[ 0 г, ----------------------, ~Г2 7-\2 ,

ёх х - х-1 ёх2 (х,- - хг-1 )2

поскольку хг е [хг -1 , хг ]. С другой стороны

у(х ) = б ёу(хг + 0) = буг1 г ёу(хг + 0) = буг 2 г

Дх^=^уг 0 г, т 72,

^ х+1 - х ёх2 (х,-+1 - хг )2

поскольку х{ е [хг , х{+1 ]. Потребуем, чтобы постоянные Qyгji и 0У^1 удовлетворяли следующим условиям:

буг 0,=бу,0 г. , ( буг 2 г )2 = ( бу12 г )2 , г = 1,2,* , т -1. (22)

х+1 - хг хг - хг-1 (хг+1 - хг )2 (хг - хг-1 )2

Эти условия гарантируют однозначность и непрерывность функции у(х) = Гу (гх ( 1)( х)) и её производных до второго порядка включительно на [х0 , хт ] и выполнение следующих равенств.

/ ) = б ёу(Х0) = буг 1 0 а2у{х0) = Оуг2 0

У\х0)=буг 0 0, ; , ~2 7 \2

ах (х1 - Х0)

йх х1 - Х0

у(х ) = О йУ(хт ) = буЛ т й2У(хт ) = бу/2

у\хт) = 0у/0 т , , ~Т2 Т

йх хт - хт-1

йх 2 (хт - хт-1 )2 ’

у(хі ) = буг 0 і = бу/0 і , й2 у (хі ) _ буг 2 і = <2у/2 і

йу(хі)_ буг1 і = Єу/1 і

ах

(хі+1- хі)2 (хі- хі-1)2’

йх хі+1 - хі хі- хі-1’

і _ 1,2,* , т -1. (23)

і _ 1,2,* , т -1;

Кривую, заданную формулами (3), (8), (9) с коэффициентами, удовлетворяющими условиям (10) и (22), назовем составной функциональной кинематической кривой 5 порядка. Она представляет собой график однозначной и дважды непрерывно дифференцируемой на

[х0 , хт ] функции у(х)_ гу (гх( 1)(х)),

їх) , неявно задаваемой с помощью

функций гх (і) и гу (і).

2. Задача Коши и представление её приближенного решения

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

йу _ /(х, у), х є [а, Ь\ (24) ах

у(а)_ у, (25)

где а, Ь (а<Ь), у - заданные постоянные, / (х, у) - заданная функция такая, что существуют ее непрерывные частные производные до второго порядка включительно, а задача Коши (24), (25) имеет единственное решение, которое мы в дальнейшем будем называть точным решением и обозначать у(х).

Приближенное решение ут (х) задачи Коши (24), (25) будем задавать либо неявно, формулами вида (8), (9)

х _гтх(і) = хі-1 +(хі -хі-1 )'(і-і +1) при іє [і-1,і], і _1,2,к ,т , (26)

у _ гту (і ) = X (і - і + 1)' °угІ і-1 + X (і - і + 1)' °у/ 5-

І_0 і_3

при і є [і -1, і], і _ 1,2, * , т ,

либо явно, формулой вида (14)

] і

ут (х) гту (гтх (х)) Т)

І _ 0

х-х

і-1

хі - хі-1 V ' ' 1 У

■ °уг]і-1 + Х ТІ

І _3

х-х

і-1

хі - хі-1 V ' ' 1

■ бу/5-

І і

при х є \хі-1 , х}

[хі-1, хі], і _1,2, * , т.

(27)

(28)

Здесь х{ (г = 0,1,* ,т ) - заданные узлы сетки точек на [а ,Ь], удовлетворяющие условиям:

а = х0 < х1 < ...хт = Ь, (29)

т

а Qyгj г -1 и I - заданные постоянные, удовлетворяющие условиям (22), которые гаран-

тируют однозначность и непрерывность приближенного решения ут (х) = гту (гтх( 1)(х)) и его

производных до второго порядка включительно на [а , Ь] , а также выполнение равенств, аналогичных (23):

Фт (х0 ) = буг 1 0 ё 2 ут (х0 ) = буг 2 0

ут (х0 ) = буг 0 0 ■

ёх х1 - х0

ёх

(х1 - х0 )2

ут (хт ) бу10 т-

Фт (хт ) = бу11 т ё 2 ут (хт ) = 2

ёх хт - хт-1

ёх

у т (хг) буг 0 г бу10 г

ёут (хг ) = буг 1 г = бу11 г

ёх хг+1 - х х - х-1

(хт - хт -1 )2

г = 1,2,* ,т-1;

ё ут (хг )= буг 2 г = бу12,

г' = 1,2,* ,т-1.

ёх (хг+1 - х) (хг - хг-1У

Значения приближенного решения ут (х ) = гту (гтх ( 1)(х)) и его производных выражаются

рез значения гтх (г) и гту (г) .

ут (х ] х=

х=гпх (г) 'ту

(г),

ёут (х)

= ёгту(г) / ёгтх (г)

ту\ / I ^‘тх\

,=т () Л ' Л

ё 2ут (х )

С

с=гтх ()

ёх

^ ё2гту (г) . Фтх () - ёгту (г) _ ё2гтх ()

ёг2 ёг ёг ёг2

/ ёгтх ()

ч ёг

3

(30) че-

(31)

(32)

(33)

Для определения приближенного решения остается вычислить значения величин

Qугj 0 , Qу/jm ( j = 0,1,2) Qугji, Qу/ji ( г = 1,2,* , т -1, j = 0,1,2), которые удовлетворяют

условиям (22) и определяют функцию гту (г) и приближенное решение ут (х) = гту (гтх( 1)(х)).

4. Получение приближенного решения задачи Коши

Для определения неизвестных постоянных мы, так же как в работе [5], введем и будем использовать понятие невязки. Невязкой дифференциального уравнения (24) на его прибли-

женном решении ут (х) = гту (г,

-1)(х ))

ту тх мы назовем величину

Кт (х )

ёУт (х) - I(^ ут (х)), х е [a, Ь] •

ёх

(34)

Производная невязки

ёКт (х ) = ё 2 ут (х ) Э/ (Х ут (х )) Э/ (Х ут (х )) Фт (х )

ёх

Эх

Эу

, хе [а,Ь]. ёх

Учитывая (30), отсюда получим значение невязки и её производной в узлах хг:

Кт (х0

(х0 ) = ^у1° - 1 ^ ^г00 ^ Кт (хт ) = ----У (хт, ^0т )

х1 - х0

хт хт-1

т

Кт (х ) = -°у-~-/(х,Оуг0, ) = -/(х,Qуl0г), г = 1,2,- , т-1, (36)

хг+1 - хг хг - хг-1

ёКт (х0 ) = Qуг 2 0 - Э/(Хo_’Qут0o) Э/ (х0, ^ут 00 ) Qуг 1 0 ёх (х1 - х0 )2 Эх Эу х1 - х0 ’

ёКт (хт ) = Qуl 2 т Э/ (хт, ^у10 т ) Э/ (х0, Qу/ 0 т ) Qу/1 т

ёх (хт - хт-1 )2 Эх Эу хт - хт-1

ёЯт (х, ) = Qуг 2 г -Э/ІХi,Qуr0i]-Э/ІХi,Q^0i] ^г1 г =

ёх (х,+1 - х, )2 Эх Эу х,+1- х,

= Qуl 2 г -Э/(х^бу/0г)-Э/(iQd0i) Qуl 1 г

(х, - х,-1 )2 Эх эу хг - хг-1

г = 1,2,* ,т-1. (37)

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

Потребуем, чтобы для приближения ут (х) = гту (гтх^ ^(х)) в точках х, невязка и её

производная обращались в ноль. Исходя из этих требований и (37) получаются выражения

для Qуг 1 г, Qуl 1 г, Qуг 2 г и Qуl 2 г через Qуг 0 г и Qуl 0 г:

^ 10 =(х1 - х0)-у (x0, ^ 00 ),

^ 1 г = (хг -хг-1Уу(хг,^0г^ ^ 1 г = (хг+1 -хгУу(хг,^0г^ г = 1,2,* ,т-1,

Qуl 1 т = (хт хт-1 )-1 (хт, Qуl 0 т ) . (38)

( )2 Э/(х0,^°уг00 ) ( )Э/(х0, ^°уг00 ) ~

Qуг 20 =(х1 - х0 ) -----Э^----- ^ - х0 )-------Эу-------Qуг10,

^ ( )2 Э/ (х,, ^г 0,),( )Э/ (х,, ^г 0, )Q , 12 1

^г2, =(х,+1- х,)---------эг— (х,'+1- х,)----Эу— ,, г =1,2,1* , т -1,

^ ( )2 Э/(хг,Qуl0,) . ( )Э/(хг,Qуl0г) ^ 1 2 1

^ 2 г = (хг - х,-1) ----э^— ^хг - х,-1)----Эу------^ г, г = '* , т - ^

Q =( )2 Э/ (хт, Qуl 0 т )+( )Э/ (хт, Qуl 0 т ) Q (39)

Qуl 2 т = чхт - хт-1) + ^т - хт-1) Эу ^^у1 1 т . (39)

Зная Qуl 0, (, = 1,2,* ,т ) и Qуг 0, (, = 0,1,* ,т -1), по формулам (38) и (39) можно найти

^ 1, (г = 1,2,* , т X ^г 1 г (г = 0,1,* ,т -1), а затем ^ 2 г (г = 1,2,* , т X Qуг 2 г

(, = 0,1,* ,т -1). Таким образом, построение приближенного решения свелось к определению ^ 0 г (г = 1,2,* , т ) и ^г 0 г (г = °Л* ,т -1).

Неизвестное значение QуГ 0 0 мы найдем из начального условия (25) и условия (30):

^ 0 0 = ут (х0 ) = ут (а) = у . (40)

Остальные значения Пуг 0 і и Пу/ 0 і будем искать последовательно минимизируя квадраты

величин невязки в серединах отрезков [хі-1 , хі ]. Координаты середин хі _ (хі-1 + хі)/2. В результате возникает цепочка задач минимизации

(41)

Эти задачи минимизации заменяют требования Ят (х,) = 0 . Значение (Ят (х, ))2 зависит от Qуг 0,-1 и Qуl 0,. При г = 1 первая из этих величин известна: Qуг 0,-1 = Qуг 00 = у . Таким

образом, (Ят (х ))2 зависит только от одной неизвестной Qуl 01. Решив первую задачу минимизации (41) для (Ят (х))2 как функции одной переменной, мы найдем Qуl 01. Используя связи (22), мы найдем Qуг 01 = Qуl 01. Тогда (Ят (щ ))2 также будет зависеть только от одной неизвестной величины Qуl 02. Решив вторую задачу минимизации (41) для (Ят(щ))2 как функции одной переменной, мы найдем Qуl 0 2. Используя связи (22), мы найдем Qуг 02 = Qуl 02. И так далее. Продолжая этот процесс мы последовательно найдем все неизвестные величины.

На г -ом шаге описанного процесса решения задач минимизации вычисляется точка минимума (Ят (х, ))2 как функции одной переменной т = Qуl 0,. Представление для Ят (х,), как функции т = Qуl 0,, мы получим, подставив в выражение для невязки (34) представление

Т (12)

(28) и вычислив значения Т^ (12) и

I_____

йі

60

+ ■

32

(хі - хі

(хі хі-1 ) ^т (хі ) 22 Оуг 0 і-1 32 буг 1 і-1 32 Пуг 2 і-1 +

, )2 №■£). у (х„ т )"

Эх Эу

14 60

— (хі- хі-1> / (хі, т)+—■ т

/

-(хі- хі-1)■ /

V

1

+ — 64

хі + хі-1

(хі- хі -1)2

32 п 10 п 1 п

64 Пуг 01-1 + 64 Пуг 1і-1 + 64 Пуг 2г'-1

+

Эх

-+-

Эу

/ (хі, т)

10 32

-—(х- хі-1> / (хі, тН—■ т

64 64

(42)

Величины (Ят (х,)) при г = 1,2,* ,т представляют собой функции т . Обозначим их угт(т). Для определения очередного значения Qуl 0, необходимо найти точку минимума

функций у,т (т).

Уіт (т)® тіп _ °.

(43)

1

2

После чего значение Qyi 0 i полагается равным координате этой точки. Для решения этой задачи можно использовать различные итерационные численные методы минимизации [6]. Бу-

дем использовать для решения описанной одномерной задачи минимизации функции

Vim (m) разностный метод парабол. Важно подобрать начальное приближение moi искомой

точки минимума функции yim (m), достаточно близкое к этой точке, чтобы обеспечивалась

сходимость применяемой последовательности итераций. Для этого мы используем разложение функции ym (x) по формуле Тейлора с центром в точке x = xi_j. Запишем это разложение при х = х{:

Qyl 0 i = Ут (xi) » Ут (xi_1)+ dymd*l_ ) • (xi _ xi_1 ) +1 d Ут (2X?_ ) • (xi _ xi_1 )2 =

dx 2 dx2

= Qyr 0 i_1 + Qyr 1 i_1 + ^ Qyr 2 i_1 . (44)

Полученное приближенное значение Qyi o i выберем в качестве начального приближения

m0i = Qyr 0 i_1 + Qyr 1 i_1 + 2 Qyr 2 i_1. (45)

Последовательность приближений к точке минимума msi строится с помощью рекуррентной формулы разностного метода парабол:

Л Vim (Msi + Л)_ Vim (msi _Л) ni ^

ms+1 i = mSi _- —(———(—V+—(----------------n, s = 0,1v- (46)

2 Vim (ms i + Л)_ 2Vim т i )+ Vim Ws i _ Л)

Здесь Л - заданное фиксированное маленькое положительное число. Перед каждым вычислением очередного члена последовательности приближений по формуле (46) необходимо проверять условие

Vim (ms i + Л)_ 2Vim (ms i )+ Vim (Ws i _ Л) ^ d > 0, (47)

где d - заданное фиксированное маленькое положительное число. Если это условие не выполняется, то это означает, что вторая производная функции Vim (m) в точке msi либо отрицательна либо близка к нулю. Чаще всего такая ситуация возникает, когда начальное приближение слишком грубое. В этом случае следует прекратить вычисления. Улучшить качество начального приближения m0 i можно, выбрав большее значение величины m .

Если в качестве приближений для точек минимума выбрать m0 i, не проводя дальнейших итераций по формуле (46) (Qyi 0 i положить равным m0i), то в результате мы получим приближенное сеточное решения задачи Коши (24), (25): m0i = Qyl 0 i = ym (xi) (i = 1,2,* , m ).

Таким образом, формула (45) порождает вычислительную схему получения приближенного сеточного решения m0i задачи Коши (24), (25), которая в книге [3] названа исправленным

методом Эйлера. Шаговая погрешность этого метода составляет величину o(h3) при h ® 0, если точное решение задачи Коши, y(x), существует, единственно и трижды непрерывно дифференцируемо на [a, b], а функция f (x, y) и ее частные производные первого и второго порядка непрерывны и ограничены. При выполнении этих условий компоненты приближенного сеточного решения m0 i сходятся к компонентам точного сеточного решения при

m ® ¥.

Сходимость начальных приближений m0 i к точному сеточному решению при m ® ¥ и при x{ _ x{_1 ® 0 позволяет добиваться высокого качества начальных приближений за счет увеличения m (уменьшения xi _ xi_1). Итерации, проводимые по формуле (46) при выполнении условия (47), уменьшают невязку и, следовательно, уточняют значения Qyi 0 i. В пределе при m ® ¥ значения Qyi 0 i сходятся к компонентам точного сеточного решения. Поэтому невязка должна стремиться к нулю при m ® ¥ и в качестве условия окончания итераций можно использовать неравенство:

Vim (ms+1 i )£ 1 , (48)

где 1 - заданное фиксированное маленькое положительное число.

В любом случае количество итераций следует ограничить. Иначе количество вычислительных операций для получения результата может стать неоправданно большим. Поэтому мы введем величину S максимального количества итераций. Даже если условие (48) не будет выполнено, вычисления по формуле (46) прекратятся при s > S . Заметим, что если выбрать отрицательное значение 1, то условие (48) никогда не будет выполнено и количество итераций будет фиксированным и равным S, если не возникнет ситуация, когда будет нарушено условие (47).

Итак, мы получили в общих чертах метод вычисления приближенного решения задачи Коши (24), (25). Перечислим его основные этапы.

Компоненты приближенного решения задачи Коши (24), (25), функции x = rmx (t) и y = rmy (t) задаются формулами (26) и (27), а само приближенное решение, функция

ym(x) = rmy(rmx^_1^(x)), - формулой (28). В эти формулы входят величины xt (i = 0,1,* ,m),

Qyi 0 i, Qyi 1 i, Qyi 2 i (i = 1,2,K ,m) Qyr 0 i , Qyr 1 i , Qyr 2 i (i = 0,1 * ,m _ 1). Величины xi за-

даются произвольно, но так, чтобы (a = x0 < x1 < ...xm = b) . Иными словами, точки x{ образуют сетку на [a, b], вообще говоря, неравномерную. Величины Qyi j i и Qyr j i будем вычислять в цикле. На i-ом шаге цикла будут вычисляться Qyi 0 i, Qyi 1 i, Qyi 2 i. При этом

чт° значения Qyi 0 i_l, Qyi 1 i_b Qyi 2 i_1, Qyr 0 i_1, Qyr 1 i_Ъ Qyr 2 i_1 вычислены на предыдущем шаге этого цикла. Процедура вычисления величин Qyi 0 i, Qyi 1 i, Qyi 2 i следующая. Задаем m0i по формуле (45) и проводим итерации по формуле (46). По окончании итерационного процесса получаем Qyi 0 i = Qyi 0 i = mSmax i (здесь smax - номер последней итерации). Далее по формулам (38) и (39) получаем Qyi 1 i, Qyi 2 i, Qyr 1 i, Qyr 2 i и определяем тем самым приближенное решение ym (x) на [xi_1, x{ ]. До начала описанного цикла значение Qyr 0 0 определяется по формуле (40), а значения Qyr 1 0 и Qyr 2 0 - по формулам

(38) и (39). По завершении описанного цикла приближенное решение задачи Коши (24), (25) будет построено на всем [a, b] .

Описанная алгоритмическая схема дает нам принципиальную возможность получения алгоритма с переменным шагом сетки h = x{ _ x{_1, который будет подбираться таким образом, чтобы в результате было получено приближенное решение задачи Коши (24), (25) с заданной точностью e . Погрешность приближенного решения определим как разность между приближенным решением и точным:

em (x) = ym (x)_ y(x) . (49)

Легко видеть, что

ет(*Лх.гтх(,) = (ут (х)-у(*))*(,)= Гту (')-гу (<) ■ (50)

Поэтому

ет (х, ) = (ут (х, ) - у(х, )) = гту (г) - гу (г) . (51)

В дальнейшем нам понадобятся только приближенные значения е т (х,) и приближенные оценки тах |ет(х)|. Обозначим их 5,(, = 0,1,* ,т ) и А, (, = 1,2,* ,т). Их можно

хе[х,-1,х, ]

получить многими разными способами. Один из способов получения оценок 5г и Аг описан ниже.

Запишем алгоритм получения приближенного решения задачи Коши (24), (25) с заданной точностью е . Исходными данными являются функции и величины: /(х, у), Э^ у),

Эх

Э/ (х у)

—Э~—-, а, Ь , у . Результатами являются значения величин: т , х, (, = 0,1,* ,т), QуГ j,

( j = 0,1,2; г = 0,1,* , т -1), Qуl j, ( j = 0,1,2; , = 1,2,* , т ), которые, в свою очередь, позволяют найти приближенное решение по формулам (26) - (28).

В начале задается небольшое начальное значение величины т , равное т0 , на отрезке [а,Ь] строится начальная равномерная сетка точек х{ = а + И•, (, = 1,2,* ,т0) с постоянным крупным шагом И = (Ь - а) / т{). Затем, как описано выше, на этой сетке определяются неизвестные постоянные QуГ j,, Qуl j, и компоненты оценки погрешности 5,. Далее строится цикл по пх от 1 до Ых тах, в котором производится измельчение отрезков разбиения [х,-1, х{ ]

начальной сетки точек. Измельчение отрезков разбиения производится путём введения новых узлов сетки в серединах отдельных отрезков разбиения (при этом соответствующим образом меняется нумерация узлов сетки и увеличивается значение величины т ). Измельчение отрезков производится до тех пор, пока все значения Аг не станут меньше е или пока не

кончится цикл. В описываемом цикле используются два способа измельчения отрезков разбиения. При измельчении отрезков первым способом измельчаются только те отрезки [х,-1, х{ ] на которых значения оценок погрешности А, превышают е . При измельчении отрезков вторым способом, кроме отмеченных выше отрезков, измельчаются также все предыдущие отрезки (с меньшими номерами г ). После каждого измельчения отрезков и введения новых узлов вычисляются QуГ j,, Qуl j, и оценки погрешности 5, и А,. На каждом шаге

описываемого цикла, в начале производится N измельчений первым способом, а затем одно измельчение вторым способом.

Описанный метод решения задачи Коши мы будем называть модифицированным методом Эйлера с итерационным уточнением и переменным шагом. Нам осталось уточнить способ получения оценок погрешности 5г и Аг .

5. Оценка погрешности приближенного решения задачи Коши

Для получения оценки погрешности приближенного решения можно воспользоваться невязкой Ят (х )= С^т (х) - / (х, ут (х)), которая становится известной после вычисления

Сх

приближенного решения ут (х) = гту (гтх ( 1)(х)). Это обстоятельство позволяет получить

оценку погрешности приближенного решения, связав погрешность ет (х) = ут (х) - у (х) с невязкой. Из определения невязки (34) и требования (40) следует, что приближенное решение удовлетворяет задаче Коши, аналогичной задаче (24), (25):

йУт}Х) = I(Х ут (х)) + Ят (х), х е к Ь].

Сх

ут (а ) = у .

Отсюда несложно понять, что и погрешность приближенного решения также будет удовлетворять аналогичной задаче Коши:

Ст (х) = 1 (x, ут (х))- 1 (x, ут (х)-ет (х))+ Ят (х) = 1(x, ет (х)), х е к 4

(52)

(53)

дх

е т (а ) = У - У = 0 .

(54)

(55)

Поскольку приближенной решение ут (х) задается на разных отрезках [хг--ь х{ ] разными формулами, мы сведем задачу (54), (55) к цепочке задач Коши. Обозначим е ті (х) сужения

функции ет (х) на [хі-ь х{ ] и, при каждом значении і = 1,2,* , т, на каждом из этих отрезков будем последовательно решать задачи Коши для уравнения (54)

де ті(х) /Л Г і

—;-= їкх е ті(х Д х є [хі-1, хі].

ах

(56)

с начальными условиями вида:

0,

і = 1,

Єт‘(хм У "| Є ті—1 (х,‘—1), і > 1.

(57)

Для приближенного решения задач Коши (56), (57) можно использовать различные методы. Используем, например, схему Рунге-Кутта 4 порядка. Для этого введем равномерную сетку точек т,■ = хг-1 + х—-г-1 • j (j = 0,1,* ,Ы ) на [хг-1,х,]. Обозначим компоненты

N

приближенного сеточного решения задачи Коши (56), (57) Zj »еmiTj). Тогда для получения оценок 5 и А получается следующий алгоритм:

50 = 0,

z0 = 5,-ь А г = Ы.

К1 = /(т 1 ■ Zj )1

К 2 = / (т l + (хг - х, )/(2 N)

Кз = .^(т, + (х, - х,_1 )/(2 N)

К4 = 7 (т/ + (хг- хг-1)/N

+

( - х-1 ^/(2М)) (хг - х-1)2/(2М) (хі- х-1 )К з/М}

2і+1 = 2і +

(хі - х-1)

•[К1 + 2 • К2 + 2 • К3 + К4 ],

6 N

А,- := тах+^Аг} у = ^..^М-1, і = 1,2,* ,т.

N ’

При выборе значения N желательно добиваться, чтобы относительная погрешность сеточного решения ZN - етг )/1^| не превышала 30%. Оценить ее можно по правилу Рунге.

Но можно работать и при фиксированном значении N .

6. Алг°ри™ вычисления бу 0 1, 0у[ 1 1, 0у[ 2 I .

Для завершения алгоритма исправленного метода Эйлера с итерационным уточнением и переменным шагом нам осталось записать алгоритм вычисления величин бу/ о г, бу11 г,

бу/2 г по известным значениям буг о г-1, буг 1 г -1, буг 2 г -1. Исходными данными являются

функЦии и величины: / (x, у) , Э/ У) , ^ У) , Г/ (*) (1 Хг -1, Хг , буг 0 г -1,

0уг 1 г-1, буг 2г-1, А, 8 , Я, £ . Результатами являются значения величин: 0^0г, бу/1 г,

бу/ 2 г. Введем также дополнительный результат 1 = 1,2,...,т). Величина полагается равной 0, если на [хг-1, х, ] итерации завершаются при выполнении условия (48). Если итерации завершаются при ^ = £, то = 1. А если итерации завершаются при нарушении условия (47) или если ^ < £ , то = 2 . На основе исходных данных производятся следующие вычисления:

т := буг 0 г -1 + 0-уг 1 г -1 + ^ буг 2 г -1

если £ < 1

то бу/ 0 г := т ; бу/ 1 г := Ь-/ (хг, V )

б := ь2 э/ (хг, т). ь б э/ (хг, т)

бу/ 2 г := Ь---5----+ Ь-бу/1 г

эх эу

все

для ^ от 1 до £ с шагом 1 начало цикла по ^

*1, := Ь-/(х,,т); К12 := Ь2 -+Ь-Кп

эх эу

У1 := 32буг 0 г-1 + 10буг 1 г-1 + буг 2 г-1 ; У2 := -60буг 0 г-1 - 14буг 1 г-1 - буг 2 г-1

*13 := Ь-/

х-1 + х , (у1 + К12 - 10КП + 32т V64

2

V У

/ := (у2 + К12 -14*11 + 60т V32 -*13; ^ := ^ ^ если 1,1 < Я

то бу/ 0 г := т ; бу/ 1 г := К11; бу/ 2 г := К12 ; е := 0 конец цикла по ^

все

К21 = Ь-/(х,,т--Д); к22 := ь2-э/(х,эт~А)+Ь-К21 э/(х,э^~А)

эх эу

K23 := h'f

xi-l + xi 2

(v1 + K22 - 10K21 + 32(m -А))/64

V /

L2 := (v2 + K22 - 14K21 + 60(m -А))/32 -K23 ; L2 := L2 ' L

K31 := h'f (xi,I + А); K32 := h

2' э/(xi, m+А)+h ' K f(xi, m+А)

Эx

31

ЭУ

K33 := h'f

xi-1 + xi 2

, (v1 + K32 - 10K31 + 32(m + А))/64

V

L3 := (v2 + K32 - 14K31 + 60(m + А))/32 -K33; L3 := L3 'L3 L4 : = L3 - 2 L1 + L2 ; L5 : = L3 - L2 ; если L4 < <5

то Qyi 01 := m; Qyi 11 := Kii; Qyi 21 := Ki2; e. := 2 конец цикла по

все

если s = S

то Qyi01 := m; Qyi 11 := h'f(xi,m)

Qyi 2/:= h+ h-Qy, 1 і Э/(Э^ ; e. := ї

ЭУ

все

конец цикла по s

Подбор параметров и некоторые численные результаты

Описанный численный метод был реализован в виде программы на языке Visual Basic и исследован на трёх модельных задачах Коши.

Первая модельная задача Коши

ddX = У, X є [0,8],

dx

У(0)= I

(59)

(60)

имеет известное точное решение y(x) = ex, экспоненциально растущее на [0,8] и принимающее большие значения вместе со своими производными в окрестности точки 8.

Вторая модельная задача

dy = -100' У +100, x є [o,l] dx

y(o) = 2,

(б1)

(б2)

имеет известное точное решение у(х) = 1 + е 100 - х. Это решение имеет большие по модулю значения производных на левом конце отрезка интегрирования (при х = 0 ), где ставится начальное условие. Поэтому данную задачу Коши можно назвать жесткой.

s

Третья модельная задача

^ = -2 • х • е-у, х є [- 0.9 , 0.9], дх

у (- 0.9 ) = 1п (0.19), (64)

имеет известное точное решение у (х ) = 1п(1 - х2). Это решение имеет большие по модулю значения производные на концах отрезка интегрирования. Поэтому и третью модельную задачу можно назвать жесткой. Рассмотрим некоторые результаты этого исследования.

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

[5]: А = 10-6, 8 = 10 -11

-24

21

£ = 5 . Во всех случаях, когда заданная точность

е > 10 , можно рекомендовать следующие значения параметров, регулирующих алгоритм

дробления шага сетки и получения оценок погрешности: Ых Мх = 6 = 1, N = 2. Во всех

таких случаях заданная точность достигается быстро. Если е < 10-11, то время счета начинает быть заметным, а заданная точность может и не достигаться. В таких случаях иногда удаётся добиться увеличения заданной точности за счет увеличения значений Nx Мх и N . Начальное число отрезков разбиения т0 > 1. Но возможны случаи, когда при маленьких значениях т0 возникает переполнение разрядной сетки. Причиной его являются некоторые очень большие значения квадрата невязки, которые реально возникают при решении некоторых задач Коши, когда значения т0 малы. В таких случаях надо увеличить значение т0. Во всех расчетах, описанных ниже, основные параметры алгоритма имели следующие значения:

6 Ns = 1, N = 2, т0 = 8.

А = 10-6, 8 = 10-24, 1 = 10-21, £ = 5, Nxmax

На рисунках 1 - 3 изображены график погрешности е т (х) и точки графика приближенной оценки погрешности (хі,8І), і = 0,1,...,т для первой, второй и третьей модельной задач, соответственно.

I --11 1

10

а-;;:--" * С '■ Ч I ч ; - г-. ■■ Л і і ? Г~*‘ >' Л , м ■ *

V - - 1' 1,5 И.З ’ /,? а) б) Рис. 1. График погрешности ет(х) и точки (хг-,81), 1 = 0,1,...,т для первой модельной

7 задачи при е = 10 ; значение т = 256; а) х є [1.5 , 2.5]; б) х є [6.5 , 7.5].

Первая модельная задача имеет экспоненциально растущее точное решение, имеющее большую по модулю производную в окрестности правого конца интегрирования. На рисунке

1 можно заметить, что с ростом значений аргумента шаг сетки точек х{ уменьшается.

Модуль производной точного решения второй модельной задачи имеет большие значения в окрестности левого конца отрезка интегрирования. На рисунке 2 можно увидеть, что именно в окрестности левого конца интегрирования сетка точек х{ - густеет.

Производная точного решения третьей модельной задачи велика по модулю как в окрестности левого конца, так и в окрестности правого конца отрезка интегрирования. На рисунке 3 можно увидеть, что именно в окрестности левого и правого концов отрезка интегрирования сетка точек х{ - самая густая.

10'1(1

-7 Ю L J

Г I

в m's

I '

■ " о л' о::

а) б)

Рис. 2. График погрешности ет(х) и точки (хг-,81), I = 0,1,...,т для второй модельной

—7

задачи при е = 10 ; значение т = 52;

а) хе [0 , 0.2]; б) хе [6.5 , 7.5].

Во всех рассмотренных случаях сгущение сетки приводит к локальному уменьшению погрешности.

Нетрудно заметить, что точки (хг-, 8t) не всегда лежат на графике погрешности. Добиться этого можно увеличив значение переменной N, но это приведет к дополнительным вычислительным затратам.

Заключение

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

The new numerical method of the solution of Cauchy problems for the ordinary differential equations is proposed.

The key words: numerical method, the Cauchy problem, ordinary differential equations, variable step.

Список литературы

1. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. [Текст] / Н.С.Бахвалов, Н.П. Жидков, Г.М. Кобельков. М.: Наука, 1987.

2. Березин И.С., Жидков Н.П. Методы вычислений. Том 2. [Текст] / И.С.Березин, Н.П. Жидков. М.: Наука, 1960.

3. Вержбицкий В.М. Основы численных методов: Учебник для вузов. [Текст] /

B.М. Вержбицкий. М.: Высш. шк., 2002.

4. Трубников С.В. Кинематические кривые (текст) / С.В.Трубников // Вестник Брянского государственного университета. №4 (2004) Естественные и точные науки. Брянск: РИО БГУ, 2004. С. 122-128.

5. Трубников С.В. О новом подходе к построению численных методов решения одномерных задач Коши на основе эрмитовой кусочно-многочленной интерполяции (текст) /

C.В.Трубников // Вестник Брянского государственного университета. №4 (2006) Естественные и точные науки. Брянск: РИО БГУ, 2006. С. 199-217.

6. Васильев Ф.П. Численные методы решения экстремальных задач. [Текст] / Ф.П.Васильев. М.: Наука, 1980 Васильев Ф.П. Численные методы решения экстремальных задач. [Текст] / Ф.П. Васильев. М.: Наука, 1980.

Об авторах

С.В. Трубников - канд. физ-мат. наук, доц., Брянский государственный университет им. академика И.Г. Петровского, ЬгуапБк§и@ mail.ru.