num-meth

где n = 0, 1, 2, . . . и y0(x) ≡ y0

Качественный результат: если в некоторой односвязной области G, содержащей такую точку (x0, y0), что:

|f(x, y)| ≤ C, |fy′ (x, y)| ≤ C1,

то найдется такая постоянная h > 0, что на отрезке [x0, x0+h] последовательность функций yn(x), определяемая методом 4.5, равномерно сходится к решению y(x) задачи Коши 4.1- 4.2, и справедлива оценка погрешности

Отметим две характеристики метода последовательных приближений Пикара, которые можно отнести к негативным. Во-первых, в силу известных проблем с нахождением первообразных, в чистом виде метод 4.5 редко реализуем. Во-вторых, как видно из вышеприведенного утверждения, этот метод следует считать локальным, пригодным для приближения решения в малой правой окрестности начальной точки.

4.1.2Метод Эйлера – разные подходы к построению

Рассмотрим несколько способов вывода метода Эйлера. При этом будем считать, что расчеты проводятся с расчетным шагом h = b−nx0 , расчетными точками (узлами)

служат точки xi = x0 + ih(i = 0, 1, . . . , n) промежутка [x0, b] и целью является построение таблицы

приближенных значений yi решения y = y(x) задачи 4.1-4.2 в расчетных точках xi.

Геометрический способ. Пользуясь тем, что в точке x0 известно и значение решения y(x0) = y0, и значение его производной y′(x0) = f(x0, y0), можно записать уравнение касательной к графику искомой функции y = y(x) в точке (x0, y0) :

по непрерывности должна мало отличаться от ординаты y(x1) решения y(x) задачи 4.1- 4.2. Следовательно, точка (x1, y1) пересечения касательной 4.6 с прямой x = x1 может быть приближенно принята за новую начальную точку. Через эту точку снова проведем прямую

y = y1 + f(x1, y1)(x − x1),

94 Глава 4. Численные методы решения обыкновенных дифференциальных уравнений

которая уже приближенно отражает поведение касательной к y = y(x) в точке (x1, y(x1)). Подставляя сюда x = x2(= x1 + h), иначе, пересекая эту "касательную"прямой x = x2, получим приближение значения y(x2) значением

y2 = y1 + hf(x1, y1),

и т.д. В итоге этого процесса, определяемого формулой

и называемого методом Эйлера, график решения y = y(x) данной задачи Коши 4.1-4.2 приближенно представляется ломаной, составленной из отрезков приближенных касательных, откуда происходит другое название метода 4.8 -метод ломанных.

Применение формулы Тейлора. Линеаризуя решение в окрестности начальной точки по формуле Тейлора, имеем

y(x) = y(x0) + y′(x0)(x − x0) + y′′(ξ)(x − x0)2. 2

Отсюда при x = x1 получаем

Точное равенство 4.9, переписанное в виде

y(x1) = y1 + r1(h),

говорит о том, что здесь мы имеем одновременно как саму формулу Эйлера для вычисления значения y1 ≈ y(x1), так и ее остаточный член

где ξ1− некоторая точка интервала (x0, x1)

Остаточный член 4.10 характеризует локальную(иначе шаговую) ошибку метода Эйлера, т.е. ошибку, совершаемую на одном шаге. За n шагов, т.е. в точке b, образуется

глобальная ошибка. Порядок глобальной ошибки (относительно шага h) на единицу ниже, чем порядок локальной ошибки, а порядком глобальной ошибки и определяется порядок соответстующего численного процесса решения задачи Коши. Таким образом, метод Эйлера относится к методам первого порядка.

Разностный способ. Рассматривая уравнение 4.1 в точке x0, с учетом 4.2 имеем равенство

Применяя к его левой части аппроксимацию производной правым разностным отношением первого порядка точности

что идентично равенству 4.9, поставляющему формулу для вычисления y1 вида 4.7 и локальный остаточный член 4.10. Порядок получающегося таким способом метода численного интегрирования дифференциальной задачи 4.1-4.2 совпадает с порядком аппроксимации производной в левой части уравнения 4.1.

Квадратурный способ. Начальную задачу для ОДУ 4.1-4.2 можно заменить эквивалентным интегральным уравнением 4.3. При x = x1 из него получается равенство

Применение к интегралу в правой части равенства 4.12 простейшей (одноточечной) формулы левых прямоугольников 4.6 дает приближенную формулу

y(x1) ≈ y0 + f(x0, y(x0))(x1 − x0),

правая часть которой, очевидно, совпадает с выражением 4.7 для подсчета значения y1. В общем случае расчетная формула 4.8 метода Эйлера получается численным интегрированием посредством простейшей формулы левых прямоугольников в равенстве

∫xi+1

xi

в предположении, что на каждом i-ом шаге в роли начальной точки выступает точка

(xi, yi).

4.1.2.1Несколько простых модификаций метода Эйлера

Если в 4.13 использовать простейшую квадратурную формулу правых прямоугольников 4.7, придем к методу

Этот метод называют неявным(или обратным)методом Эйлера.

Применение к интегралу в 4.13 простейшей квадратурной формулы трапеции приводит

который будем называть методом трапеций. Квадратурная формула трапеций на порядок точнее формул левых и правых прямоугольников. Отсюда вытекает более высокий (на единицу) порядок точности метода трапеций 4.15 по сравнению с явным и с неявным методом Эйлера, т.е. метод трапеций - это метод второго порядка.

Некоторый интерес представляет совместное применение явного метода Эйлера и неявного метода трапеций.

По форме равенство 4.15 представляет собой скалярную задачу о неподвижной точке относительно неизвестного значения yi+1. Поэтому если в правую часть 4.15 подставить хорошее начальное приближение yi0+1, подсчитываемое по формуле 4.14, то тогда само это равенство можно считать шагом метода простых итераций. Получаем гибридный метод

96 Глава 4. Численные методы решения обыкновенных дифференциальных уравнений

который называют методом Хойна.

Ясно, что можно достичь большей точности, если, исходя из того же начального приближения

yi0+1 = yi + hf(xi, yi),

Такой вариант совместного применения метода Эйлера и метода трапеций называют усовершенствованным методом Эйлера-Коши с итерационной обработкой

(или методом Эйлера с пересчетом). Делать много итераций по формуле 4.17 не рекомендуется (обычно их выполняют не более трех-четырех). Совпадение определенного числа разрядов в полученных числах yi(+1k) и yi(+1k−1) говорит о точности, с которой решено методом простых итераций уравнение 4.15 относительно yi+1, а вовсе не о том, что с такой точностью найдено значение y(xi+1).

Чтобы получить следующую модификацию метода Эйлера, проинтегрируем уравнение 4.1 по отрезку [xi−1, xi+1]. Имеем

Применяя к посленему интегралу одноточечную квадратурную формулу средних прямоугольников 4.10 и заменяя значения y(xi−1) и y(xi) известными приближенными значениями yi−1 и yi соответственно, из 4.18 выводим формулу для подсчета приближенного значения y(xi+1)

которую будем называть уточненным методом Эйлера.

Квадратурная формула прямоугольников (средней точки) имеет тот же порядок точности, что и квадратурная формула трапеций, так что уточненный метод Эйлера тоже является методом второго порядка. Подтверждением этого факта может служить вывод метода 4.19 на разностной основе. Применив к равенству 4.11 формулу симметричной аппроксимации y′(xi) второго порядка точности, получим

y(xi+1) − y(xi−1) ≈ f(xi, y(xi)),

2h

откуда после приближенной замены y(xi−1) ≈ yi−1, y(xi) ≈ yi, y(xi+1) ≈ yi+1 следует 4.19. Принципиальное отличие метода 4.19 от всех других рассмотренных до этого момента методов: метод 4.19 является двухшаговым. Здесь для вычисления значения yi+1 привлекаются два предыдущих значения yi и yi−1. Двухшаговость накладывает определенные ограничения, по крайней мере, на начало численного процесса: значение y1 ≈ y(x1) не может быть найдено непосредственно этим методом с тем же шагом h. Поэтому недостающую вторую начальную для процесса 4.19 точку приходится получать другим путем, например, явным методом Эйлера, а чтобы не сделать сразу большой ошибки, применяя на

старте метод более низкого порядка точности, рекомендуется осуществлять постепенное вхождение в процесс 4.19. Так, “разгон” можно выполнить по формулам

а далее уже переключаться на счет по формуле 4.19.

4.1.2.2Исправленный метод Эйлера

Пусть найдено приближенное значение yi ≈ y(xi) решения y = y(x) задачи 4.1 и требуется вычислить yi+1 ≈ y(xi+1), где xi+1 = xi + h. Запишем разложение решения по формуле Тейлора p-го порядка, принимая за базовую точку xi. Имеем

Если ограничится двумя слагаемыми в правой части разложения 4.21, то получим обычный метод Эйлера. Учет третьего слагаемого дает:

Подставляя приближенные выражения y(xi), y′(xi), y′′(xi) в равенство 4.22, получаем

Определямый ею метод будем называть исправленным методом Эйлера. Так как при i = 0 формулы 4.23 и 4.24 точны, а y0 = y(x0), то на первом шаге вычислений по формуле 4.25 будет совершаться ошибка, связанная только с усечением ряда Тейлора. Следовательно, локальная ошибка или, иначе, шаговая погрешность метода 4.25 составляет величину (h3), а это означает, что исправленный метод Эйлера относится к методам второго порядка.

4.2О семействе методов Рунге-Кутты. Методы второго порядка

Недостатком исправленного метода Эйлера 4.25 является необходимость вычисления на каждом шаге частных производных функции f(x, y). Идея построения явных методов Рунге-Кутты р-ого порядка заключается в получении приближений к значениям f(xi+1) по формуле вида

98 Глава 4. Численные методы решения обыкновенных дифференциальных уравнений

где φ(x, y, h) - некоторая функция, приближающая отрезок ряда Тейлора 4.21 до р-ого порядка и не содержащая частных производных функции f(x, y).

Полагая в 4.26 φ(x, y, h) ≡ f(x, y), приходим к методу Эйлера 4.28, т.е. метод Эйлера можно считать простейшим примером методов Рунге-Кутты, соответствующим случаю p = 1. Для построения методов Рунге-Кутты порядка, выше первого, функцию φ(x, y, h) берут многопараметрической и подбирают ее параметры сравнением выражения 4.26 с многчленом Тейлора для y(x) соответствующей желаемому порядку степени. Рассмотрим случай p = 2. Возьмем функцию φ в 4.26 следующей структуры:

φ(x, y, h) := c1f(x, y) + c2f(x + ah, y + bhf(x, y)).

Ее параметры c1, c2, a, b будем подбирать так, чтобы записанная, согласно 4.26, формула

определяла метод второго порядка, т.е. чтобы максимальная локальная ошибка составляла величину (h3). Разложим функцию f(x + ah, y + bhf(x, y)) по формуле Тейлора, ограничиваясь линейными членами:

f(x + ah, y + bhf(x, y)) = f(x, y) + fx′ (x, y)ah + fy′ (x, y)bhf(x, y) + (h2).

Ее подстановка в 4.27 дает

шения y(x) 4.28 с точностью до (h3) равносильно сравнению его с выражением yi+1 по формуле 4.25, т.е. с исправленным методом Эйлера.

В результате подстановки этих значений параметров в формулу 4.27 приходим к однопараметрическому семейству методов Рунге-Кутты второго порядка.

h

yi+1 = yi + 2 [f(xi, yi) + f(xi + h, yi + hf(xi, yi))] ,

При α = 1 из 4.30 новый простой метод, который назовем методом средней точки:

4.2.1Методы Рунге-Кутты произвольного и четвертого порядков

Любой метод из семейчства методов Рунге-Кутты второго порядка 4.30 реализуют по следующей схеме. На каждом шаге, т.е. при каждом i = 0, 1, 2, . . . , вычисляют значения функции:

а затем находим шаговую поправку

[ ]

yi := h (1 − α)η1i + αη2i ,

Прибавление которой к результату предыдущего шага дает приближенное значение решения y(x) в точке xi+1 = xi + h :

yi+1 = yi + yi.

Метод такой структуры называют двухэтапным по количеству вычислений значений функции на одном шаге. По аналогии с предыдущим для семейства методов Рунге-Кутты р-ого порядка используется запись, состоящая из следующей совокупности формул:

где k = 2, 3, . . . , p (для p-этапного метода). Многочисленные параметры ck, ak, bkj, фигурирующие в формулах 4.32, подбираются так, чтобы получаемое методом 4.32 значение yi+1 совпадало со значением разложения y(xi+1) по формуле Тейлора с погрешностью(hp+1) (без учета погрешностей, совершаемых на предыдущих шагах). Наиболее употребительным частным случаем семейства методов 4.32 является метод Рунге-Кутты четвертого порядка, относящийся к четырехэтапным и определяемый следующими частными формулами:

Заметим, что если первый этап, как уже упоминалось, соответствует применению явного метода Эйлера, то четвертый-неявного, а второй и третий-уточненного методов Эйлера.