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

uk+1 = 2uk

−uk−1 + 2r2(uk

−uk −hf (t )),

0

0

0

1

0

k

k=1, 2, ….

Исследуем свойства схемы «крест». Используя последнее соотношение из (11.11),

получим следующее выражение

∂2u

(x

,t

) −c2 ∂2u (x

,t ) =

∂t2

n

k

∂x2

n

k

.

4

4

1

2

∂ u

2 2 ∂ u

4

4

+O (h

)

−τ

4

(xn

,tk ) + c h

4 (xn

,tk )

+ τ

12

∂t

∂x

Это значит, что разностная схема (11.52) аппроксимирует дифференциальную задачу (11.51) со вторым порядком по h и τ. Далее, подставляя в разностную схему решение в виде (11.14), получим уравнение для собственных значений оператора перехода

λ2 −2λ + 1 = r2 (exp(iα) −2 + exp(−iα)) = −4r2 sin2(21 α) .

Решение этого уравнения имеет вид

			2			(			2		2
	2	2					2	2			)
λ(α) = 1 −2r		sin	(		α) ±	1 −2r		sin	(		α) −1 ,

и, согласно критерию фон Неймана, схема «крест» будет устойчива при r≤1.

________________________________________________________________________________

Пример 11.4 (схема «крест»)

Рассмотрим задачу (11.51) с нулевыми начальными условиями и граничными условиями

		f t(1 −t)n , 0 ≤t ≤1
			0				,
u(0,t) =				>	1
		0 , t
			u(l,t) = 0 ,				(11.53)
где
f	= (1		1	+ n n		, n = 10 .
		+ n)
0				n

Решение, полученное с помощью схемы (11.52) показано на рисунке 11.13. Как показывают эти результаты, данная схема сохраняет амплитуду, но искажает форму волны. Более того, появляются высокочастотные осцилляции, отсутствующие в точном решении.

41

Рисунок 11.13. Решение уравнения (11.51) с граничными условиями (11.53) в момент времени t=1 (l=1, c=1); - точное решение, - приближённое решение (h=0.05,

τ=0.0333, r=2/3).

________________________________________________________________________________

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

11.3.2. Диссипация и дисперсия сеточного волнового решения.

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

u(x,t)=u0exp(i(ωt-kx)),

где ω=2πν - круговая частота, k=2π/λ - волновое число. Если мы подставим это решение в какое-нибудь линейное волновое уравнение, то получим зависимость ω=ω(k), которое называется дисперсионным соотношением. Если ω принимает вещественные значения, то амплитуды гармоник не изменяются со временем. В противном случае, когда ω принимает

42

комплексные значения, волна будет затухать (диссипировать) как exp(-(Imω)t)=exp(-γt). Отношение

	Re ω	= cf
	k

называется фазовой скоростью, то есть скоростью, с которой движется фаза или узел отдельной гармоники. Величина

d

dk Re ω = cg

называется групповой скоростью. Она может интерпретироваться как скорость волнового пакета, состоящего из гармонических волн с близкими волновыми числами. Передача энергии в волне осуществляется с групповой скоростью. Если фазовая (групповая) скорость зависит от k, то гармоники с разными волновыми числами распространяются с разными скоростями. Это явление называется дисперсией. Для уравнения (11.51) нетрудно получить, что ω=ck и, следовательно, cf=cg=c. Это значит, что точное решение уравнения (11.51) представляет собой волну, которая распространяется без затухания и дисперсии. Другой пример: для линеаризованного уравнения Бюргерса

∂u	+ c	∂u	= μ	∂2u	,
∂t		∂x		∂x2

зависимость ω(k) имеет следующий вид

ω(k)=ck+iμk2.

Это соотношение описывает затухающую недиспергирующую волну. В частности, при μ=0 мы получаем уравнение переноса (11.1), для которого справедливо соотношение cf=cg=c и коэффициент затухания равен нулю.

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

Анализ диссипации и дисперсии сеточного решения мы будем проводить на решении unk = u exp(i(ωtk −kxn )).

43

Если подставить это решение в разностную схему, то мы получим некоторое соотношение вида ω=ω(k,τ,h). Сравнивая это соотношение с соответствующим дисперсионным соотношением для дифференциального уравнения, которое эта схема аппроксимирует, мы можем анализировать свойства данной схемы. В качестве примера рассмотрим свойства схемы «крест». Подставим приведённое выше решение в схему (11.52) и получим следующее

выражение exp(iωτ)-2+exp(-iωτ)=r2(exp(-ikh)-2+exp(ikh)) из которого следует, что sin(12 ωτ) = r sin(21 kh).

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

ω = τ2 arcsin(r sin(21 kh)).

Тогда групповая скорость волны распространяющейся на сетке, определяется как

c	= dω = c	cos(	1	kh)			,
			2
m	dk	1 −r2 sin2 (			1	kh)
					2

то есть, в отличии от уравнения (11.1), возникает дисперсия. Только когда r=1, групповая скорость приближённого решения совпадает с групповой скоростью точного решения.

Минимальная длина волны, которая может быть представлена на сетке, составляет λ=4h (kh=π/2). Понятно, что гармоники с такой длиной волны сильнее всего подвержены влиянию дисперсии. Для этого случая, зависимость безразмерного параметра a=сm/c от параметра Куранта показана на рис. 11.14.

44

Рисунок 11.14. Зависимость параметра a=сm/c от параметра Куранта для схемы «крест»;

- безразмерная групповая скорость волны cg/c=1 для дифференциального уравнения;

- параметр a при kh=π/2.

Зависимость, показанная на этом рисунке, представляет наихудшую ситуацию. Для длин волн много больших пространственного шага, дисперсия приближённого решения уменьшается.

11.3.3. Схема Лакса-Вендроффа.

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

		∂s	+ A		∂s	= o ,				(11.54)
		∂t
					∂x
		0 ≤ x ≤ l, 0 ≤ t ≤ T,
где
u				0 1/					0
							ρ
		, A =
s =				2				, o =		.
p				ρc					0
					0

Здесь u – скорость частиц среды, p – давление, ρ - плотность среды, и c – скорость звука.

Уравнения (11.54) дополняются начальными условиями

u(x, 0) = g1(x) , , p(x, 0) = g2(x)

и граничными условиями

α1u(0,t) + α2p(0,t) = f1(t) , . α3u(l,t) + α4p(l,t) = f2(t).

В дальнейшем, мы будем рассматривать построение разностных схем для системы (11.54). Предположим, что решение в момент времени t известно. Тогда решение в следующий

момент времени t+ t можно получить с помощью разложения функции s(x,t+ t) в ряд Тейлора в точке (x,t) по степеням t, то есть

s(x,t +	t) = s(x,t) +		∂s	(x,t) +	1		2	∂2s	(x,t) + ... .
		t	∂t		2	(	t)	∂t2

Производные по t заменяются производными по x путём использования системы (11.4):

45

	∂s		= −A				∂s		,
	∂t					∂x					.
∂2s	= −		∂	A	∂s			= A2		∂2s
∂t2		∂t			∂x
										∂x2

Подставим эти выражения в наше разложение и получим

s(x,t +

t) = s(x,t) −

tA

∂s

(x,t) +

1

2

A

2 ∂2s

(x,t) + ...

(11.55)

2

(

t)

∂x

∂x2

Теперь введём разностную сетку с узлами (xn,tk) и шагами h и τ. Затем построим разностную аппроксимацию выражения (11.5) на этой сетке, используя приближения второго порядка для производных по x (смотри (3.2)) и полагая τ= t. Это приводит к схеме Лакса-Вендроффа

k+1	k τ		k	k	1	τ 2	2	k	k k
sn	= sn −		A(sn+1 −sn−1 )+				A	(sn+1	−2sn + sn−1 ),	(11.56)
		2h
					2 h

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

или в развёрнутом виде

unk+1

pnk+1

	= unk	−	r	(pnk+1 − pnk−1 )
			2cρ
	= pnk	− rc2ρ(unk+1 −unk−1 )

n=1, …, N-1;

+ 1 r2 (uk+ −2uk + uk− ),

2 n 1 n n 1

+	1 r2	(pnk+1 −2pnk	+ pnk−1 ).
	2

k=0, 1, … .

Из процедуры построения следует, что схема (11.56) имеет второй порядок по h и τ. Рассмотренная нами схема удобна, когда матрица A постоянна. В случае переменной матрицы, проще использовать другой подход, который исключает вычисление производных от матрицы A. Сначала в центрах прямоугольных ячеек на плоскости (x,t) вычисляются промежуточные значения (рисунок 11.15):

	snk++1/1/22	=	1	(snk +1 + snk )−	τ	(Ank +1snk +1 −Ank snk ),	(11.57)
			2		2h

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

46

Рисунок 11.15. Вспомогательные узлы «×» для вычисления промежуточного решения.

Здесь

			0
	k					1/ ρ(xn ,tk )
A
	n	=		2			.
					(xn ,tk )	0
		ρ(xn ,tk )c

Окончательные значения получается затем из уравнения

snk+1 = snk −	τ (Ank++1/1/22snk++1/1/22 −Ank+−1/1/	22snk+−1/1/	22 ),	(11.58)
	h

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

Выражения (11.57) и (11.58) называются двухшаговой схемой Лакса-Вендроффа. Легко проверить, что в частном случае, когда матрица A постоянна, эта схема переходит в схему

(11.56).

Перейдём к анализу устойчивости схемы (11.6). Только теперь мы имеем систему

уравнений, поэтому вместо (11.14) зададим решение в виде skn = λk (α)exp(iαn)s ,

где s* - некоторый постоянный вектор. Подстановка этого решения в разностную схему (11.56) даёт следующее выражение:

	λs = s −				τ	(exp(iα) −exp(−iα))As +
					2h
								.
	1		τ 2			2	s
	2			(exp(iα) −2 + exp(−iα))A
		h

После некоторых преобразований, мы получим стандартную задачу на собственное значение

47

	τ		τ 2		2			= Bs			.
I −(i		sin(α))A −			(1 −cos(α))A	s				= λs
	h		h

Решение характеристического уравнения det(B-λI)=0 даёт следующее выражение для

собственного значения оператора перехода

λ(α) = 1 −r2(1 −cos(α)) ± ir sin(α) .

Тогда из критерия фон Неймана следует, что разностная схема (11.6) будет устойчива при r≤1.

Анализ диссипативных и дисперсионных свойств разностных схем для системы

уравнений (11.4) мы будем проводить на решении

skn = s exp(i(ωtk −kxn )) .

Если подставить это решение в разностную схему (11.6),то мы получим систему уравнений

τ

τ

2

2

= Ds

= o .

(exp(iωτ) −1)I −i

sin(kh) A +

(1 −cos(kh))A

s

h

h

Из условия det(D)=0 следует дисперсионное соотношение

exp(iωτ) = 1 −r2 (1 −cos(kh))+ ir sin(kh)

и, следовательно, ω принимает комплексные значения, что говорит о диссипации сеточного решения. Введём обозначение ω=ω0+iγ, и приравнивая вещественные и мнимые части в

предыдущем выражении, получим

exp(−γτ)cos(ω0τ) = 1 −r2 (1 −cos(kh)), . exp(−γτ)sin(ω0τ) = r sin(kh).

Отсюда нетрудно получить коэффициент затухания

γτ = −	1		−r	2	+ r	2	2	+ r	2	2
	2	ln (1					cos(kh))			sin	(kh)

и групповую скорость волны на сетке

c	=	dω		0	= c				r2 + (1 −r2 )cos(kh)								.
						(						2
m
			d k								)
							1 −r	2	+ r	2	cos(kh)		+ r	2	2	(kh)
															sin

На рисунках 11.16 и 11.17 представлены зависимости параметра диссипации d=γτ и групповой скорости от параметра Куранта.

48

Рисунок 11.16. Зависимость параметра диссипации d=γτ от параметра Куранта для схемы Лакса-Вендроффа (kh=π/2).

Рисунок 11.17. Зависимость параметра a=сm/c от параметра Куранта для схемы ЛаксаВендроффа; - безразмерная групповая скорость волны cg/c=1 для дифференциального уравнения; - параметр a при kh=π/2.

49

u0k+1

Для аппроксимации граничных условий со вторым порядком по h и τ, можно использовать метод фиктивных областей. Пусть, например, задано граничное условие при x=0: u(0,t)=f(t). В разностном виде это граничное условие записывается очень просто:

= f (tk+1) . Для вычислений по схеме (11.56) нам необходимо также значение давления в граничном узле. Для определения этого значения введём фиктивный слой с координатами

(x-1,tk) и запишем второе уравнение (11.56) в виде

p0k+1 = p0k − rc2ρ(u1k −u−k 1 )+ 21 r2 (p1k −2p0k + p−k 1 ).

Для исключения фиктивных значений u−k1 и p−k1 используем очевидное приближение

0.5(u1k + u−k 1) = u0k = f (tk ) .

Тогда первое уравнение из (11.56) примет вид

2rρc (p1k − p−k 1 ) = u0k −u0k+1 = f (tk ) − f (tk+1) .

После несложных преобразований, мы получим формулу для вычисления второго граничного условия

pk+1 = pk −ρc r (uk − f (t + ))+ r2 (pk − pk ).

0 0 1 k 1 1 0

В случае, когда на границе задано давление, аналогичная процедура применяется для определения скорости в граничной точке.

________________________________________________________________________________

Пример 11.5 (схема Лакса-Вендроффа)

Рассмотрим задачу (11.54) (которая эквивалентна задаче (11.51)) с нулевыми начальными условиями и граничными условиями (11.53). Приближённое решение, полученное с помощью схемы Лакса-Вендроффа (11.56) показано на рисунке 11.18.

50