Разностные аппроксимации

1.
Примеры разностных аппроксимаций.

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

w
_h
={x_i=ih, i=0,
±
1,
±
2,…}.

Пусть u(x) – достаточно гладкая функция, заданная на отрезке [x_i-1, x_i+1]. Обозначим
Разностные отношения
называются соответственно правой, левой и центральной разностными производными функции u(x) в точке x_i , т.е. при фиксированном x_i  и при h®0 (тем самым при i®¥) пределом этих отношений является u’(x_i). Проводя разложение по формуле Тейлора, получим
u_x,i – u’(x_i) = 0,5hu’’(x_i) + O(h²),

u_x,i – u’(x_i) = -0,5hu’’(x_i) + O(h²),

u_x,i – u’(x_i) = O(h²),


Отсюда видно, что левая и правая разностные производные аппроксимируют u’(x) с первым порядком по h, а центральная разностная производная – со вторым порядком. Нетрудно показать, что вторая разностная производная




аппроксимирует u’’(x_i) со вторым порядком по h, причем справедливо разложение




Рассмотрим дифференциальное выражение




(1)
с переменным коэффициентом k(x). Заменим выражение (1) разностным отношением




                                                                                                                 (2)
где a=a(x) – функция, определенная на сетке w_h. Найдем условия, которым должна удовлетворять функция a(x) для того, чтобы отношение (au_x)_x,i аппроксимировало (ku’)’ в точке x_i со вторым порядком по h. Подставляя в (2) разложения





где u_i’ = u’(x_i), получим
С другой стороны,  Lu = (ku’)’ = ku’’ + k’u’,


т.е.
Отсюда видно, что L_hu–Lu = O(h²), если выполнены условия




(3)
Условия (3) называются достаточными условиями второго порядка аппроксимации. При их выводе предполагалось, что функция u(x) имеет непрерывную четвертую производную и k(x) – дифференцируемая функция. Нетрудно показать, что условиям (3) удовлетворяют, например, следующие функции:




Заметим, что если положить a_i = k(x_i), то получим только первый порядок аппроксимации.

В качестве следующего примера рассмотрим разностную аппроксимацию оператора Лапласа




(4)
Введем на плоскости (x₁, x₂) прямоугольную сетку с шагом h₁ по направлению x₁ и с шагом h₂ по направлению x₂, т.е. множество точек
w
_h
= {(xⁱ₁, x^j₂) | xⁱ₁ = ih₁, x^j₂ = jh₂; i, j = 0,
±
1,
±
2,…},
и обозначим




Из предыдущих рассуждений следует, что разностное выражение




(5)
аппроксимирует дифференциальное выражение (4) со вторым порядком, т.е. L_hu_ij – Lu(xⁱ₁, x^j₂) = O(h²₁) + O(h²₂). Более того, для функций u(x₁, x₂), обладающих непрерывными шестыми производными, справедливо разложение





Разностное выражение (5) называется пятиточечным разностным оператором Лапласа, так как оно содержит значения функции u(x₁, x₂) в пяти точках сетки, а именно в точках (x₁ⁱ, x₂^j), (x₁^i±1, x₂^j), (x₁ⁱ, x₂ ^j±1). Указанное множество точек называется шаблоном разностного оператора. Возможны разностные аппроксимации оператора Лапласа и на шаблонах, содержащих большее число точек.
2. Исследование аппроксимации и сходимости

2.1. Аппроксимация дифференциального уравнения. Ранее рассматривалась краевая задача

 (k(x) u’(x))’ – q(x) u(x) + f(x) = 0,   0 < x < l,                                       (1)

– k(0) u’(0) +
b
u(0) =
m
₁
,   u(l) =
m
₂
,                                                      (2)

k(x)
³
c₁ > 0,
b

³
0,

для которой интегро-интерполяционным методом была построена разностная схема




(3)

(4)
где




(5)




(6)




Обозначим через Lu(x) левую часть уравнения (1) и через L_hy_i – левую часть уравнения (3), т.е.






Пусть u
(x) – достаточно гладкая функция и u
(x_i) – ее значение в точке x_i сетки

w
_h
= {x_i = ih, i = 0, 1, …,N, hN = l}                                                        (7)
Говорят, что разностный оператор L_h аппроксимирует дифференциальный оператор L в точке x=x_i, если разность L_hu
_i
– L_h
u
(x_i) стремится к нулю при h®0. В этом случае говорят также, что разностное уравнение (3) аппроксимирует дифференциальное уравнение (1).

Чтобы установить наличие аппроксимации, достаточно разложить по формуле Тейлора в точке x=x_i значения u
_i
±
1
=
u
(x_i
±
h), входящие в разностное выражение L_hu
_i. Большая часть этой работы проделана в предыдущей главе, где показано, что при условиях
(8)
выполняется соотношение






Если кроме того, докажем, что
d_i = q(x_i) + O(h²),   
j
_i
= f(x_i) + O(h²)                                                      (9)
то тем самым будет установлено, что оператор L_h аппроксимирует L со вторым порядком по h, т.е.
L_h
u
_i
– L
u
(x_i) = O(h²),   i = 1, 2,…, N–1                                                   (10)
Итак, доказательство второго порядка аппроксимации сводится к проверке сводится к проверке условий (8), (9) для коэффициентов (5), (6). Проверим сначала выполнение условий (8). Обозначая p(x) = k^-1(x), получим





следовательно,





Аналогично





Отсюда получим





т.е. условия (8) выполнены. Условия (9) выполнены в силу того, что замена интегралов (6) значениями q_i, f_i соответствует приближенному вычислению этих интегралов по формуле прямоугольников с узлом в середине отрезка интегрирования.

2.2. Аппроксимация граничного условия. Исследуем погрешность аппроксимации разностного граничного условия (4). Обозначим l_hu
(0) = –a₁
u
_{x, 0}
+
b
u
₀. Если u
(x) – произвольная достаточно гладкая функция, то очевидно

l_h
u
(0) = –k(0)
u
’(0) +
b
u
(0) + O(h),
т.е. имеет место аппроксимация первого порядка по h. Однако если u
=u(x) – решение задачи (1), (2), то разностное граничное условие (4) имеет второй порядок аппроксимации, т.е.






Докажем последнее утверждение. Используя разложение
u_{x, 0} = (u₁ – u₀)/h = u’(x_1/2) + O(h²),   x_1/2 = 0,5h,

a₁ = k_1/2 + O(h²)

получим





Отсюда имеем





Учитывая граничное условие (2), получаем
l_hu(0) = 0,5h [– (ku’)’(0) + d₀u₀ –
j
₀
] + O(h²).
Выражение, стоящее в квадратных скобках, преобразуем, учитывая уравнение (1), к виду
– (ku’)’(0) + d₀u₀ –
j
₀
= – (ku’)’(0) + q(0)u(0) – f(0) +

                                  + (d₀ – q(0))u₀ – (f(0) –
j
₀
) = (d₀ – q(0))u₀ – (f(0) –
j
₀
).
Из соотношений





получаем






что и требовалось доказать.

Таким образом, при достаточной гладкости коэффициентов k(x), q(x), f(x) и решения u(x) разностная схема (10) аппроксимирует исходную задачу (2) со вторым порядком по h.

При практическом использовании разностной схемы для нахождения ее коэффициентов не обязательно вычислять интегралы (4), (6) точно. Можно воспользоваться коэффициентами, полученными путем замены этих интегралов квадратурными формулами, имеющими точность O(h²) и выше. Например, в результате применения формулы прямоугольников получим следующие коэффициенты: a_i = k(x_i – 0,5h), d_i = q(x_i), j
_i
= f(x_i).

Применяя формулу трапеций, получим





Представление коэффициентов разностной схемы в виде интегралов (4), (6) оказывается полезным при исследовании сходимости в случае разрывных функций k(x), q(x), f(x).
2.3. Уравнение для погрешности. Решение y_i = y(x_i) разностной задачи (3), (4) зависит от шага h сетки, y(x_i) = y_h(x_i). По существу, мы имеем семейство решений {y_h(x_i)}, зависящее от параметра h. Говорят, что решение y_h(x) разностной задачи сходится к решению u(x) исходной дифференциальной задачи, если при h®0 погрешность y_h(x_i) – u(x_i), i = 0, 1,…, N, стремится к нулю в некоторой норме. В настоящем параграфе в качестве такой нормы будем брать норму в сеточном пространстве C(w
_h
), т.е. положим






Говорят, что разностная схема имеет m-й порядок точности (или сходится с порядком m), если






где m>0, M>0 – константы, не зависящие от h.

Выше было установлено, что схема (3), (4) имеет второй порядок аппроксимации. Докажем теперь, что эта схема имеет и второй порядок точности. Для этого прежде всего выпишем уравнение, которому удовлетворяет погрешность z_i = y_i – u(x_i). Поставим y_i = z_i + u(x_i) в уравнения (3), (4). Тогда получим уравнения
(11)

(12)

где обозначено





Функция y
_i, входящая в правую часть уравнения (11), называется погрешностью аппроксимации дифференциального уравнения (1) разностным уравнением (3) на решении задачи (1), (2). В п.1 было доказано, что y
_i
= O(h²) при h®0, i=1, 2,…, N–1. Аналогично, величина n₁ является по определению погрешностью аппроксимации краевого условия (2) разностным краевым условием (4) на решении задачи (1), (2), причем n
₁
=O(h²). Таким образом, структура уравнений для погрешности (11), (12) та же, что и у разностной схемы (3), (4), отличаются только правые части.

Чтобы доказать сходимость разностной схемы, оценим решение задачи (11), (12) через правые части y
_i
,
n
₁, т.е. получим неравенство вида
                                                                                                                 (13)
с константой M₁, не зависящей от h. Из этого неравенства и будет следовать, что
Отметим, что неравенства вида (13), называемые априорными оценками, нашли широкое применение в теории разностных схем. Поскольку структура для погрешности (11), (12) та же, что и у разностной схемы (3), (4), а отличаются только правые части, то оценка (13) выполняется одновременно с аналогичной оценкой





для разностной схемы (3), (4) при m
₂
= 0. Последняя оценка выражает устойчивость решения разностной задачи по правым частям j и m
₁.

3.3. Неявные схемы. Чисто неявной разностной схемой для уравнения теплопроводности теплопроводности (схемой с опережением) называется разностная схема, использующая шаблон (x_i, t_n), (x_i±
1
, t_n+1), (x_i, t_n+1) и имеющая вид




(12)
Здесь j
ⁿ
_i
= f(x_i, t_n+1) + O(
t
+ h²). Схема имеет первый порядок аппроксимации по t и второй – по h. Решение системы (12) находится, как и в случае явной схемы, по слоям, начиная с n = 1. Однако, теперь, в отличие от явной схемы, для нахождения yⁱ_n+1 по известным y_iⁿ требуется решить систему уравнений




(13)
где g
=
t
/h², F_iⁿ = y_iⁿ +
t
j
_i
ⁿ. Эту систему можно решать методом прогонки, так как условия устойчивости прогонки выполнены.

Для исследования устойчивости разностной схемы (12) будем искать частные решения уравнения





имеющие вид (10). Тогда получим





следовательно, |q| £
1 при любых j
,
t
, h. Таким образом, схема (12) абсолютно устойчива, т.е. устойчива при любых шагах t и h. Абсолютная устойчивость является основным условием неявных схем. Теперь уже не надо брать шаг t слишком малым, можно взять, например, t
= h = 10^-2. Величина шагов сетки t
, h определяются теперь необходимой точностью расчета, а не соображениями устойчивости.

Шеститочечной симметричной схемой называется разностная схема

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

Обобщением трех рассмотренных схем является однопараметрическое семейство схем с весами. Зададим произвольный действительный параметр s и определим разностную схему




(15)
При s = 0 получим отсюда явную схему, при s = 1 – чисто неявную схему и при s = 0,5 – симметричную схему (14). Исследуем погрешность аппроксимации схемы (15) на решении исходной задачи (1) – (3). Представим решение задачи (15) в виде y_iⁿ = u(x_i, t_n) + z_iⁿ, где u(x_i, t_n) – точное решение дифференциальной задачи (1) – (3). Тогда для погрешности получим систему уравнений




(16)
i = 1, 2,…, N – 1,   n = 0, 1,…, K – 1,

z₀ⁿ⁺¹ = z_Nⁿ⁺¹ = 0,   n = 0, 1,…, K – 1, z_i⁰ = 0, i = 0, 1,…, N.
Сеточная функция y_iⁿ, входящая в правую часть уравнения (16) и равная
(17)

называется погрешностью аппроксимации схемы (15) на решении задачи (1) – (3). Получим первые члены разложения функции y_iⁿ по степеням h и t. Будем разлагать все функции, входящие в выражение для y_iⁿ, по формуле Тейлора в точке (x_i, t_n + 0,5t). Учитывая разложения
где






получим





Отсюда, проводя разложение в точке (x_i, t_n+1/2) и обозначая u = u (x_i, t_n+1/2), будем иметь





и, перегруппировывая слагаемые, получим, что





Учитывая уравнение (1) u’’ – u = – f и следствие из него u^IV – u’’ = –f’’, окончательно можно записать, что




(18)
Из формулы (18) можно сделать следующие выводы. Если
то схема (15) имеет второй порядок аппроксимации по t и четвертый – по h. Такая схема называется схемой повышенного порядка аппроксимации. Если






то схема (15) имеет второй порядок аппроксимации по t и по h. При остальных значениях s и при j
_i
ⁿ

º
0 в виде (10), то получим





и |q| £
1 при всех j, если




(19)
Отсюда видно, в частности, что все схемы с s ³ 0,5 абсолютно устойчивы. Схема повышенного порядка аппроксимации (s = s_*) также абсолютно устойчива, что проверяется непосредственно.

При s ¹ 0 разностная схема (15) является неявной схемой. Для нахождения решения y_iⁿ⁺¹ по заданным y_iⁿ требуется решать систему уравнений




(20)
где
Система (20) решается методом прогонки. Условия устойчивости прогонки при s ¹ 0 сводятся к неравенству
|1 + 2
s
g
|
³
2 |
s
|
g

и выполнены при s ³ – 1/(4g). Последнее неравенство следует из условия устойчивости (19) разностной схемы.
3.4. Уравнения с переменными коэффициентами и линейные уравнения. Рассмотрим первую краевую задачу для уравнения теплопроводности с переменными коэффициентами




(21)
где r
(x, t), k(x, t), f(x, t) – достаточно гладкие функции, удовлетворяющие условиям
0 < c₁
£
k(x, t)
£
c₂,  
r
(x, t)
³
c₃ > 0.                                                     (22)




Дифференциальное выражение                                                               при каждом

фиксированном t аппроксимируем в точке (x_i, t) так же, как и в стационарном случае, разностным отношением




(23)
где разностный коэффициент теплопроводности a(x_i, t) должен удовлетворять условиям второго порядка аппроксимации





Наиболее употребительны следующие выражения для a(x_i, t):





Разностная схема с весами для задачи (21) имеет вид




(24)
Здесь в качестве t можно взять любое значение t Î
[t_n, t_n+1], например t = t_n + 0,5t. Если в уравнении (24) t = t_n + 0,5t
,
s
= 0,5, то схема (24) имеет второй порядок аппроксимации по t и по h. При остальных значениях s и t выполняется первый порядок аппроксимации по t и второй – по h.

При исследовании устойчивости разностных схем с переменными коэффициентами иногда применяется принцип замороженных коэффициентов, сводящий задачу к уравнению с постоянными коэффициентами. Рассмотрим явную схему, соответствующую уравнению (24) с s
= 0 и f(x_i, t) º
0, т.е. схему




(25)
Предположим, что коэффициенты r
(x_i, t), a(x_i, t) – постоянные, r
(x_i, t)
º

r
= const,   a(x_i, t)
º
a = const. Тогда уравнение (25) можно записать в виде






или
Из п.2 известно, что последнее уравнение устойчиво при t
’
£
0,5h², т.е. при




(26)
Принцип замороженных коэффициентов утверждает, что схема (25) устойчива, если условие (26) выполнено при всех допустимых значениях a(x_i, t), r
(x_i, t), т.е. если при всех x, t выполнены неравенства




(27)
Если известно, что 0 < c₁ £
a(x_i, t)
£
c₂,  
r
(x_i, t)
³
c₃ > 0, то неравенство (27) будет выполнено при





Строгое обоснование устойчивости схемы (25) будет дано в примере 2 из главы 2.

Если параметр s ³ 0,5, то из принципа замороженных коэффициентов следует абсолютная устойчивость схемы (24).

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




(28)
В случае нелинейных уравнений, когда заранее неизвестны пределы изменения функции k(u), избегают пользоваться явными схемами. Чисто неявная схема, линейная относительно y_iⁿ⁺², i = 1, 2,…, N – 1, имеет вид




(29)
где a_i = 0,5 (k(yⁿ_i) + k(yⁿ_i-1)). Эта схема абсолютно устойчива, имеет первый порядок аппроксимации по t и второй – по h. Решение y_iⁿ⁺¹, i = 1, 2,…, N – 1, находится методом прогонки. Заметим, что схему (29) можно записать в виде





где k_i = k(y_iⁿ).

Часто используется нелинейная схема






(30)
Для реализации этой схемы необходимо применить тот или иной итерационный метод. Например такой:




(31)
Здесь s – номер итерации. Как видим, нелинейные коэффициенты берутся с предыдущей итерации, а в качестве начального приближения для y_iⁿ⁺¹ выбирается y_iⁿ. Это начальное приближение тем лучше, чем меньше шаг t. Число итераций M задается из соображений точности. В задачах с гладкими коэффициентами при k(u) ³
c₁ > 0 часто бывает достаточно провести две – три итерации. Значения y_i^(S+1) на новой итерации находятся из системы (31) методом прогонки. При M = 1 итерационный метод (31) совпадает с разностной схемой (29).

Для приближенного решения нелинейного уравнения (28) применяются также схемы предиктор – корректор второго порядка точности, аналогичные методу Рунге – Кутта для обыкновенных дифференциальных уравнений. Здесь переход со слоя n на слой n+1 осуществляется в два этапа. Приведем пример такой схемы. На первом этапе решается неявная линейная система уравнений





из которой находятся промежуточные значения y_i^n+1/2, i = 0, 1,…, N. Затем на втором этапе используется симметричная шеститочечная схема для уравнения (28), в которой нелинейные коэффициенты a(y), f(y) вычисляются при y = y_i^n+1/2, т.е. схема