Г.М. Кобельков - Курс лекций по численным методам (1160467), страница 12
Текст из файла (страница 12)
Поскольку нам нужно сделать h1 шагов,то в итоге точность будет довольно низкой: O(h).Всё зависит от того, насколько точно выбирается квадратурная формула. Модифицируем этот метод, выбравдругую формулу (формулу прямоугольников):hhyn+1 = yn + hf xn + , y xn +.(17)22Эта формуладаёт точность на каждом шаге O(h3 ), но есть одна неприятность: мы ничего не знаем проhy xn + 2 . Поэтому её нельзя использовать без дополнительных ухищрений. Поступим так: вычислим yn+ 12по обычному методу Эйлера:hyn+ 21 = yn + f (xn , yn ).(18)2Она даёт погрешность O(h2 ), но поскольку перед f в основной формуле уже есть один множитель h, в итогеточность O(h3 ) не изменится.
Таким образом, схема вычислений в этом случае имеет вид y 1 = y + h f (x , y ),nn nn+ 22(19)yn+1 = yn + hf (xn+ 21 yn+ 12 ).Другая модификация заключается в использовании формулы трапеций, которая выбирается в качестве основной:hyn+1 = yn + f (xn , yn ) + f (xn+1 , yn+1 ) .(20)2Тут та же трудность: мы не знаем ничего про yn+1 , поэтому его приходится приближать опять же с помощьюформулы Эйлераyn+1 = yn + hf (xn , yn )(21)и потом подставлять в предыдущую формулу. Фактически, здесь мы в каком-то смысле делаем один шаг методапростой итерации: сначала вычисляем начальное приближение, а потом уточняем его.Естественно, за повышение точности нам приходится платить: мы дважды вычисляем значение f .5.2.2. Метод Рунге – КуттаЭтот метод является естественным обобщением для модификаций метода Эйлера. Именно, рассмотрим систему поправокk1 (h) := hf (x, y),k2 (h) := hf (x + α2 h, y + β21 k1 ),...q−1Xkq (h) := hf x + αq h, y +βqj kj .j=143(22)После этого вычисление ведётся так:y(x + h) ≈ z(h) := y(x) +qXpj kj (h).(23)j=1Здесь имеется куча параметров, а именно q штук параметров αj и ещё столько же pj , а кроме того, ещё целаяштук).строго нижнетреугольная матрица параметров βkj (их будет q(q−1)2Определение.
Такой метод называется явным методом Рунге – Кутта.Положим ϕ(h) := y(x + h) − z(h).Определение. Будем говорить, что метод имеет порядок s, еслиϕ(0) = ϕ′ (0) = · · · = ϕ(s) (0) = 0,ϕ(s+1) (0) 6= 0.(24)Если метод имеет порядок s, то, по формуле Тейлора получаем такое выражение для ϕ:ϕ(h) =ϕ(s+1) (0) s+1 ϕ(s+2) (θh) s+2h+h .(s + 1)!(s + 2)!(25)Отсюда мораль: нужно стараться выбирать параметры так, чтобы s было максимальным.Пример 2.1. Пусть q = 1.
Тогда k1 (h) = hf (x, y), z(h) = y + p1 hf (x, y), отсюдаϕ(h) = y(x + h) − y − p1 hf (x, y).(26)′Очевидно, что ϕ(0) = 0. Вычислим p1 из условия ϕ (0) = 0:0 = ϕ′ (0) = y ′ (x) −p1 f (x, y) = 0,| {z }(27)=f (x,y)откуда p1 = 1. Таким образом, при q = 1 получаем обычный метод Эйлера.Пример 2.2. Пусть q = 2. Как и ранее, k1 (h) = hf (x, y).
Обозначим для краткостиx := x + α2 h,y := y + β21 k1 ,(28)тогда k2 (h) = hf (x, y). Функция ϕ имеет видϕ(h) = y(x + h) − y − p1 hf (x, y) − p2 hf (x, y).(29)0 = ϕ′ (0) = y ′ − p1 f − p2 f,(30)Имеем ϕ(0) = 0. Далее, из условия ϕ′ (0) = 0 получаемоткуда p1 + p2 = 1. Продифференцируем второй раз:0 = ϕ′′ (0) = y ′′ (x) − 2p2 fx (x, y)α2 + fy (x, y)β21 f == fx (x, y) + fy (x, y)f − 2p2 fx (x, y)α2 + fy (x, y)β21 f = fx (1 − 2p2 α2 ) + fy f (1 − p2 β21 ) .
(31)Чтобы это было так, нужно, чтобы коэффициенты при fx и fy f были нулевые. Значит, получается ещё двауравнения:2p2 α2 = 1, 2p2 β21 = 1.(32)Однако заметим, что одного уравнения, вообще говоря, не хватает: у нас 3 уравнения и 4 неизвестных. Такчто далее можно рассматривать различные дополнительные условия и получать разные формулы.Например, возьмём p1 = 0. Тогда p2 = 1, α2 = 21 и β21 = 12 . В этом случае получаются в точности формулыдля метода прямоугольников.А если взять p1 = p2 = 12 , то получим α2 = 1 и β21 = 1, поэтому возникнут в точности формулы, полученныеиз квадратурной формулы трапеций.При этом возникает естественное желание: раз у нас есть ещё одна степень свободы, то, быть может, можновыбрать константы так, чтобы метод имел порядок 3, а не 2.
Однако нас ждёт неудача: простейшее уравнениеy ′ = y даёт, как несложно показать, ϕ′′′ (0) 6= 0.Теперь приведём ещё немного результатов (без всяких обоснований). Если q = 3, то можно добиться порядкааппроксимации s = 3. Более того, если q = 4 или q = 5, то, как ни странно, можно получить всего лишь s = 4.Для реальных вычислений были построены методы 8-го порядка2 (Mérson). Кроме того, для специальныхсистем построены методы 14-го порядка. Очень недавно какой-то хмырь (история не сохранила его фамилии.
. . )построил схему, которая позволяет для полиномиальных уравнений получать методы сколь угодно высокогопорядка.2 Остаётсятолько догадываться, сколько там коэффициентов :)445.2.3. Метод Рунге априорной оценки погрешностиВозникает вопрос о том, а как оценить погрешность, если мы совсем ничего не знаем про точное решение?Мы уже отмечали, что если метод имеет порядок s, то ϕ(h) = M hs+1 + o(hs+1 ). Здесь M — главный членпогрешности. Это число, конечно, зависит от точки (x, y), но будем считать, что оно не сильно меняется.Пусть мы уже знаем y(x). Вычислим по методу Рунге – Кутта значение y(x + h).
Тогда мы при этом огребёмпогрешность порядка M hs+1 . Сделав ещё один шаг длины h, мы вычислим y(x + 2h), при этом получим числоyh , а суммарная погрешность составит 2M hs+1 .С другой стороны, можно сразу вычислить y(x + 2h), используя для этого удвоенный шаг. При этом получимнекоторое число y2h с погрешностью M (2h)s+1 . Итого получаем такую систему:(yh = y(x + 2h) + 2M hs+1 ,(33)y2h = y(x + 2h) + M (2h)s+1 .В этой системе нам неизвестны M и y(x + 2h). Решив её, мы получим хоть какую-то оценку на величинуконстанты M .5.2.4. Обобщение метода Рунге – КуттаКак мы видели из примеров, сколько ни увеличивай число q, ощутимых результатов (в смысле повышенияпорядка метода) мы не добьёмся.
Значит, нужно что-то ещё.Модификация метода Рунге – Кутта приводит к так называемому неявному методу Рунге – Кутта. Суть его втом, что мы дополняем матрицу βkj главной диагональю и несколькими (сколькими именно — вопрос отдельный,и мы его обсуждать не будем) строками над ней.5.3. Разностные схемы для решения дифференциальных уравненийРассмотрим уравнение y ′ = f (x, y(x)) и проинтегрируем его по отрезку I := [−nh, 0]:ZZ′y (x + t) dt = f (x + t, y(x + t)) dt.I(34)IЗаменим интеграл квадратурной формулой, получимqXa−i yn−ii=0h=qXi=0b−i f (xn−i , yn−i ).(35)Получили разностную схему нашего уравнения.Погрешность такой схемы — это величинаr(x) :=qXa−i y(x − ih)i=0h−qXi=0b−i f x − ih, y(x − ih) .(36)Определение.1◦ .
Если a0 6= 0, а b0 = 0, то схема называется явной.2◦ . Если a0 6= 0 и b0 6= 0, то получается, вообще говоря, нелинейное уравнение на yn видаa0yn − b0 f (xn , yn ) = . . . ,h(37)поэтому этот случай безнадёжен.3◦ . Если всё-таки a0 = 0, а b0 6= 0, то это неявная схема, или, как ещё говорят, «схема с забеганием вперёд».Теперь нужно найти условия на коэффициенты a−i и b−i , чтобы погрешность была как можно более высокогопорядка по h.Определение. Будем говорить, что схема имеет порядок s, еслиr(x) = E0 yh−1 + E1 y ′ + E2 y ′′ h + . . .
+ Es y (s) hs−1 + Es+1 hs + . . .(38)причём E0 = E1 = · · · = Es = 0, а Es+1 6= 0.Разложим решение в ряд Тейлора и подставим его в формулу для погрешности:y(x − ih) = y(x) − (ih)y ′ +45(ih)2 ′′y + ....2(39)Тогда сумма первых двух членов разложения первого слагаемого погрешности будут иметь видqqX1Xa−i y −a−i iy ′ .h i=0i=0Отсюда следует, чтоqXi=0a−i = 0,qXi=0ia−i = −1,(40)(41)потому что это слагаемое должно (при h → 0) сходиться к y ′ . Аналогично поступая со вторым слагаемымпогрешности, получаем, чтоqXb−i = 1,(42)i=0потому что это слагаемое в пределе должно давать f (x, y(x)).Если мы хотим получать схемы аппроксимации более высокого порядка, то нужно, чтобы занулялись коэффициенты при y ′′ и так далее.
Для второго порядка, например, появится уравнениеа для третьего порядка — ещё и уравнение−Xa−ii2 X+b−i i = 0,2(43)Xa−ii2i3 X−b−i = 0.62(44)5.3.1. Устойчивость схем в определениях и примерахРассмотрим уравнение y ′ = 0 и соответствующую схемуkXi=0a−i yn−i = 0.(45)Это некоторое рекуррентное соотношение, и, как мы знаем, его решение нужно искать в виде yn = µn . Подставляем в схему, получаем уравнение на µ:kXa−i µk−i = 0.(46)i=0Найдём решения этого уравнения, получим k корней (с учётом кратности). Если найдётся корень µj , длякоторого |µj | > 1, то существует экспоненциально растущее решение (при этом неважно, вещественный этоткорень или комплексный).Если |µi | 6 1 для всех i, но имеются кратные корни, то будет полиномиальный рост порядка на единицуменьше кратности корня.Если |µi | 6 1 для всех i, и на окружности {|µ| = 1} нет кратных корней, то уже всё хорошо: решениябудут ограниченными (кратные корни строго внутри круга не страшны, потому что их задавит экспонента сотрицательным показателем).Сейчас мы рассмотрим несколько примеров схем для решения одного и того же уравнения и посмотрим, какони себя ведут.
Будем решать уравнение y ′ + M y = 0 с начальным условием y(0) = y0 и M > 0.Пример 3.1. Рассмотрим самую простую схему (метод Эйлера). Она имеет видyn+1 = yn − M yn h = (1 − M h)yn .(47)yn = y0 (1 − M h)n ,(48)Отсюда получаемзначит, по крайней мере нужно, чтобы 0 6 1 − M h 6 1, тогда решение, вычисленное по этой схеме, будет похоже1на убывающую экспоненту. Это задаёт условие h 6 M.Возникает вопрос: а что, если взять схему с более высоким порядком аппроксимации.
Вот пример, которыйпоказывает, что «больше» — не всегда значит «лучше».Пример 3.2. Рассмотрим такую схему:yn+1 − yn−1+ M yn = 0.2h46(49)Преобразуя уравнение, получаемyn+1 + 2M hyn − yn−1 = 0.(50)Соответствующее характеристическое уравнение имеет видµ2 + 2M hµ − 1 = 0,откудаµ1,2 = −M h ±У этого уравнения всегда есть «плохой» кореньµ1 = −M h −(51)pM 2 h2 + 1.pM 2 h2 + 1,|µ1 | > 1,(52)(53)поэтому наша схема никогда не будет устойчивой.Чтобы читателю не показалось, что всё совсем плохо, мы приведём пример хорошей схемы с аппроксимациейпроизводной второго порядка.Пример 3.3.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
