Разностные схемы

общем же случае он представляет собой некоторую матрицу). Когда α изменяется от 0 до 2π, значения λ(α) образуют спектр оператора перехода. Тогда необходимое условие устойчивости можно сформулировать следующим образом: спектр оператора перехода на один шаг по времени должен лежать внутри единичного круга на комплексной плоскости. Если оператор перехода S является нормальной матрицей (SSH-SHS=0, где символ “H” обозначает последовательное применение двух операций: комплексного сопряжения и транспонирования), то критерий фон Неймана (11.15) является также и достаточным условием устойчивости. В случае одного уравнения, выражение для λ(α) можно получить, если подставить (11.14) непосредственно в разностную схему; это даст некоторое уравнение на λ(α).

Рассмотрим теперь устойчивость построенных выше схем. Параметр λ(α) для схемы

(11.2) мы уже определили:

λ(α) = 1 −r (1 −exp(−iα)) = 1 −r + r cos(α) −ir sin(α),

где r=cτ/h. Тогда

| λ(α) | = (1 −r)2 + 2r(1 −r)cos(α) +r2 .

Когда α изменяется от 0 до 2π, |λ(α)| пробегает значения от 1 до |1-2r|. Следовательно, по критерию устойчивости (11.18) эта схема будет устойчива, если

r≤1 или chτ ≤1.

Схема (11.2) является примером условно устойчивой схемы.

Это условие совпадает с условием Куранта-Фридрикса-Леви h ≥cτ , записанного для случая гиперболических уравнений. Это условие имеет простое физическое объяснение. Предположим, что начальное условие u(x,0) равно 0 во всех точках пространственной сетки кроме одной. Тогда решение в узле сетки удалённом от этой точки на n пространственных шагов должно оставаться нулевым, по крайней мере, в течении n временных шагов. Это эквивалентно условию, что шаг по времени τ должен быть меньше того промежутка времени, в течении которого волна проходит расстояние равное пространственному шагу h или, что скорость распространения информации на сетке, h/τ, должна быть больше скорости

волны.

Если мы подставим разностное решение вида (11.14) в схему (11.4), то получим

λk+1 exp(iαn) −λk exp(iαn)

τ

+ cλk+1 exp(iαn) −exp(iα(n −1)) = 0 . h

11

Далее, разделим это выражение на λkexp(iαn) и получим следующее уравнение для λ:

λ−1 + cλ1 −exp(−iα) = 0 ,

τ h

из которого нетрудно получить, что

| λ(α) | =			1
		+ r)2 −2r(1	+ r)cos(α)	+r2
(1

Так как r≥0, то |λ(α)| будет всегда меньше единицы. Это говорит о том, что разностная схема (11.4) устойчива при любых значениях параметров h и τ (безусловно устойчива). Обратимся теперь к схеме (11.5). Применяя использованный выше подход, получим следующее уравнение

λ−1 + cλ exp(iα) −1 = 0 ,

τ h

из которого следует, что

| λ(α) | = (1 + r)2 −2r(1 + r)cos(α) +r2 .

Когда α изменяется от 0 до 2π, |λ(α)| пробегает значения от 1 до 1+2r. Так как r≥0, то |λ(α)| будет всегда больше единицы. Это говорит о том, что разностная схема (11.5) неустойчива. Неустойчивые разностные схемы обычно дают осциллирующие решения возрастающей амплитуды, подобные тому, показанному на рисунке 11.4, и, естественно, такие схемы не имеют никакого практического значения.

11.1.4. Принцип замороженных коэффициентов

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

∂∂ut + u ∂∂ux = 0 , -∞≤x≤∞, t≥0,

u(x, 0) = g(x), -∞≤x≤∞,

По аналогии с (11.2) мы можем построить следующую схему:

uk+1	−uk	uk	−uk
n	n + unk	n	n−1	= 0 ,	(11.16)
τk			h

12

n=0, 1, …; k=0, 1, …

un0 = g(xn ), n=0, 1, ….

Для проведения вычислений по этой формуле, нам необходимо выбрать шаг по времени, а это требует анализа устойчивости. Непосредственно применить критерий фон Неймана мы не можем. Поэтому сначала применим следующий подход: «заморозим» переменный коэффициент в некоторой точке (x*,t*) области определения задачи. В результате, мы получим набор разностных схем

uk+1	−uk	uk	−uk
n	n + u(x*,t* )	n	n−1	= 0 ,	(11.17)
τk			h

n=0,1,…; k=0,1,… , un0 = g(xn ), n=0,1,…. ,

где u(x*,t*) некоторый постоянный коэффициент. Принцип замороженных коэффициентов можно сформулировать следующим образом: разностная схема (11.16) будет устойчива, если схема (11.17) будет устойчива для всех значений x*,t* из области определения задачи. Для c=u(x*,t*)=const условие устойчивости нам известно: cτ/h ≤1 . Учтём теперь, что

1)в исходной задаче скорость распространения волны зависит от самого решения и,

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

В результате мы приходим к следующему условию:

τ	k	max uk			h
		n	n	≤1 или τk ≤		.
		h			max unk
					n

11.1.5. Шаблон разностной схемы

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

– горизонтально вправо. Если два узла связаны вертикальной или горизонтальной прямой, то

13

они используются для аппроксимации первой производной по времени или по пространству. Если три узла связаны вертикальной или горизонтальной прямой, то они используются для аппроксимации второй производной по времени или по пространству. Например, разностные схемы (11.2), (11.4) и (11.5) имеют следующие шаблоны, показанные на рис. 11.5.

Рисунок 11.5. Шаблоны для схем (11.2), (11.4) и (11.5), соответственно.

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

uk+1	−uk		uk+1	−uk+1	uk	−uk
n	n	+ c β	n	n−1 + (1 −β)	n	n−1		= 0 ,	(11.18)
	τ			h		h

0≤β≤1,

которая переходит в схему (11.2) при β=0 и в схему (11.4) при β=1.

Рисунок 11.6. Шаблон для схемы (11.18).

11.2.Распространение тепла (диффузия)

Вкачестве примера диффузионного процесса мы рассмотрим распространение тепла в неподвижной среде. Эта задача формулируется следующим образом: заданы начальное

14

распределение температуры, граничные условия, источники тепла и требуется найти распределение температуры в последующие моменты времени.

Уравнение теплопроводности в декартовой системе координат имеет вид

∂u

∂

∂u

∂

∂u

∂

∂u

c(x,t)ρ(x,t)

=

a(x,t)

+

a(x,t)

+

a(x,t)

+ fs (x,t) ,

∂t

∂x1

∂x1

∂x2

∂x2

∂x3

∂x3

где u(x,t) температура среды в точке x=(x1, x2, x3) в момент времени t, ρ(x,t) – плотность, c(x,t) – теплоёмкость, a(x,t) – коэффициент теплопроводности и fs(x,t) – источник тепла. Вначале мы рассмотрим основные идеи построения разностных схем для задач теплопроводности на примере простого одномерного уравнения с постоянными коэффициентами.

11.2.1.Разностная схема для одномерного уравнения теплопроводности с постоянными коэффициентами.

Вслучае одномерной однородной среды предыдущее уравнение принимает вид

∂u	= κ	∂2u	, 0≤t≤T, 0≤x≤l,	(11.19)
∂t		∂x2

где κ=a/cρ - коэффициент температуропроводности. Конкретная физическая задача определяется начальным условием (начальным распределением температуры) u(x,0)=g(x) и граничными условиями, которые мы запишем в общем виде:

αu(0,t) + α		∂u	(0,t) = f (t),
1	2	∂x	1
α u(l,t) + α		∂u	(l,t) = f (t) ,
3	4	∂x	2

где g(x), f1(t), f2(t) – заданные функции.

Введём разностную сетку с узлами (xn,tk) и шагами h и τ. Аппроксимируем производную по времени разностью вперёд относительно точки (xn,tk), а производную ∂2u/∂x2 трёхточечным разностным оператором в момент времени tk. В результате

дифференциальное уравнение (11.19) заменяется разностной схемой
uk+1	−uk	= κ	uk	−2uk	+ uk
n	n		n+1	n	n−1	,	(11.20)
				h2
	τ

k=0,1,…; n=1,…,N-1,

un0 = g(xn ), n=0,…,N,

+разностная форма граничных условий.

15

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

uk+1 = uk + κτ (uk+ −2uk + uk− ),

n n h2 n 1 n n 1

k=0,1,…; n=1,…,N-1.

Так как распределение температуры в момент времени tk всегда известно, то мы получили формулу для расчёта распределения температуры в момент времени tk+1. Рассмотренный нами пример является частным случаем более общей схемы вида

k+1

k

k+1

k+1

k+1

k

k

k

un

−un

un+1

−2un

+ un−1

un+1

−2un

+ un−1

= κ γ

+ (1 − γ)

,

(11.21)

τ

h

2

h

2

0≤γ≤1, k=0,1,…; n=1,…,N-1,

которая включает в себя все возможные двухслойные схемы. Шаблон этой схемы показан на рис. 11.7.

Рисунок 11.7. Шаблон схемы (11.21).

Введение дополнительного параметра позволяет управлять свойствами разностной схемы, но это приводит к неявным разностным уравнениям. Решение таких уравнений мы обсудим ниже.

Рассмотрим погрешность аппроксимации схемы (11.21). Подставляя точное решение в (11.21) и используя соотношения (11.14), получим

∂u

1

2

2

2

1

4

u

∂ u

) = κ

∂ u

2 ∂

(xn ,tk ) +

τ

2

(xn ,tk ) +O(τ

γ

2

(xn ,tk+1) +

h

4

(xn ,tk+1)

+

∂t

2

∂t

∂x

12 ∂x

	2		1 2			4	u			4
	∂ u					∂
(1 − γ)		2 (xn ,tk ) +			h			4		+O(h ).
									(xn ,tk )
	∂x		12			∂x
			16

Проведём исследование аппроксимации в точке (xn,tk). Тогда величины определённые в точке (xn,tk+1) необходимо разложить в ряд Тейлора в точке tk по степеням τ. Окончательно мы приходим к следующему дифференциальному приближению

∂u

(x

,t ) +

1

τ

∂2u

(x

,t ) +O(τ2 ) = κ

∂2u

(x

,t ) + κγτ

∂3u

(x

,t ) +

∂t

2

∂t2

∂x2

∂x2∂t

n

k

n

k

n

k

n

k

κh2 ∂4u (xn ,tk ) +O(h4 ), 12 ∂x 4

из которого следует, что в общем случае схема (11.3) аппроксимирует уравнение (11.1) с точностью до членов O(τ+h2). Выбирая параметр γ определённым образом, можно повысить порядок аппроксимации схемы (11.11). Действительно, из уравнения (11.19) следует, что

∂2u

= κ2 ∂4u

и

∂3u

= κ

∂4u .

∂t2

∂x2∂t

∂x 4

∂x 4

Тогда дифференциальное приближение можно записать в виде

∂u

2

1

h

2

4

(xn ,tk ) −κ

∂ u

κτ

+ κγτ

+

∂ u

(xn ,tk ) +O(τ

2

+ h

4

) .

2

(xn ,tk ) = κ −

4

∂t

∂x

2

12

∂x

Отсюда следует, что при

γ =

1

−

h2

,

(11.22)

2

12κτ

выражение в квадратных скобках равно нулю и разностная схема (11.18) аппроксимирует уравнение (11.19) с точностью до величин O(τ2+h4).

Для анализа устойчивости схемы (11.21), подставим в неё решение вида (11.14) и

получим
λ −1	= κ(γλ + (1 − γ))	exp(iα) −2 + exp(−iα)		= −4	κ	(γλ + 1 − γ)sin2(	1	α)
τ		h2			h2		2

и окончательное выражение для собственного значения оператора перехода принимает вид

	λ(α) =	1 −4(1 − γ)β sin2(			1	α)		,
					2
		1 + 4γβ sin2(	1	α)
			2

0 ≤ α ≤ 2π,	β =	κτ .
		h2

В данном случае λ(α) принимает вещественные значения в интервале

1 −4(1 − γ)β ≤ λ(α) ≤1 . 1 + 4γβ

Из спектрального критерия устойчивости необходимо чтобы

17

1 −4(1 − γ)β ≥ −1 1 + 4γβ

и мы приходим к условию

	κτ	= β ≤	1		, 0 ≤ γ < 1	2 ,	(11.23)
	h2		2 −	4γ

β- любое>0, 12 ≤ γ ≤1.

Вчастности, для явной схемы (11.20) (γ=0) условие устойчивости приобретает вид β≤0.5. Как следует из выражения (11.22), эта схема будет иметь четвёртый порядок аппроксимации при ещё более строгом ограничении β=1/6. Напротив, схема (11.21) будет устойчива при любом значении β, если γ определяется выражением (11.22).

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

особого труда. Например, заданное при x=0 условие u(0,t)=f(t) записывается как
uk+1	= f (t	k+1	)	(11.24)
0

Граничные условия второго и третьего рода определяют значения теплового потока в граничных точках. В качестве примера можно привести следующие физические процессы, приводящие к условиям такого рода:

1)	воздействие теплового потока на вещество
		a ∂u (0,t) = f (t),	(11.25)
		∂x
2)	конвективный теплообмен
	a	∂u (0,t) = ko (u(0,t) −uc ),
		∂x
где uc – температура окружающей среды, ko – коэффициент теплоотдачи,
3)	тепловое излучение	∂u
			4
		a ∂x (0,t) = σ(u(0,t)) ,

где σ - коэффициент излучения.

Для простоты, рассмотрим граничное условие (11.25). Очевидное представление этого условия в разностной форме

	uk +1	− uk +1
a	1	0	= f (tk +1)
		h
		18

непригодно, так как оно имеет первый порядок аппроксимации по h, что значительно снижает точность приближённого решения (смотри пример 7.11). Для построения более точной аппроксимации, второго порядка по h, применим метод фиктивных областей. Для этого введём фиктивный слой с узлами (x-1,tk). Тогда, используя центральную разность, граничное условие (11.25) можно записать в виде

uk −uk

a 1 2h −1 = f (tk ), k=0, 1, … .

Значения температуры в граничном узле вычисляется из разностного уравнения (11.3) при n=0:

uk+1 −uk			uk+1 −2uk+1 + uk+1				+ (1 − γ)	uk −2uk + uk
0	0	= κ γ	1		2	−1		1		2		.
				0						0	−1
	τ			h					h

Исключение из этих двух выражений, фиктивных значений с номером n=-1 даёт разностную форму граничного условия (11.7), которую мы запишем в виде

	(1 + 2βγ)u0k+1 −(2βγ)u1k+1 =
	(11.26)
(1 −2β(1 − γ))u0k	+ (2β(1 − γ))u1k − 2βah (γf (tk+1) + (1 − γ)f (tk )).

В случае явной схемы (γ=0), это выражение непосредственно даёт формулу для вычисления температуры в граничном узле.

В общем случае, когда γ≠0, разностная схема (11.3) является неявной, так как при каждом значении n разностное уравнение включает в себя неизвестные в момент времени tk+1 значения температуры в трёх узлах сетки. Для того чтобы вычислить решение в момент времени tk+1, необходимо рассмотреть совокупность всех разностных уравнений. Перепишем выражение (11.21), группируя неизвестные величины в левой части, а известные в правой:

	−(βγ)unk−+11 + (1 + 2βγ)unk+1 −(βγ)unk++11 =	(11.27)
	(1 −2β(1 − γ))unk + β(1 − γ)(unk−1 + unk+1 )
		= dn ,
	n=1,…,N-1.

Введём обозначения

w = (u0k+1,...,uNk+1 )T , d = (d0,...,dN )T .

Тогда предыдущее выражение представляет собой систему линейных уравнений вида

Aw=d,

(11.28)

с трёхдиагональной матрицей

19

s	0	r
		0
		s1	r1		0
p1
		.	.	.
A =							.
			pn	sn	rn
		0		.	.	.
					pN
						sN

Коэффициенты этой матрицы следуют из (11.27) и равны

pn = −βγ , sn = 1 + 2βγ , rn = −βγ,

n=1,…,N-1,

а коэффициенты первого и последнего уравнений определяются граничными условиями. Например, для граничного условия (11.24) s0=1, r0=0 и d0=f(tk+1), а для граничного условия

(11.26)

s0 = 1 + 2βγ , r0 = −2βγ ,

d = (1 −2β(1 − γ))uk + (2β(1 − γ))uk − 2βh (γf (t + ) + (1 − γ)f (t )).

0 0 1 a k 1 k

Так как элементы матрицы А удовлетворяют условию |sn|>|pn|+|rn|, то для решения системы уравнений можно применить метод прогонки (смотри п. 3.6.1). Таким образом, для выполнения одного шага по времени в схеме (11.21) требуется большее количество операций по сравнению с явной схемой (11.20). С другой стороны, в неявной схеме мы можем выбрать больший шаг по времени, поэтому при заданной длительности процесса потребуется меньшее число шагов и общее количество операций будет сравнимо или даже меньше, чем для явной схемы. Для вычисления одного временного шага по явной схеме (11.20) необходимо 4N операций. Для вычислений по неявной схеме (11.21) сначала необходимо определить вектор d (4N операций), а затем решить систему уравнений ( 8N операций), что в сумме даёт приблизительно 12N операций на один шаг по времени. Поэтому, если в неявной схеме шаг по времени выбрать из соотношения β≈1.5, то общее количество операций будет сравнимо с количеством операций требуемых для вычисления по явной схеме.

________________________________________________________________________________

Пример 11.1 (остывание бесконечной пластины)

20