Решение дифференциальных уравнений второго порядка методом адамса. Многошаговые методы (методы Адамса). Применение дифференциальных уравнений с малым параметром для решения нелинейных трансцендентных и алгебраических уравнений

Адамс – английский астроном и математик XIXвека, который много занимался небесной механикой. При изучении траекторий планет ему постоянно приходилось численно интегрировать уравнения их движения. Желая минимизировать объем вычислений, Адамс разработал один из наиболее экономичных методов численного решения дифференциальных уравнений, к обсуждению которого мы теперь переходим.

Пусть - решение дифференциального уравнения. Для производной этой функции имеет место равенство

Интегрируя его между двумя точками сетки, получим соотношение

Мы не можем использовать это соотношение непосредственно для перехода в процессе решения задачи от -ой точки сетки к
-ой, поскольку функция
нам не известна. Чтобы сделать следующий шаг, нужно приближенно заменить эту функцию на такую функцию, которую можно вычислить. Опишем, как эта проблема решается в методе Адамса.

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

,
.

Его можно записать по формуле Лагранжа

где
специальные многочлены вида

которые мы уже рассматривали в третьей главе.

Главная идея метода Адамса заключается в том, чтобы для расчета
использовать формулу типа, приближенно заменяя в ней функцию
на интерполяционный многочлен
, составленный согласно по результатам предыдущих вычислений. Это приводит к рекуррентной формуле

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

В этом случае формула переходит в рекуррентную формулу метода Эйлера

обеспечивающую первый порядок точности. Такой результат сам по себе тривиален. Мы привели его только для того, чтобы показать, что для метода Адамса, как и для метода Рунге-Кутта, исходной точкой является схема Эйлера.

Перейдем к исследованию варианта
. В этом случае для аппроксимации функции
используется полином первой степени, построенный по значениям функциив двух точках
и
:

Подставляя его в формулу и проводя интегрирование, получим

Отметим следующую особенность рекуррентной формулы. Для расчета очередного значения сеточной функции
нужно знать ее значения в двух предыдущих точкахи
. Таким образом, формула начинает работать только со второй точки. Вычислить по нейнельзя. Это значение решения разностной задачи приходится вычислять каким-нибудь другим методом, например, методом Рунге-Кутта.

Рекуррентную формулу можно записать в виде разностного уравнения

Подсчитаем для него погрешность аппроксимации дифференциального уравнения

Предположим, что функция
имеет в интересующей нас области изменения аргументов непрерывные вторые производные, так что решение задачи
трижды непрерывно дифференцируемо. Запишем разложения Тейлора

Подставляя их в формулу, получим

Отсюда можно написать оценку

где
- постоянная, мажорирующая третью производную функции
:

,
.

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

Процесс построения более точных схем можно продолжить за счет увеличения
. При
получается схема третьего порядка точности, при
- четвертого и т.д. Схема четвертого порядка, как и в методе Рунге-Кутта, является наиболее употребительной, поэтому мы коротко остановимся на ее выводе и обсуждении.

Если написать интерполяционный полином третьей степени
на сетке из четырех точек,
,
,
и провести интегрирование, то рекуррентная формула примет вид:

Приведем еще одно форму записи этой формулы через так называемые конечные разности

Первая, вторая и третья разности приближенно соответствуют первой, второй и третьей производной функции
. Эквивалентность формул и легко проверить непосредственно. Формула иногда более удобна для организации вычислительного процесса и контроля точности.

Особенность метода Адамса проявляется в формуле еще сильнее, чем в формуле. Здесь для расчета очередного значения
нужно знать значенияв четырех предыдущих точках -,
,
,
. Таким образом, формула начинает работать только с четвертой точки. Вычислить по ней,,нельзя. Эти значения решения разностной задачи приходится рассчитывать другим методом, например, методом Рунге-Кутта.

Перейдем к обсуждению точности схемы. Если функция
имеет непрерывные четвертые производные по своим аргументам в интересующей нас области их изменения, так что решение задачи
пять раз непрерывно дифференцируемо, то разностное уравнение аппроксимирует дифференциальное уравнение с четвертым порядком точности относительно. Доказательство этого утверждения проводится также, как и для схемы второго порядка, только теперь в разложениях типа нужно удерживать больше членов. Четвертый порядок точности при аппроксимации уравнения обеспечивает четвертый порядок точности для погрешности решения
при предположении, что начальные значения для метода Адамса,,вычислены с такой же точностью. Они рассчитываются независимо и при этом важно, чтобы начальный этап вычислительного процесса не внес такую погрешность, которая исказит все последующие результаты.

Задача 5.

Построить решение задачи Коши, на отрезке
с шагом
по схеме Адамса второго и четвертого порядка. Сравнить результаты расчетов между собой, с результатами расчетов по схеме Рунге-Кутта и с аналитическим решением задачи.

Результаты расчетов приведены в четвертом и пятом столбцах таблицы 2. В соответствии с заданием, нужно сравнивать четвертый столбец со вторым и шестым, а пятый – с третьим и шестым. Напомним, что в шестом столбце приведено аналитическое решение (53) рассматриваемой задачи, так что сравнение с ним позволяет судить о точности приближенного решения по схеме Рунге-Кутта и схеме Адамса.

Расчет по схеме Адамса второго порядка точности начинается с , четвертого - с . Значение в четвертом столбце,,,в пятом столбце рассчитывались по схеме Рунге-Кутта соответствующего порядка, поэтому в таблице они оказываются одинаковыми с соответствующими данными второго и третьего столбцов. Сравнение результатов проведенных расчетов двумя методами с аналитическим решением задачи показывает, что их точность примерно одинакова.

Сравним схемы четвертого порядка точности в методе Рунге-Кутта и Адамса с точки зрения организации вычислительного процесса. Чтобы сделать один шаг по методу Рунге-Кутта, необходимо вычислить функцию
четыре раза, а в методе Адамса только один раз. В трех предшествующих точках функция
была уже вычислена на предыдущих шагах и вычислять ее снова нет необходимости. В этом заключается главное достоинство метода Адамса, которое особенно высоко ценилось в докомпьютерную эру.

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

С этой особенностью метода Адамса связана еще одна проблема. При численном интегрировании дифференциального уравнения часто приходится менять шаг . В методе Рунге-Кутта это не составляет труда, поскольку каждый шаг делается независимо от предыдущего. В методе Адамса ситуация иная. Здесь нужно либо изначально программировать весьма сложные формулы расчета с переменным шагом, либо после каждой смены шага заново проводить расчет первых трех точек по методу Рунге-Кутта. Только после этого можно переходить на стандартный счет по методу Адамса. Эти недостатки приводят к тому, что сегодня при компьютернах расчетах предпочтение часто отдается более удобному методу Рунге-Кутта.

Перед нами все та же задача Коши

f (1) (t )=F (t , f (t )), a £t £b , f (a )=f a .

В одношаговых методах значение f (t k +1) определялось только информацией в предыдущей точке t k . Представляется возможным повысить точность решения, если использовать информацию в нескольких предыдущих точках при ее наличии. Так и поступают в методах, которые называются многошаговыми. С первого взгляда на постановку задачи становится очевидным, что в момент старта t =t a есть только одно начальное условие и, если мы собираемся работать с двумя, тремя или четырьмя предыдущими точками, то не видно, как получить вторую, кроме использования одношаговых методов. Так и поступают; «комплексный» алгоритм решения может выглядеть так:

на первом шаге одношаговым методом получают вторую точку, на втором получают третью с помощью двухшагового метода, на третьем – четвертую с помощью трехшагового метода и т.д., пока для основного метода, который предполагается использовать, не наберется достаточно предыдущих точек.

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

Многошаговые алгоритмы можно создать так. Учитывая, что

f (t k +1)=f (t k )+ ,

можно численно проинтегрировать стоящую под знаком интеграла правую часть ОДУ. Если использовать метод прямоугольников (интерполяционный полином для интегрируемой функции – константа), получим обычный метод Эйлера. Если использовать 2 точки и интерполяционный полином первого порядка

p (x )= ,

то интегрирование по методу трапеций от t k до t k +1 даст следующий алгоритм:

f (t k +1)=f (t k )+0.5h (3F k -F k -1).

Аналогично для трех точек будем иметь квадратичный интерполирующий полином по данным (t k -2 , F k -2), (t k -1 , F k -1), (t k , F k ) и интегрирование по методу Симпсона даст алгоритм:

f (t k +1)=f (t k )+ (23F k –16F k -1 +5F k -2).

Для 4-х точек полином будет кубическим и его интегрирование даст:

f (t k +1)=f (t k )+ (55F k –59F k -1 +37F k -2 –9F k -3).

В принципе мы могли бы продолжать так сколь угодно долго.

Приведенные алгоритмы носят название методов Адамса-Башфорта второго, третьего и четвертого порядков.

Формально мы можем при построении интерполяционного полинома помимо N уже просчитанных точек использовать и еще R будущих t k +1 , t k +2 ; в простейшем случае набор

t k +1 , t k , t k -1 ,…, t k -N .

При этом порождается класс так называемых методов Адамса-Моултона. В четырехшаговом варианте он оперирует с данными (t k +1 , F k +1), (t k , F k ), (t k -1 , F k -1), (t k -2 , F k -2) и его алгоритм:

f (t k +1)=f (t k )+ (9F k +1 +19F k –5F k -1 +F k -2).

Нельзя, разумеется, вести расчет по отсутствующим данным, поэтому алгоритмы Адамса объединяют в последовательность алгоритмов Адамса-Башфорта и Адамса-Моултона, получая при этом так называемые методы прогноза и коррекции. Например, метод прогноза и коррекции четвертого порядка выглядит так: вначале прогнозируем по алгоритму Адамса-Башфорта с использованием «прошлых» точек

f (t k +1)=f (t k )+ (55F k –59F k -1 +37F k -2 –9F k -3).

Затем по вычисляем приближенное значение правой части уравнения

F k +1 =F (t k +1 , f (t k +1).

И, наконец, корректируем f (t k +1) с использованием его же приближенного значения

f (t k +1)=f (t k )+ (9F k +1 +19F k –5F k -1 +F k -2).

Наиболее эффективные из имеющихся компьютерных программ, позволяющих пользователю менять величину шага и порядок метода, основаны на методах Адамса высокого порядка (свыше 10). Опыт эксплуатации этих программ показывает, что различия в их реализации могут оказывать более существенное влияние на точность, чем различия во внутренних свойствах самих методов.

Явная схема Адамса.

Рассмотренные выше методы являются явными одношаговыми (для нахождения последующего приближения используется лишь одно предыдущее). Приведённый ниже метод является многошаговым.

Пусть задана задача Коши:

Для точного решения (которое нам не известно) выполнено:

Предположим, нам известны приближенные значения функции u(x) в k точках (стартовые k точек, в частности, можно найти методом Эйлера или методом Рунге-Кутта того или иного порядка), тогда функцию f (x, u(x)) в (2.4.2) для приближенного вычисления интеграла можно заменить на интерполяционный полином порядка k-1, построенный по k точкам, интеграл от которого считается явно и представляет собой линейную комбинацию значений c некоторыми множителями. Таким образом, мы получаем следующую рекуррентную процедуру вычисления приближенных значений функции u(x) (являющимся точным решением задачи Коши) в точках:

Описанная схема является k-шаговой явной формулой Адамса.

Неявная схема Адамса.

Пусть - интерполяционный полином порядка k, построенный по k+1 значению б одно из которых, именно, мы будем считать неизвестным. Модифицируем (2.4.3), заменив в нём на полином более высокой степени, интеграл от которого выражается в виде линейной комбинации значений с некоторыми новыми коэффициентами:

Формула (2.4.4) представляет собой неявную схему Адамса и является уравнением на, которое можно решать методом последовательных приближений. Естественно, что начальное приближение, должно быть разумно выбрано. Для этого удобно объединить явную и неявную схемы Адамса в одну, называемую «методом коррекции». Именно с помощью явной схемы определяется начальное приближение (прогноз), а затем по неявной схеме оно необходимое число раз (обычно один или два) корректируется методом последовательных приближений до достижения заданной точности (коррекция).

При S = 1 формула (6.16) примет вид

Если Q = 2, получим следующее вычислительное правило:

Обычно на практике используют экстраполяционную формулу (6.18), а затем корректируют полученное значение по формуле (6.23). И если результат уточненного значения не превышает допустимую погрешность расчета, то шаг H считается допустимым .

Для расчетов на компьютере формулы (6.18) и (6.23) в конечно-разностном виде неудобны. С учетом (6.21) их можно представить в виде

(6.24)

Приведенные формулы имеют достаточно большую точность. Они дают погрешность порядка ~ О (H4 ), но сами формулы оценки погрешности достаточно сложны. Приближенно погрешность можно оценить по правилу Рунге.

Пример 6.2. Решить дифференциальное уравнение на отрезке c начальным условием Y (X = 0) = 1. Найти методом Адамса (с коррекцией) в точке X 4 , в трех первых точках найти методом Рунге- Кутта, приняв шаг .

Решение. Значения функции в четырех первых точках возьмем из табл. 6.1 (см. пример в предыдущем разделе). Теперь стало понятно, зачем мы сохраняли значения первой производной в этих точках (см. формулы (6.24)).

X4 = X3 + H = 0.15 + 0.05 = 0.2;

Для того чтобы скорректировать полученный результат, необходимо вычислить значение производной в этой точке:

Теперь уточним значение по интерполяционной формуле (а можно этого и не делать, тогда погрешность метода будет больше):

Так как в качестве нового значения функции принято скорректированное, то Обязательно Следует пересчитать значение производной. В нашем случае модуль разности экстраполяционной и интерполяционной формул меньше ε, Что позволяет продолжить вычисления с тем же шагом.

Вопросы для самопроверки

· Сформулируйте задачу Коши для обыкновенных дифференциальных уравнений первого порядка.

· Что является решением дифференциального уравнения: а) в высшей математике, б) в прикладной математике?

· Какие методы дифференциальных уравнений называются одношаговыми, многошаговыми? Приведите примеры.

· Сравните , полученные на первом, втором шаге методами Эйлера, Рунге-Кутта и разложением в ряд Тейлора (трудоемкость, погрешность…).

· Как оценить погрешность применяемого метода? Как ее уменьшить?

· Сравните одношаговые и многошаговые методы решения дифференциальных уравнений, указав достоинства и недостатки первых и вторых.

· Что такое экстраполяционные и интерполяционные методы (формулы) Адамса?

· Можно ли применять: а) только экстраполяционные методы Адамса,
б) только интерполяционные?

· Можно ли использовать: а) многошаговые методы без одношаговых;
б) одношаговые методы без многошаговых?

· При решении дифференциального уравнения методом Адамса на 27-м шаге необходимо сменить шаг. Как это сделать?

В настоящее время методы Адамса являются одними из перспективных численных методов интегрирования для решения задачи Коши. Доказано, что при применении многошаговых численных методов Адамса для решения задачи Коши до 12 порядка область устойчивости уменьшается. При дальнейшем увеличении порядка область устойчивости, а также точность метода возрастает. Кроме того, при одинаковой точности для многошаговых методов на одном шаге интегрирования требуется меньше вычислений правых частей дифференциальных уравнений, чем в методах Рунге-Кутты. К достоинствам методов Адамса относится и то обстоятельство, что в них легко меняется шаг интегрирования и порядок метода.

На практике широко используются два типа методов Адамса - явные и неявные. Явные методы известны как методы Адамса-Бэшфорта, неявные - как методы Адамса-Мултона.

Рассмотрим применение численных методов для решения задачи Коши

При решении задачи (2. 1) с помощью одношаговых методов значение yn+1 зависит только от информации в предыдущей точке xn. Можно предположить, что можно добиться большей точности, если использовать информацию о нескольких предыдущих точках xn, xn-1... xn-k. На этой идее основаны многошаговые методы.

Большинство многошаговых методов возникает на основе следующего подхода. Если подставить в уравнение (2. 1) точное решение y (x) и проинтегрировать уравнение на отрезке , то получим:

Заменяя в формуле (2. 2) функцию f (x, y (x)) интерполяционным полиномом P (x), получим приближенный метод

Для того, чтобы построить полином P (x), предположим, что yn, yn-1... yn-k - приближения к решению в точках xn, xn-1... xn-k. Полагаем, что узлы xi расположены равномерно с шагом h. Тогда fi=f (xi, yi), (i=n, n-1.. n-k) - есть приближения к f (x, y (x)) в точках xn, xn-1... xn-k.

В качестве P (x) возьмем интерполяционный полином степени, k удовлетворяющий условиям

Если проинтегрировать этот полином явно, то получим следующий метод:

При k=0 полином P (x) - есть константа, равная fn, и формула (2. 4) превращается в обычный метод Эйлера.

При k=1 полином P (x) является линейной функцией, проходящей через точки (xn-1, fn-1) и (xn, fn), т. е.

Интегрируя этот полином от xn до xn+1, получим двухшаговый метод

который использует информацию в двух точках xn и xn+1.

Если k=2, то P (x) представляет собой квадратичный полином, интерполирующий данные (xn-2, fn-2), (xn-1, fn-1) и (xn, fn). Можно показать, что соответствующий метод имеет вид

Если k=3, то соответствующий метод определяется формулой

При k=4 имеем

Отметим, что метод (2. 7) является трехшаговым, (2. 8) - четырехшаговым и (2. 9) - пятишаговым. Формулы (2. 6) - (2. 9) известны как методы Адамса-Бэшфорта. Метод (2. 6) имеет второй порядок точности, поэтому его называют методом Адамса-Бэшфорта второго порядка. Аналогично, методы (2. 7), (2. 8) и (2. 9) называются соответственно методами Адамса-Бэшфорта третьего, четвертого и пятого порядков.

Продолжая этот процесс, используя все большее число предыдущих точек, а также интерполяционный полином более высокой степени, получим методы Адамса-Бэшфорта сколь угодно высокого порядка.

Многошаговые методы порождают трудности, которых не возникает при использовании одношаговых методов. Эти трудности становятся понятными, если, например, обратиться к методам Адамса-Бэшфорта пятого порядка (2. 9).

В задаче (2. 1) задано начальное значение y0 но при n=0 для счета по формуле (2. 9) необходима информация в точках x-1, x-2, x-3, x-4, которая естественно отсутствует. Обычный выход из данной ситуации заключается в использовании какого-либо одношагового метода того же порядка точности, например метода Рунге-Кутты, до тех пор, пока не будет получено достаточно значений для работы многошагового метода. Или же можно на первом шаге использовать одношаговый метод, на втором - двухшаговый и так далее, пока не будут получены все стартовые значения. При этом существенно, чтобы эти стартовые значения были вычислены с той же степенью точности, с какой будет работать окончательный метод. Поскольку стартовые методы имеют более низкий порядок точности, вначале приходится считать с меньшим шагом и использовать больше промежуточных точек.

Вывод методов (2. 6) - (2. 9) основан на замене функции f (x, y) интерполяционным полиномом P (x). Известно, что имеет место теорема, доказывающая существование и единственность интерполяцион ного полинома. Если узлы x0, x1… xn различны, то для любых f0, f1… fn существует единственный полином P (x) степени не выше n такой, что P (xi) =fi, i=0, 1,.. n.

Хотя интерполяционный полином является единственным, имеется несколько способов представления этого полинома. Чаще всего используются полиномы Лагранжа, но и они оказываются неудобными, если к набору данных нужно добавить (или удалить из него) какой-либо узел. В этом случае имеется другое представление интерполяционного полинома. Это представление Ньютона

Полином Pn+1 (x) можно записать в виде

Представление интерполяционного полинома в виде (2. 11) в ряде случаев бывает особенно полезным для практики.

Методы Адамса-Бэшфорта используют уже известные значения в точках xn, xn-1... xn-k. При построении интерполяционного полинома можно использовать и точки xn, xn, xn-1... xn-k. При этом возникает класс неявных m -шаговых методов, известных как методы Адамса-Мултона.

Если k=0, то P (x) - линейная функция, проходящая через точки (xn, fn) и (xn+1, fn+1), и соответствующий метод

является методом Адамса-Мултона второго порядка.

При k=1, 2, 3 получаем соответствующие методы

третьего, четвертого и пятого порядков аппроксимации. Соотношения (2. 12) - (2. 15) содержат искомые значения yn+1 неявно, поэтому для их реализации необходимо применять итерационные методы.

На практике обычно не решают непосредственно уравнений (2. 12) - (2. 15), а используют совместно явную и неявную формы, что приводит к методу прогноза и коррекции.

Например, для метода Адамса второго порядка, используя обозначения, где г - номер итерации, имеем для г=1 следующую схему вычислений:

Этот процесс называют методом PECE (P означает применение предсказывающей формулы, С - применение исправляющей формулы, Е - вычисление функции f). Можно сократить процесс вычисления, отбросив последнюю формулу. Это приводит к так называемому методу PEC.

Рассмотрим второй метод решения уравнений (2. 12) - (2. 15). Формулы (2. 12) - (2. 15) можно переписать в виде

где gn содержит известные величины. Доказано, что если, где L - константа Липшица, то существует единственное решение уравнения (2. 17), которое можно получить с помощью итерационного процесса

где - произвольно.

Итерации в выражении (2. 18) продолжаются до тех пор, пока не будет достигнута сходимость. При этом число вычислений функции f меняется от точки к точке и может быть достаточно велико.

С другой стороны, если уменьшить величину h, то сходимость может быть достигнута за фиксированное число итераций. Этот метод называется исправлением до сходимости.

На первый взгляд может показаться, что явный многошаговый метод является самым простым методом с точки зрения вычислений. Однако на практике явные методы используются очень редко. Неявный метод Адамса-Мултона является более точным, чем явный метод Адамса-Бэшфорта. Например, вычислительная схема для метода Адамса-Мултона 5-го порядка имеет следующий вид:

Методы Адамса до пятого порядка включительно могут быть использованы для решения обыкновенных дифференциальных уравнений, не требующих высокой степени точности.

Как и в случае с методом Адамса-Бэшфорта, при использовании метода Адамса-Мултона важным вопросом является вопрос выбора оптимального соотношения шага интегрирования и порядка метода. Следует отметить, что при создании эффективных алгоритмов и программ увеличение порядка метода является более предпочтительным по сравнению с уменьшением шага интегрирования.

Для решения более сложных задач необходимо применять методы Адамса более высокого порядка. В таблице 2. 1 приведены значения коэффициентов для методов Адамса. В первой строке указан порядок метода; во второй - значения коэффициентов Ck для соответствующего порядка k; в последующих строках - пары коэффициентов Bkj и Mkj для методов Адамса-Бэшфорта и Адамса-Мултона соответственно. Тогда, с учетом данных таблицы 2. 14, коэффициенты вj в выражении

для метода Адамса-Бэшфорта k-гo порядка могут быть найдены из cоотношения

а для метода Адамса-Мултона k-гo порядка по аналогичной формуле

Формулы для предикторно-корректорных методов Адамса с 6-го по по 14-ый порядок имеют следующий вид:

6 порядок:
7 порядок:
8 порядок:
9 порядок:
10 порядок:
11 порядок:
12 порядок:
13 порядок:
14 порядок:
15 порядок:
16 порядок:

Формулы, приведенные выше, предпочтительнее использовать для практического применения решения обыкновенных дифференциальных уравнений или систем дифференциальных уравнений первого порядка с постоянным шагом интегрирования. Если в процессе решения уравнения шаг интегрирования переменный, то для методов Адамса существуют специальные приемы для закладки новых начальных данных при смене шага интегрирования.