Конспект лекций - 3ий поток, лектор - Ионкин (1113828), страница 14

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 14)

Метод Рунге-Кутта (двухэтапный) или схема предиктор-корректор.Поставим в соответствие задаче (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) получаем:1yn+1 = yn + τ f (tn + , 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 (tn + a2 τ, un + b21 τ f (tn , un ))τРазложим un+1 в ряд Тейлора в окрестности точки (tn , un ):un+1 − unτ= u0n + u00n + O(τ 2 )τ2(6)99Далее разложим 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 с учетом проведенных преобразований:∂fn∂fnτ 0020+ b21 τ fn+ O(τ 2 ).ψn = − un + un + O(τ ) + σ1 fn + σ2 fn + a2 τ2∂t∂nЗаметим, чтоu00nd=dtdudt=d∂fn∂fn(fn ) =+ fn.dt∂t∂tТогда погрешность аппроксимации приобретет вид:∂fn∂fn∂fn∂fn0ψn = − un + 0.5τ+ fn+ (σ1 + σ2 ) fn + σ2 a2 τ+ σ2 b21 τ+ O(τ 2 ) =∂t∂n∂t∂u∂fn∂fn+ τ (σ2 b21 − 0.5)fn+ O(τ 2 ).= −u0n + (σ1 + σ2 )fn + (σ2 a2 − 0.5)τ∂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)),dtu(0) = u0t>0(1)100Поговорим об 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 τ 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 ).Пример. Схема Рунге–Кутта четвертого порядка.yn+1 − yn1= (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)1011Выше было показано, что если σa2 = σb2 = , то будет второй порядок аппрок2симации.yn+1 − yn= (1 − σ)fn + σf (tn + aτ, yn + aτ fn ),τ(2)y0 = u0 , n = 0, 1, .

. .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) в эквивалентном виде, только сформировав погрешность аппроксимации:un+1 − unzn+1 − zn=−+ (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 ;2102Тогда(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 τ.103§3 Многошаговые разностные методы решения задачи Коши du= f (t, u(t)), t > 0dtu(0) = u u = u(t )0nn(1)Определение. М-шаговым разностным методом решения задачи (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=0104Переходим к вычислению оценки погрешности аппроксимации на решение:ψ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(tn )(l+1) + O(τ p )Подставим все это в погрешность аппроксимации и проведем очевидные преобразования:pp−1mmXXXak X (−kτ )l (l)(−kτ )l (l+1)u (tn ) +bku(tn ) + O(τ p ) =ψn = −τl!l!l=0k=0l=0k=0pmXXak (−kτ )lpmXX(−kτ )l−1 (l)u (tn ) + O(τ p ) =τl!l(l−1)!l=1 k=0l=0 k=0!pmmmXXXakak (−kτ )l (l) X (−kτ )l−1 l=−lbku(tn ) +−u +u (tn ) + O(τ p ) =ττl!l!k=0k=0k=1l=1!pmmXXXak(−kτ )l−1 (l)=−u(tn ) −u (tn )(kak + lbk ) + O(τ p ).τl!k=0l=1k=0=−(l)u(tn )+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.105Среди всех многошаговых методов широко известен метод Адамса:myn − yn−1 X=bn fn−kτk=0a0 = 1,a1 = −1,ak = 0,k>2Выбирая коэффициенты bk в правой части, можно достичь любую погрешность.§4 Понятие устойчивости многошаговых разностныхметодовНикакой разностный метод с помощью компьютеров точно не реализуется, потому что:1.

входные данные могут быть не точными2. будет происходить машинное округлениеРассмотрим для примера следующую схему:yn+1 = qy(1)n = 0, 1, . . . ,где y0 задано, а q — любое число, возможно комплексное.Если |q| > 1, то процесс вычисление по этой формуле будет неустойчивым, таккак на каждом шаге мы находим приближенное значение ỹn = yn +δn . Следовательноỹn+1 = qyn + qδn = yn+1 + δn+1 , откуда видно, что δn+1 = qδn . А так как |q| > 1, тоδn+1 → ∞.Если же |q| 6 1, то δn+1 не возрастает, и мы получим устойчивый метод.Рассмотрим модельную задачу:(U 0 (t) + λU (t) = 0, t > 0, λ > 0(3)U (0) = U0Её решение имеет вид U (t) = U0 e−λt .

Если λ > 0, то |U (t)| 6 U0 , т.е. имеет местоустойчивость по начальному условию.Рассмотрим следующую задачу: dU= f (t, U (t)), t > 0(2)dtU (0) = U0Явная схема Эйлера для задачи (2) представляется в видеyn − yn−1= f (tn , yn )τ106А применительно к модельной задаче, она будет выглядеть следующим образом: yn − yn − 1 + λyn−1 = 0τy = U00Мы можем разрешить относительно yn :yn = yn−1 − τ λyn−1Следовательно, вводя обозначение q = 1 − τ λ, получимyn+1 = (1 − τ λ)yn = qynЕсли |q| 6 1, то эта разностная схема устойчива.Получим, что−1 6 1 − τ λ 6 1.Правая часть автоматически выполняется, следовательно τ λ 6 2 и0<τ 62λ(4)Условие (4) означает устойчивость разностной схемы.Рассмотрим неявную схему Эйлера для этой задачи:yn+1 − yn= f (tn+1 , yn+1 )τА применительно к модельной задаче, она будет выглядеть следующим образом:yn+1 − yn+ λyn+1 = 0τВыразим yn+1 через yn .yn+1 − τ λyn+1 = yn(1 + τ λ)yn+1 = yn11yn , q =1 + τλ1 + τλВидно, что 0 < q < 1.

Следовательно, вне зависимости от шагов метод сталабсолютно устойчивым.yn+1 =Общий m-шаговый разностный метод dU= f (t, U (t)), t > 0dtU (0) = U0(1)107Модельная задача примет следующий вид: dU+ λU (t) = 0, t > 0dtU (0) = U0(2)Тогда разностный метод для задачи (1) примет вид:mXakk=0τyn−k =mXbk fn−k(3)k=0Где U0 = y0 , y1 , . . . , ym−1 заданы и y0 6= 0, а ak , bk не зависят от τ .Этот же разностный метод, применительно к модельной задаче запишется ввиде:mXak(4)( + bk λ)yn−k = 0τk=0Перепишем (4) какmX(ak + τ λbk )yn−k = 0(5)k=0Решение данного уравнения ищется в виде yj = q j . Если эту формулу подставитьв уравнение (5), то в силу однородности сократим на q n−m и получим уравнеиеmX(ak + τ λbk )q m−k = 0.F (q, τ ) =(6)k=0Уравнение (6) называется характеристическим для разностной схемы (4).

Характеристики

Список файлов лекций