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

Глава 6 ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ КОШИ ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ Дифференциальные уравнения и системы таких уравнений являются математическими моделями множества научно­ технических задач. В ряде случаев систему можно решить в конечном виде, т. е. получить решение в виде элементарных функций или системы первых интегралов. Большое количество таких уравнений можно найти в [14, 15, 18]. Однако большин­ ство таких задач в конечном виде не решается, а если такое ре­ шение получить можно, то это часто сопряжено с трудностями, которые не оправдывают затрат труда на его получение. В явном виде представляется решение уравнения у"+у=О. Его общим решением, очевидно, яв-ляется функция у= =Acosx + Bsinx. Мы имеем явное выражение для решения. Решением уравнения , х+у у=-­ х-у будет общий интеграл 2 arctg Jf_- ln( х 2 + у 2 ) =С. Выразить из него х у в явном виде невозможно. Решение уравнения у'= х 2 +у2 в виде общего интеграла за­ писать не удается. Таким образом, в большинстве случаев задача в конечном виде не решается, а если такое решение получить можно, то это часто сопряжено с трудностями, которые не оправдывают затрат на его получение. В настоящее время разработаны приближенные методы решения таких задач. Одна из основных математических моделей для дифферен­ циальных уравнений - Постановка задачи задача Коши. 1. Требуется на отрезке [а,Ь] решить уравнение у' =f(x,y) 146 (6.1) при начальных условиях (6.2) Будем считать, что решение задачи 1 существует и един­ ственно. Для этого достаточно потребовать выполнения условий теоремы 6.1 [37]. Теорема 6.1 (теорема Коши). Если в полосе Т = ={а :5 х :5 Ь,-оо <у< оо} функция f(x,y) удовлетворяет следую­ щим условиям: 1) она непрерывна по совокупности переменных; 2) удовлетворяет условию Липшица по у: с положительным L для любых точек из Т, то на отрезке [а, Ь] существует единственное решение задачи 1. Условие Липшица проверить достаточно трудно. Поэтому его обычно заменяют на более сильное: потребуем, чтобы в полосе Т производная fy была непрерывной и ограниченной, ltyl(x ) L kl k=O n-1 = Lу (k)(x ) k=O kl (x-x 0 )k + y(~) п! (х-х 0 )п = (x-xo)k +О[(х-х0 )п]. (6.3) Здесь точка ~ лежит между точками х и х0 и зависит от их положения. Для того чтобы представить решение задачи (6.1) в таком виде, нужно найти значения производных y(x0 ). Будем обозначать значения производных функции f(x,y) в начальной точке (х 0 ,у 0 ) индексом О справа. Теорема 6.2. Если f(x,y) имеет непрерывные вторые произ­ водные в окрестности точки (х 0 ,у 0 ), то справедлива формула у(х) =Уо + fo(X-Xo)+~(fxo + fyofo)(X-Xo) 2 + + 3 +О[(Х-Хо) 4 ]. (6.4) Величина у(х 0 ) нам известна, у(х 0 )=у 0 • Подставим в уравнение (6.1)-(6.2). (6.1) точное решение задачи Коши При этом оно превращается в тождество у'(х) =f(x,y(x)), которое мы имеем право дифференцировать. 148 (6.5) Подставим в (6.1) х = х 0 , у= у0 : у'(хо) = f(xo.Yo) = fo· Продифференцируем тождество у"= Подставим в (6. 7) у"(хо) = (6.5) (6.6) по х: fx(x,y(x)) + fy(x,y(x))y'. х = х 0 , у= Уо• у'(хо) = (6.7) fo: fx(Xo,y(xo ))+ fy(Xo ,у(хо))у'(хо) = = fхо + fuofo • Продифференцируем у"'(х) = (6.5) (6.8) еще раз: fxx (х,у(х)) + 2fxy (х,у(х))у'(х) + +fyy(x,y(x))[y'(x)] 2 + fy(x,y(x))y"(x) = = fxx (х,у(х)) + 2fxy (х,у(х))у'(х) + fuu (х,у(х))[у'(х)] 2 + +fy (x,y(x))[fx (х,у(х)) + fu (х,у(х))у']. Подставим в последнее равенство х 0 : у"'(хо) = fxxo +2fxuofo + fuuoM + fxofyo + fu2ofo· (6.9) Таким образом, получена требуемая формула: у(х) = Уо + fo(X-Xo)+l. k2(h) = kз(h) = hf( Хо +~,у 0 + ki ~h)). hf( Хо+~ ,у 0 + k2~h)). (6.22) k 4 (h)=hf(x 0 +h,y 0 +k3 (h)). При r= 5 формулы пятого порядка построить невозможно. Формулы более высокого порядка громоздки и используются в специальных случаях. К достоинствам метода Рунге- Кутты относят следующие: в 1) формулы Рунге- Кутты имеют высокую точность; 2) формулы являются явными, т. е. для получения решения очередной точке вычисления ведутся по формулам (6.13); 3) метод допускает расчет с переменным шагом, т. е. на участ­ ках быстрого изменения решения требуется применять мелкий шаг, а на пологих участках - большой; 4) метод «самостартующий•, т. е. для расчета очередной точки Yk+l необходимо знать только Yk и h. К недостаткам метода следует отнести то, что на одном шаге ,приходится несколько раз вычислять промежуточные значения функции f(x,y), которые в дальнейшем не используются. 6.1.3. Метод Адамса Как упоминалось ранее, существенным недостатком метода Рунге- Кутты являете.я то, что на одном шаге численного ин­ тегрирования приходится многократно вычислять правую часть уравнения (6.1). Если она имеет сложный вид, то это сопряже­ но со значительной вычислительной работой. Для устранения данного недостатка были разработаны разностные методы, тре­ бующие лишь одно вычисление правой части на шаге, в част­ ности метод Адамса. Предположим, что нам удалось вычислить приближенное решение в точках х0 =а, х 1 Yi =y(xi) - 154 = х0 + h, ... , Хт =Xm-l +h. Обозначим точное решение задачи (6.1) в узлах. Обозначим также приближенное решение в этих же точках, как У;= У(х;) и fi =f(xi•Yi), Fi =f(xi,Yi). Подставим в уравнение (6.1) точное решение и проинтегри­ руем его по отрезку [хт-;•Хт+1]: Xm+l f Xm+l y'dx = Xm-j f f(x,y(x))dx, Xm-j Ут+1 -Ут-; (6.23) Xm+l f = f(x,y(x))dx. Xm-j Поскольку точное решение нам неизвестно, заменим правую часть (6.23) по какой-либо приближенной формуле, например формуле левых прямоугольников, подставив в нее полученные приближенные значения. В результате получим формулу i (6.24) Ym+1-Ym-j =hLFm-k• k=O Так как величины, стоящие в правой части, нам известны, можно вычислить решение в точке часть (6.23) Xm+l· Заменим теперь правую по формуле правых прямоугольников и получим Ут+1 - Ym-j i =h L (6.25) Fm-k. k=-1 Эта формула связывает значения Ym+l и Fm+l• которые пока неизвестны. Следовательно, мы получили уравнение для опреде­ ления этого решения. ·Рассмотрим формулы более общие, чем k (6.24) и (6.25): k (6.26) La.iYm+i =hL/3iFm+i• i=O где a.i, /3i - i=O коэффициенты, не зависящие от У и h. Получим уравнения для определения коэффициентов фор­ мулы (6.26). Будем считать, что решение задачи (6.1) имеет р непрерывных производных на отрезке [а,Ь]. Подставим в (6.26) это точное решение и разложим полученное выражение в точке х,,,, по степеням h, учтя тот факт, что f; =у;: 155 +~hi yУт = =h[P21m+2 +Р11т+1 +(Р1 +3Р2 -2)fт]• (6.29) В случае .явной формулы (Р 2 =О) получаем семейство формул, зависящих от одного параметра р: Ут+2 +(2Р-4)Ут+1 +{3-2Р)Ут =h[Pfт+1 +(Р-2)fт]; (6.30) в случае Р=-1 получим формулу Ут+2 -6Ут+1 +5ут =-h(fт+1 +3fт); при (6.31) P=l (6.32) 157 при Р=3/2 (6.33) Неявное уравнение можно получить, полагая, например, Р2 =1, Р1 =-1. При этом а. 0 =1, а. 1 =-2, Ро =0 и формула при­ нимает вид Ут+2 -2Ут+1 Положим в формуле + Ут = h(f m+2-f т+д• (6.29) т = О: (6.34) У2 +(2Р1 +4Р2 -4)У1 +(3-2Р1 -4Р2>Уо = (6.35) =h[P2f2 +P1ft +<Р1 +3Р2 -2)fo]· Нетрудно видеть, что для начала расчета по формуле нам не хватает одного значения - (6.35) У 1 • Способ его получения обсудим далее, а теперь, считая, что оно известно, протестируем полученные формулы на задаче { у'= у, ХЕ [0,1], y(O)=l, решение которой известно у= ех. Точное (у - = 2, 71828) приближенное (У) решения будем сравнивать в точке х = Результаты сведем в табл. и 1. 6.2. Рассмотрим теперь пример применения неявной формулы, рассчитав решение задачи Коши с точностью Е +у2, { у'=х2 у(О) = l в точке х = 0,1. Воспользуемся формулой = 0,001: (6.36) (6.28): Таблица Номер формулы (6.31) (6.32) (6.33) 158 h Y(l) 8 = IY(l) - y(l)I 0,1 100508,9 100 506 0,01 7,38168. 1066 7 ,38168. 1066 0,1 2,676154 0,042128 0,01 2,713366 0,004915 0,1 2,708814 0,009468 0,01 2, 71817 0,000111 6.2 У1 -уо = 0·:} (f1 + fo). У1 (6.37) -1=0,05(0,0l+yf +1), У1 -1=0,05(1,0l+yf). Подставим в правую часть уравнения (6.37) yi > = 1 и получим yi > = 1+0,05(1,01+1) = 1,1005. Подставим теперь в правую часть (6.37) yi >. Имеем 1 1 У1 2 ) = 1+0,05(1,0l+(l,1005) 2 ) = 1,111055. Разность приближений У1 2 ) -ур> =1,111055-1,1005=0,01055. Поскольку разность велика, рассчитываем очередное при­ ближение: УiЗ) = 1+0,05(1,01+(1,111055) 2 ) = 1,11222; У1З) _у1 2 ) =1,11222-1,111055=0,001165. Четвертое приближение дает требуемую точность yf 4) = 1+0,05(1,01+(1,11222) 2 ) = 1,11235; yf 4 >-yi 3 >=1,11235-1,1122=0,00015. Принимаем полученное значение в качестве у 1 и переходим к вычислению значения у 2 • Таким образом, мы видели особенности применения формулы Адамса. Разберемся теперь с результатами счета по формулам (6.31)-(6.33). 6.2. Понятие устойчивости Определение 6.6. Будем говорить, что уравнение (6.26) устойчиво по начальным данным, если существует постоянная М, не зависящая от h, если для любых начальных данных Уо,У1, ." ,Yk-1 выполняется неравенство IYml~M 0:5.j:5,k-1 max JYjJ• m=k,k+1, ... (6.38) f Понятие устойчивости означает, что решение уравне­ • ния (6.26) равномерно ограничено. 159 Сформулируем условие устойчивости. Определение 6. 7. Многочлен k p(z)= ,Laizi (6.39) j=O будем называть характеристическим мн.огочяен.ом для урав­ нения (6.26), а уравнение p(z) =0 - (6.40) характеристическим уравн.ен.ием для уравнения Теорема (6.26). 6.3. Для того чтобы решение уравнения (6.26) было устойчиво по начальным данным, необходимо и достаточно, чтобы уравнение имело корни, лежащие внутри еди­ (6.40) ничного круга на комплексной плоскости. Корни, лежащие на единичной окружности, должны быть простыми. f • Из первого уравнения (6.27) следует, что 1 всегда явля­ ется корнем характеристического уравнения (6.40). Рассмотрим с этой точки зрения уравнения (6.31) - (6.34). (6.31) будет Характеристическим уравнением для уравнения уравнение z 2 -6z+5=0. (6.41) (6.41) являются z = 1 и z = 5. Следова­ (6.31) неустойчиво. Результаты его приме­ нения видны в табл. 6.2. Для уравнения (6.32) корнями характеристического уравне­ ния будут z 1,2 = 1. Уравнение (6.32) также неустойчиво. Несмо­ Корнями уравнения тельно, уравнение тря на то что результаты достаточно приемлемы, может найтись уравнение, для которого неустойчивость проявится. Для уравнения (6.33) корнями характеристического уравне­ = 1 и z = О. Уравнение (6.33) устойчиво. ния будут z Уравнение (6.34), как и (6.32), неустойчиво. 6.3. Задача Коши ДJIЯ сисrемы уравнений первого порядка Постановка задачи 2. Требуется на отрезке [а,Ь] решить си­ стему п уравнений y'=f(x,y) (6.42) y(a)='ifo, (6.43) при начальных условиях где у и 160 f(x,y) - п-мерные векторы. Будем считать, что решение задачи существует и един­ 2 ственно. Для этого достаточно потребовать выполнения условий теоремы 6.4. Теорем~ 6.4 (теорема Коши). Если в слое Т ={а::;; х::;; Ь, lliJll < оо} функция f (х,У): 1) непрерывна по совокупности переменных; 2) каждая компонента вектора 7 удовлетворяет условию Лип­ шица по у: n it; 2 . j=l В дальнейшем будем считать, что правая часть уравнения (6.42) в области Т имеет нужное по ходу изложения количество производных. 6.3.1. Сведение уравнения вьи:окоzо порядка к системе уравнений первого порядка Пусть дана задача Коши для уравнения порядка п, разре­ шенного относительно старшей производной y(n) = f(x,y,y', ... ,y - бук­ вой Е. Тогда применение метода Хемминга можно символиче­ ски записать как РС. На практике между стадиями Р и С часто делают оценку причем итерации осуществляют несколько раз. В этом случае процесс получения уточненного решения можно символически записать в виде Р(ЕС)п. Достоинствами методов прогноза-коррекции являются сле­ f, дующие: 1) разность между прогнозированным и скорректированным значениями дает оценку ошибки, сделанной на шаге, и может быть использована для контроля величины шага интегрирова­ ния; 2) на одном шаге интегрирования правые части уравнения вычисляются только два раза по сравнению с четырьмя вычис­ лениями для стандартного метода Рунге- Кутты. 163 К недостаткам методов относятся сложность их программи­ рования и необходимость вычисления некоторого количества начальных точек. Пример 6.5. Рассмотрим применение формул (6.50) к задаче взяв в качестве шага интегрирования h = 0,1 и вычислив (6.36), по формуле Эйлера значение У1 =1+0,1(0 2 +1 2 )=1,1. Формулы примут следующий вид ур> = У1 +0,05(3F1 -F0 ) = У1 +0,05(3(xf + У12 )-(х~ + Уа2)), YJ 2 > = У1 +О,05((х~ +(YJ1» 2)+(xf + Yi2}). При наших данных система примет вид y~l) = 1, 1+о,05(3(0, 01+(1,1) 2)-1), У~2 ) = 1,1 +0,05((0,01 +(У~1 » 2 > +(0,01 +(1,1) 2 )). Шаг коррекции дает значение y~l) = 1,233, у~2 > = 1,1 +О,05(0,02+(1,233) 2 +(1,1) 2 )=1,2375. Разность спрогнозированного и скорректированного значе­ ний 1,2375-1, 233 =О, 0045 слишком велика, поэтому применим итерационное уточнение. Следующая итерация имеет вид YJ 3) =1,1+0,05(0,02+(1,2375) 2 +(1,1) 2 )=1,238. Разность второй и первой итераций 1,238-1,2375=0,0005 достаточно мала, поэтому принимаем в качестве У2 = YJ 3> = 1,238 и переходим к вычислению значения У3 • 6.5. Выбор mara и.нтеrрирования Выбор шага интегрирования - сложная проблема. Сначала пред­ положим, что все вычисления проводятся точно. В этом случае рас­ смотренные методы имеют погрешность, т. е. норма разности между дифференциальным уравнением у'= f (х, у) и формулой численного интегрирования L(y) имеет вид llY' - f(x,y)-L(y)ll = O(hP ). Если точность решения _~авнения равна&, то в качестве на­ чального шага возьмем h = f./&. Для дальнейшего уточнения шага можно использовать правило Рунге. Вычислим решение в точке Х; сначала с шагом h, У?>, затем с шагом h/2, Y/hl 2>. Из (6.29) следует, что главный член погрешности O(hP) имеет вид у(х;)-У; =KhP, где константа К нам не известна. 164 Тогда получим y(xi)-Y/h> =KhP, y(xi)-Y/hl 2>=K(h/2)P. Исключив из системы неизвестное y(xi), можно оценить по­ стоянную К: y_(h/2) - y_(h) к- i i - (h/2)P(2P -1) • Значит, критерием nрекращения деления шага может быть выполнение условия y_(h/2) _ y_(h) 1 l i 2Р -1 (6.51) <Е. i Итак, схема вычислений такова: 1) вычисляется начальный шаг h = efi; 2) вычисляются решения с шагом h и h/2; 3) проверяется выполнение условия (6.51); 4) в случае выполнения условия (6.51) аргументу х задают приращение h; 5) в противном случае шаг h делят пополам и повторяют п. 2. Подобная схема решения задачи заложена в большинство современных стандартных программ. При этом в процедуру решения заложено автоматическое прекращение решения за­ дачи, если после некоторого количества делений шага точность решения не достигается. 6.6. Оценка поrрешносrи вычислений Теперь рассмотрим случай, более приближенный к реаль­ ности, т. е. учтем, что существующие ЭВМ имеют конечную разрядную сетку. Покажем, что в этом случае неограниченно уменьшать шаг нельзя. Обозначим через у(х,х 0 ,у 0 ) решение задачи (6.1) с начальными условиями у(х 0 ) у 0 , вычисленное в = точке х, через У решение в точке xk результат вычисления. Тогда приближенное можно рассматривать как точное решение с начальным условием Ek y(x,xk,Yk). Вычислим погрешность = y(xk,x 0,y0)-Yk =y(xk,x 0,yo)-y(xk,Xk, Yk). В силу условий, наложенных на уравнение (6.1), (6.52) оно должно иметь единственное решение. Потребуем также выполнения не­ прерывной зависимости решения от начальных данных. В при­ нятых обозначениях это равносильно условию 165 lду{~,11>1~м. Преобразуем разность (6.53) (6.52) Ek =: y(xk,Xo,Yo)-y(xk,Xk, Yk) = =y(xk ,хо ·Уо )- y(xk ,хо' Уо) + y(xk ,хо' Уо )-y(xk ,х1 'У1) + k-1 = y(xk,xo.Yo)-y(xk,Xo, Уо)+ L[y(xk,xi, Yi )-y(xk,xi+l, Yi+1)]. j=O Так как то k-1 +L[y(xk,xi+1 •Y(Xi+1 •xi, Yj))-y(xk ,xi+1 • Yi+1)]. j=O По теореме Лагранжа имеем _ Ek - 0 ,у0 )( -У. )+~дy(xk,xi,Yi)[ ( . . У·)-У. ] ду(хk,х ~ Уо О ~ ~ У Х1+1•Х1, J 1+1 • }=0 VI( где у0 , Yi - VI( промежуточные точки. Оценим модуль погрешности -У.1+ ~ дy(xk,xi,Yi) 1У(Х1+1•Х1, . . У·)-У· 1< 1Ek l -