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

= J
p
(
√
λx). Подставив решение в краевое условие и учитывая, что dJ
p
(
√
λx)
dx
=
√
λJ
0
p
(
√
λx),
получим
R
√
λJ
0
p
(
√
λT ) + SJ
p
(
√
λT ) = 0.
(1.26)
Обозначим
√
λT =
γ
и перепишем уравнение (1.26) в виде
R
γ
J
0
p
(
γ
) + ST J
p
(
γ
) = 0.
(1.27)
Как известно из 1.5, это уравнение имеет множество простых корней
γ
k
, k = 1, 2, ... . В случае, если R = 0 или S = 0, уравнение (1.27) превраща- ется в J
p
(
γ
) = 0 или J
0
p
(
γ
) = 0 с аналогичными свойствами. Следовательно,
λ будет собственным числом задачи, если
√
λT =
γ
k
. Таким образом, по- лучено множество собственных чисел
λ
k
=

γ
k
T

2
,
где
γ
k
– корень уравнения (1.27), k = 1, 2, ..., и множество собственных функций y
k
(x) = J
p

γ
k
T
x

Система собственных функций n
J
p

γ
k
T
x
o
+∞
k=1
образует полную ортого- нальную систему в пространстве L
2
[0, T ; x].
Для того чтобы разложить функцию в ряд по системе ортогональных функций Бесселя, потребуется несколько интегральных тождеств.
26

Утверждение 1.5. Для любых
α
,
β
∈ R,
α
6=
β
, выполняется
T
Z
0
J
p

α
T
x

J
p

β
T
x

x dx = T
2
β
J
p
(
α
)J
0
p
(
β
) −
α
J
0
p
(
α
)J
p
(
β
)
α
2
−
β
2
(1.28)
Доказательство. Заметим, что функция y(x) = J
p

α
T
x

является решением уравнения с λ =
α
2
T
2
, т. е. справедливо равенство
1
x

x

J
p

α
T
x

0

0
+

α
2
T
2
−
p
2
x
2

J
p

α
T
x

= 0.
Аналогично, если λ =
β
2
T
2
, то
1
x

x

J
p

β
T
x

0

0
+

β
2
T
2
−
p
2
x
2

J
p

β
T
x

= 0.
Умножим первое равенство на J
p

β
T
x

x, а второе – на J
p

α
T
x

x и вы- чтем из первого второе:

x

J
p

α
T
x

0

0
J
p

β
T
x

−

x

J
p

β
T
x

0

0
J
p

α
T
x

+
+
α
2
−
β
2
T
2
J
p

α
T
x

J
p

β
T
x

x = 0.
Проинтегрировав полученное выражение, получим
α
2
−
β
2
T
2
T
Z
0
J
p

α
T
x

J
p

β
T
x

x dx =
=
T
Z
0

x

J
p

β
T
x

0

0
J
p

α
T
x

dx −
T
Z
0

x

J
p

α
T
x

0

0
J
p

β
T
x

dx.
Интегралы в правой части выражения возьмем по частям:
α
2
−
β
2
T
2
T
Z
0
J
p

α
T
x

J
p

β
T
x

x dx =
27

= x

J
p

β
T
x

0
J
p

α
T
x

T
0
−
T
Z
0
x

J
p

β
T
x

0

J
p

α
T
x

0
dx −
− x

J
p

α
T
x

0
J
p

β
T
x

T
0
−
T
Z
0
x

J
p

α
x
T

0

J
p

β
T
x

0
dx =
=
β
J
0
p
(
β
)J
p
(
α
) −
α
J
0
p
(
α
)J
p
(
β
).
Поделив на
α
2
−
β
2
T
2
, получим формулу (1.28).
Утверждение 1.6. Для любого
α
6= 0 справедлива формула
T
Z
0
J
2
p

α
T
x

x dx =
T
2 2

J
0
p
(
α
)

2
+

1 −
p
2
α
2

J
2
p
(
α
)

(1.29)
Доказательство. В выражении (1.28) сделаем предельный переход при
β
→
α
и применим правило Лопиталя:
T
Z
0
J
2
p

α
T
x

x dx = T
2
lim
β
→α
β
J
p
(
α
)J
0
p
(
β
) −
α
J
0
p
(
α
)J
p
(
β
)
α
2
−
β
2
=
= T
2
lim
β
→α
J
p
(
α
)J
0
p
(
β
) +
β
J
p
(
α
)J
00
p
(
β
) −
α
J
0
p
(
α
)J
0
p
(
β
)
−2
β
Из уравнения Бесселя (1.19) следует, что J
00
p
(
β
) = −
1
β
J
0
p
(
β
)−

1 −
p
2
β
2

J
p
(
β
),
поэтому
T
Z
0
J
2
p

α
T
x

dx =
= T
2
lim
β
→α
J
p
(
α
)J
0
p
(
β
) − J
p
(
α
)J
0
p
(
β
) −
β

1 −
p
2
β
2

J
p
(
α
)J
p
(
β
) −
α
J
0
p
(
α
)J
0
p
(
β
)
−2
β
=
=
T
2 2

J
0
p
(
α
)

2
+

1 −
p
2
α
2

J
2
p
(
α
)

28

1.7. Оператор Лежандра. Многочлены Лежандра
Оператором Лежандра называется дифференциальный оператор вида
L(y) = − (1 − x
2
)y
0

0
За область определения оператора возьмем функции, дважды дифферен- цируемые на (−1; 1). Так как на концах интервала функция p(x) = 1 − x
2
обращается в нуль, то в качестве однородных краевых условий потребуем ограниченность y(x) при x → ±1.
Поставим для оператора Лежандра задачу на собственные значения
(
− (1 − x
2
)y
0

0
= λy,
y(x) − ограничена при x → ±1.
(1.30)
Уравнение задачи можно переписать в виде
(x
2
− 1)y
00
+ 2xy
0
− λy = 0.
(1.31)
Рассмотрим теперь выражения вида
P
n
(x) =
1 2
n n!
(x
2
− 1)
n

(n)
Достаточно очевидно, что P
n
(x) является многочленом степени n,
P
0
(x) ≡ 1, P
1
(x) =
1 2
(x
2
− 1)
0
= x, P
2
(x) =
1 4 · 2
(x
2
− 1)
2

00
=
1 2
(3x
2
− 1).
Функции P
n
(x) называются многочленами Лежандра.
Утверждение 1.7. Функция P
n
(x) является собственной функци- ей оператора Лежандра, соответствующей собственному числу λ
n
=
= n(n + 1).
Доказательство. При доказательстве используется формула диффе- ренциального бинома
(uv)
(n)
= uv
(n)
+ C
1
n u
0
v
(n−1)
+ C
2
n u
00
v
(n−2)
+ ... + u
(n)
v.
(1.32)
Здесь C
k n
=
n(n − 1)...(n − k + 1)
k!
– биномиальные коэффициенты.
Запишем очевидное равенство
(x
2
− 1)
n

0
= n(x
2
− 1)
n−1 2x.
Умножим левую и правую части этого равенства на x
2
− 1 и n + 1 раз продифференцируем:
h
(x
2
− 1) (x
2
− 1)
n

0
i
(n+1)
=
n(x
2
− 1)
n
2x

(n+1)
29

Последнее равенство распишем по дифференциальному биному, учитывая,
что (m + 1)-я производная от многочлена степени m равна нулю:
(x
2
−1)
(x
2
− 1)
n

(n+2)
+(n+1)2x
(x
2
− 1)
n

(n+1)
+
(n + 1)n
2!
2
(x
2
− 1)
n

(n)
=
= 2nx
(x
2
− 1)
n

(n+1)
+ 2n(n + 1)
(x
2
− 1)
n

(n)
,
или после приведения подобных
(x
2
− 1)
(x
2
− 1)
n

(n+2)
+ 2x
(x
2
− 1)
n

(n+1)
− n(n + 1)
(x
2
− 1)
n

(n)
= 0.
Умножив последнее соотношение на
1 2
n n!
, получим
(x
2
− 1)P
00
n
(x) + 2xP
0
n
(x) − n(n + 1)P
n
(x) = 0.
Отсюда видно, что P
n
(x) удовлетворяет уравнению (1.31) с λ =
= n(n + 1).
Можно доказать [2], что при других значениях λ задача (1.30) имеет только тривиальное решение y ≡ 0.
Таким образом, {λ
n
= n(n + 1)}
+∞
n=0
образуют систему собственных чи- сел оператора Лежандра, а многочлены {P
n
(x)}
+∞
n=0
– систему собственных функций. Система {P
n
(x)}
+∞
n=0
полна в пространстве L
2
[−1; 1].
Перед тем как вычислять норму P
n
(x), докажем вспомогательное утвер- ждение.
Утверждение 1.8. Многочлен Q(x) =
(x
2
− 1)
n

(m)
, где m < n, об- ращается в нуль на концах промежутка [−1; 1].
Доказательство. Распишем Q(x), используя формулу дифференци- ального бинома (1.32):
Q(x) = [(x − 1)
n
(x + 1)
n
]
(m)
=
= ((x − 1)
n
)
(m)
(x+1)
n
+C
1
n
((x − 1)
n
)
(m−1)
((x+1)
n
)
0
+...+(x−1)
n
((x + 1)
n
)
(m)
Так как m < n, все слагаемые этого выражения содержат множители (x−1)
и (x + 1) и, следовательно, Q(−1) = Q(1) = 0.
Найдем теперь норму собственной функции P
n
(x):
kP
n k
2
=
1
Z
−1
P
2
n
(x)dx =
1
(2
n n!)
2 1
Z
−1
(x
2
− 1)
n

(n)
(x
2
− 1)
n

(n)
dx.
Применим формулу интегрирования по частям:
kP
n k
2
=
1
(2
n n!)
(x
2
− 1)
n

(n)
(x
2
− 1)
n

(n−1)
1
−1
−
30

−
1
(2
n n!)
2 1
Z
−1
(x
2
− 1)
n

(n+1)
(x
2
− 1)
n

(n−1)
dx.
Согласно утверждению 1.8 первое слагаемое равно нулю, а интеграл опять преобразуем по частям и т. д., пока не получим kP
n k
2
=
1
(2
n n!)
2
(−1)
n
1
Z
−1
(x
2
− 1)
n

(2n)
(x
2
− 1)
n dx =
=
(−1)
n
(2n)!
(2
n n!)
2 1
Z
−1
(x
2
− 1)
n dx.
Последний интеграл снова возьмем по частям:
1
Z
−1
(x
2
− 1)
n dx =
1
Z
−1
(x − 1)
n
(x + 1)
n dx =
= (x − 1)
n
(x + 1)
n+1
n + 1 1
−1
−
n n + 1 1
Z
−1
(x − 1)
n−1
(x + 1)
n+1
dx =
=
(−1)
n n!
(n + 1)(n + 2)...2n
1
Z
−1
(x + 1)
2n dx =
(−1)
n n!n!
(2n)!
2 2n+1 2n + 1
В итоге получим kP
n k
2
=
(−1)
n
(2n)!
(2
n n!)
2
(−1)
n
(n!)
2
(2n)!
2 2n+1
(2n + 1)
=
2 2n + 1
(1.33)
Таким образом, любую функцию f ∈ L
2
[−1; 1] можно разложить в ряд
Фурье по многочленам Лежандра f (x) =
∞
X
n=0
c n
P
n
(x),
где коэффициент Фурье вычисляется по формуле c
n
=
(f, P
n
)
kP
n k
2
=
2n + 1 2
1
Z
−1
f (x)P
n
(x)dx.
Отметим еще несколько свойств многочленов Лежандра, которые яв- ляются общими для любых систем ортогональных многочленов.
31

Утверждение 1.9. Пусть Q(x) – произвольный многочлен степени n, тогда его можно представить в виде Q
n
(x) =
n
X
k=0
α
k
P
k
(x), при этом
α
k
=
(Q
n
, P
k
)
kP
k k
2
Доказательство. Этот факт очевидным образом следует из того, что система P
0
(x), ..., P
n
(x) ортогональна и, следовательно, линейно независи- ма, а потому ее можно принять за базис линейного пространства много- членов степени не выше n. Поскольку многочлены ортогональны, то коэф- фициенты разложения
α
k вычисляются по формулам
α
k
=
(Q
n
, P
k
)
kP
k k
2
Утверждение 1.10. Пусть P
n
(x) – многочлен Лежандра, Q
m
(x) –
произвольный многочлен степени m < n. Тогда P
n и Q
m ортогональны.
Доказательство. Запишем Q
m
(x) =
m
X
k=0
α
k
P
k
(x). Скалярное произ- ведение
(Q
m
, P
n
) =
m
X
k=0
α
k
P
k
, P
m
!
=
m
X
k=0
α
k
(P
k
, P
n
) = 0.
Утверждение 1.11. Все корни многочлена Лежандра простые и ле- жат в интервале (−1; 1).
Доказательство. Обозначим через x
1
, x
2
, ..., x m
∈ (−1; 1) точки, где многочлен P
n
(x) меняет знак, и предположим, что m < n. Тогда много- член Q
m
(x) = (x − x
1
)(x − x
2
)...(x − x m
) меняет знак в тех же точках и произведение Q
m
(x)P
n
(x) на промежутке [−1; 1] знак не меняет. А тогда
(Q
m
, P
n
) =
1
Z
−1
Q
m
(x)P
n
(x)dx 6= 0,
что противоречит утверждению 1.8 .
Для многочленов Лежандра выполняются рекуррентные соотноше- ния. Приведем без доказательства две рекуррентные формулы:
(n + 1)P
n+1
(x) − (2n + 1)xP
n
(x) + P
n−1
(x) = 0,
(1.34)
(2n + 1)P
n
(x) = P
0
n+1
(x) − P
0
n−1
(x).
(1.35)
Вывод этих формул можно найти, например, в [7] или [8].
32

1.8. Присоединенные функции Лежандра
Присоединенной функцией Лежандра степени n порядка k называется функция вида
P
k n
(x) = (1 − x
2
)
k/2
P
(k)
n
(x),
k = 0, 1, ...,
n = k, k + 1, ...,
где P
n
(x) – многочлен Лежандра. Заметим, что для четных k функция
P
k n
(x) будет многочленом степени n, а для нечетных – многочленом степени n − 1, умноженным на
√
1 − x
2
. Например,
P
0 0
(x) = 1,
P
0 1
(x) = P
1
(x) = x,
P
1 1
(x) = (1 − x
2
)
1/2
P
0 1
(x) =
√
1 − x
2
,
P
0 2
(x) = P
2
(x) =
3 2
x
2
−
1 2
,
P
1 2
(x) = (1 − x
2
)
1/2
P
0 2
(x) =
√
1 − x
2 3x,
P
2 2
(x) = (1 − x
2
)P
00 2
(x) = (1 − x
2
)3,
P
2 3
(x) = (1 − x
2
)P
00 3
(x) = (1 − x
2
)15x и т. д.
Перед тем как исследовать свойства присоединенных функций Ле- жандра, получим уравнение для производных многочленов Лежандра.
Утверждение 1.12. Производные многочленов Лежандра P
(k)
n
(x) удо- влетворяют уравнению
−(1 − x
2
)z
00
+ 2(k + 1)xz
0
+ k(k + 1)z = n(n + 1)z.
(1.36)
Доказательство. Многочлен Лежандра удовлетворяет уравнению
(1.31)
(x
2
− 1)P
00
n
(x) + 2xP
0
n
(x) = n(n + 1)P
n
(x).
Продифференцируем данное соотношение k раз:
(x
2
− 1)P
00
n
(x)

(k)
+ (2xP
0
n
(x))
(k)
= n(n + 1)P
(k)
n
(x)
и воспользуемся формулой дифференциального бинома (1.32):
(x
2
−1)P
(k+2)
n
(x) + k2xP
(k+1)
n
(x) +
k(k − 1)
2 2P
(k)
n
(x) + 2xP
(k+1)
n
(x) + k2P
(k)
n
=
= n(n + 1)P
(k)
n
После приведения подобных получаем
(x
2
− 1)(P
(k)
n
)
00
+ 2(k + 1)x(P
(k)
n
)
0
+ k(k + 1)P
(k)
n
= n(n + 1)P
(k)
n
Утверждение 1.13. Присоединенная функция Лежандра P
k n
(x) удо- влетворяет уравнению
−(1 − x
2
)y
00
+ 2xy
0
+
k
2 1 − x
2
y = n(n + 1)y.
(1.37)
33

Доказательство. Сделаем в уравнении замену y = (1 − x
2
)
k/2
z. Най- дем из этого равенства y
0
и y
00
и, подставив в (1.37), после преобразований получим
−(1 − x
2
)z
00
+ 2(k + 1)xz
0
+ k(k + 1)z = n(n + 1)z,
т. е. уравнение (1.36).
Уравнение (1.37) можно записать в симметричной форме
− (1 − x
2
)y
0

0
+
k
2 1 − x
2
y = n(n + 1)y.
(1.38)
Если к уравнению (1.38) добавить однородные краевые условия: y(x) огра- ничена при x → ±1, то получим задачу на собственные значения для опе- ратора L
k
(y) = − (1 − x
2
)y
0

0
+
k
2 1 − x
2
y. Таким образом, для оператора
L
k
(y) найден набор собственных чисел {λ
n
= n(n + 1)}
+∞
n=k и набор соб- ственных функций {P
k n
(x)}
+∞
n=k
. Присоединенные функции Лежандра P
k n
и
P
k m
(n 6= m) ортогональны, так как они являются собственными функция- ми симметричного оператора.
Утверждение 1.14. Для присоединенных функций Лежандра спра- ведлива формула
(P
k m
, P
k n
) = (n − k + 1)(n + k)(P
k−1
m
, P
k−1
n
).
(1.39)
Доказательство. Найдем скалярное произведение
(P
k m
, P
k n
) =
1
Z
−1
P
(k)
m
(x)P
(k)
n
(x)(1 − x
2
)
k dx.
Данный интеграл возьмем по частям, считая, что v
0
(x) = P
(k)
m
(x), а u(x) =
= P
(k)
n
(x)(1 − x
2
)
k
:
(P
k n
, P
k m
) = P
(k−1)
m
(x)P
(k)
n
(x)(1 − x
2
)
k
1
−1
−
−
1
Z
−1
h
P
(k+1)
n
(x)(1 − x
2
)
k
− 2kxP
(k)
n
(x)(1 − x
2
)
k−1
i
P
(k−1)
m
(x)dx.
Внеинтегральное слагаемое очевидно равно нулю, поэтому
(P
k n
, P
k m
) =
1
Z
−1
h
−(1 − x
2
)P
(k+1)
n
+ 2kxP
(k)
n i
P
(k−1)
m
(1 − x
2
)
k−1
dx.
(1.40)
34

Согласно формуле (1.36) P
(k−1)
n будет удовлетворять уравнению
−(1 − x
2
)

P
(k−1)
n
(x)

00
+ 2kx

P
(k−1)
n
(x)

0
+ (k − 1)kP
(k−1)
n
(x) =
= n(n + 1)P
(k−1)
n
(x).
Отсюда
−(1 − x
2
)P
(k+1)
n
(x) + 2kxP
(k)
n
(x) = (n
2
+ n − k
2
+ k)P
(k−1)
n
(x).
Подставив полученное соотношение в (1.40), получаем
(P
k n
, P
k m
) = (n
2
+ n − k
2
+ k)
1
Z
−1
P
(k−1)
n
(x)P
(k−1)
m
(x)(1 − x
2
)
k−1
dx =
= (n − k + 1)(n + k)(P
k−1
n
, P
k−1
m
).
Заметим, что в формуле (1.39) коэффициент не зависит от m. Это совершенно естественно, если вспомнить, что при m 6= n скалярное произ- ведение (P
k m
, P
k n
) = 0. При n = m имеем kP
k n
k
2
= (P
k n
, P
k n
) = (n − k + 1)(n + k)(P
k−1
n
, P
k−1
n
) =
= (n − k + 1)(n + k)(n − k + 2)(n + k − 1)(P
k−2
n
, P
k−2
n
) = ... =
= (n + k)(n + k − 1)...(n + 1)(n − k + 1)(n − k + 2)...n(P
0
n
, P
0
n
).
Учитывая, что P
0
n
(x) = P
n
(x), и формулу (1.33), получим kP
k n
k
2
=
(n + k)!
(n − k)!
2 2n + 1
Для фиксированного значения k набор {P
k n
}
+∞
n=k образует полную орто- гональную систему в пространстве L
2
[−1; 1]. Любую функцию f ∈ L
2
[−1; 1]
можно разложить в ряд Фурье по присоединенным функциям Лежандра f (x) =
∞
X
n=k c
k
P
k n
(x),
где c k
– коэффициенты Фурье:
c k
=
(f, P
k n
)
kP
k n
k
2
=
(n − k)!
(n + k)!
(2n + 1)
2 1
Z
−1
f (x)P
k n
(x)dx.
35

1.9. Решение краевой задачи методом конечных разностей
В предыдущих параграфах было получено решение краевой задачи в виде аналитического выражения или в виде ряда. Теперь рассмотрим численный метод решения краевых задач, который называется методом конечных разностей, или методом сеток.
Пусть отрезок [a; b] разбит на n частей. Введем шаг сетки h =
b − a n
,
точки x
0
= a, x
1
= x
0
+ h, x
2
= x
0
+ 2h, ..., x n
= x
0
+ nh = b назовем узлами сетки, а само множество узлов
ω
h
= {x k
} (k = 0, ..., n) – сеткой. Если на отрезке задана непрерывная функция y(x), то функцию {y k
= y(x k
)}
n k=0
естественно называть сеточным аналогом функции y(x). Для аппроксима- ции производных функции в узлах сетки будем рассматривать разностные отношения (линейные комбинации значений сеточной функции в несколь- ких узлах сетки). Из геометрических соображений понятно, что производ- ную функции в узле x k
можно аппроксимировать следующим образом:
y
0
(x k
) ≈
y k+1
− y k
h или y
0
(x k
) ≈
y k
− y k−1
h
Первое из этих выражений принято называть разностным отношени- ем “вперед”, и его можно использовать в узлах x
0
, x
1
, ..., x n−1
, второе –
разностное отношение “назад” подходит для узлов x
1
, x
2
, ..., x n
. Найдем связь между производными дифференцируемой функции y(x) и ее сеточ- ным аналогом.
Утверждение 1.15. Пусть y ∈ C
2
[a,b]
, тогда справедливы оценки y
0
(x k
) −
y(x k+1
) − y(x k
)
h
≤
M
2
h
2
,
y
0
(x k
) −
y(x k
) − y(x k−1
)
h
≤
M
2
h
2
,
где M
2
= max x∈[a,b]
|y
00
(x)|.
Доказательство. В узле x k
запишем формулу Тейлора первого по- рядка с остатком в форме Лагранжа y(x) = y(x k
) + y
0
(x k
)(x − x k
) +
1 2
y
00
(˜
x)(x − x k
)
2
,
где ˜
x ∈ [x k
, x k+1
]. Положим x = x k+1
, тогда x k+1
− x k
= h и y(x k+1
) = y(x k
) + y
0
(x k
)h +
1 2
y
00
(˜
x)h
2
,
36

откуда y(x k+1
) − y(x k
)
h
− y
0
(x k
) =
1 2
y
00
(˜
x)h.
(1.41)
Заменив y
00
(˜
x) на максимум второй производной, получим оценку y
0
(x k
) −
y(x k+1
) − y(x k
)
h
≤
M
2
h
2
Для разностного отношения “назад” доказательство аналогично.
Отметим принципиальный момент. Для того чтобы аппроксимировать первую производную функции разностным отношением, нужно, чтобы функ- ция y(x) была дважды дифференцируемой. Для аппроксимации первой производной y
0
(x k
) можно использовать и другие выражения, например симметричную разность (y k+1
− y k−1
)/2h.
Утверждение 1.16. Пусть y ∈ C
3
[a,b]
, тогда справедлива оценка y
0
(x k
) −
y(x k+1
) − y(x k−1
)
2h
≤
M
3
h
2 6
,
где M
3
= max x∈[a,b]
|y
000
(x)|, k = 1, 2, ..., n − 1.
Доказательство. Запишем формулу Тейлора 2-го порядка с остатком в форме Лагранжа y(x) = y(x k
) + y
0
(x k
)(x − x k
) +
1 2
y
00
(x k
)(x − x k
)
2
+
1 6
y
000
(˜
x)(x − x k
)
3
Для точек x k+1
и x k−1
получим:
y(x k+1
) = y(x k
) + y
0
(x k
)h +
1 2
y
00
(x k
)h
2
+
1 6
y
000
(˜
x)h
3
,
y(x k−1
) = y(x k
) − y
0
(x k
)h +
1 2
y
00
(x k
)h
2
−
1 6
y
000
(
≈
x)h
3
,
где ˜
x ∈ [x k
, x k+1
],
≈
x
∈ [x k−1
, x k
]. Вычтем из первой формулы вторую:
y(x k+1
) − y(x k−1
) = 2y
0
(x k
)h +
1 6

y
000
(˜
x) + y
000
(
≈
x)

h
3
,
откуда y(x k+1
) − y(x k−1
)
2h
− y
0
(x k
) =
1 12

y
000
(˜
x) + y
000
(
≈
x)

h
2
(1.42)
Далее, заменив y
000
(˜
x) и y
000
(
≈
x)) на M
3
, получим искомую оценку.
Для аппроксимации второй производной обычно используется выра- жение y
00
(x k
) ≈
y k−1
− 2y k
− y k+1
h
2 37

Утверждение 1.17. Пусть y ∈ C
4
[a,b]
, тогда справедлива оценка y
00
(x k
) −
y(x k−1
) − 2y(x k
) + y(x k+1
)
h
2
≤
M
4
h
2 12
,
где M
4
= max x∈[a,b]
|y
(4)
(x)|.
Доказательство. Запишем формулу Тейлора 3-го порядка с остатком в форме Лагранжа y(x) = y(x k
) + y
0
(x k
)(x − x k
) +
1 2
y
00
(x k
)(x − x k
)
2
+
1 6
y
000
(x k
)(x − x k
)
3
+
+
1 24
y
(4)
(˜
x)(x − x k
)
4
Подставив x k+1
и x k−1
вместо x, соответственно, получим:
y(x k+1
) = y(x k
) + y
0
(x k
)h +
1 2
y
00
(x k
)h
2
+
1 6
y
000
(x k
)h
3
+
1 24
y
(4)
(˜
x)h
4
,
y(x k−1
) = y(x k
) − y
0
(x k
)h +
1 2
y
00
(x k
)h
2
−
1 6
y
000
(x k
)h
3
+
1 24
y
(4)
(
≈
x)h
4
,
где ˜
x ∈ [x k
, x k+1
],
≈
x
∈ [x k−1
, x k
]. Сложим полученные формулы:
y(x k+1
) + y(x k−1
) = 2y(x k
) + y
00
(x k
)h
2
+
1 24

y
(4)
(˜
x) + y
(4)
(
≈
x)

h
4
Отсюда y(x k−1
) − 2y(x k
) + y(x k+1
)
h
2
− y
00
(x k
) =
1 24

y
(4)
(˜
x) + y
(4)
(
≈
x)

h
2
(1.43)
Заменив в полученном выражении значения y
(4)
(˜
x) и y
(4)
(
≈
x) на их макси- мум M
4
, получим нужную оценку.
Из полученных равенств (1.41)–(1.43) следует, что соответствующее разностное отношение аппроксимирует производную в узле x k
с погреш- ностью, которая ведет себя как произведение ограниченной функции на h
m
(m зависит от вида разностного отношения). Для описания таких по- грешностей удобно использовать символ O(h m
). Опираясь на определение,
известное из курса математического анализа, поясним смысл этого симво- ла. Символом O(h m
) в точке x k
обозначается любая функция, которая в окрестности точки x k
ведет себя как произведение ограниченной функции на h m
, где h = x k+1
− x k
или h = x k
− x k−1
. Использование этого сим- вола позволит однотипно описывать погрешности аппроксимации во всех точках сетки
ω
h
. Кроме того, при решении сеточных задач важно видеть,
38

как зависит погрешность аппроксимации от шага сетки. Из полученных оценок ясно, что эта погрешность зависит только от h (множитель при h m
с уменьшением шага существенно не влияет на изменение погрешности).
Таким образом, в каждой внутренней точке сетки
ω
h справедливы равенства:
y
0
k
(x k
) =
y k+1
− y k
h
+ O(h),
(1.44)
y
0
k
(x k
) =
y k
− y k−1
h
+ O(h),
(1.45)
y
0
k
(x k
) =
y k+1
− y k−1 2h
+ O(h
2
),
y
00
(x k
) =
y k−1
− 2y k
+ y k+1
h
2
+ O(h
2
).
(1.46)
Основные идеи метода конечных разностей рассмотрим на примере краевой задачи





−y
00
+ q(x)y = f (x),
y(a) = 0,
y(b) = 0,
(1.47)
где q(x) ≥ q
0
> 0. Введем на [a, b] сетку {x
0
, x
1
, ..., x n
} с шагом h. Во внут- ренних узлах сетки x k
запишем уравнение, заменяя производную y
00
(x k
)
соответствующим разностным отношением (1.46), в граничных узлах сет- ки x
0
и x n
– краевые условия:
−
y k−1
− 2y k
+ y k+1
h
2
+ O(h
2
) + q(x k
)y k
= f (x k
), k = 1, ..., n − 1,
(1.48)
y
0
= 0,
y n
= 0.
Получилась система уравнений относительно y k
= y(x k
), k = 0, ..., n. Пред- полагая, что погрешности O(h
2
) малы, отбросим их. Тогда придем к систе- ме уравнений относительно сеточной функции ˜
y
−
˜
y k−1
− 2˜
y k
+ ˜
y k
h
2
+ q(x k
)˜
y k
= f (x k
),
k = 1, ..., n − 1,
(1.49)
˜
y
0
= 0,
˜
y n
= 0.
Эта система называется разностной схемой для краевой задачи (1.47).Числа
˜
y k
– это приближенные значения функции y в узлах сетки.
Учитывая граничные условия, систему можно записать только для внутренних узлов сетки. Преобразуем k-е уравнение к виду
−˜
y k−1
+ 2 + h
2
q(x k
)
˜
y k
− ˜
y k+1
= h
2
f (x k
).
39

Последовательно подставим k = 1, 2, ..., n − 1. С учетом того, что ˜
y
0
= ˜
y n
=
= 0, получим









(2 + q(x
1
)h
2
)˜
y
1
− ˜
y
2
= f (x
1
)h
2
,
−˜
y
1
+ (2 + q(x
2
)h
2
)˜
y
2
− ˜
y
3
= f (x
2
)h
2
,
−˜
y
2
+ (2 + q(x
3
)h
2
)˜
y
3
− ˜
y
4
= f (x
3
)h
2
,
−˜
y n−2
+ (2 + q(x n−1
)h
2
)˜
y n−1
= f (x n−1
)h
2
(1.50)
Обозначим
˜y =




˜
y
1
˜
y
2
˜
y n−1




;

f =




f (x
1
)
f (x
2
)
f (x n−1
)




;
A =


2 + q(x
1
)h
2
−1 0
0 0
−1 2 + q(x
2
)h
2
−1 ...
0 0
0 0
0
... −1 2 + q(x n−1
)h
2


и запишем систему (1.50) в матричном виде
A˜
y =
f h
2
Матрица A – трехдиагональная с диагональным преобладанием (модуль диагонального элемента больше, чем сумма модулей остальных элементов строки или столбца). Справедливо следующее утверждение.
Утверждение 1.18. Если A = [a ij
], i = 1, ..., n, j = 1, ..., n – матри- ца со строгим диагональным преобладанием
(|a ij
|
>
X
j6=i
|a ij
|,
i = 1, ..., n), тогда эта матрица невырожденна.
Доказательство. Предположим, что матрица A вырожденна, тогда существует вектор
x 6= 0, для которого A
x = 0. Пусть x p
– максимальный по модулю элемент вектора
x. Очевидно, что |x p
| > 0. Из уравнения A
x = 0
следует, что
|a pp x
p
| =
X
j6=p a
pj x
j
≤ |x p
|
X
j6=p
|a pj
|.
Учитывая, что |a pp x
p
| = |a pp
| |x p
|, получим
|a pp
| ≤
X
j6=p
|a pj
|.
Это неравенство противоречит тому, что A – матрица со строгим диаго- нальным преобладанием. Значит, A – невырожденная матрица.
40

Из этого утверждения следует, что система уравнений (1.50) однознач- но разрешима.
Покажем теперь, что с уменьшением шага сетки h значения реше- ния системы (1.50) будут стремиться к значениям решения краевой задачи
(1.47) в узлах сетки. Для этого введем понятие сеточной нормы.
Сеточной нормой функции v, определенной в узлах сетки
ω
h
= {x k
},
k = 0, ..., n, назовем наибольшее по модулю значение этой функции на сетке kvk h
= max k=0,...,n
|v k
|.
Решение ˜
y разностной краевой задачи при измельчении сетки сходится к решению y краевой задачи (1.47), если ky − ˜
yk h
→ 0
при h → 0.
Если, кроме того, выполняется неравенство ky − ˜
yk h
≤ Ch k
,
где C не зависит от h, тогда разностная схема имеет k-й порядок точности.
Сходимость разностной схемы связана с двумя понятиями: устойчиво- стью и точностью аппроксимации.
Разностная схема называется устойчивой по правой части, если норма ее решения не превосходит нормы функции f , заданной в задаче, умножен- ной на число, не зависящее от шага сетки:
k˜
yk h
≤ Ckf k h
(1.51)
Утверждение 1.19. Разностная схема (1.49) устойчива по правой части, т. е. справедлива оценка (1.51), где C не зависит от шага h.
Доказательство. Пусть наибольшее среди чисел |˜
y k
|, k = 0, ..., n,
есть число |˜
y p
|. Если p = 0 или p = n, то неравенство (1.51) выполняется
(˜
y
0
= 0, ˜
y n
= 0). Предположим, что 0 < p < n. Для всех k = 0, ..., n
|˜
y k
| ≤ |˜
y p
|. Рассмотрим разностную схему (1.49) для узла сетки x p
:
(2 + h
2
q(x p
))˜
y p
= ˜
y p−1
+ ˜
y p
+ h
2
f (x p
).
Так как 2 + h
2
q(x p
) > 0, то
(2 + h
2
q(x p
)) |˜
y p
| = |˜
y p−1
+ ˜
y p
+ h
2
f (x p
)| ≤ |˜
y p−1
| + |˜
y p
| + h
2
|f (x p
)|.
Учитывая, что |˜
y k
| ≤ |˜
y p
| для всех k = 0, ..., n, получим
(2 + h
2
q(x p
))|˜
y p
| ≤ 2|˜
y p
| + h
2
|f (x p
)|,
или
|˜
y p
| ≤
|f (x p
)|
q(x p
)
41

Из этого равенства следует оценка k˜
yk h
≤
1
q
0
kf k h
,
означающая устойчивость разностной схемы (1.49).
Подставим в разностные уравнения (1.49) вместо значений сеточной функции ˜
y k
значения точного решения дифференциальной задачи y(x k
).
Для сохранения равенства в правую часть (1.49) необходимо ввести допол- нительное слагаемое. Если это слагаемое с уменьшением h изменяется как
O(h m
), то говорят, что разностная схема аппроксимирует дифференциаль- ную задачу с погрешностью O(h m
).
Утверждение 1.20. Разностная схема (1.49) аппроксимирует кра- евую задачу (1.47) с погрешностью O(h
2
).
Доказательство. Это утверждение следует из равенств (1.48). Во внутренних узлах сетки функция y удовлетворяет разностной схеме с по- грешностью O(h
2
). Граничные условия задаются точно.
Из устойчивости по правой части и аппроксимации с погрешностью
O(h
2
) следует сходимость разностной схемы.
Утверждение 1.21. Разностная схема (1.49) сходится к решению задачи (1.47) в сеточной норме при h → 0.
Доказательство. Обозначим w k
= y k
− ˜
y k
. Из уравнений (1.48), (1.49)
относительно функций y и ˜
y следует, что сеточная функция w удовлетво- ряет разностной схеме
−
w k−1
− 2w k
+ w k+1
h
2
+ q(x k
)w k
= O(h
2
),
k = 1, ..., n − 1,
(1.52)
w
0
= 0,
w n
= 0.
Таким образом, разность между точным и приближенными решениями краевой задачи удовлетворяет таким же сеточным уравнениям, что и функ- ция ˜
y, только в правой части уравнений стоит погрешность аппроксимации второй производной.
Согласно утверждению 1.19 об устойчивости разностная схема (1.52)
устойчива и справедлива оценка kwk h
≤ C kO(h
2
)k h
,
где C не зависит от h.
Из этого следует, что kwk h
→ 0 при h → 0. Так как справедлива оцен- ка |O(h
2
)| ≤ M h
2
, то построенная разностная схема (1.49) имеет второй порядок точности.
Если краевые условия в задаче (1.47) более сложные, например:
y
0
(a) − S
1
y(a) = t
1
,
(1.53)
42

y
0
(b) + S
2
y(b) = t
2
,
(1.54)
то их тоже надо аппроксимировать разностными уравнениями.
Так, для условия (1.53) используют разностное отношение (1.44)
y
1
− y
0
h
+ O(h) − S
1
y
0
= t
1
,
а для (1.54) – соотношение (1.45)
y n
− y n−1
h
+ O(h) + S
2
y n
= t
2
Тогда, если отбросить погрешности, получается
(1 + S
1
h)˜
y
0
− ˜
y
1
= −t
1
h,
−˜
y n−1
+ (1 + S
2
h)˜
y n
= t
2
h.
Эти уравнения нужно добавить к системе (1.50). Количество уравнений увеличится, однако матрица системы по-прежнему будет трехдиагональ- ной с диагональным преобладанием. Полученная разностная схема будет аппроксимировать краевую задачу с погрешностью O(h).
2. УРАВНЕНИЕ ТЕПЛОПРОВОДНОСТИ
2.1. Вывод одномерного уравнения теплопроводности
Рассмотрим однородный стержень длины l, боковая поверхность кото- рого теплоизолирована. Стержень предполагается настолько тонким, что в каждый момент времени температуру в любом его поперечном сечении можно считать постоянной. Если этот стержень неравномерно нагрет, то- гда в нем будет происходить процесс передачи тепла от более нагретых участков к менее нагретым. Для вывода уравнения теплопроводности бу- дем использовать следующую модель: ось 0x направим вдоль оси стержня,
совместив начало координат с его левым концом. Тогда процесс распро- странения тепла в стержне может быть описан функцией u(x, t), представ- ляющей температуру в сечении стержня x в момент времени t. Частная производная
∂u(x, t)
∂x выражает скорость изменения температуры в направ- лении оси 0x. Если температура в направлении оси 0x растет, тогда
∂u
∂x
> 0,
если температура уменьшается, то
∂u
∂x
< 0.
Вывод дифференциального уравнения теплопроводности базируется на следующих экспериментальных положениях:
43

1. Количество тепла ∆Q, которое необходимо сообщить однородному телу, чтобы повысить его температуру на ∆u, равно:
∆Q = c
ρ
V ∆u,
(2.1)
где c – удельная теплоемкость;
ρ
– плотность; V – объем тела.
2. Количество тепла q(x), протекающее через поперечное сечение x стержня за промежуток времени [t, t + ∆t], пропорционально скорости из- менения температуры в направлении, перпендикулярном к сечению, т. е.
∂u
∂x
, площади сечения S и времени ∆t:
q(x) = −k
∂u
∂x
S ∆t,
(2.2)
где k – коэффициент теплопроводности.
Величину −k
∂u
∂x называют удельным тепловым потоком. Знак минус объясняется тем, что величина теплового потока считается положитель- ной, когда тепло передается в сторону возрастания x. Если
∂u
∂x
> 0, то в направлении оси 0x температура увеличивается, а так как тепло перехо- дит от более нагретых участков к менее нагретым, то тепловой поток будет направлен в сторону уменьшения x.
Рассмотрим малый элемент стержня [x; x + ∆x]. Выведем уравнение,
которому будет удовлетворять функция u(x, t). Для этого составим уравне- ние теплового баланса для выделенного элемента стержня за промежуток времени ∆t.
Количество тепла q(x), входящее через поперечное сечение с абсцис- сой x за промежуток времени ∆t, определяется соотношением (2.2). Если воспользоваться формулой Тейлора f (x + ∆x) = f (x) + f
0
(x)∆x + o(∆x),
величина теплового потока, выходящего через сечение x + ∆x, будет сле- дующей:
q(x + ∆x) = −k
∂u
∂x
S ∆t −
∂
∂x

k
∂u
∂x
S ∆t

∆x + o(∆x).
Считаем приращение ∆x настолько малым, что величиной o(∆x) можно пренебречь.
Если найти разность значений входящего и выходящего тепловых по- токов на участке [x; x + ∆x] за время ∆t, получится количество тепла ∆Q,
сообщенное выбранному участку за указанное время:
∆Q = −k
∂u
∂x
S ∆t +

k
∂u
∂x
S ∆t +
∂
∂x

k
∂u
∂x
S

∆x∆t

=
44

=
∂
∂x

k
∂u
∂x
S

∆x∆t.
При этом за тот же промежуток времени температура в каждой точке выделенного участка изменилась на величину
∆u = u(x, t + ∆t) − u(x, t) =
∂u
∂t
∆t + o(∆t).
Считаем,что ∆t мало. Пренебрегая величиной o(∆t), используя соотноше- ние (2.1) и учитывая, что V = S∆x, получим
∆Q = c
ρ
S∆x
∂u
∂t
∆t.
Приравнивая полученные выражения для ∆Q, составим уравнение тепло- вого баланса:
ρ
cS∆x∆t
∂u
∂t
=
∂
∂x

k
∂u
∂x

S∆x∆t.
(2.3)
После сокращения на общий множитель S∆x∆t получим уравнение c
ρ
∂u
∂t
=
∂
∂x

k
∂u
∂x

,
которое называется линейным уравнением теплопроводности без тепловых источников. Если k, c,
ρ
– постоянные величины, то полученное уравнение можно записать в виде
∂u
∂t
= a
2
∂
2
u
∂x
2
(2.4)
Число a
2
=
k c
ρ
называется коэффициентом температуропроводности.
Предположим теперь, что в некоторых частях стержня находятся ис- точники тепла (выделение или поглощение тепла может происходить в ре- зультате прохождения электрического тока, вследствие химической реак- ции и т. п.). Пусть g(x, t) – функция, описывающая плотность тепловых ис- точников, т. е. такая функция, что на малом участке стержня [x; x + ∆x] за малый промежуток времени ∆t выделяется или поглощается тепло, равное g(x, t)S∆x∆t. Если тепло выделяется, тогда g(x, t) > 0, если поглощается,
то g(x, t) < 0.
Если в уравнении (2.3) учесть тепло, выделяемое (поглощаемое) ис- точниками на участке стержня [x; x + ∆x] за время ∆t, тогда уравнение теплового баланса после сокращения на S∆x∆t примет вид
ρ
c
∂u
∂t
=
∂
∂x

k
∂u
∂x

+ g(x, t).
45

Пусть f (x, t) =
1
ρ
c g(x, t), тогда в случае постоянных k, c,
ρ
полученное уравнение можно записать в виде
∂u
∂t
= a
2
∂
2
u
∂x
2
+ f (x, t).
(2.5)
Это линейное уравнение теплопроводности с учетом источников тепла.
Уравнение теплопроводности (2.4) является линейным однородным, а уравнение (2.5) – линейным неоднородным дифференциальным уравнени- ем.
Очевидно, что однородное уравнение (2.4) имеет тривиальное решение u = 0. Кроме того, любая функция, не зависящая от t и линейная относи- тельно x, удовлетворяет этому уравнению. Из этого следует, что уравнение теплопроводности имеет бесконечное множество решений.
Для того чтобы получить единственное решение, к уравнению тепло- проводности добавляют обычно начальное и краевые условия.
2.2. Постановка начально-краевых задач
Начальное условие для уравнения теплопроводности состоит в задании значений функции u(x, t) в начальный момент времени t = 0:
u(x, 0) =
ϕ
(x),
где
ϕ
(x) – функция, описывающая начальную температуру стержня. Кра- евые (граничные) условия задаются на торцевых концах стержня в точках x = 0 и x = l в соответствии с теплообменом стержня с окружающей средой.
Рассмотрим различные случаи краевых условий.
1. На концах стержня поддерживается заданная температура (в точке x = 0 – температура u
0
(t), а в точке x = l – температура u l
(t)):
u(0, t) = u
0
(t),
u(l, t) = u l
(t).
Это условия первого рода, или условия Дирихле.
2. Концы стержня теплоизолированы. Это означает, что поток тепла
−k
∂u
∂x через единицу поверхности соответствующего конца стержня будет равен нулю. Тогда краевые условия на концах стержня можно задать сле- дующим образом:
∂u(0, t)
∂x
= 0;
∂u(l, t)
∂x
= 0.
Полученные условия – это однородные условия второго рода, или условия
Неймана.
46

В общем случае могут быть известны значения теплового потока на границе: q
0
(t) (на конце x = 0) и q l
(t) (на конце x = l). Тогда краевые условия второго рода записывают в виде k
∂u(0, t)
∂x
= q
0
(t),
−k
∂u(l, t)
∂x
= q l
(t).
(Если тепло уходит в окружающую среду, то q
0
(t), q l
(t) > 0.)
3. На концах стержня происходит теплообмен с окружающей средой по закону Ньютона. В этом случае поток тепла пропорционален разности температур тела и окружающей среды и краевые условия имеют вид k
∂u(0, t)
∂x
= h(u(0, t) − u
0
(t));
−k
∂u(l, t)
∂x
= h(u(l, t) − u l
(t)).
Здесь h – коэффициент теплоотдачи (h > 0); u
0
(t), u l
(t) – температура среды на концах x = 0 и x = l соответственно. Такие краевые условия называются условиями третьего рода.
Все рассмотренные краевые условия являются линейными. Они могут быть описаны одним соотношением на соответствующей границе:
R
1
∂u(0, t)
∂x
− S
1
u(0, t) = g
1
(t),
R
2
∂u(l, t)
∂x
+ S
2
(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) –
условие третьего рода.
4. Если тело находится в вакууме, тогда изменение температуры на границе может происходить вследствие теплоизлучения по закону Стефана–
Больцмана. Граничное условие в этом случае, например на конце x = 0,
имеет вид k
∂u(0, t)
∂x
=
σ
u
4
(0, t) − u
4 0
(t)
,
где
σ
– постоянная Стефана–Больцмана; u
0
– температура окружающей среды. Такое краевое условие, в отличие от рассмотренных ранее, является нелинейным.
При решении прикладных задач на концах x = 0 и x = l могут быть поставлены краевые условия разного рода.
5. Для бесконечного стержня (−∞ < x < +∞) предполагается, что температура в его бесконечно удаленных точках ограничена. Подобные за- дачи возникают при изучении процесса теплопроводности в очень длинном стержне. Температурный режим на концах такого стержня слабо влияет на температуру в его центральной части. Существенное значение имеет только начальное распределение температуры в стержне.
47

В общем трехмерном случае, если изучается нагрев тела Ω с кусочно- гладкой границей Γ, уравнение теплопроводности с учетом источников или стоков тепла для функции u(M, t) (M ∈ Ω и t > 0) записывается в виде
ρ
c
∂u
∂t
= div(k grad u) + g,
где k – коэффициент теплопроводности; c – удельная теплоемкость;
ρ
–
плотность тела; g – функция, описывающая объемную плотность источни- ков тепла. Если тело однородное и изотропное, т. е. можно считать, что
ρ
,
c, k – постоянные величины, то уравнение теплопроводности обычно пре- образуют к виду
∂u
∂t
= a
2
∆u + f,
где a
2
=
k
ρ
c
; f =
g
ρ
c
; ∆ = div(grad u) – оператор Лапласа. Вид оператора зависит от выбранной системы координат:
∆u =
∂
2
u
∂x
2
+
∂
2
u
∂y
2
+
∂
2
u
∂z
2
(в декартовой),
∆u =
1
ρ
∂
∂
ρ

ρ
∂u
∂
ρ

+
1
ρ
2
∂
2
u
∂
ϕ
2
+
∂
2
u
∂z
2
(в цилиндрической),
∆u =
1
ρ
2
∂
∂
ρ

ρ
2
∂u
∂
ρ

+
1
ρ
2
sin
θ
∂
∂
θ

sin
θ
∂u
∂
θ

+
1
ρ
2
sin
2
θ
∂
2
u
∂
ϕ
2
(в сфери- ческой).
Вывод трехмерного уравнения подробно описан в [2].
К этому уравнению добавляется начальное условие при t = 0
u(M, 0) =
ϕ
(M ),
где
ϕ
(M ) – функция, описывающая начальную температуру тела.
На границе области Γ ставятся краевые условия в соответствии с усло- виями теплообмена тела с окружающей средой. Это может быть одно из трех условий вида:
u
Γ
=
µ
(M, t),
(2.6)
∂u
∂
n
Γ
=
ν
(M, t),
(2.7)
∂u
∂
n
+ hu
Γ
=
χ
(M, t).
(2.8)
Здесь
µ
,
ν
,
χ
, h – заданные и непрерывные на границе Γ функции; M
– точка границы области Γ;
∂u
∂
n
– производная по направлению внешней нормали к границе Γ.
48

Условие (2.6) – это условие первого рода, которое ставится в том слу- чае, когда известна температура тела на границе. Условие (2.7) – это усло- вие второго рода, ставится тогда, когда на границе задан тепловой поток.
Условие (2.8) – условие третьего рода, ставится, если теплообмен тела со средой происходит по закону Ньютона.
Начально-краевая задача для уравнения теплопроводности с краевым условием первого рода на всей границе называется задачей Дирихле, а с краевым условием второго рода – задачей Неймана.
На разных частях границы Γ могут быть заданы краевые условия раз- ного рода. В этом случае задача называется смешанной.
2.3. Решение краевой задачи для уравнения теплопроводности методом Фурье
Метод Фурье является одним из наиболее распространенных методов решения линейных уравнений в частных производных.
Рассмотрим начально-краевую задачу для одномерного уравнения теп- лопроводности. Пусть требуется найти функцию u(x, t), удовлетворяющую уравнению
∂u
∂t
= a
2
∂
2
u
∂x
2
+ f (x, t)
(0 < x < l,
t > 0),
начальному условию u(x, 0) =
ϕ
(x),
краевым условиям
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
| + |S
1
| 6= 0, |R
2
| + |S
2
| 6= 0). Решением поставленной задачи назовем функцию u(x, t), обладающую следующими свойствами:
а) u(x, t) определена и непрерывна в области Ω = [0, l] × [0, T ];
б) u(x, t) удовлетворяет уравнению теплопроводности при 0 < x < l,
t > 0;
в) u(x, t) удовлетворяет начальному условию и краевым условиям.
Метод Фурье непосредственно применяется к задачам, в которых ис- комая функция удовлетворяет однородным краевым условиям (g
1
(t) ≡ 0 и g
2
(t) ≡ 0). Если краевые условия неоднородные, тогда сначала задачу сле- дует свести к задаче с однородными условиями. Для этого функцию u(x, t)
представляют в виде u(x, t) = v(x, t) + w(x, t),
49

где w(x, t) – дифференцируемая функция, удовлетворяющая тем же крае- вым условиям, что и функция u(x, t):
R
1
∂w(0, t)
∂x
− S
1
w(0, t) = g
1
(t),
R
2
∂w(l, t)
∂x
+ S
2
w(l, t) = g
2
(t).
В этом случае функция v(x, t) будет удовлетворять однородным условиям
(см. 1.1). Выбор функции w(x, t) зависит от типа граничных условий. Для многих краевых задач функция w(x, t) может быть задана как линейная по переменной x w(x, t) =
α
(t)x +
β
(t).
Сумму функций v(x, t) + w(x, t) подставляют в дифференциальное уравнение и начальные условия, и всю краевую задачу записывают для функции v(x, t), которая удовлетворяет однородным краевым условиям.
Задача сводится к нахождению функции v(x, t), удовлетворяющей уравне- нию
∂v
∂t
= a
2
∂
2
v
∂x
2
+ f v
(x, t)
(0 < x < l,
t > 0),
(2.9)
начальному условию v(x, 0) =
ϕ
v
(x),
(2.10)
краевым условиям
R
1
∂v(0, t)
∂x
− S
1
v(0, t) = 0,
R
2
∂v(l, t)
∂x
+ S
2
v(l, t) = 0,
(2.11)
где f v
(x, t) = f (x, t) −
∂w(x, t)
∂t
+ a
2
∂
2
w
∂x
2
,
ϕ
v
(x) =
ϕ
(x) − w(x, 0).
Далее для задачи с однородными краевыми условиями метод Фурье применяется по следующей схеме:
1. Для линейного дифференциального оператора второго порядка
L
x
(u) = −
∂
2
u
∂x
2
решают соответствующую задачу Штурма–Лиувилля:
−y
00
(x) = λy(x),
0 < x < l
(y(x) 6≡ 0),
R
1
y
0
(0) − S
1
y(0) = 0,
R
2
y
0
(l) + S
2
y(l) = 0.
В результате находят собственные числа {λ
k
}
+∞
k=1
и систему собственных функций {y k
(x)}
+∞
k=1
оператора краевой задачи.
2. Функцию v(x, t) представляют в виде ряда Фурье по собственным функциям {y k
(x)}
+∞
k=1
v(x, t) =
+∞
X
k=1
c k
(t)y k
(x).
Однородные краевые условия (2.11) для v(x, t) при этом автоматически выполняются.Этот ряд подставляют в уравнение (2.9) и начальное условие
50

(2.10), предварительно разложив функции f v
(x, t) и
ϕ
v
(x) в ряды Фурье по той же системе функций {y k
(x)}
+∞
k=1
. Далее, для коэффициентов c k
(t)
получают и решают задачи Коши.
3. Искомую функцию u(x, t) записывают в виде u(x, t) = w(x, t) +
+∞
X
k=1
c k
(t)y k
(x).
Покажем, как применяется этот метод, на примерах.
Пример 2.1. Найти температуру тонкого однородного стержня дли- ны l с теплоизолированной поверхностью, если его начальная температура u(x, 0) =
T
0
x l
, конец x = 0 теплоизолирован, а конец x = l поддерживается при постоянной температуре T
0
Функция u(x, t), описывающая распределение температуры в стержне,
удовлетворяет уравнению
∂u
∂t
= a
2
∂
2
u
∂x
2
(0 < x < l,
t > 0),
начальному условию u(x, 0) =
T
0
x l
и краевым условиям
∂u(0, t)
∂x
= 0,
u(l, t) = T
0
Метод Фурье в случае данных краевых условий непосредственно непри- меним, так как условие на конце стержня x = l неоднородное. Поэто- му сведем сначала поставленную задачу к задаче с однородными краевы- ми условиями. Представим функцию u(x, t) в виде суммы двух функций u(x, t) = v(x, t) + w(x, t).
Пусть w(x, t) =
α
x+
β
. Подберем
α
и
β
так, чтобы для функции w(x, t)
выполнялись краевые условия



∂w(0, t)
∂x
=
α
= 0,
w(l, t) =
α
l +
β
= T
0
Тогда
α
= 0,
β
= T
0
, w(x, t) = T
0
и, следовательно, u(x, t) = v(x, t) + T
0
В этом случае справедливы равенства
∂u
∂t
=
∂v
∂t
,
∂
2
u
∂x
2
=
∂
2
v
∂x
2
,
u(x, 0) = v(x, 0) + T
0 51

Тогда для функции v(x, t) получится следующая задача с однородны- ми краевыми условиями:
∂v
∂t
= a
2
∂
2
v
∂x
2
(0 < x < l,
t > 0),
v(x, 0) =
T
0
x l
− T
0
,
∂v(0, t)
∂x
= 0,
v(l, t) = 0.
Эту задачу решим методом Фурье.
1. Функцию v(x, t) будем искать в виде ряда Фурье по собственным функциям линейного дифференциального оператора
L
x
(u) = −
∂
2
u
∂x
2
второго порядка краевой задачи. Эти функции найдем, ре- шая соответствующую задачу Штурма–Лиувилля
−y
00
(x) = λy(x),
0 < x < l
(y(x) 6≡ 0),
y
0
(0) = 0,
y(l) = 0.
Собственные числа оператора −y
00
удовлетворяют условию λ ≥ 0.
При λ = 0 общее решение уравнения y(x) = C
1
+ C
2
x. Подставив y(x)
в краевые условия, получим C
1
= 0, C
2
= 0, тогда y(x) ≡ 0 и λ = 0 не является собственным числом оператора.
При λ > 0 обозначим λ =
µ
2
. Общим решением уравнения будет функ- ция y(x) = C
1
cos(
µ
x) + C
2
sin(
µ
x).
Функция y(x) удовлетворяет однородным краевым условиям:
(
y
0
(0) = C
2
µ
= 0,
y(l) = C
1
cos(
µ
l) + C
2
sin(
µ
l) = 0.
Из этих равенств следует, что при
µ
6= 0
C
2
= 0
и
C
1
cos(
µ
l) = 0.
Так как C
1 6= 0 (функция y(x) не может быть нулевой), то cos(
µ
l) = 0.
Полученное равенство справедливо для чисел
µ
k l =
π
2
+
π
k, т. е.
µ
k
=
(2k + 1)
π
2l
,
k = 0, 1, 2, ... .
Таким образом, для данной задачи Штурма–Лиувилля собственными числами и собственными функциями являются:
λ
k
=
µ
2
k
=
(2k + 1)
π
2l

2
,
y k
(x) = cos(
µ
k x),
k = 0, 1, 2, ... .
52

2. Функцию v(x, t) будем искать в виде ряда по найденным собствен- ным функциям дифференциального оператора:
v(x, t) =
+∞
X
k=0
c k
(t)y k
(x).
Подставим ряд в уравнение теплопроводности и в начальное условие:
+∞
X
k=0
c
0
k
(t)y k
(x) = a
2
+∞
X
k=0
c k
(t)y
00
k
(x),
+∞
X
k=0
c k
(0)y n
(x) =
+∞
X
k=0
ϕ
k y
k
(x).
Здесь
+∞
X
k=0
ϕ
k y
k
(x) – ряд Фурье для функции
ϕ
v
(x) =
T
0
x l
− T
0
, коэффици- енты Фурье которого вычисляются по правилу
ϕ
k
=
(
ϕ
v
, y k
)
ky k
k
2
Используя равенство y
00
k
(x) = −
µ
2
k y
k
(x) и свойство единственности раз- ложения функции в ряд Фурье, получим задачи Коши для коэффициентов c
k
(t):
(
c
0
k
(t) = −a
2
µ
2
k c
k
(t),
c k
(0) =
ϕ
k
(2.12)
Найдем сначала коэффициенты
ϕ
k
:
ky k
k
2
=
l
Z
0
cos
2
(
µ
k x)dx =
l
2
,
(
ϕ
v
, y k
) =
l
Z
0
T
0
x l
− T
0

cos(
µ
k x)dx = −
T
0
l
µ
2
k
,
тогда
ϕ
k
= −
2T
0
l
2
µ
2
k
, k = 0, 1, 2, ... .
Общее решение уравнения (2.12) c k
(t) = A
k e
−a
2
µ
2
k t
. Используя равен- ство c k
(0) =
ϕ
k
, получим A
k
=
ϕ
k
, а значит,
c k
(t) = −
2T
0
l
2
µ
2
k e
−a
2
µ
2
k t
и v(x, t) = −
2T
0
l
2
+∞
X
k=0
e
−a
2
µ
2
k t
µ
2
k cos(
µ
k x).
53

3. Окончательно решение исходной задачи получим как сумму функ- ций w(x, t) и v(x, t):
u(x, t) = T
0
−
8T
0
π
2
+∞
X
k=0
e
−a
2
µ
2
k t
(1 + 2k)
2
cos
(1 + 2k)
π
x
2l

Пример 2.2. В цилиндрическом проводнике радиуса R вследствие прохождения постоянного тока в соответствии с законом Джоуля–Ленца выделяется тепловая энергия с объемной плотностью Q
0
. Теплоотдача с поверхности проводника происходит по закону Ньютона. Найти распреде- ление температуры по сечению проводника, если его начальная темпера- тура и температура внешней среды равны T
0
В примере изучается нагрев трехмерного цилиндрического тела. Вы- берем произвольное поперечное сечение проводника. Введем цилиндриче- скую систему координат, поместив ее начало в центре выбранного сечения и направив ось 0z вдоль оси проводника. Тогда процесс нагрева проводника описывается функцией u(
ρ
,
ϕ
, z, t), удовлетворяющей уравнению теплопро- водности
∂u
∂t
= a
2
∆u + Q,
где Q =
Q
0
ρ
c
,
ρ
– плотность, c – удельная теплоемкость проводника.
Запишем уравнение в цилиндрической системе координат
∂u
∂t
= a
2
1
ρ
∂
∂
ρ

ρ
∂u
∂
ρ

+
1
ρ
2
∂
2
u
∂
ϕ
2
+
∂
2
u
∂z
2

+ Q,
0 <
ρ
< R,
0 ≤
ϕ
< 2
π
,
−∞ < z < +∞,
t > 0.
Распределение температуры в каждом сечении проводника можно счи- тать одинаковым (функция u от z не зависит), т. е. в уравнении
∂
2
u
∂z
2
= 0.
Можно заметить, что распределение температуры в выбранном сечении зависит только от времени и от расстояния до центра проводника, т. е.
∂
2
u
∂
ϕ
2
= 0. Тогда функцию u можно рассматривать как функцию двух пе- ременных u(
ρ
, t).
Функция u(
ρ
, t) удовлетворяет уравнению теплопроводности:
∂u
∂t
= a
2 1
ρ
∂
∂
ρ

ρ
∂u
∂
ρ

+ Q,
0 <
ρ
< R,
t > 0,
начальному условию: u(
ρ
, 0) = T
0
и краевым условиям:
u(
ρ
, t) ограничена при
ρ
→ 0 + 0,
∂u(R, t)
∂
ρ
+ h(u(R, t) − T
0
) = 0 54

(h =
h
0
k
, где h
0
– коэффициент теплоотдачи, k – коэффициент теплопро- водности материала проводника).
Условие ограниченности функции u при
ρ
→ 0 + 0 является естествен- ным однородным краевым условием. Условие при
ρ
= R неоднородное.
Заменой u(
ρ
, t) = v(
ρ
, t) + w(
ρ
, t) сведем данную краевую задачу к краевой задаче относительно функции v(
ρ
, t) с однородными краевыми условиями.
Положим w(
ρ
, t) = C = const. Первое краевое условие будет выполняться.
Подставив w во второе краевое условие, убеждаемся в том, что оно будет выполняться при C = T
0
. Таким образом, w(
ρ
, t) = T
0
. Для функции v(
ρ
, t)
получим следующую начально-краевую задачу:
∂v
∂t
= a
2 1
ρ
∂
∂
ρ

ρ
∂v
∂
ρ

+ Q
(0 <
ρ
< R,
t > 0),
v(
ρ
, 0) = 0,
v(
ρ
, t) ограничена при
ρ
→ 0 + 0,
∂v(R, t)
∂
ρ
+ hv(R, t) = 0.
Краевые условия для функции v однородные, и метод Фурье непосред- ственно применим к этой задаче.
1. Сначала найдем собственные функции и собственные числа опера- тора B
0
(v) = −
1
ρ
∂
∂
ρ

ρ
∂v
∂
ρ

– оператора Бесселя нулевого порядка. Для этого решим соответствующую задачу Штурма–Лиувилля:
−
1
ρ
(
ρ
y
0
(
ρ
))
0
= λy(
ρ
),
0 <
ρ
< R
(y(
ρ
) 6≡ 0),
y(
ρ
) ограничена при
ρ
→ 0 + 0, y
0
(R) + hy(R) = 0.
Дифференциальное уравнение можно записать в виде y
00
+
1
ρ
y
0
+ λy = 0.
Это уравнение Бесселя нулевого порядка. Из общей теории следует, что собственные числа рассматриваемого оператора неотрицательны.
При λ = 0 уравнение Бесселя становится уравнением Эйлера:
ρ
2
y
00
+
ρ
y
0
= 0.
Выполнив подстановку
ρ
= e t
, найдем y(
ρ
) = C
1
+ C
2
ln
ρ
. Учитывая кра- евые условия,получим y(
ρ
) ≡ 0. Значит, λ = 0 не является собственным числом оператора.
При λ > 0 решением уравнения Бесселя является функция y(
ρ
) = C
1
J
0
(
√
λ
ρ
) + C
2
N
0
(
√
λ
ρ
),
55

где J
0
и N
0
– функции Бесселя и Неймана соответственно.
Функция N
0
(
√
λ
ρ
) не ограничена при
ρ
→ 0 + 0 (см. 1.5). Ограни- ченным в нуле решением уравнения будет функция y(
ρ
) = C
1
J
0
(
√
λ
ρ
).
Подставим ее во второе краевое условие
(J
0
(
√
λ
ρ
))
0
+ hJ
0
(
√
λ
ρ
)
ρ
=R
= 0.
Так как J
0 0
(
ρ
) = −J
1
(
ρ
), справедливо равенство
−
√
λJ
1
(
√
λR) + hJ
0
(
√
λR) = 0.
Обозначим
√
λR =
γ
и получим уравнение относительно
γ
−
γ
J
1
(
γ
) + hRJ
0
(
γ
) = 0.
(2.13)
Покажем, что это уравнение имеет корни. Для этого запишем его в виде
J
1
(
γ
)
J
0
(
γ
)
=
hR
γ
Корни полученного уравнения являются абсциссами точек пересечения гра- фиков функций y =
J
1
(
γ
)
J
0
(
γ
)
и y =
hR
γ
при
γ
> 0 (рис. 2.1).
γ
γ
γ
γ
γ
2 3
4 5
6

1   2   3   4   5   6   7   8   9   10