НЕЯВНЫЕ РАЗНОСТНЫЕ ФОРМУЛЫ

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

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

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

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

.

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

y_n+1 = y_n + h(f_n+1 – Df_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 – Df_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 – Df_n – D²f_n–1)

= y_n + (f_n+1 + f_n) – ( f_n+1 – 2f_n + f_n–1) =

= y_n + (5f_n+1 + 8f_n – f_n–1). (6.18)

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

4) y_n+1 = y_n + h(f_n+1 – Df_n – D²f_n–1 – D³f_n–2) =

= y_n + (9f_n+1 + 19f_n – 5f_n–1 + f_n–2). (6.19)

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

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

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

y_n+1 = y_n–1 + (f_n+1 + 4f_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 + (3ly_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)² – 2hl = 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 + 3lh| £ 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 + (5f_n+1 + 8f_n – f_n–1)

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

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

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

ЖЕСТКИЕ ЗАДАЧИ

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

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

Рассмотрим систему дифференциальных линейных уравнений вида:

,

где y и f – векторы, A = {a_i,j} – постоянная матрица, у которой все собственные значения l_i имеют отрицательные вещественные части (i, j = 1, …, M). Таким свойством обладает, например, симметричная матрица. Общее решение записанной системы уравнений может быть представлено в виде суммы частных решений

,

где - собственный вектор, соответствующий собственному значению l_i.

Эту систему называют жесткой, если отношение

,

где S – коэффициент жесткости.

На практике система считается жесткой, если S > 10, однако в задачах химической кинетики, в задачах управления, в расчетах электрических цепей коэффициент S достигает величин ~10⁶ и более.

Посмотрим на поведение решения для жесткой системы из двух уравнений:

.

Пусть и l – вещественные.

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

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

Для решения жестких дифференциальных уравнений применяются неявные схемы. Хотя при этом приходится на каждом шаге решать в общем случае систему нелинейных уравнений, но за счет возможности значительно увеличить шаг h общий объем вычислений в случае использования неявных схем может быть существенно меньше, чем явных [11].

КРАЕВЫЕ ЗАДАЧИ

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

Методы решения краевых задач можно в целом разделить на три основных класса:

1) конечно-разностные методы;

2) методы стрельбы (пристрелки);

3) приближенные методы (решение представляется в виде линейной комбинации некоторых базисных функций): метод коллокаций, метод Галеркина, метод Релея-Ритца.