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
).
63
Кратко эту систему можно записать в видеA˜ u j = F j , где A − матрица коэффициентов; ˜ u j − вектор значений функции ˜ u; F j − вектор значений правой части уравнений для j-го временного слоя. Вектор F j меняется от слоя к слою. Матрица коэффициентов системы A трехдиагональная. Для строк этой матрицы выполняется условие строго- го диагонального преобладания |B i | > |A i | + |C i | (i = 0, ..., n). Это условие обеспечивает однозначную разрешимость системы (см. [ 14 ], 1.9). Кроме то- го, для решения системы уравнений с подобной матрицей коэффициентов можно применить метод Гаусса без перестановки строк и столбцов матри- цы. Такой метод решения системы уравнений с трехдиагональной матри- цей коэффициентов называется методом прогонки. Сначала выполняется прямой ход метода прогонки: вычисляются прогоночные коэффициенты ˜ C 0 = C 0 B 0 , ˜ F 0 = F 0 B 0 , (7.9) ˜ 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. (7.10) При этом матрица коэффициентов системы преобразуется к верхней тре- угольной матрице, на главной диагонали которой стоят единицы; ˜ C i – нену- левые элементы матрицы, стоящие выше главной диагонали; ˜ F i – числа, стоящие в правой части системы уравнений. Затем выполняется обратный ход (решается система уравнений с верх- ней треугольной матрицей коэффициентов): ˜ u j n = ˜ F n , ˜ u j i = ˜ F i − ˜ C i ˜ u j i+1 , i = n − 1, ..., 0. (7.11) Метод прогонки можно применять только в том случае, когда в фор- мулах для вычисления прогоночных коэффициентов ˜ C i , ˜ F i знаменатели не обращаются в нуль. Это требование будет выполняться, если для матри- цы коэффициентов справедливо условие строгого диагонального преобла- дания. Таким образом, при применении неявной разностной схемы для на- хождения приближенного решения поставленной задачи значения ˜ u j i будут определяться последовательно по временным слоям. Для каждого времен- ного слоя необходимо будет решать систему уравнений с трехдиагональной матрицей коэффициентов. Преимущество неявной схемы заключается в том, что в отличие от явной разностной схемы построенная неявная схема является безусловно устойчивой, т. е. она устойчива при любых τ и h. 64 Справедливо следующее утверждение. Утверждение 7.4. Неявная разностная схема (7.8) устойчива по начальным данным и правой части, т. е. k˜ uk h,τ ≤ C(k ϕ k h + kf k h,τ ), где C не зависит от h и τ Доказательство данного утверждения (в частном случае) можно найти в пособии [ 12 ]. 7.4. Типовой расчет по теме „Решение уравнения теплопроводности методом Фурье и разностными методами“ (ТР 3.2) Студентам выдаются индивидуальные задания следующего вида: ТР 3.2. Вар. 0. Решить начально-краевую задачу для уравнения тепло- проводности методом Фурье и сеточными методами: ∂u(x, t) ∂t = ∂ 2 u(x, t) ∂x 2 + (−3x + 3)t + (−0.2 sin t)x + 0.1 t + 1 , x ∈ [ 0, 3.0 ], t ∈ [ 0, 0.4 ], u(0, t) = 0.1 ln(t + 1), ∂u(3.0, t) ∂x = 0.2 cos t+0.1, u(x, 0) = x 2 − 5.70x+0.0. Требуется: 1. В области x ∈ [0, l] × [0, T ] получить решение краевой задачи для уравнения теплопроводности в виде ряда Фурье; используя первые 6 слага- емых ряда Фурье, найти приближенное решение задачи при t = T в точках x 1 , . . . , x 6 , разбив промежуток [0, l] на 5 равных частей. 2. Построить явную схему разностной аппроксимации заданной крае- вой задачи; проверить выполнение условия устойчивости явной разностной схемы; разбив каждый промежуток [0, l], [0, T ] на 5 равных частей найти приближенное решение краевой задачи в узлах сетки. 3. Построить неявную схему разностной аппроксимации заданной кра- евой задачи; разбив каждый промежуток [0, l], [0, T ] на 5 равных частей найти приближенное решение краевой задачи в узлах сетки. Пример выполнения типового расчета ТР 3.2. Вар. 0. 1. Решение уравнения теплопроводности методом Фурье. 1) Сведем исходную задачу к задаче с однородными краевыми усло- виями. Для этого выполним замену u(x, t) = v(x, t) + w(x, t). Зададим 65 w(x, t) = α (t)x + β (t) и подберем α (t) и β (t) так, чтобы функция w(x, t) удовлетворяла тем же краевым условиям, что и функция u(x, t): w(0, t) = 0.1 ln(t + 1), ∂w(3, t) ∂x = 0.2 cos t + 0.1. Тогда v(x, t) будет удовлетворять однородным условиям: v(0, t) = 0, ∂v(3, t) ∂x = 0. Подставляя w(x, t) в краевые условия, получим β (t) = 0.1 ln(t + 1), α (t) = 0.2 cos t + 0.1. Значит, w(x, t) = (0.2 cos t + 0.1)x + 0.1 ln(t + 1). Запишем краевую задачу относительно функции v(x, t): ∂v(x, t) ∂t = ∂ 2 v(x, t) ∂x 2 + (−3x + 3) t, x ∈ (0, 3), t ∈ (0, 0.4], v(0, t) = 0, ∂v(3, t) ∂x = 0, v(x, 0) = x 2 − 6x. (7.12) Эту задачу будем решать методом Фурье. Функцию v(x, t) будем ис- кать в виде ряда по собственным функциям оператора L x (v) = − ∂ 2 v ∂x 2 2) Решим задачу Штурма–Лиувилля: ( −y 00 = λ y, x ∈ ( 0, 3 ), y(x) 6= 0, y(0) = 0, y 0 (3) = 0. Краевые условия соответствуют условиям решаемой задачи. Известно, что λ ≥ 0 (см. [12]). Пусть λ = 0. Тогда −y 00 = 0, значит, y(x) = C 1 + C 2 x. Подставим эту функцию в краевые условия( y(0) = C 1 = 0, y 0 (3) = C 2 = 0 =⇒ y(x) = 0 – не является собственной функцией. Пусть λ > 0. Обозначим λ = µ 2 . Тогда −y 00 = µ 2 y. Запишем характе- ристическое уравнение −k 2 = µ 2 . Его корни: k 1,2 = ± µ i, следовательно, y(x) = C 1 cos( µ x) + C 2 sin( µ x). Подставим найденное решение в краевые условия ( y(0) = C 1 = 0, y 0 (3) = C 2 µ cos(3 µ ) = 0. Тогда y(x) = C 2 sin( µ x), C 2 6= 0, при этом cos(3 µ ) = 0 ⇔ 3 µ k = π 2 + π k, λ k = π + 2 π k 6 2 , k = 0, 1, 2, . . . . 66
Найдены собственные числа:
λ
k
=
π
+ 2
π
k
6
2
и собственные функции:
y k
(x) = sin
π
+ 2
π
k
6
x
, k = 0, 1, 2, . . . оператора Штурма–Лиувилля рассматриваемой задачи.
3) Функцию v(x, t) будем искать в виде ряда Фурье по найденным собственным функциям дифференциального оператора y k
(x):
v(x, t) =
+∞
X
k=0
c k
(t) y k
(x).
Функции f (x, t) = (−3x + 3)t и
ϕ
(x) = x
2
− 6 x из дифференциального уравнения и начального условия также разложим в ряды Фурье:
f (x, t) = (−3x + 3)t =
+∞
X
k=0
f k
(t) y k
(x),
f k
(t) =
(f (x, t), y k
(x))
ky k
(x)k
2
,
ϕ
(x) = x
2
− 6 x =
+∞
X
k=0
ϕ
k y
k
(x),
ϕ
k
=
(
ϕ
(x), y k
(x))
ky k
(x)k
2
Вычислим ky k
(x)k
2
и скалярные произведения:
ky k
(x)k
2
=
3
Z
0
sin
2
π
+ 2
π
k
6
x
dx =
3 2
,
(f (x, t), y k
(x)) =
3
Z
0
(−3x + 3)t sin
π
+ 2
π
k
6
x dx =
=
18
π
+ 2
π
k
−
108 (−1)
k
(
π
+ 2
π
k)
2
t ,
(
ϕ
(x), y k
(x)) =
3
Z
0
(x
2
− 6 x) sin
π
+ 2
π
k
6
x dx =
−432
(
π
+ 2
π
k)
3
Найдем коэффициенты Фурье:
f k
(t) =
12
π
+ 2
π
k
−
72 (−1)
k
(
π
+ 2
π
k)
2
t,
ϕ
k
=
−288
(
π
+ 2
π
k)
3
Подставим полученные ряды в дифференциальное уравнение и на- чальные условия задачи (7.12).
67
Выполним подстановку в уравнение ∂v(x, t) ∂t = ∂ 2 v(x, t) ∂x 2 + (−3x + 3) t: +∞ X k=0 c 0 k (t) y k (x) = +∞ X k=0 c k (t) y 00 k (x) + +∞ X k=0 f k (t) y k (x) ⇔ ⇔ +∞ X k=0 c 0 k (t) y k (x) = +∞ X k=0 (− λ c k (t) + f k (t)) y k (x). Подставим ряды в начальное условие v(x, 0) = x 2 − 6x: +∞ X k=0 c k (0) y k (x) = +∞ X k=0 ϕ k y k (x). Поскольку разложение в ряд Фурье единственно ([ 14 ]), то для коэф- фициентов c k (t) выполняются равенства: c 0 k (t) y k (x) = − λ c k (t) + f k (t), c k (0) = ϕ k , (7.13) т. е. коэффициенты c k (t) являются решениями задачи Коши. Уравнение в задаче (7.13) – линейное неоднородное дифференциальное уравнение пер- вого порядка с постоянными коэффициентами и с правой частью специ- ального вида f k (t) = 12 π + 2 π k − 72 (−1) k ( π + 2 π k) 2 t Его можно решить методом неопределенных коэффициентов, или с помощью преобразования Лапласа, или методом вариации постоянной (см. [ 12 ]). Для решения уравнения задачи (7.13) воспользуемся методом неопре- деленных коэффициентов. Сначала найдем общее решение однородного уравнения c 0 k,o (t) y k (x) = − λ c k,o (t). Очевидно: c 0 k,o (t) = C e −λ k t Частное решение ˜ c k (t) неоднородного уравнения будем искать в виде˜ c k (t) = A t + B. Подставив ˜ c k (t) в уравнение, получим A = − λ (A t + B) + G t, где G = 12 π + 2 π k − 72 (−1) k ( π + 2 π k) 2 . Сравнивая коэффициенты при одинаковых степенях t найдем A = 36 G ( π + 2 π k) 2 , B = − 1296 G ( π + 2 π k) 4 68
Запишем общее решение неоднородного уравнения:
c k
(t) = C e
−λ
k t
+ A t + B.
Коэффициент C определим используя начальное условие
C = −
288
(
π
+ 2
π
k)
2
+
1296 G
(
π
+ 2
π
k)
4
Окончательное решение задачи Коши (7.13) имеет вид c
k
(t) =
432 t
(
π
+ 2
π
k)
3
−
15552
(
π
+ 2
π
k)
5
−
2592 t (−1)
k
(
π
+ 2
π
k)
4
+
93312 (−1)
k
(
π
+ 2
π
k)
6
+
+
−
288
(
π
+ 2
π
k)
3
+
15552
(
π
+ 2
π
k)
5
−
93312 (−1)
k
(
π
+ 2
π
k)
6
e
−λ
k t
При найденных значениях c k
(t) искомая функция v(x, t) – решение краевой задачи (7.12) – представляется в виде ряда:
v(x, t) =
+∞
X
k=0
c k
(t) sin
π
+ 2
π
k
6
x
4) Получим решение исходной задачи u(x, t):
u(x, t) = w(x, t) + v(x, t) =
= (0.2 cos t + 0.1) x + 0.1 ln(t + 1) +
∞
X
k=0
c k
(t) sin
π
+ 2
π
k
6
x
5) Найдем приближенное решение исходной задачи, используя N -ю частичную сумму ряда Фурье при N = 6, в точках x
1
, . . . , x
6
, разбив про- межуток [0, l] на 5 равных частей при t = T = 0.4:
u
6
(x i
, T ) = (0.2 cos T + 0.1) x i
+ 0.1 ln(T + 1) +
5
X
k=0
c k
(T ) sin
π
+ 2
π
k
6
x i
Вычислим значения первых шести коэффициентов c k
(T ) ряда Фурье
(табл. 7.1).
Таблица 7.1
k
0 1
2 3
4 5
c k
(T )
−8.59186
−0.00511 0.01336 0.01672 0.00534 0.00454
Вычислим значения функции u(x, T ) в точках x
1
, . . . , x
6
(табл. 7.2).
69
Таблица 7.2 x 0 0.6 1.2 1.8 2.4 3.0 u(x, T ) 0.03365 −2.42786 −4.69670 −6.41486 −7.44352 −7.70303 2. Решение уравнения теплопроводности по явной разностной схеме. Обозначим f (x, t) = (−3x + 3) t + (−0.2 x sin(t)) + 0.1 t + 1 , ϕ (x) = x 2 − 5.7 x – функции из уравнения и начального условия, g 1 (t) = 0.1 ln(t + 1), g 2 (t) = 0.2 cos t + 0.1 – функции из краевых условий заданной задачи для уравнения теплопроводности. Разобьем точками x 0 , . . . , x n промежуток [0, l] на n равных частей и точками t 0 , . . . , t m промежуток [0, T ] на m равных частей (n = m = 5). Получилась сетка, покрывающая область Ω = [0, l] × [0, T ]. Обозначим (x i , t j ) узлы сетки. Пусть h = l n = 3 5 = 0.6, τ = T m = 0.4 5 = 0.08 – шаги сетки по осям Ox и Ot соответственно. Обозначим u j i = u(x i , t j ) – значение искомой функции в узле (x i , t j ). Для нулевого временного слоя (x i , t 0 ), i = 0, ..., n, используем началь- ное условие u 0 i = ϕ (x i ), i = 0, ..., n. В узлах сетки (x i , t j ) i = 1, . . . , n − 1; j = 0, . . . , m − 1 запишем диф- ференциальное уравнение, заменив производные соответствующими раз- ностными отношениями ∂u(x i , t j ) ∂t = u j+1 i − u j i τ + O( τ ) (производная по вре- мени аппроксимируется разностным отношением „вперед“) и ∂ 2 u(x i , t j ) ∂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( τ ) = 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. В граничных узлах сетки (x 0 , t j+1 ) и (x n , t j+1 ), j = 0, ..., m−1 запишем заданные краевые условия: u j+1 0 = g 1 (t j+1 ), u j+1 n − u j+1 n−1 h + O(h) = g 2 (t j+1 ). Считая, что погрешности аппроксимации O( τ ), O(h 2 ), O(h) в уравнениях и граничных условиях малы, отбросим их. В результате получим явную 70 разностную схему для уравнения теплопроводности – систему уравнений относительно неизвестных ˜
u j
i
(приближенных значений функции u(x, t) в узлах сетки):
˜
u
0
i
=
ϕ
(x i
), i = 0, ..., n,
˜
u j+1
i
− ˜
u j
i
τ
=
˜
u j
i−1
− 2˜
u j
i
+ ˜
u j
i+1
h
2
+f (x i
, t j
), i = 1, ..., n−1,
˜
u j+1 0
= g
1
(t j+1
),
˜
u j+1
n
− ˜
u j+1
n−1
h
= g
2
(t j+1
),
j = 0, ..., m−1.
Из полученных уравнений видно, что значения искомой функции на
(j + 1)-м временном слое выражаются через значения на предыдущем слое.
Перепишем полученные уравнения в удобном для расчетов виде, выразив
˜
u j+1
i через ˜
u j
i
Сначала вычисляются значения функции на нулевом временном слое:
˜
u
0
i
=
ϕ
(x i
), i = 0, ..., n.
Затем определяются значения ˜
u j+1
i на следующих временных слоях во внут- ренних точках:
˜
u j+1
i
= ˜
u j
i
+
τ
h
2
(˜
u j
i−1
− 2˜
u j
i
+ ˜
u j
i+1
) +
τ
f (x i
, t j
),
i = 1, . . . , n − 1,
и граничных точках:
˜
u j+1 0
= g
1
(t j+1
),
u j+1
n
= u j+1
i−1
+ hg
2
(t j+1
),
j = 0, . . . , m − 1.
Проверим выполнение условия устойчивости явной разностной схемы.
Условие устойчивости:
τ
≤
h
2 2a
2
. В этой задаче: 0.08 ≤
0.6 2
2 · 1
⇔ 0.08 ≤
≤ 0.18. Условие выполнено, следовательно, применение явной разностной схемы возможно.
Воспользуемся полученной разностной схемой.
Вычислим значения функции на нулевом временном слое при t = 0:
˜
u
0
i
= x
2
i
− 5.7 x i
,
i = 0, . . . , n (табл. 7.3).
Таблица 7.3
i
0 1
2 3
4 5
x
0 0.6 1.2 1.8 2.4 3.0
˜
u
0
i
0.00
−3.06
−5.40
−7.02
−7.92
−8.10
Для краевого условия при x = 0 найдем значения функции
˜
u j+1 0
= 0.1 ln(t j+1
+ 1),
j = 0, . . . , m − 1 (табл. 7.4).
71
Таблица 7.4 i 0 1 2 3 4 t j+1 0.08 0.16 0.24 0.32 0.40 ˜ u j+1 0 0.00770 0.0184 0.02151 0.02776 0.03365 Согласно полученным формулам вычислим значения ˜ u j+1 i : ˜ u j+1 i = ˜ u j i + τ h 2 (˜ u j i−1 − 2˜ u j i + ˜ u j i+1 ) + τ (−3x i + 3) t j+1 + + τ −0.2 x i sin(t j+1 ) + 0.1 t j+1 + 1 , i = 1, . . . , n − 1, u j+1 n = u j+1 i−1 + h(0.2 cos t j+1 + 0.1), x = x n = 3, j = 0, . . . , m − 1. Получим значения функции u j i (i = 0, . . . , 5) на временных слоях j = 0, . . . , 5 (табл. 7.5). Таблица 7.5 i j 0 1 2 3 4 5 0 0.00 −3.06 −5.40 −7.02 −7.92 −8.10 1 0.00770 −2.88397 −5.23797 −6.87026 −7.78254 −7.60292 2 0.01484 −2.74221 −5.08147 −6.73869 −7.59286 −7.41439 3 0.02151 −2.62217 −4.93956 −6.60673 −7.44673 −7.27017 4 0.02776 −2.51593 −4.81043 −6.48734 |
Скачать 0.49 Mb.