Численные методы. Соснин (2005) (1160462), страница 15
Текст из файла (страница 15)
Будем предполагать, что функция f (t, u) Липшиц-непрерывна по второму аргументу, то есть ∀ u1 , u2 |f (t, u1 ) − f (t, u2 )| 6 L|u1 − u2 |, где L — некоторая константа.Рассмотрим функцию погрешности: zn = yn − un . Выразим из нее yn = zn + un и подставим влевую часть (5.2). Тогда получим:mmXzn+1 − znun+1 − un X=−+σi Ki (u) +σi (Ki (y) − Ki (u)).ττi=1i=1mmXun+1 − un X+σi Ki (u) = ψn , и обозначим ψn =σi (Ki (y) − Ki (u)).τi=1i=1Оценим для разных i выражение |Ki (y) − Ki (u)| :Заметим, что −i=1:i=2:|K1 (y) − K1 (u)| = |f (tn , yn ) − f (tn , un )| 6 L|yn − un | = L|zn |.|K2 (y) − K2 (u)| = |f (tn + a2 τ, yn + b21 τ K1 (y)) − f (tn + a2 τ, un + b21 τ K1 (u))| 66 L|yn − un + b21 τ (K1 (y) − K1 (u))| 6 L(|zn | + τ |b21 | · |K1 (y) − K1 (u)|) 66 L(|zn | + τ |b21 |L|zn |) 6 {обозначим b = max bij } 6 L|zn |(1 + τ bL).i=2, nj=1, m−176Глава 5.
ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧi=3: Если тоже самое проделать для i = 3, то получим такую оценку:|K3 (y) − K3 (u)| 6 L|zn |(1 + τ bL)2 .Теперь перейдем к общей оценке. Докажем, что|Kl (y) − Kl (u)| 6 L|zn |(1 + τ bL)l−1 ,l = 1, m.Пусть эта оценка верна для некоторого i :|Ki (y) − Ki (u)| 6 L|zn |(1 + τ bL)i−1 ,докажем ее для i + 1 : iXb(i+1)j Kj (y) −|Ki+1 (y) − Ki+1 (u)| = f tn + ai+1 τ, yn + τj=1 X iiX −f tn + ai+1 τ, un + τb(i+1)j Kj (u) 6 L yn − un + τb(i+1)j (Kj (y) − Kj (u)) 6 j=1j=16 {воспользуемся неравенством bij 6 b и подсчитаем сумму геометрической прогрессии} 6iiXXj−1(1 + τ bL)=|Kj (y) − Kj (u)| 6 L |zn | + τ bL|zn |6 L |yn − un | + τ bj=1j=1bL)i= L|zn | 1 + τ bL 1−(1+τ= L|zn |(1 + τ bL)i .1−(1+τ bL)Таким образом, оценка (5.3) действительно имеет место.
Получим теперь оценку на |ψn | :|ψn | 6mX|σi | · |Ki (y) − Ki (u)| 6 {обозначим σ = max |σi |} 6i=1, mi=16σmXL|zn |(1 + τ bL)i−1 = σL|zn |(1 + τ bL)m−1 m.i=1Оценим два последних множителя:m(1 + τ bL)m−1 6 {(1 + y)α 6 eαy } 6 meτ bL(m−1) 6 {τ 6 T} 6 meT bL(m−1) — обозначим за α,тогда |ψn | 6 |zn |σLα.Откуда, учитывая, чтоzn+1 − zn= ψn + ψn , получаем оценку на погрешность:τ|zn+1 | 6 |zn | + τ |ψn | + τ |ψn | 6 |zn | + τ |ψn | + τ |zn |σLα == |zn |(1 + τ σLα) + τ |ψn | 6 |zn−1 |(1 + τ σLα)2 + τ |ψn−1 |(1 + τ σLα) + τ |ψn |.Применив эту же операцию n − 1 раз, получим:|zn+1 | 6 |z0 |(1 + τ σLα)n + τnX|ψj |(1 + τ σLα)n−j .j=0Если мы обозначим ψ = max |ψj | и учтем, что z0 = y0 − u0 = 0, то получим:j=0, n|zn+1 | 6 ψτnXj=0(1 + τ σLα)j 6 ψτ max (1 + τ σLα)j (n + 1).j=0, n(5.3)5.2.
МЕТОДЫ РУНГЕ-КУТТА ВТОРОГО ПОРЯДКА АППРОКСИМАЦИИ77В силу того, что τ n = T,|zn+1 | 6 ψT eT σLα ,откуда следует, что если наша схема аппроксимирует на всей сетке ОДУ, то есть |ψi | → 0, то имеетместо сходимость (|zn | → 0 ). Кроме того, если наша разностная схема аппроксимирует исходное ОДУс p-м порядком аппроксимации (ψn = O(τ p )), то погрешность имеет соответственно p-й порядок:|zn | = O(τ p ).Таким образом, теорема полностью доказана.Теперь свяжем требование на порядок аппроксимации с количеством промежуточных точек, в которых требуется вычислять значение f (t, y).5.2Методы Рунге-Кутта второго порядка аппроксимацииРассмотрим семейство методов Рунге-Кутта при m = 2. Схема для вычисления приближенного значения ( yn+1 ) будет выглядеть так:yn+1 − ynτ= σ1 K1 (y) + σ2 K2 (y);K1 (y) = f (tn , yn );K2 (y) = f (tn + a2 τ, yn + τ b21 K1 (y)).Найдем, что является достаточным условием для достижения 2-го порядка точности.
Для этогорассмотрим выражение для невязки:ψn = −un+1 − un+ σ1 K1 (u) + σ2 K2 (u),τ(5.4)и потребуем, чтобы ψn = O(τ 2 ).Для начала распишем дробь в правой части (5.4), применив формулу Тейлора:un+1 − un11τ2τ= (u(tn + τ ) − un ) = (un + u0n τ + u00n+ O(τ 3 ) − un ) = u0n + u00n + O(τ 2 ).τττ22Согласно постановке, K1 (u) = f (tn , yn ).
Обозначим это число за fn . Используя это обозначение,разложим выражение для K2 (u) по формуле Тейлора для функции двух переменных:K2 (u) = f (tn + a2 τ, un + τ b21 fn ) = fn + ft0 (tn , yn )a2 τ + fu0 (tn , un )b21 fn τ + O(τ 2 ).Подставив это выражение в (5.4), получим такое выражение для невязки:ψn = −u0n − u00nτ+ σ1 fn + σ2 (fn + ft0 (tn , yn )a2 τ + fu0 (tn , un )b21 fn τ ) + O(τ 2 ).2Так как u — точное решение, то для него справедливы формулы: 0u = f (t, u);u00 = ft0 + fu0 u0 = ft0 + fu0 f.С их использованием выражение для невязки перепишется так:ψn = fn (−1 + σ1 + σ2 ) + ft0 (tn , yn )(−ττ+ σ2 a2 τ ) + fu0 (tn , yn )fn (− + σ2 b21 τ ) + O(τ 2 ).2278Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧНетрудно подобрать такие σ1 , σ2 , a2 , b21 , чтобы первые три слагаемых обратились в ноль — это идаст требуемую оценку для невязки.
Коэффициенты σ1 , σ2 , a2 , b21 должны быть такими, чтоσ1 + σ2 = 1;1σ2 a2 = ;2 σ b = 1.2 212(5.5)Из последних двух уравнений следует, что a2 и b21 равны друг другу. Обозначим это число за a.Теперь обозначим σ2 = σ, и тогда из первого уравнения будет следовать, что σ1 = 1−σ.
Это позволяетпереписать систему (5.5) так: σ1 = 1 − σ; σa = 1 .2Таким образом, методы с расчетной схемой следующего вида:yn+1 − yn = τ ((1 − σ)f (tn , yn ) + σf (tn + aτ, yn + af (tn , yn )τ ))1имеют 2-й порядок аппроксимации, если выполняется условие σa = . Согласно доказанной теореме,2построенное решение будет иметь 2-й порядок точности.1Можно привести примеры таких методов.
При σ = 1, a = получается такая схема:2yn+1 − yn = τ f (tn +А при σ =ττ f (tn , yn ), yn +).221, a = 1 — вот такая:2yn+1 − yn =τ[f (tn , yn ) + f (tn+1 , yn + τ f (tn , yn )] .2Несмотря на бесконечное количество схем заданного порядка точности, нам может не подойти ниодна. Это связано с тем, что методы Рунге-Кутта не являются устойчивыми, и при их использованиинакапливается машинная погрешность, в конце вычислений сравнимая с полученными величинами.Подробнее о вычислительной устойчивости речь пойдет в следующих разделах.Методы Рунге-Кутта четвертого порядка точностиПриведем без вывода расчетную схему 4-го порядка точности в методе Рунге-Кутта с параметромm=4:1yn+1 − yn=(K1 (y) + 2K2 (y) + 2K3 (y) + K4 (y));τ6K1 (y) = f (tn , yn );ττK2 (y) = f (tn + , yn + K1 (y));22ττK3 (y) = f (tn + , yn + K2 (y));22K4 (y) = f (tn + τ, yn + τ K3 (y)).Все методы Рунге-Кутта требуют возможность вычислять значение функции в произвольной точке.Теперь рассмотрим методы, в которых этого делать не надо.5.3.
ОПИСАНИЕ МНОГОШАГОВЫХ МЕТОДОВ5.379Описание многошаговых методовОпределение. Линейным m-шаговым методом называется метод с расчетной схемой следующего вида:a0 yn + a1 yn−1 + a2 yn−2 + . . . + am yn−m= b0 fn + b1 fn−1 + . . . + bm fn−m .(5.6)τгде ai , bj — параметры метода, а yn−k и fn−i означают следующее:yn−k = y(tn−k );fn−i = f (tn−i , yn−i ).Таким образом, для реализации m-шагового метода на первом шаге требуется знать значенияy0 , y1 , . . . , ym−1 . Значение y0 можно взять равным u0 — начальному условию, а вот для вычисленияy1 , .
. . , ym−1 применяют методы Рунге-Кутта соответствующего порядка точности.Заметим также, что в методе используются только табличные данные о f (x) — то есть уметьвычислять функцию f на промежуточных точках в общем случае не требуется, а может понадобитьсятолько для получения y1 , . . . , ym−1 .Если в схеме (5.6) коэффициент b0 равен нулю, то в правой части fn не присутствует, и соответствующий метод называется явным (по тем же причинам, что и раньше). Если же b0 6= 0, то методназывается неявным (как ни странно, тоже по тем же причинам, что и раньше); возникает нелинейноеуравнение относительно yn , которое в общем случае решается методом Ньютона.Очевидно, что метод не изменится, если выражение (5.6) домножить на какую-нибудь ненулевуюmXконстанту. Поэтому устраним неоднозначность, введя условие нормировки:bi = 1.
Покажем, что вi=0этом случае правая часть уравнения (5.6) будет аппроксимировать правую часть дифференциальногоуравнения исходной задачи:f (tn , un ) −mXbi f (tn − iτ, u(tn − iτ )) = {разлагая слагаемые в сумме в ряд Тейлора} =i=0= f (tn , un ) −mXbi [f (tn , un ) + O(τ )]gg = fn (1 −i=0mXbi ) + O(τ ) = O(τ )i=0— это и означает аппроксимацию.Теперь выведем достаточные условия для достижения k-го порядка аппроксимации исходной функции. Для этого рассмотрим выражение для невязки:mXψn = −mXai u(tn − iτ )i=0+τmXbi f (tn − iτ, u(tn − iτ )) = −ai un−ii=0i=0+τmXbi f (tn−i , un−i ).i=0Теперь разложим составляющие равенства в ряд Тейлора:u(tn − iτ )=k(j)Xunj=0f (tn−i , un−i )j!(−iτ )j + O(τ k+1 );= {u0 = f(t, u)} = u0n−i = u0 (tn − iτ ) =k−1Xj=0(j+1)unj!(−iτ )j + O(τ k ).Подставив эти формулы в выражение для невязки, получим:mkmk−1(j)XXXX u(j+1)1unnψn = −ai (−iτ )j +bi(−iτ )j + O(τ k ).τ i=0j!j!i=0j=0j=080Глава 5.
ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧПоменяем порядки суммирования в двойных суммах, при этом из первой суммы вынесем отдельносоставляющую при j = 0, а во второй сделаем замену j = j + 1. Тогда выражение для невязкиперепишется так:ψn= −mkmk(j) m(j)XX1X1 X un Xunun ai −ai (−iτ )j +bi (−iτ )j−1 + O(τ k ) =τ i=0τ j=1 j! i=0(j−1)!i=0j=1= −mkm(j)Xun Xun j−1 Xai +τ(−i)j−1 (iai + jbi ) + O(τ k ).τ i=0j!j=1i=0Заметим, что если все суммы подбором коэффициентов ai , bj обратить в нуль, то для невязки будетсправедлива оценка:ψn = O(τ k ).Таким образом, достаточным условием k-го порядка аппроксимации будет выполнение системыравенств: Xmai = 0;i=0(5.7)mXj−1(−i) (iai + jbi ) = 0, j = 1, k.i=0Рассмотрим отдельно последнее условие при j = 1 :mXi=0Согласно условию нормировки,mXiai +mXbi = 0.(5.8)i=0bi = 1, поэтому (5.8) перепишется так:i=0mXi=0iai = −1 ⇐⇒mXiai = −1.i=1Добавив это уравнение и условие нормировки в систему (5.7), получим окончательный вариантдостаточного условия k-го порядка аппроксимации: mXbi = 1;i=0m Xai = 0;i=0(5.9)mXiai = −1;i=1mXij−1 (iai + jbi ) = 0, j = 2, k.i=1Мы получили систему из k + 2 линейных уравнений, решив которую, мы получим параметры,определяющие метод k-го порядка аппроксимации.
Система содержит 2m + 2 неизвестных. Чтобы онане была переопределенной, потребуем, чтобы k + 2 6 2m + 2.Таким образом, порядок аппроксимации m -шагового линейного метода не может превышать 2m —для неявного метода. Если же метод явный, то одним неизвестным в системе становится меньше, имаксимально возможный порядок аппроксимации будет равен 2m − 1.Перейдем к практическим примерам.5.4. МЕТОДЫ АДАМСА И ГИРА5.481Методы Адамса и ГираОпределение.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
