Численные методы. Ионкин (2012) (неоффициальные) (косяки есть) (1160437), страница 15
Текст из файла (страница 15)
будет происходить машинное округлениеРассмотрим для примера следующую схему:yn+1 = qyn(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(2)U (0) = U0Её решение имеет вид U (t) = U0 e−λt . Если λ > 0, то |U (t)| 6 U0 , т.е.
имеет местоустойчивость по начальному условию.Рассмотрим следующую задачу: dU= f (t, U (t)), t > 0(3)dtU (0) = U0Явная схема Эйлера для задачи (3) представляется в видеyn − yn−1= f (tn , yn )τ107А применительно к модельной задаче, она будет выглядеть следующим образом: 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.
Следовательно, вне зависимости от шагов метод сталабсолютно устойчивым. Таким образом показано, что для устойчивости дифференциальной задачи (2) разностные методы могут быть как устойчивыми, так и неустойчивыми.yn+1 =108Общий m-шаговый разностный метод dU= f (t, U (t)), t > 0dtU (0) = U0Модельная задача примет следующий вид: dU+ λU (t) = 0, t > 0dtU (0) = U0(1)(2)Тогда разностный метод для задачи (1) примет вид:mXakk=0τyn−k =mXbk fn−k(3)k=0где y0 , y1 , . . . , ym−1 заданы, а ak , bk не зависят от τ .Этот же разностный метод, применительно к модельной задаче запишется ввиде:mXak( + bk λ)yn−k = 0(4)τk=0Перепишем (4) какmX(ak + τ λbk )yn−k = 0(5)k=0Решение данного уравнения ищется в виде yj = q j .
Если эту формулу подставитьв уравнение (5), то в силу однородности сократим на q n−m и получим уравнеиеmXF (q, τ ) =(ak + τ λbk )q m−k = 0.(6)k=0Уравнение (6) называется характеристическим для разностной схемы (4). Дляустойчивости необходимо, чтобы его корни по модулю не превосходили 1. Поисккорней уравнения (6) представляет собой, как правило, сложнейшую задачу.
Поэтомубудем считать, что τ — мало и положим τ = 0. Тогда получимF (q, 0) =mXak q m−k = 0(7).k=0Уравнение (7) также называется характеристическим.Определение. Говорят, что разностная схема (3) удовлетворяет условию (α),если все корни характеристического уравнения (7) лежат внутри или на границеединичного круга комплексной плоскости, причем на границе нет кратных корней.109Теорема 9. Пусть разностная схема (3) удовлетворяет условию (α) и |fn0 | 6 L для0 6 t 6 T, при 0 6 tn = nτ 6 T . Тогда для всех достаточно малых τ выполняетсяоценка!nX|y(tn ) − U (tn )| 6 Mτ |ψ| + max |yi − U (ti )| ,06i6m−1j=mгде M не зависит от τ (M = M(L,T)).Замечание (1).
Метод Адамса удовлетворяет условию (α):myn − yn−1 X=bk fn−kτk=0y0 = U0a0 = 1,a1 = −1Тогда характеристическое уравнение имеет вид:q n − q n−1 = 0,оно имеет корни q0 = 0 и q1 = 1, причем q1 — некратный корень.Замечание (2). Говорить об условной или безусловной устойчивости не имеетсмысла. Она всегда условная, т.к. рассматриваются малые τ .Замечание (3). Для неявных схем наивысший порядок погрешности аппроксимации p 6 2m.
Для явных схем p 6 2m − 1.Однако, схемы высокого порядка не удовлетворяют условию (α), т.е. не являются устойчивыми. Наивысший порядок аппроксимации для схем, удовлетворяющих условию (α), следующий:1. Для явных схем:(a) Если m — четно, то p 6 m + 2(b) Если m — нечетно, то p 6 m + 12. Для неявных схем p 6 mЗадача. Доказать, что для разностной схемы2fn−1 + fn−2yn + 4yn+2 − 5yn2=6τ3погрешность аппроксимации O(τ 3 ).Доказательство.ψn = −Un + 4Un−1 − 5Un−2 2fn−1 + fn−2+6τ3110Запишем условия, налагаемые на многошаговый разностный метод для того, чтобыпогрешность аппроксимации имела третий порядок:mPb0 = 1 −bkk=1mPa=−ak 0k=1mPak k = −1k=1mP k l−1 (ak k + bk ) = 0, l = 2, 3.k=012521, a1 = , a2 = − , b0 = 0, b1 = , b2 = .63633Выписанные условия легко проверяются, следовательно, ψn = O(τ 3 ).В данном случае m = 2, a0 =§5 Жесткие системы ОДУРассмотрим систему ОДУ:dU1 (t)+ a1 U (t) = 0, t > 0dtU1 (0) = U1,0dU2 (t)+ a2 U (t) = 0, t > 0dtU2 (0) = U2,0a1 > 0, a2 > 0, a2 a1 ,то есть a2 много больше a1(1)Решение имеет вид:U1 (t) = U1,0 e−a1 t ,U2 (t) = U2,0 e−a2 t .В приведенном случае одна функция убывает медленно, а другая быстро.
Тогда получим,что при t = t∗ вторая компонента решения U2 (t) близка к нулю, и если еенаходить численно, то шаг сетки можно брать не обязательно маленьким:111u(t)u1,0u2,0tt∗1) Рассмотрим явную схему Эйлера:y1n+1 − y1n+ a1 y1n = 0τ(2)y2n+1 − y2n+ a2 y2n = 0τ(3)22, а для (3) 0 < τ <.a1a2А чтобы обеспечить устойчивость системы (1) мы должны выбрать минимальноезначение τ . Если учесть, что a2 a1 , то минимальный шаг будет у (3).2) Рассмотрим неявную схему Эйлера:Устойчивость схемы (2) достигается при 0 < τ <y1n+1 − y1n+ a1 y1n+1 = 0τy2n+1 − y2n+ a2 y2n+1 = 0τЭти схемы абсолютно устойчивые. А шаг ограничен только условием точности.Введем понятие жесткости системы ОДУЗадача Коши для линейных уравнений dU+ AU (t) = 0,dtU (0) = U 0Где A = A(m × m) не зависит от t.t>0(5)112Определение. Система дифференциальных уравнений (5) называется жесткой,если выполняются два условия:k = 1, m1) ReλAk > 0,max16k6m |ReλAk|2) s =>> 1Amin16k6m |Reλk |Задача Коши для нелинейных уравненийdU= f (t, U ),dtt > 0,U (0) = U 0(6)где U (t) = (U1 (t), U2 (t), .
. . , Um (t)), f (t, U (t)) = (f1 (t, U (t)), f2 (t, U (t)), . . . , fm (t, U (t))).Пусть V (t) — известное решение. Проведем процесс линеаризации в окрестностиизвестного решения V (t) : z(t) = U (t) − V (t). Получаем k уравнений вида:dz k= fk (t, z(t) + V (t) − fk (t, V (t))),dtk = 1, m.Раскладывая в окрестности известного решения относительно точки (t, V (t)), получим∂fk (t, V (t))∂fk (t, V (t))∂fk (t, V (t))dz kz1 (t)+z2 (t)+· · ·+zm (t)−fk (t, V (t))+o(|z|).= fk (t, V (t))+dt∂U1∂U2∂UmОбозначимJ(t, V (t)) = aij =∂fi (t, V (t)),∂Uji, j = 1, m.Следовательно:dz k= J(t, V (t))z(7)dtЛинейная система (7) — система первого приближения. Теперь введем число жесткости, как отношениеmaxk |ReλJk |s(t) =(8)mink |ReλJk |Определение.
Система (6) называется жесткой на решение V (t), 0 6 t 6 T , есливыполнено 2 условия:1)ReλJk < 0, k = 1, m2)s(t) 1113§6 Дальнейшее определение устойчивости и примеры разностных схем, интегрирования жестких системдифференциальных уравненийЗапишем задачу Коши: dU= f (t, U (t)),dtU (0) = U0t>0(1)Запишем модельную задачу: dU= λU (t)dtU (0) = U0(2)где λ — собственные значения матрицы J.Введем комплексное число µ = τ λ (µ = µ0 + iµ1 ).Определение. Областью устойчивости разностного метода называется множество точек комплексной плоскости µ, для которых данный метод, примененный куравнению (2), устойчив.Рассмотрим явную схему Эйлера :yn+1 − yn= f (tn , yn )τПрименительно к модельной задачи (2) разностная схема примет вид:yn+1 − yn= λynτВыразим yn+1 :yn+1 = yn + λτ yn = (1 + µ)ynМетод является устойчивым, если |q| 6 1 или, в нашей задаче, |µ + 1| 6 1.
Тогдаполучаем, что|1 + µ0 + iµ1 | 6 1,(1 + µ0 )2 + µ21 6 1.В данном случае область устойчивости представляет собой внутренностькруга с центром в точке (0, -1) и радиусом 1 в системе координат (µ0 , µ1 ):µ11µ01114Рассмотрим неявную схему Эйлера :yn+1 − yn= f (tn+1 , yn+1 )τПрименительно к модельной задачи (2) разностная схема примет вид:yn+1 − yn= λyn+1τРазрешим её относительно yn+1 :yn+1 = yn + λτ yn+1 ,(1 − µ)yn+1 = yn ,1yn .1−µ 1 6 1.Для устойчивости необходимо, чтобы q = 1−µПолучаем, что|1 − µ| > 1,yn+1 =|1 − µ0 − iµ1 | > 1,(1 − µ0 )2 + µ21 > 1.Областью устойчивости неявной схемы Эйлера является внешность круга радиуса 1 с центром в точке (1, 0).µ1µ0Определение.
Разностный метод А-устойчив, если область его устойчивостисодержит всю левую полуплоскость комплексной плоскости, т.е. Reµ < 0.Утверждение. Если разностный метод А - устойчив, то он абсолютно устойчив,т.е. ∀τ > 0.Утверждение. Явных А-устойчивых методов нет.Утверждение. Среди неявных А-устойчивых методов существуют разностныеметоды не выше второго порядка.115Рассмотрим симметричную схему:yn−1 − yn= 0.5(f (tn , yn ) + f (tn+1 , yn+1 ))τПрименительно к модельной задачи (2) разностная схема примет вид:yn−1 − yn= 0.5λ(yn + yn+1 ).τВыразим yn+1 :yn+1 = yn + 0.5µ(yn + yn+1 ),(1 − 0.5µ)yn+1 = (1 + 0.5µ)yn ,yn+1 =Получаем, что q =1 + 0.5µ.1 − 0.5µ1 + 0.5µ.