Неявные разностные формулы

До сих пор мы записывали явные разностные схемы, т.е. необходимое значение функции определялось явным образом по известным точкам.

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

Аппроксимируем подынтегральную функцию интерполяционным полиномом Ньютона N_{n+ 1}(x), т.е. начнем с точки x_{n+ 1}, а не с точки x_n, как это мы делали в п.6.3.

Тогда для x = x_n+th имеем

.

Проведя интегрирование, получаем неявную (интерполяционную) формулу Адамса (в разностном виде):

y_{n+ 1} = y_n + h (f_{n+ 1} – D f_n – D² f_{n– 1} – D³ f_{n– 2} –

– D⁴ f_{n– 3} – D⁵ f_{n– 4} – …). (6.15)

Получим конкретные интерполяционные формулы Адамса в ординатном виде:

1) y_{n+ 1} = y_n + h f_{n+ 1} (неявная формула Эйлера). (6.16)

Формула первого порядка. Локальная погрешность ;

2) y_{n+ 1} = y_n + h f_{n+ 1} – D f_n =

= y_n + h f_{n+ 1} – (f_{n+ 1} – f_n) = y_n + (f_{n+ 1} + f_n) (6.17)

Формула второго порядка. Локальная погрешность ;

3) y_{n+ 1} = y_n + h (f_{n+ 1} – D f_n – D² f_{n– 1})

= y_n + (f_{n+ 1} + f_n) – (f_{n+ 1} – 2 f_n + f_{n– 1}) =

= y_n + (5 f_{n+ 1} + 8 f_n – f_{n –1}). (6.18)

Формула третьего порядка. Локальная погрешность ;

4) y_{n+ 1} = y_n + h (f_{n+ 1} – D f_n – D² f_{n –1} – D³ f_{n– 2}) =

= y_n + (9 f_{n +1} + 19 f_n – 5 f_{n –1} + f_{n– 2}). (6.19)

Формула четвертого порядка. Локальная погрешность .

В принципе, целый класс неявных разностных формул можно получить интегрированием интерполяционного многочлена Ньютона по отрезку [ x_n–j, x_{n+ 1}], j ³1 (при j = 0 – неявные формулы Адамса). Это приводит к неявным разностным формулам Милна.

	C 0	C 1	C 2	C 3
j = 1		–4
j = 2		–9		–9
j = 3		–16		–64

Например, в ординатном виде для j =1 получаем неявную формула Милна четвертого порядка (только по двум точкам!):

y_{n+ 1} = y_{n– 1} + (f_{n+ 1} + 4 f_n + f_{n– 1}).

Локальная погрешность

Знакомая формула, знакомая погрешность (см. формулу Симпсона). Т.е. здесь просто интеграл заменяется известной квадратурной формулой.

В чем недостаток неявных разностных формул по сравнению с явными?

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

В чем преимущество неявных разностных формул?

1. Используют меньше точек, чем соответствующие явные схемы.

2. Их точность выше, чем у соответствующих явных схем (сравните локальные погрешности формул (6.16) и (6.11), (6.17) и (6.12), (6.18) и (6.13)).

3. Шире диапазон устойчивости по сравнению с явными схемами.

С последним достоинством неявных разностных схем следует разобраться более подробно.

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

Проиллюстрируем это на примере.

Пример. Рассмотрим линейное уравнение вида:

y′ = ly, y (0) =y ₀, l < 0 – некоторая константа.

Точное решение этой задачи y (x) = e^lxy ₀. При возрастании x это решение убывает, сохраняя знак.

1. Для явного метода Эйлера имеем

y_{n+ 1} = y_n + hly_n = (1 + hl) y_n.

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

0 < 1 +lh < 1 Þ | l|h < 1.

Таким образом, есть ограничение на шаг:

h < 1/| l|.

2. Для явного метода Рунге-Кутта второго порядка:

y_{n+ 1} = y_n + (ly_n + l (y_n + hl y_n)) = y_n (1 + hl + ).

Чтобы решение этой разностной задачи монотонно убывало, необходимо выполнение условия

0 < 1 + hl + < 1 или – 2 < lh (lh + 2) < 0;

– 2 < lh < 0;| l|h < 2.

Ограничение на шаг: h < 2/| l|.

3. Для явного метода Адамса второго порядка:

y_{n+ 1} = y_n + (3 ly_n – ly_{n– 1}) = y_n (1 + hl) – hl y_{n –1};

y_{n+ 1} – (1 + hl) y_n + hl y_{n– 1} = 0.

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

Запишем для удобства это уравнение в виде

y_{n+ 2} – (1 + hl) y_{n+ 1} + hl y_n = 0.

Общее решение этого уравнения: y_n = C ₁ y_n ⁽¹⁾ + C ₂ y_n ⁽²⁾ (сумма двух линейно независимых решений). Это следует из теории дифференциальных уравнений. Будем искать линейно независимые решения в виде: y_n = qⁿ,где q ¹ 0 – некоторое число.

После подстановки получаем квадратное уравнение

q ² – (1 + hl) q + hl = 0.

Его дискриминант

D = (1 + hl)² – 2 hl = 1 + hl + (hl)² > 0при любых hl.

Следовательно, q ₁ и q ₂ действительны и различны.

Устойчивость означает |q ₁ |< 1, |q ₂ |< 1.

Воспользуемся следующим свойством квадратного уравнения q ² +bq+c= 0: корни уравнения не превосходят по модулю единицу, если

|b| £1 + c, c £ 1Þ | 1 + h| £1 + h; h £1;

| 2 + 3 lh| £ 2 + lh; |lh| £1.

Ограничение на шаг: h < 1/| l|.

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

y_{n+ 1} = y_n + hly_{n+ 1}; ;

, что для l< 0 выполняется всегда.

5. Неявный метод Адамса второго порядка.

y_{n+ 1} = y_n + (f_{n+ 1} + f_n) = y_n + (y_{n+ 1} + y_n)

,

что для всех l<0 выполняется всегда.

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

y_{n+ 1} = y_n + (5 f_{n+ 1} + 8 f_n – f_{n– 1})

устойчива при | lh|< 6. Сравните, явная схема Адамса третьего порядка имеет область устойчивости | lh|< 0.545.

Таким образом, мы видим, что при использовании явных схем существуют довольно жесткие ограничения на шаг. Например, для модельной задачи они имеют вид h £ C/|l|.

При больших |l| это условие ограничивает выбор шага. Отсюда и возник термин жесткие уравнения. Подчеркнем, что жесткость – это свойство уравнения, а не численной схемы. Следовательно, для обеспечения устойчивости явной схемы необходимо выбирать очень маленький шаг интегрирования, что приводит к большому числу шагов на больших промежутках интегрирования и, как следствие, чрезмерному возрастанию времени решения задач на ЭВМ и увеличению вычислительной погрешности.