Файл: Учебное пособие СанктПетербург 2016.pdf

x
x
x
x
x
l
x
t
t
t
t
t
T
t
i
n
j
m
0 0
0 1
2 2
1
Рис. 2.2
Множество точек плоскости
ω
hτ
= {(x i
, t j
)}, i = 0, ..., n, j = 0, ..., m,
называется сеткой, покрывающей область Ω. Точки плоскости (x i
, t j
) –
это узлы сетки. При фиксированном j совокупность узлов сетки (x i
, t j
),
i = 0, ..., n, называется j-м временным слоем. Будем считать, что функ- ции f (x, t),
ϕ
(x), g
1
(t), g
2
(t) из условия задачи определены только в узлах сетки.
Непрерывной в области Ω функции u(x, t) поставим в соответствие сеточную функцию u(x i
, t j
) = u j
i
Заданное дифференциальное уравнение, начальное и краевые условия запишем в узлах сетки. При этом каждую производную заменим разност- ным отношением, связывающим значения сеточной функции в несколь- ких узлах сетки. Производную по времени
∂u
∂t в узле сетки (x i
, t j
) можно заменить разностным выражением различными способами. Простейшими являются замены разностным отношением “вперед”
u j+1
i
− u j
i
τ
+ O(
τ
)
!
или “назад”
u j−1
i
− u j
i
τ
+ O(
τ
)
!
. В зависимости от способа аппроксимации
67

производных разностными отношениями получаются различные разност- ные схемы для уравнения теплопроводности.
Явная разностная схема.
Для узлов сетки
(x i
, t j
),
i = 1, ..., n − 1, j = 0, ..., m − 1, запишем уравнение теплопроводности,
заменив
∂u
∂t разностным отношением “вперед”
u j+1
i
− u j
i
τ
+ O(
τ
), а
∂
2
u
∂x
2
–
разностным отношением u
j i−1
− 2u j
i
+ u j
i+1
h
2
+ O(h
2
):
u j+1
i
− u j
i
τ
+ O(
τ
) = a
2
u j
i−1
− 2u j
i
+ u j
i+1
h
2
+ O(h
2
) + f (x i
, t j
),
i = 1, ..., n − 1,
j = 0, ..., m − 1.
(2.23)
Для нулевого временного слоя (x i
, t
0
), i = 0, ..., n, запишем начальное усло- вие u
0
i
=
ϕ
(x i
),
i = 0, ..., n.
В граничных узлах сетки (x
0
, t j+1
) и (x n
, t j+1
), j = 0, ..., m − 1, запишем краевые условия, заменив производные соответствующими разностными отношениями:
R
1
u j+1 1
− u j+1 0
h
+ O(h) − S
1
u j+1 0
= g
1
(t j+1
),
R
2
u j+1
n
− u j+1
n−1
h
+ O(h) + S
2
u j+1
m
= g
2
(t j+1
).
Считая, что погрешности аппроксимации O(
τ
), O(h
2
), O(h) в урав- нениях и граничных условиях малы, отбросим их. В результате получим систему уравнений относительно неизвестных ˜
u j
i
(приближенных значений функции u(x, t) в узлах сетки), которая называется явной разностной схе- мой для уравнения теплопроводности:
˜
u j+1
i
− ˜
u j
i
τ
= a
2
˜
u j
i−1
− 2˜
u j
i
+ ˜
u j
i+1
h
2
+f (x i
, t j
), i = 1, ..., n−1, j = 0, ..., m−1,
˜
u
0
i
=
ϕ
(x i
), i = 0, ..., n,
R
1
˜
u j+1 1
− ˜
u j+1 0
h
− S
1
˜
u j+1 0
= g
1
(t j+1
),
R
2
˜
u j+1
n
− ˜
u j+1
n−1
h
+ S
2
˜
u j+1
n
= g
2
(t j+1
),
j = 0, ..., m−1.
На рис. 2.3 указаны узлы сетки, которые используются при аппрокси- мации производных, входящих в уравнение. Такой рисунок принято назы- вать шаблоном схемы. Данная разностная схема является четырехточечной двухслойной. Схема называется явной, поскольку значения искомой функ- ции ˜
u j+1
i на (j + 1)-м временном слое можно последовательно в явном виде
68

определить, если известны значения ˜
u j
i на j-м временном слое. Для этого полученные уравнения преобразуем к виду, удобному для вычислений:

























˜
u j+1
i
= ˜
u j
i
+
τ
a
2
˜
u j
i−1
− 2˜
u j
i
+ ˜
u j
i+1
h
2
+ f (x i
, t j
)
!
,
i = 1, 2, ..., n − 1,
˜
u j+1 0
=
R
1
R
1
+ S
1
h
˜
u j+1 1
−
h
R
1
+ S
1
h g
1
(t j+1
),
˜
u j+1
n
=
R
2
R
2
+ S
2
h
˜
u j+1
n−1
+
h
R
2
+ S
2
h g
2
(t j+1
),
˜
u
0
i
=
ϕ
(x i
),
i = 0, ..., n.
(2.24)
(x
i −1
, t
j
( x
i
, t
j
(x
i
,
t
j
)
)
)
+1
( x
i
,
t
j+1
)
Рис. 2.3
Сначала находятся значения ˜
u
0
i
=
= u
0
i
=
ϕ
(x i
), i = 0, ..., n, на нулевом слое. Затем последовательно определяют- ся значения ˜
u j+1
i на всех временных сло- ях.
Исследуем построенную разностную схему и найдем условия, при выполнении которых решение сеточной задачи схо- дится к решению исходной задачи для уравнения теплопроводности. Исследование сходимости проведем для за- дачи (2.20),(2.21), когда искомая функция u(x, t) удовлетворяет краевым условиям первого рода u(0, t) = g
1
(t),
u(l, t) = g
2
(t).
Разностная схема в этом случае примет вид





















˜
u j+1 0
= g
1
(t j+1
),
˜
u j+1
i
= ˜
u j
i
+
τ
a
2
˜
u j
i−1
− 2˜
u j
i
+ ˜
u j
i+1
h
2
+ f (x i
, t j
)
!
,
i = 1, 2, ..., n − 1,
˜
u j+1
n
= g
2
(t j+1
),
˜
u
0
i
=
ϕ
(x i
),
i = 0, ..., n.
(2.25)
Сеточной нормой функции называется максимальное значение модуля
69

функции на той сетке, на которой она задана:
kuk h,τ
= max i=0,...,n;
j=0,...,m
|u(x i
, t j
)|,
k
ϕ
k h
= max i=0,...,n
|
ϕ
(x i
)|.
Разностная схема называется сходящейся к решению дифференциаль- ной задачи в сеточной норме, если ku − ˜
uk h,τ
→ 0 при h → 0 и
τ
→ 0.
Сходимость разностной схемы связана с двумя важными понятиями:
устойчивостью разностной схемы и точностью аппроксимации (см. 1.9).
Разностная схема называется устойчивой по начальным данным и пра- вой части, если норма ее решения не превосходит суммы норм заданных в задаче функций, умноженных на число, не зависящее от шагов сетки:
k˜
uk h,τ
≤ C(k
ϕ
k h
+ kf k h,τ
+ kg
1
k
τ
+ kg
2
k
τ
).
(2.26)
Это условие показывает, что при малом изменении исходных данных или при уменьшении шагов сетки решение разностной задачи будет меняться мало. Выполнение этого условия означает непрерывную зависимость реше- ния, полученного с помощью разностной схемы, от входных данных. При решении разностных уравнений, например (2.25), естественно использовать вычислительные средства. А значит, расчеты будут выполняться не точ- но, а с погрешностью, связанной с конечным представлением чисел (длина разрядной сетки ограничена).
Если малые ошибки округления чисел приводят к большим искаже- ниям решения, то это явление называется вычислительной неустойчиво- стью разностной схемы. Разностная схема, обладающая таким свойством,
не пригодна для вычислений. При исследовании разностной схемы обычно используется следующий подход. Считается, что ошибки, возникающие при выполнении арифметических операций, можно рассматривать как возму- щения, которые вносятся в начальные данные и в правую часть уравнения.
От разностной схемы требуют, чтобы для нее выполнялось условие устой- чивости по начальным данным и правой части. Выполнение этого условия обеспечивает ее вычислительную устойчивость, при этом разностную схему называют устойчивой, а условие (2.26) – условием устойчивости разностной схемы.
Докажем устойчивость явной разностной схемы при условии, что функ- ция u(x, t) удовлетворяет однородным краевым условиям первого рода:
u(0, t) = 0,
u(l, t) = 0
(g
1
(t) = 0, g
2
(t) = 0).
Утверждение 2.3. Если
τ
≤
h
2 2a
2
, то явная разностная схема (2.25)
устойчива по начальным данным и правой части, т. е.
k˜
uk h,τ
≤ C(k
ϕ
k h
+ kf k h,τ
),
70

где C не зависит от h и
τ
Доказательство. Обозначим k˜
u j
k h
, kf j
k h
наибольшие по модулю зна- чения сеточных функций ˜
u и f на j-м временном слое:
k˜
u j
k h
= max i=0,...,n
|˜
u j
i
|,
kf j
k h
= max i=0,...,n
|f (x i
, t j
)|.
Из равенств (2.24) следует, что
|˜
u j+1
i
| =
τ
a
2
h
2
˜
u j
i−1
+

1 −
2
τ
a
2
h
2

˜
u j
i
+
τ
a
2
h
2
˜
u j
i+1
+
τ
f (x i
, t j
)
Тогда
|˜
u j+1
i
| ≤
τ
a
2
h
2
k˜
u j
k h
+
1 −
2
τ
a
2
h
2
k˜
u j
k h
+
τ
a
2
h
2
k˜
u j
k h
+
τ
kf j
k h
Так как ˜
u j+1 0
= ˜
u j+1
n
= 0, то неравенство справедливо для всех i = 0, ..., n.
Пусть 1−
2
τ
a
2
h
2
≥ 0 (т. е.
τ
≤
h
2 2a
2
), тогда, заменив
1 −
2
τ
a
2
h
2
на 1−
2
τ
a
2
h
2
,
получим
|˜
u j+1
i
| ≤ k˜
u j
k h
+
τ
kf j
k h
при всех i = 0, ..., n. Значит,
k˜
u j+1
k h
≤ k˜
u j
k h
+
τ
kf j
k h
Используя оценку k˜
u j
k h
≤ k˜
u j−1
k h
+
τ
kf j−1
k h
и т. д., придем к цепочке неравенств k˜
u j+1
k h
≤ k˜
u j−1
k h
+
τ
kf j−1
k h
+ kf j
k h
≤
≤ k˜
u j−1
k h
+
τ
kf j−2
k h
+ kf j−1
k h
+ kf j
k h
≤ ... ≤
≤ k˜
u
0
k h
+
τ
kf
0
k h
+ kf
1
k h
+ ... + kf j
k h
.
Поскольку kf j
k h
≤ kf k h,τ
для всех j = 0, 1, ..., m − 1 и k˜
u
0
k h
= k
ϕ
k h
,
получим k˜
u j+1
k h
≤ k
ϕ
k h
+
τ
(j + 1)kf k h,τ
,
где
τ
(j + 1) ≤ T (j = 0, ..., m − 1). Тогда для всех j = 0, ..., m − 1
k˜
u j+1
k h
≤ k
ϕ
k h
+ T kf k h,τ
,
а значит,
k˜
uk h,τ
≤ k
ϕ
k h
+ T kf k h,τ
Явная разностная схема будет устойчивой только при выполнении условия
τ
≤
h
2 2a
2
. В противном случае явная разностная схема не пригодна
71

для вычислений. В случае a
2
= 1 условие устойчивости записывается в виде
τ
≤
h
2 2
Если в систему разностных уравнений вместо значений сеточной функ- ции ˜
u j
i подставить значения точного решения дифференциальной задачи u(x i
, t j
), то для сохранения равенства в правую часть разностных урав- нений необходимо ввести дополнительное слагаемое. Если с уменьшением шагов сетки это слагаемое меняется как O(h k
) + O(
τ
m
), то говорят, что разностная схема аппроксимирует задачу с погрешностью O(h k
) + O(
τ
m
).
Утверждение 2.4. Явная разностная схема (2.25) аппроксимирует задачу для уравнения теплопроводности (2.20), (2.21) с краевыми услови- ями u(0, t) = g
1
(t), u(l, t) = g
2
(t) с погрешностью O(h
2
) + O(
τ
).
Доказательство. Запишем разностную схему для указанных крае- вых условий:
˜
u j+1
i
− ˜
u j
i
τ
= a
2
˜
u j
i−1
− 2˜
u j
i
+ ˜
u j
i+1
h
2
+ f (x i
, t j
),
i = 1, ..., n − 1,
j = 1, ..., m − 1,
˜
u
0
i
=
ϕ
(x i
),
i = 0, ..., n,
˜
u j+1 0
= g
1
(t j+1
),
˜
u j+1
n
= g
2
(t j+1
),
j = 0, ..., m − 1.
Функции u(x, t) и ˜
u j
i в граничных узлах сетки (при t = 0, x = 0, x = l)
удовлетворяют одинаковым условиям. Поэтому погрешность аппроксима- ции в данной задаче связана только с заменой производных разностными отношениями в уравнении теплопроводности. Из равенств (2.23) следует,
что погрешность аппроксимации для рассматриваемой разностной схемы есть O(h
2
) + O(
τ
).
Заметим, что разностная схема (2.24) аппроксимирует начально-крае- вую задачу для уравнения теплопроводности (2.20)–(2.22) с краевыми усло- виями третьего рода с погрешностью O(h) + O(
τ
). Это связано с тем, что появляется погрешность O(h) при аппроксимации производных в гранич- ных условиях.
Докажем теперь сходимость явной разностной схемы (2.25).
Утверждение 2.5. При
τ
≤
h
2 2a
2
явная разностная схема (2.25) схо- дится в сеточной норме к решению задачи (2.20), (2.21) с краевыми усло- виями u(0, t) = g
1
(t), u(l, t) = g
2
(t).
Доказательство. Доказательство проведем, используя точность ап- проксимации и устойчивость разностной схемы. Обозначим w = u − ˜
u раз- ность между точным и приближенным решениями задачи. Очевидно, что
72

во внутренних узлах сетки функция w будет удовлетворять уравнениям w
j+1
i
− w j
i
τ
= a
2
w j
i−1
− 2w j
i
+ w j
i+1
h
2
+ O(h
2
) + O(
τ
),
i = 1, ..., n − 1,
j = 0, ..., m − 1;
на нулевом временном слое – начальному условию w
0
i
= 0,
i = 0, ..., n;
в граничных узлах (x
0
, t j+1
), (x n
, t j+1
) – краевым условиям w
j+1 0
= 0,
w j+1
n
= 0,
j = 0, ..., m − 1.
В данном случае f (x, t) = O(h
2
)+O(
τ
) и
ϕ
(x) = 0. Применив утверждение
2.3 к разностной схеме для функции w, при выполнении условия
τ
≤
h
2 2a
2
получим kwk h,τ
≤ T kO(h
2
) + O(
τ
)k h,τ
≤ T kO(h
2
)k h,τ
+ kO(
τ
)k h,τ
.
Поскольку справедливы оценки kO(h
2
)k h,τ
≤ M
1
h
2
,
kO(
τ
)k h,τ
≤ M
2
τ
(см. 1.9), то kwk h,τ
≤ M (h
2
+
τ
),
где M = T max(M
1
, M
2
), а значит, kwk h,τ
→ 0 при h → 0 и
τ
→ 0.
Сходимость доказана только при условии, что шаги сетки по времени
τ
достаточно малы (
τ
≤
h
2 2a
2
). Если условие устойчивости нарушается, то явная схема не будет сходиться к решению поставленной задачи. От слоя к слою значения ˜
u j+1
i начинают неограниченно расти. Схема “разваливается”.
Неявная разностная схема. Построим для задачи (2.20)–(2.22) неяв- ную разностную схему. Для внутренних узлов сетки запишем уравнение теплопроводности, заменив
∂u(x i
, t j
)
∂t разностным отношением “назад”
u j
i
− u j−1
i
τ
+ O(
τ
), а
∂
2
u(x i
, t j
)
h
2
τ
∂x
2
, как и для явной схемы, разностным от- ношением u
j i−1
− 2u j
i
+ u j
i+1
h
2
+ O(h
2
):
u j
i
− u j−1
i
τ
+ O(
τ
) = a
2
u j
i−1
− 2u j
i
+ u j
i+1
h
2
+ O(h
2
) + f (x i
, t j
),
i = 1, ..., n − 1,
j = 1, ..., m.
Для слоя (x i
, t
0
) запишем начальное условие:
u
0
i
=
ϕ
(x i
),
i = 0, ..., n,
73

а для граничных узлов – краевые условия:
R
1
u j
1
− u j
0
h
+ O(h) − S
1
u j
0
= g
1
(t j
),
R
2
u j
n
− u j
n−1
h
+ O(h) + S
2
u j
n
= g
2
(t j
),
j = 1, ..., m.
Отбросив в уравнениях неизвестные погрешности O(h), O(
τ
), O(h
2
)
получим систему относительно приближенных значений ˜
u j
i искомой функ- ции в узлах сетки:
˜
u j
i
− ˜
u j−1
i
τ
= a
2
˜
u j
i−1
− 2˜
u j
i
+ ˜
u j
i+1
h
2
+ f (x i
, t j
),
i = 1, ..., n − 1,
j = 1, ..., m;
начальное условие:
˜
u
0
i
= u
0
i
=
ϕ
(x i
),
i = 0, ..., n;
краевые условия:
R
1
˜
u j
1
− ˜
u j
0
h
− S
1
˜
u j
0
= g
1
(t j
),
R
2
˜
u j
n
− ˜
u j
n−1
h
+ S
2
˜
u j
n
= g
2
(t j
),
j = 1, ..., m.
(x
i −1
, t
j
( x
i
, t
j
(x
i +1
,
t
j
)
(x
i
,
t
j −1
)
)
)
Рис. 2.4
Построенная разностная схема, как и явная, является четырехточечной двух- слойной разностной схемой. Ее шаблон изображен на рис. 2.4.
Схема называется неявной, посколь- ку значения функции ˜
u j
i в данном случае не могут быть вычислены последователь- но через значения функции на (j − 1)-м временном слое. Для их нахождения на каждом временном слое следует решать систему уравнений











−(R
1
+ S
1
h)˜
u j
0
+ R
1
˜
u j
1
= hg
1
(t j
),
˜
u j
i−1
− (2 +
h
2
τ
a
2
)˜
u j
i
+ ˜
u j
i+1
= −
h
2
τ
a
2
˜
u j−1
i
− h
2
f (x i
, t j
),
i = 1, ..., n − 1,
−R
2
˜
u j
n−1
+ (R
2
+ S
2
h)˜
u j
n
= hg
2
(t j
).
Значение ˜
u
0
i для нулевого временного слоя определяется из начального условия ˜
u
0
i
=
ϕ
(x i
) (i = 0, ..., n).
74

Систему уравнений для j-го временного слоя удобно записать в мат- ричном виде









B
0
C
0 0
0 0
A
1
B
1
C
1 0
0 0
A
2
B
2
C
2 0
0 0
... A
n−1
B
n−1
C
n−1 0
0 0
A
n
B
n


















˜
u j
0
˜
u j
1
˜
u j
2
˜
u j
n−1
˜
u j
n









=








F
0
F
1
F
2
F
n−1
F
n








Здесь
B
0
= R
1
+ S
1
h, C
0
= R
1
, F
0
= hg
1
(t j
),
A
i
= 1,
B
i
= −2 +
h
2
τ
a
2
, C
i
= 1,
F
i
= −
h
2
τ
a
2
˜
u j−1
i
− h
2
f (x i
, t j
),
i = 1, ..., n − 1
A
n
= −R
2
, B
n
= R
2
+ S
2
h,
F
n
= hg(t j
).
Кратко эту систему можно записать в виде
A˜
u j
=
F
j
,
где A – матрица коэффициентов; ˜
u j
– вектор значений функции ˜
u;
F
j
–
вектор значений правой части уравнений для j-го временного слоя. Вектор

F
j меняется от слоя к слою. Матрица коэффициентов системы A трехдиа- гональная. Для строк этой матрицы выполняется условие строгого диаго- нального преобладания |B
i
| > |A
i
| + |C
i
| (i = 0, ..., n). Это условие обес- печивает однозначную разрешимость системы (см. 1.9). Кроме того, для решения системы уравнений с подобной матрицей коэффициентов можно применить метод Гаусса без перестановки строк и столбцов матрицы. Такой метод решения системы уравнений с трехдиагональной матрицей коэффи- циентов называется методом прогонки. Сначала выполняется прямой ход метода прогонки: вычисляются прогоночные коэффициенты
˜
C
0
=
C
0
B
0
,
˜
F
0
=
F
0
B
0
,
˜
C
i
=
C
i
B
i
− A
i
˜
C
i−1
,
i = 1, ..., n − 1,
˜
F
i
=
F
i
− A
i
˜
F
i−1
B
i
− A
i
˜
C
i−1
,
i = 1, ..., n.
При этом матрица коэффициентов системы преобразуется к верхней тре- угольной матрице, на главной диагонали которой стоят единицы; ˜
C
i
– нену- левые элементы матрицы, стоящие выше главной диагонали; ˜
F
i
– числа,
стоящие в правой части системы уравнений.
75

Затем выполняется обратный ход (решается система уравнений с верх- ней треугольной матрицей коэффициентов):
˜
u j
n
= ˜
F
n
,
˜
u j
i
= ˜
F
i
− ˜
C
i
˜
u j
i+1
,
i = n − 1, ..., 0.
Метод прогонки можно применять только в том случае, когда в фор- мулах для вычисления прогоночных коэффициентов ˜
C
i
, ˜
F
i знаменатели не обращаются в нуль. Это требование будет выполняться, если для матри- цы коэффициентов справедливо условие строгого диагонального преобла- дания. Таким образом, при применении неявной разностной схемы для на- хождения приближенного решения поставленной задачи значения ˜
u j
i будут определяться последовательно по временным слоям. Для каждого времен- ного слоя необходимо будет решать систему уравнений с трехдиагональной матрицей коэффициентов.
Преимущество неявной схемы заключается в том, что в отличие от явной разностной схемы построенная неявная схема является безусловно устойчивой, т. е. она устойчива при любых
τ
и h. Устойчивость и сходи- мость неявной разностной схемы рассмотрим для задачи (2.20),(2.21) при условии,что функция u(x, t) удовлетворяет краевым условиям первого ро- да:
u(0, t) = g
1
(t),
u(l, t) = g
2
(t).
Неявная разностная схема в этом случае имеет вид



















˜
u j
0
= g
1
(t j
),
˜
u j
i−1
− (2 +
h
2
τ
a
2
)˜
u j
i
+ ˜
u j
i+1
= −
h
2
τ
a
2
˜
u j−1
i
− h
2
f (x i
, t j
),
i = 1, ..., n − 1,
˜
u j
n
= g
2
(t j
),
˜
u
0
i
=
ϕ
(x i
),
i = 0, ..., n.
(2.27)
Пусть функция u(x, t) удовлетворяет однородным краевым условиям g
1
(t) = 0 и g
2
(t) = 0. Справедливо следующее утверждение.
Утверждение 2.6. Неявная разностная схема (2.27) устойчива по начальным данным и правой части, т. е.
k˜
uk h,τ
≤ C(k
ϕ
k h
+ kf k h,τ
),
где C не зависит от h и
τ
76

Доказательство. Для всех внутренних узлов j-го временного слоя выполняется равенство
˜
u j
i−1
−

2 +
h
2 2
τ

˜
u j
i
+ ˜
u j
i+1
= −
h
2
τ
˜
u j−1
i
− h
2
f (x i
, t j
)
или

2 +
h
2
τ

˜
u j
i
=
˜
u j
i−1
+ ˜
u j
i+1
+
h
2
τ
˜
u j−1
i
+ h
2
f (x i
, t j
)
Тогда справедлива оценка

2 +
h
2
τ

|˜
u j
i
| ≤ |˜
u j
i−1
| + |˜
u j
i+1
| +
h
2
τ
|˜
u j−1
i
| + h
2
|f (x i
, t j
)|.
Пусть k˜
u j
k h
и kf j
k h
– наибольшие по модулю значения сеточных функций
˜
u и f на j-м временном слое, тогда

2 +
h
2
τ

|˜
u j
i
| ≤ 2k˜
u j
i k
h
+
h
2
τ
k˜
u j−1
i k
h
+ h
2
kf j
k h
Функция ˜
u удовлетворяет однородным краевым условиям ˜
u j
0
= ˜
u j
n
= 0,
поэтому последнее неравенство выполняется для всех i = 0, ..., n, и тогда

2 +
h
2
τ

k˜
u j
k h
≤ 2k˜
u j
i k
h
+
h
2
τ
k˜
u j−1
i k
h
+ h
2
kf j
k h
Преобразуем это неравенство:
k˜
u j
k h
≤ k˜
u j−1
k h
+
τ
kf j
k h
Последовательно применяя его к (j−1)-му, (j−2)-му и т. д. слоям, получим:
k˜
u j
k h
≤ k˜
u j−2
k h
+
τ
kf j−1
k h
+ kf j
k h
≤ ... ≤
≤ k˜
u
0
k h
+
τ
kf
1
k h
+ kf
2
k h
+ ... + kf j
k h
.
Учитывая равенство k˜
u
0
k h
= k
ϕ
k h
и оценки kf j
k h
≤ kf k h,τ
,
τ
j ≤ T для j = 1, ..., m, получим неравенство k˜
u j
k h
≤ k
ϕ
k h
+ T kf k h,τ
,
которое будет выполняться для всех j = 1, ..., m, а значит,
k˜
uk h,τ
≤ k
ϕ
k h
+ T kf k h,τ
Очевидно, что неявная разностная схема аппроксимирует задачу для уравнения теплопроводности с той же погрешностью, что и соответствую- щая явная схема. Если функция u(x, t) удовлетворяет краевым условиям
77

первого рода u(0, t) = g
1
(t) и u(l, t) = g
2
(t), то погрешность аппроксимации разностной схемы (2.27) будет O(h
2
) + O(
τ
).
Из устойчивости и аппроксимации следует сходимость неявной раз- ностной схемы. Доказательство сходимости абсолютно аналогично доказа- тельству сходимости явной разностной схемы.
Схема Кранка–Николсона. Неявная разностная схема хотя и не имеет ограничения на шаги по времени
τ
, но поскольку ее погрешность аппроксимации равна O(
τ
) + O(h
2
), то разумно выбрать шаги
τ
достаточ- но малыми (порядка O(h
2
)). Более высокую точность аппроксимации по времени имеет схема Кранка–Николсона. Это шеститочечная двухслойная неявная разностная схема (рис. 2.5).

1   2   3   4   5   6   7   8   9   10

j
i
−1
+1
( x
i
j
)
( x
i
,
t
j
)
(x
i +1
, t
j
)
( x
i −1
, t
j −1
)
( x
i
, t
)
(x
,
t
j −1
−1
, t
)
Рис. 2.5
Построим эту схему для уравнения теплопроводности (2.20) при условии, что функция u(x, t) удовлетворяет начально- му условию (2.21) и краевым условиям первого рода u(0, t) = g
1
(t), u(l, t) = g
2
(t).
Между узлами сетки (x i
, t j
) и (x i
, t j+1
)
выбирается точка (x i
, t j
+
τ
2
) (рис. 2.5).
Для этой точки записывается уравнение теплопроводности. При этом
∂u
∂t заменя- ется центральным разностным отношением u
j i
− u j−1
i
2(
τ
/2)
+O


τ
2

2

(см. 1.9),
а
∂
2
u
∂x
2
в этой точке находится как полусумма u
j i−1
− 2u j
i
+ u j
i+1
h
2
+ O(h
2
) и u
j−1
i−1
− 2u j−1
i
+ u j−1
i+1
h
2
+ O(h
2
). После отбрасывания погрешностей аппрокси- мации производных получается неявная разностная схема Кранка–Никол- сона:
˜
u j
i
− ˜
u j−1
i
τ
=
=
˜
u j
i−1
− 2˜
u j
i
+ ˜
u j
i+1
+ ˜
u j−1
i−1
− 2˜
u j−1
i
+ ˜
u j−1
i+1 2h
2
+ f (x i
, t j−1
+
τ
2
),
i = 1, ..., n − 1,
j = 1, ..., m.
Начальное условие:
˜
u
0
i
=
ϕ
(x i
),
i = 0, ..., n.
Граничные условия:
˜
u j
0
= g
1
(t j
),
˜
u j
n
= g
2
(t j
).
Погрешность аппроксимации такой разностной схемы равна O(h
2
) + O(
τ
2
).
78

Иногда значение f (x i
, t j−1
+
τ
2
) заменяют полусуммой f (x i
, t j
) + f (x i
, t j−1
)
2
Тогда для внутренних узлов сетки (x i
, t j
) схема Кранка–Николсона полу- чается как полусумма неявной схемы, записанной для узлов (x i
, t j
), и явной схемы, построенной для узлов (x i
, t j−1
). Погрешность аппроксимации раз- ностной схемы в этом случае также O(h
2
) + O(
τ
2
).
3. УРАВНЕНИЕ КОЛЕБАНИЙ СТРУНЫ.
ВОЛНОВОЕ УРАВНЕНИЕ
3.1. Уравнение колебаний струны
Рассмотрим туго натянутую струну, т. е. тонкую упругую (подчиняю- щуюся закону Гука) гибкую нить. Ось 0x направим вдоль струны. Предпо- ложим, что струну вывели из положения равновесия, например оттянули или ударили по ней. Будем считать, что возникающие при этом поперечные колебания струны происходят в одной плоскости 0xu и все точки струны движутся перпендикулярно оси 0x. Обозначим через u(x, t) отклонение то- чек струны от прямолинейной формы. Функция u(x, t) описывает попереч- ные колебания струны при t ≥ 0. График функции u(x, t) (рис. 3.1) зависит от момента времени t и представляет собой форму колеблющейся струны в данный момент времени t.
x
u
B
A
α
T
( x , t)
(
∆
,
x
t
)
+
x
α
∆
x
T
+
F
0
x
x
1 2
Рис. 3.1
При фиксированном x функция u(x, t) описывает закон движения, про- изводная
∂u(x, t)
∂t
– скорость движения, а вторая производная
∂
2
u(x, t)
∂t
2
–
ускорение выбранной точки струны.
Для вывода уравнения поперечных колебаний струны сделаем следу- ющие предположения:
1. Будем считать, что струна совершает только малые колебания, т. е.
ее форма незначительно отличается от прямой u = 0. Касательная в каж- дой точке x, проведенная к графику u(x, t), почти параллельна оси абсцисс.
Пусть
α
(x, t) – угол, который образует касательная с осью 0x. Условие ма-
79

лости колебаний означает, что величиной
α
2
можно пренебречь. Разложим функции sin
α
, tg
α
, cos
α
по формуле Тейлора первого порядка:
sin
α
=
α
+ o(
α
),
tg
α
=
α
+ o(
α
),
cos
α
= 1 + o(
α
).
При сделанном предположении
α
≈ sin
α
≈ tg
α
=
∂u
∂x и cos
α
≈ 1.
2. На любой выбранный участок струны действуют упругие силы на- тяжения. Выделим на струне произвольный малый участок [g
AB](рис. 3.1).
Его проекция на ось 0x есть [x, x + ∆x]. Будем считать, что упругие силы,
приложенные к концам выбранного участка, направлены по касательной к графику u(x, t), их модули |
T (x, t)| = |
T (x + ∆x, t)| = T и не зависят от x и t.
При изучении малых поперечных колебаний струны такое предполо- жение допустимо. Обозначим
α
1
и
α
2
углы, которые образуют касательные,
проведенные в точках A и B соответственно к графику струны в момент времени t, с осью 0x. По условию все точки струны движутся параллельно оси 0u. Значит, силы, действующие на участок [g
AB], таковы, что сумма их проекций на ось 0x должна быть равна нулю:
−|
T (x, t)| cos
α
1
+ |
T (x + ∆x, t)| cos
α
2
= 0.
Поскольку cos
α
1
≈ cos
α
2
≈ 1, то можно считать, что |
T (x, t)| = T = const.
3. В положении равновесия масса участка струны [x, x + ∆x] равна
ρ
∆x, где
ρ
– линейная плотность струны. Считаем, что в процессе колеба- ний масса этого участка не меняется и масса [g
AB] тоже равна
ρ
∆x.
Такое предположение также допустимо, так как длина участка струны
[g
AB] в момент времени t
|g
AB| =
x+∆x
Z
x s
1 +
∂u
∂x

2
dx.
При этом
∂u
∂x
≈
α
и величиной
α
2
можно пренебречь. Значит, |g
AB| ≈ ∆x.
4. Предположим,что на струну в плоскости колебаний действуют непре- рывно распределенные внешние силы, перпендикулярные оси 0x. Плот- ность распределения этих сил, рассчитанную на единицу длины, обозначим g(x, t). Если сила
F , приложенная к участку длины ∆x, направлена вверх,
то
F (x, t) = g(x, t)∆xj, если – вниз, то
F (x, t) = −g(x, t)∆xj.
Рассмотрим выбранный малый участок [g
AB]. Удалим части струны,
расположенные справа и слева от него. Воздействие отброшенных частей струны заменим соответствующими силами натяжения. Тогда выделенный
80

участок можно рассматривать как материальную точку, которая находится под воздействием трех сил
T (x, t),
T (x + ∆x, t) и
F (x, t) = g(x, t)∆xj.
Проецируя на ось 0u действующие на участок [g
AB] силы и применяя закон Ньютона, получим
ρ
∆x
∂
2
u
∂t
2
= T sin
α
2
− T sin
α
1
+ g(x, t)∆x.
Согласно сделанным предположениям углы
α
1
и
α
2
малы и sin
α
1
≈ tg
α
1
=
∂u(x, t)
∂x
,
sin
α
2
≈ tg
α
2
=
∂u(x + ∆x, t)
∂x
Тогда, если использовать формулу Тейлора
∂u(x + ∆x, t)
∂x
=
∂u(x, t)
∂x
+
∂
2
u(x, t)
∂x
2
∆x + o(∆x)
и отбросить бесконечно малую o(∆x), придем к равенству
ρ
∆x
∂
2
u
∂t
2
= T
∂
2
u(x, t)
∂x
2
∆x + g(x, t)∆x.
Сокращая на ∆x, получим уравнение малых вынужденных поперечных колебаний струны
ρ
∂
2
u
∂t
2
= T
∂
2
u
∂x
2
+ g(x, t).
Обычно это уравнение записывают в виде
∂
2
u
∂t
2
= a
2
∂
2
u
∂x
2
+ f (x, t),
(3.1)
где a
2
=
T
ρ
; f (x, t) =
g(x, t)
ρ
. Полученное уравнение называется также одномерным волновым уравнением, или уравнением Даламбера.
Если f (x, t) = 0, уравнение (3.1) называется однородным. Оно описы- вает свободные колебания струны без воздействия внешних сил.
В случае f (x, t) 6= 0 уравнение называется неоднородным и описыва- ет вынужденные колебания струны. В частности, вынужденные колебания могут происходить под действием силы тяжести. Если же натяжение стру- ны T велико и на нее действует только сила тяжести, действием последней обычно пренебрегают и считают, что струна совершает свободные колеба- ния.
81

Отметим, что многие физические задачи приводят к полученному вол- новому уравнению. Точно так же выглядит уравнение продольных колеба- ний тонкого упругого стержня. Функция u(x, t) при изучении таких коле- баний описывает продольные смещения точек стержня, имеющих в поло- жении равновесия абсциссу x. В уравнении (3.1) в этом случае a
2
=
E
ρ
,
f (x, t) =
g(x, t)
ρ
, где E – модуль Юнга; g(x, t) – функция, описывающая плотность сил, действующих вдоль оси стержня. Вывод уравнения подроб- но описан в [9]. Кроме того, уравнение (3.1) совпадает с уравнением кру- тильных колебаний вала. Подобные уравнения появляются при изучении электрических колебаний и во многих других случаях. Меняется только физический смысл функций и коэффициентов, входящих в уравнение.
В случае двух или трех пространственных переменных волновое урав- нение имеет вид
∂
2
u
∂t
2
= a
2
∆u + f,
(3.2)
где ∆ – оператор Лапласа; u и f – функции пространственных переменных и времени.
Двумерное волновое уравнение получается, например, при изучении поперечных колебаний мембраны. Если мембрана прямоугольная, то урав- нение записывается в декартовой системе координат:
∂
2
u
∂t
2
= a
2
∂
2
u
∂x
2
+
∂
2
u
∂y
2

+ f (x, y, t).
При описании колебаний круглой мембраны переходят к полярной си- стеме координат и уравнение принимает вид
∂
2
u
∂t
2
= a
2
1
ρ
∂
∂
ρ

ρ
∂u
∂
ρ

+
1
ρ
2
∂
2
u
∂
ϕ
2

+ f (
ρ
,
ϕ
, t).
Трехмерное волновое уравнение возникает, например, в задачах, свя- занных с изучением колебаний газа, находящегося в ограниченном объеме,
в задачах распространения акустических волн и во многих других. В декар- товой системе координат это уравнение записывается следующим образом:
∂
2
u
∂t
2
= a
2
∂
2
u
∂x
2
+
∂
2
u
∂y
2
+
∂
2
u
∂z
2

+ f (x, y, z, t).
3.2. Постановка начальных и краевых условий
Волновое уравнение имеет бесчисленное множество решений. Для од- нозначного описания колебательного процесса на функцию u следует на-
82

ложить дополнительные условия, вытекающие из физического смысла за- дачи. Дополнительные условия – это начальное и краевые условия для функции u.
Покажем как ставятся эти условия при изучении поперечных колеба- ний струны.
Начальные условия обычно задают в момент времени t = 0. Условия описывают начальное положение точек струны и ее начальную скорость:
u(x, 0) =
ϕ
(x),
∂u(x, 0)
∂t
=
ψ
(x),
где
ϕ
(x) и
ψ
(x) – заданные функции.
Граничные условия показывают, что происходит на концах струны в течение всего времени колебаний. Предположим, что струна имеет конеч- ную длину l (0 ≤ x ≤ l).
1. Если концевые точки струны движутся по определенному закону,
то краевые условия имеют вид u(0, t) =
ω
1
(t),
u(l, t) =
ω
2
(t).
В частности, если концы струны закреплены, то для любого момента вре- мени t u(0, t) = 0,
u(l, t) = 0.
2. Граничные условия можно задать следующим образом:
∂u(0, t)
∂x
= h
1
(t),
∂u(l, t)
∂x
= h
2
(t).
Это соответствует случаю, когда известен закон изменения касательных в концевых точках струны. Если h
1
(t) = 0 и h
2
(t) = 0, то в концевых точках для любого t ≥ 0 струна имеет касательные, параллельные оси 0x.
3. В случае упругого закрепления концов стержня краевые условия записывают в виде
∂u(0, t)
∂x
= h
1
(u(0, t) − u
0
(t)) ,
∂u(l, t)
∂x
= h
2
(u(l, t) − u l
(t)) ,
где u
0
(t) и u l
(t) – заданные функции.
Все рассмотренные краевые условия являются линейными и описыва- ются уравнениями
R
1
∂u(0, t)
∂x
− S
1
u(0, t) = g
1
(t),
R
2
∂u(l, t)
∂x
+ S
2
u(l, t) = g
2
(t).
При R
1
= 0, R
2
= 0 – это краевые условия первого рода, или условия
Дирихле. При S
1
= 0, S
2
= 0 – условия второго рода, или условия Неймана,
а при R
1
S
1 6= 0, R
2
S
2 6= 0 – условия третьего рода.
83

Заметим, что при решении физических задач на концах x = 0 и x = l могут задаваться краевые условия разного рода.
4. Если струна бесконечная, то для функции u(x, t) задают только на- чальные условия, краевые условия на нее обычно не накладываются, но при этом предполагается, что на бесконечности функция u(x, t) ограниче- на.
Аналогично задаются начальные и краевые условия для уравнения
(3.2). Пусть в области Ω с границей Γ функция u(M, t) (M ∈ Ω и t > 0)
удовлетворяет волновому уравнению (3.2). К этому уравнению добавляют- ся начальные условия u(M, 0) =
ϕ
(M ),
∂u(M, 0)
∂t
=
ψ
(M ).
На границе Γ ставится одно из трех условий:
u
Γ
=
µ
(M, t),
∂u
∂
n
Γ
=
ν
(M, t),
∂u
∂
n
+ hu
Γ
=
χ
(M, t),
где
µ
,
ν
,
χ
, h – заданные и непрерывные на границе Γ функции; M – точка границы области Γ;
∂u
∂
n
– производная по направлению внешней нормали к границе.
На разных частях границы Γ могут быть заданы краевые условия раз- ного рода. В таком случае краевые условия называются смешанными.
3.3. Колебания бесконечной струны. Метод Даламбера
Рассмотрим бесконечную струну, которую в начальный момент вре- мени вывели из положения равновесия. Будем считать, что внешние силы отсутствуют и струна совершает свободные колебания. Для нахождения поперечных колебаний такой струны следует решить волновое уравнение
∂
2
u
∂t
2
= a
2
∂
2
u
∂x
2
,
−∞ < x < +∞,
t > 0,
(3.3)
при начальных условиях u(x, 0) =
ϕ
(x),
∂u(x, 0)
∂t
=
ψ
(x),
где функции
ϕ
(x),
ψ
(x) описывают начальное положение и начальную ско- рость точек струны.
Сформулированная задача называется задачей Коши для бесконечной струны. Решим эту задачу методом, который называется методом Далам- бера, или методом бегущих волн. Покажем сначала, что общее решение уравнения (3.3) имеет вид u(x, t) = P (x − at) + Q(x + at),
(3.4)
84

где P и Q – произвольные дважды дифференцируемые функции. Введем новые переменные
ξ
= x − at и
η
= x + at и запишем волновое урав- нение (3.3) в новых переменных. Используя правило дифференцирования сложной функции, выразим производные функции u(x, t) по x и t через производные по
ξ
и
η
:
∂u
∂x
=
∂u
∂
ξ
+
∂u
∂
η
,
∂
2
u
∂x
2
=
∂
2
u
∂
ξ
2
+ 2
∂
2
u
∂
η
∂
ξ
+
∂
2
u
∂
η
2
,
∂u
∂t
= a
∂u
∂
η
−
∂u
∂
ξ

,
∂
2
u
∂t
2
= a
2
∂
2
u
∂
η
2
− 2
∂
2
u
∂
η
∂
ξ
+
∂
2
u
∂
ξ
2

Подставим выражения для
∂
2
u
∂x
2
и
∂
2
u
∂t
2
в уравнение (3.3), приведем подобные слагаемые и сократим на (−4a
2
). Тогда в новых переменных уравнение запишется в виде
∂
2
u
∂
η
∂
ξ
= 0,
или
∂
∂
η
∂u
∂
ξ

= 0.
Из этого следует, что
∂u
∂
ξ
= f (
ξ
),
где f (
ξ
) – некоторая функция. Интегрируя последнее равенство, получим u =
Z
f (
ξ
)d
ξ
+ Q(
η
),
где Q(
η
) – произвольная функция.
Обозначим P (
ξ
) =
Z
f (
ξ
)d
ξ
и подставим вместо
ξ
и
η
их выражения через x и t. В итоге получим равенство (3.4).
Функции P (x − at) и Q(x + at) называются волнами отклонения. На- звание функций связано с их свойствами. Построим графики этих функций при t = 0: y = P (x) и y = Q(x). График функции y = P (x − at) получается параллельным переносом графика функции y = P (x) на at единиц вправо
(a > 0). Соответственно, график функции y = Q(x + at) получается парал- лельным переносом графика y = Q(x) на at единиц влево. Таким образом,
при непрерывном изменении t происходит перемещение графика функции y = P (x) вправо (рис. 3.2), а графика функции y = Q(x) влево (рис. 3.3).
Для того чтобы решить поставленную задачу Коши, следует, поль- зуясь начальными условиями, определить неизвестные функции P и Q.
85

=
y
P ( x −
y = P ( )
x
)
at
at
x
y
0
Рис. 3.2
−at
y
x
y
y
0
= Q ( x + at )
= Q ( x )
Рис. 3.3
Продифференцируем функцию u(x, t) по t:
∂u
∂t
= −aP
0
(x − at) + aQ
0
(x + at).
Подставим выражения для u и
∂u
∂t в начальные условия, положив t = 0.
В результате получим систему уравнений для функций P (x) и Q(x):
(
P (x) + Q(x) =
ϕ
(x),
−aP
0
(x) + aQ
0
(x) =
ψ
(x).
Интегрируя второе равенство в пределах от 0 до x, получим
−a(P (x) − P (0)) + a(Q(x) − Q(0)) =
x
Z
0
ψ
(
ξ
)d
ξ
Отсюда
−P (x) + Q(x) =
1
a x
Z
0
ψ
(
ξ
)d
ξ
+ C,
где C = −P (0) + Q(0) – постоянная величина. Решая систему уравнений,
найдем
P (x) =
1 2
ϕ
(x) −
1 2a x
Z
0
ψ
(
ξ
)d
ξ
−
C
2
,
Q(x) =
1 2
ϕ
(x) +
1 2a x
Z
0
ψ
(
ξ
)d
ξ
+
C
2
Подставляя полученные выражения для P (x) и Q(x) в (3.4), найдем функ- цию u(x, t):
u(x, t) =
1 2
ϕ
(x − at) −
1 2a x−at
Z
0
ψ
(
ξ
)d
ξ
+
1 2
ϕ
(x + at) +
1 2a x+at
Z
0
ψ
(
ξ
)d
ξ
86

Заметим, что
−
x−at
Z
0
ψ
(
ξ
)d
ξ
+
x+at
Z
0
ψ
(
ξ
)d
ξ
=
0
Z
x−at
ψ
(
ξ
)d
ξ
+
x+at
Z
0
ψ
(
ξ
)d
ξ
=
x+at
Z
x−at
ψ
(
ξ
)d
ξ
Тогда функцию u(x, t), являющуюся решением поставленной задачи, мож- но представить в виде u(x, t) =
ϕ
(x − at) +
ϕ
(x + at)
2
+
1 2a x+at
Z
x−at
ψ
(
ξ
)d
ξ
Полученное равенство называется формулой Даламбера решения за- дачи Коши для уравнения колебаний бесконечной струны.
Найденное решение представляет собой сумму двух волн P (x − at) и
Q(x + at). Одна волна “бежит” вправо, другая – влево. Число a =
r T
ρ
в уравнении колебаний струны называется скоростью распространения вол- ны.
Функция u(x, t), полученная методом Даламбера, будет решением по- ставленной задачи при условии, что функция
ϕ
(x) дважды дифференци- руема, а функция
ψ
(x) дифференцируема один раз. В некоторых задачах
ϕ
(x) и
ψ
(x) не имеют нужных производных. Например, если струна в на- чальный момент времени имеет форму ломаной линии (рис. 3.4). В таких случаях считают, что формула Даламбера также дает решение задачи, хо- тя при этом функция u(x, t) не всюду дважды дифференцируема. Такое решение называют обобщенным решением задачи.
Рассмотрим несколько примеров, иллюстрирующих применение мето- да Даламбера.
Пример 3.1. Решить задачу Коши для уравнения колебаний неогра- ниченной струны:
∂
2
u
∂t
2
= a
2
∂
2
u
∂x
2
,
t > 0,
−∞ < x < +∞,
u(x, 0) = he
−x
2
,
∂u(x, 0)
∂t
= ve
−x
2
Для поставленной задачи
ϕ
(x) = he
−x
2
,
ψ
(x) = ve
−x
2
. Используя фор- мулу Даламбера, получим u(x, t) =
he
−(x−at)
2
+ he
−(x+at)
2 2
+
1 2a x+at
Z
x−at ve
−ξ
2
d
ξ
=
87

=
h
2

e
−(x−at)
2
+ e
−(x+at)
2

+
v
√
π
4a
(erf(x + at) − erf(x − at)) ,
где erf x =
2
√
π
x
Z
0
e
−ξ
2
d
ξ
– функция ошибок.
Пример 3.2. Изобразить форму бесконечной струны для моментов времени t
0
= 0, t
1
=
l
2a
, t
2
=
l a
, t
3
=
2l a
, если начальная скорость точек струны равна нулю, а начальная форма – это треугольник на отрезке [−l, l]
с максимальным отклонением h (рис. 3.4).

1   2   3   4   5   6   7   8   9   10