Численные методы. Ионкин (2012) (неоффициальные) (косяки есть) (1160437), страница 14
Текст из файла (страница 14)
. . , um (t))T ,f (t, u(t)) = (f1 (t, u(t)), . . . , fm (t, u(t))T ,|u|2 = u21 + u22 + . . . + u2m ,Рассмотрим параллелепипедR = {(t, u) : |t| 6 a,|u − u0 | 6 b}.Определение. Функция f (t, u) удовлетворяет в R условию Липшица по второмуаргументу, если ∃L > 0 и выполнено неравенство:|f (t, u) − f (t, v)| 6 L|u − v|Пусть f (t, u) из (1) удовлетворяет условию Липшица в R. Тогда решение (1) u(t)существует и единственно при 0 < t < T .
Проинтегрируем первое уравнение из (1) иучтем начальное условие:Ztu(t) = u0 +f (t, u(x))dx09798На этом представлении основан метод Пикара:Ztun+1 (t) = u0 +f (t, un (x))dx0Этот метод не может быть эффективным методом решения задачи (1), так какинтеграл не всегда можно посчитать аналитически, да и сходимость была бы медленной. Поэтому для решения систем ОДУ применяются разностные методы: перваягруппа методов - методы Рунге-Кутта, вторая - многошаговые разностные методы(например, метод Адамса).Введем сетку:wτ = {tn = nτ, τ > 0, n = 0, 1, 2, .
. . }Пример 1. Схема Эйлера.Введем обозначения:un = u(tn ) fn = f (tn , un )Тогда схема Эйлера имеет вид: yn+1 − yn = f (t , y ),n nτy = u , n = 0, 1, . . .00tn ∈ wτ(2)Получили явную схему. Явная, так какyn+1 = yn + τ fn ,n = 0, 1, . . .Запишем погрешность аппроксимации:ψn = −un+1 − un+ fnτ(3)Разложим un+1 в ряд Тейлора в точке tn :un+1 − un= u0n + O(τ )τПодставим последнее выражение в (3):ψn = −u0n + fn + O(τ )Учитывая, что −u0n + fn = 0, окончательно получим:ψn = O(τ ),Это означает, что разностная схема (2) аппроксимирует исходную задачу.Если будет доказана оценка |yn − u(tn )| 6 M τ , где M не зависит от τ , то этобудет означать сходимость.99Пример 2. Метод Рунге-Кутта (двухэтапный) или схема предиктор-корректор.Поставим в соответствие задаче (1) разностную схему, введя полуцелый слой:tn+ 1 = tn + 0.5τ2Метод является двухэтапным, так как для нахождения решения в точке tn+1используются два этапа:tn −→ tn+ 1 −→ tn+12Выполним первый шаг по схеме Эйлера:yn+ 1 − yn20.5τВторой шаг:= f (tn , yn )yn+1 − yn= f (tn+ 1 , yn+ 1 ),22τ(4)(5)где y0 = u0 , n = 0, 1, .
. .Далее видно, что из (4) следует:yn+ 1 = yn + 0.5τ fn ,2а из (5) получаем:yn+1 = yn + τ f (tn+ 1 , yn + 0.5τ fn ).2Позже будет показано, что погрешность аппроксимации этого метода имеет второй порядок по τ . Оценка погрешности общего двухэтапного метода Рунге–Кутта.Рассмотрим общий вид метода:yn+1 − yn= σ1 K1 + σ2 K2 ,τy0 = u0 , n = 0, 1, 2 . . . , σ1 , σ2 ∈ R,K1 = f (tn , yn ), K2 = f (tn + a2 τ, yn + b21 τ K1 ),где σ1 , σ2 , a2 , b21 — вещественные числа.Подставим значения K1 и K2 в уравнение:yn+1 − yn= σ1 f (tn , yn ) + σ2 f (tn + a2 τ, yn + b21 τ f (tn , yn ))τТогда можем записать погрешность аппроксимации на решение (1):ψn = −un+1 − un+ σ1 f (tn , un ) + σ2 f (tn + a2 τ, un + b21 τ f (tn , un ))τРазложим un+1 в ряд Тейлора в окрестности точки (tn , un ):τun+1 − un= u0n + u00n + O(τ 2 )τ2(6)100Далее разложим f (tn + a2 τ, un + b21 τ fn ) в окрестности точки (tn , un ):f (tn + a2 τ, un + b21 τ fn ) = f (tn , un ) + a2 τ∂fn∂fn+ b21 τ fn+ O(τ 2 )∂t∂uПерепишем ψn с учетом проведенных преобразований:τ 00∂fn∂fn02ψn = − un + un + O(τ ) + σ1 fn + σ2 fn + a2 τ+ b21 τ fn+ O(τ 2 ).2∂t∂nЗаметим, чтоu00nd=dtdudt=d∂fn∂fn(fn ) =+ fn.dt∂t∂tТогда погрешность аппроксимации приобретет вид:∂fn∂fn∂fn∂fn0+ fn+ σ2 b21 τ+ O(τ 2 ) =ψn = − un + 0.5τ+ (σ1 + σ2 ) fn + σ2 a2 τ∂t∂n∂t∂u∂fn∂fn= −u0n + (σ1 + σ2 )fn + (σ2 a2 − 0.5)τ+ τ (σ2 b21 − 0.5)fn+ O(τ 2 ).∂t∂nПотребуем, чтобы были выполнены следующие условия:1.
σ1 + σ2 = 1 (условие аппроксимации)2. σ2 a2 = σ2 b21 = 0.5 (для того, чтобы достичь второго порядка погрешностиаппроксимации)Если выполнено условие (1), то ψn = O(τ ), а если выполнены оба условия, ψn =O(τ 2 ). Положим σ2 = σ, а σ1 = 1 − σ, тогда условие 1 автоматически выполнено имы получим однопараметрические семейство:yn+1 − yn= (1 − σ)K1 + σK2 .τВ примере предиктор-корректор параметры имели следующие значения:a2 = b21 = 0.5,σ=1Если взять σ = 0.5, b21 = a2 = 1, то получим симметричную схему:yn+1 − yn= 0.5(f (tn , yn ) + f (tn+1 , yn+1 )).τ§2 Методы Рунге-КуттаРешаем задачу Коши: du= f (t, u(t)),dtu(0) = u0t>0(1)101Поговорим об m–этапном методе Рунге–Кутта. Идея заключается в переходеtn → tn+1 , используя m промежуточных этапов.yn+1 − yn= σ1 K1 + σ2 K2 + .
. . + σm Kmτy0 = u0 , n = 0, 1, . . . ,(2)гдеK1 = f (tn , yn ),K2 = f (tn + a2 τ, yn + b21 τ K1 ),K3 = f (tn + a3 τ, yn + b31 τ K1 + b32 τ K2 ),...Km = f (tn + am τ, yn + bm1 τ K1 + bm2 τ K2 + . . . + bmm−1 τ Km−1 ).Условие аппроксимации:mXσi = 1i=1На практике редко используются методы Рунге–Кутта для m > 4. Приведемпримеры разностных методов Рунге–Кутта, имеющих третий и четвертый порядокпогрешности аппроксимации.Пример.
Схема Рунге–Кутта третьего порядка.yn+1 − yn1= (K1 + 4K2 + K3 )τ6K1 = f (tn , yn ),K2 = f (tn + 0.5τ, yn + 0.5τ K1 ),K3 = f (tn + τ, yn − τ K1 + 2τ K2 ).Данная схема имеет третий порядок аппроксимации по τ : ψn = O(τ 3 ).Пример. Схема Рунге–Кутта четвертого порядка.1yn+1 − yn= (K1 + 2K2 + 2K3 + K4 )τ6K1 = f (tn , yn ),K2 = f (tn + 0.5τ, yn + 0.5τ K1 ),K3 = f (tn + 0.5τ, yn + 0.5τ K2 ),K4 = f (tn + τ, yn + τ K3 ).Данная схема имеет четвертый порядок аппроксимации по τ : ψn = O(τ 4 ).Оценка точности на примере двухэтапного метода Рунге-КуттаВсе сложности порождены нелинейностью задачиdu= f (t, u(t)), t > 0,dtu(0) = u0(1)1021Выше было показано, что если σa2 = σb2 = , то будет второй порядок аппрок2симации двухэтапного метода Рунге–Кутта.yn+1 − yn= (1 − σ)fn + σf (tn + aτ, yn + aτ fn ),τy0 = u0 , n = 0, 1, .
. .(2)1Как правило σ неотрицательно. Для удобства положим 0 6 σ 6 1 и σa = .2Введем функцию погрешности zn :zn = yn − unПолучаем задачу:zn+1 − znun+1 − un=−+ (1 − σ)f (tn , yn ) + σf (tn + aτ, yn + aτ f (tn , yn )),ττz0 = 0, n = 0, 1, 2, . . .(3)Для сходимости нужно показать, что|zn | −→ 0,n −→ ∞.Покажем, что |zn | 6 M |ψ|, где M не зависит от τ . Перепишем (3) в эквивалентном виде, только сформировав погрешность аппроксимации:zn+1 − znun+1 − un=−+ (1 − σ)f (tn , un ) + σf (tn + aτ, un + aσf (tn , un ))+ττ+(1−σ)(f (tn , yn )−f (tn , un ))+σ(f (tn +aτ, yn +aτ f (tn , yn ))−f (tn +aτ, un +aτ f (tn , un ))) =(2)= ψn + (1 − σ)ϕ(1)n + σϕn ,(1)(2)где ψn , ϕn , ϕn — введенные функции.То есть получим:zn+1 − zn(2)= ψn + (1 − σ)ϕ(1)n + σϕnτВведем допущение: функция f по второму аргументу удовлетворяет условию(1)(2)Липшица с константой L.
Оценим, исходя из этого допущения, ϕn и ϕn :|ϕ(1)n | = |f (tn , yn ) − f (tn , un )| 6 L|yn − un | = L|zn |;|ϕ(2)n | 6 |f (tn + aτ, yn + aτ f (tn , yn )) − f (tn + aτ, un + aτ f (tn , un ))| 66 L|yn +aτ f (tn , yn )−un −aτ f (tn , un )| 6 L|yn −un |+aτ |f (tn , yn )−f (tn , un )|) 6 L(|zn |+aτ L|zn |)Теперь сделаем предположение следующего характера:1σa 6 ;2103Тогда(2)zn+1 = zn + τ ψn + (1 − σ)τ ϕ(1)n + σϕnЛегко видеть, что|zn+1 | 6 τ |ψn | + (1 − σ)τ L|zn | + στ L(|zn | + aτ L|zn |) == τ L|zn | + τ L (σ + σaτ L) |zn | + τ |ψn | + |zn |;|zn+1 | 6 τ |ψn | + (1 + τ L + 0.5τ 2 L2 )|zn |;Заметим, что 1+τ L+0.5τ 2 L2 являются первыми членами разложения по Тейлоруфункции eτ L .
Следовательно:|zn+1 | 6 eτ L |zn | + τ |ψn |;обозначая eτ L = ρ > 0, получим оценку:|zn+1 | 6 ρ|zn | + τ |ψn |;Это соотношение можно рассматривать как рекуррентную формулу. Легко видеть, что:nXn|zn+1 | 6 ρ |z0 | +ρn−j τ |ψj |;j=0|zn+1 | 6nXρn−jτ |ψj | = max |ψ|06j6nj=0nXτ ρn−j = tn+1 max |ψ|ρn = tn+1 etn L kψkc06j6nj=0Окончательно получаем|zn+1 | 6 M kψkc ,где M > 0 не зависит от τ .Ясно, чтоzn → 0,при n → ∞,то есть имеет место сходимость с погрешностью аппроксимации.Напомним, что схема предиктор-корректор удовлетворяет условиям:1. σ = 1, a = 0.5,ψn = O(τ 2 ) =⇒ |zn+1 | 6 M τ 2 ;2. σ = 0.5, a = 1,ψn = O(τ 2 ) =⇒ |zn+1 | 6 M τ 2 .Если σ = 0,∀a ψn = O(τ ), то|zn+1 | 6 M τ.104§3 Многошаговые разностные методы решения задачи Коши du= f (t, u(t)), t > 0dtu(0) = u u = u(t )0nn(1)Определение.
Линейным m-шаговым разностным методом решения задачи (1)называется метод, записанный уравнением:mXakk=0τyn−k =mXbk fn−k(2)k=0yn−k = y(tn − kτ ),fn−k = f (tn − kτ, yn−k ),где ak , bk , (k = 0, m) действительные числа, причем a0 6= 0, bm 6= 0, τ > 0.Если b0 = 0, то метод явный, в противном случае назовем его неявным.Для начальных вычислений по формуле (2) необходимы значения y0 , y1 , . . .
, ym−1— так называемый "разгонный этап". Это приводит к некоторой сложности, таккак нам известно только значение y0 . Остальные, как правило, получаются другимиметодами. Поэтому будем считать, что все они заданы.Попробуем сопоставить многошаговые разностные методы с методом Рунге–Кутта.Плюсы многошагового разностного метода.1. Компактная красивая формула.2. Можно легко получить более высокий порядок погрешности аппроксимации.Минусы многошагового разностного метода.1.
Наличие разгонного этапа, так как y1 , y2 , . . . , ym нужно найти, применяя другиеметоды.2. При нахождении yk используются значения yn−1 , yn−2 , . . . , yn−m — их нужнопомнить.Если будет совсем высокий порядок, то мы позже убедимся, что это плохо —потеряется устойчивость.Условие нормировки:mXbk = 1(3)k=0105Оценка погрешности m–шагового разностного методаПереходим к вычислению оценки погрешности аппроксимации на решение:ψn = −mXakτk=0un−k +mXbk f (tn − kτ, un−k )k=0(4)un−k = u(tn − kτ )Применим формулу Тейлора в окрестности точки tn :un−k =pX(−kτ )ll!l=0f (tn − kτ, un−k ) =u0n−k=u(l) (tn ) + O(τ p+1 )p−1X(−kτ )ll=0l!u(l+1) (tn ) + O(τ p )Подставим все это в погрешность аппроксимации и проведем очевидные преобразования:p−1pmmXXXak X (−kτ )l (l)(−kτ )l (l+1)bku (tn ) +u(tn ) + O(τ p ) =ψn = −τ l=0l!l!l=0k=0k=0pmXXak (−kτ )lpmXX(−kτ )l−1 (l)u (tn ) + O(τ p ) =τl!l(l − 1)!l=0 k=0l=1 k=0!pmmmll−1XXXXakak (−kτ ) (l)(−kτ )=−u(tn ) +u +lbkul (tn ) + O(τ p ) =−ττl!l!k=0l=1k=1k=1!pmml−1XXXak(−kτ )=−u(tn ) −u(l) (tn )(kak + lbk ) + O(τ p ).τl!k=0l=1k=0=−u(tn )(l) +bkУсловие аппроксимации:mXak = 0(5)k=0Для достижения аппроксимации порядка p должно быть выполнено соотношение:mX(kak + lbk ) = 0,(6)k=0где l = 1, p.В многошаговом методе 2m + 2 неизвестных — a0 , a1 , .
. . , am , b0 , b1 , . . . , bm и p + 2уравнений. Чтобы система не была переопределенной, должно выполняться неравенствоp + 2 6 2m + 2 =⇒ p 6 2m,что означает, что порядок погрешности аппроксимации не выше, чем 2m.106Среди всех многошаговых методов широко известен метод Адамса:myn − yn−1 X=bn fn−kτk=0a0 = 1,a1 = −1,ak = 0,k>2Выбирая коэффициенты bk в правой части, можно достичь любую погрешность.§4 Понятие устойчивости многошаговых разностныхметодовНикакой разностный метод с помощью компьютеров точно не реализуется, потому что:1. входные данные могут быть не точными2.