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

c2 (un,m ,un,m−1 )un,m−1 + c1 (un,m ,un−1,m )un−1,m −(c1 (un,m ,un−1,m )+ c1 (un,m ,un+1,m )+ c2 (un,m ,un,m−1 )+ c2 (un,m ,un,m+1 ))un,m +

c1 (un,m ,un+1,m )un+1,m + c2 (un,m ,un,m+1 )un,m+1 − f (xn ,ym ) = 0,

n=1, …, N-1 ; m=1, …, M-1,

где

c (u ,u + ) = 1 a (1 (u + + u )),

1 n,m n 1,m h2 2 n 1,m n,m

x .

c (u ,u + ) = 1 a (1 (u + + u )).

2 n,m n,m 1 hy2 2 n,m 1 n,m

При каждом фиксированном значении n и m мы получаем одно уравнение из системы. Далее эта система решается каким-нибудь итерационным методом (смотри главу 5).

________________________________________________________________________________

Пример 11.8 (метод SOR)

Рассмотрим уравнение (11.84) при a=1, b=0 в квадрате {0≤x≤1, 0≤y≤1} с граничными условиями

u (0,y)=u (1,y)=0 , u (x,0) = 0			, u (x,1) = x(1 −x) ln(2).	(11.99)
Правая часть уравнения
f (x,y)= −	x(1-x)	−2	ln(1+y).
	2
(1+y)

Зададим разностную сетку с N=M=16. При начальном приближении un(0),m = 0 , (n=0,…,N ; m=0,…,M) и ε=10-5 итерационный процесс сходится за 40 итераций. Приближённое решение показано на рисунке 11.28.

81

Рисунок 11.28. Приближённое решение уравнения (11.1) с граничными условиями (11.16),

метод SOR (a=1, b=0, N=M=16, ε=10-5).

Точное решение этой задачи – u(x,y)=x(1-x) ln(1+y). Для оценки ошибки приближённого решения введём u(e)={u(xn,ym)} - точное решение вычисленное в узлах сетки и u={un,m} - приближённое решение задачи. Тогда для решения показанного на рисунке 11.28

e

=

u −u(e)

2

≈ 8 10−5 .

a

u(e)

2

______________________________________________________________________________

11.4.3. Применение быстрого преобразования Фурье.

Если область определения задачи представляет прямоугольную односвязную область, коэффициенты уравнения (11.84) постоянны и граничные условия в одном из направлений имеют специальный вид, то можно применить метод, основанный на преобразовании Фурье. Мы рассмотрим основные идеи этого метода на примере одной задачи, а затем обсудим все остальные случаи, в которых этот метод можно применять.

Пусть коэффициенты уравнения (11.92) постоянны: a(x,y)=a и b(x,y)=b, область определения – прямоугольник {0≤x≤lx, 0≤y≤ly}, а граничные условия имеют вид

82

	u(0,y) = u(lx ,y) = 0,
αu(x, 0) + α			∂u	(x, 0)=g		(x), .
1		2	∂y		2
α u(x,l		) + α	∂u	(x,l	)=g		(x).
3	y	4	∂y	y	4

Тогда разностная схема (11.11) примет вид

c1un,m−1 + c2un−1,m −c3un,m + c2un+1,m + c1un,m+1 = fn,m ,	(11.100)
n=1, …, N-1 ; m=1, …, M-1,

u0,m = uN ,m = 0 , m=0, …, M,

2c α h

2c h

1 1

y

1

y

c2un−1,0

− c3

−

un,0

+ c2un+1,0 + 2c1un,1un,1

= fn,0 +

g2

(xn ), α2 ≠ 0

α2

α2

1

=

, α = 0

u

n,0

g

(x

n

)

α1

2

2

n=1, …, N-1,

2c α h

2c h

1 3

y

1

y

2c1un,M −1 + c2un−1,M − c3

+

un,M + c2un+1,M = fn,M −

g4

(xn ), α4

≠ 0

α4

α4

1

=

, α

=

u

g

(x

n

)

0

n,M

α3

4

4

m=1, …, M-1.

Здесь

	c	=		a	,
			h2
	1
				y
	c	=		a		, ,
				h2
	2
				x

c3 = 2(c1 + c2 ) −b.

(11.101)

,(11.102)

,(11.103)

и для вывода граничных условий при y=0 и y=ly применялся метод фиктивных областей. Для этого типа граничных условий решение можно представить в виде конечного ряда

Фурье:

		1	N −1
	un,m =		∑vk,m sin(πkn N ), m=0, …, M.	(11.104)
		N k=1

Правую часть уравнения также можно представить в виде

83

		1	N −1
	fn,m =		∑dk,m sin(πkn N ), m=0, …, M	(11.105)
		N k=1
	и коэффициенты dk,m вычисляются по формуле
		N −1
	dk,m = ∑fn,m sin(πkn N ), m=0, …, M.			(11.106)
		n=1

В дальнейшем, операцию вида (11.106) мы будем обозначать как dk,m = Fn (fn,m ), где n указывает на параметр, по которому происходит суммирование.

Подстановка выражений (11.104) и (11.105) в уравнения (11.100) и граничные условия (11.102) и (11.103) даёт следующую систему уравнений

γ

(1)v

k,0

+ γ(2)v

k,1

= d

,0

,

k

k

k

c1vk,m−1

−(c3 −2c2 cos(πk N ))vk,m + c1vk,m+1 = dk,m ,

(11.107)

m=1, …, M-1,

γ(3)v

k,M −1

+ γ

(4)v

k,M

= d

,

k

k

k,M

где

1

(x

)

F

g

, α = 0

α1

n

(

2

n

)

2

dk,0

,

=

2h c

y

1

F

f

+

g

2

x

, α ≠ 0

n

n,0

α2

(

n )

2

1

(x

)

F

g

,

α

=

0

α3

n

(

4

n

)

4

dk,M

.

=

2h c

y

1

F

f

−

g

x

,

α

≠

0

n

n,M

α4

4 (

n )

4

Параметры γk(1), ... , γk(4) определяются как

1

, α = 0

2

(1)

2c1α1hy

,

γk

=

cos(πk N )+

, α2 ≠ 0

−c3 + 2c2

α

2

0

, α =

0

(2)

2

=

,

γk

,

α ≠

0

2c

1

2

84

	0		, α = 0
(3)			4
					,
γk	=		, α	≠ 0
	2c
		1	4

	1			, α = 0
					4
(4)			2c1α3hy				.
γk	=	+ 2c2 cos(πk N )−
					,	α4 ≠ 0
	−c3		α
				4

Система (11.107) представляет собой систему уравнений с трёхдиагональной матрицей. При b≤0 эта матрица будет иметь свойства диагонального преобладания и для решения этой системы можно применить метод прогонки (смотри параграф 3.6).

Эффективность данного метода определяется эффективностью вычисления преобразования Fn(…). Рассматриваемый нами метод имеет широкое применение благодаря тому, что при N=2p существует алгоритм, который позволяет вычислить преобразование Fn(…) за O(log2(N)) операций. Этот алгоритм носит название быстрое преобразование Фурье (FFT, Fast Fourier Transform). В настоящее время процедура вычисления FFT является стандартной и входит в состав всех инструментальных сред и библиотек программ.

Суммируя всё изложенное выше, процедура решения уравнения (11.100) состоит в следующем:

1)проводят дискретное преобразование Фурье (11.106) правой части, используя FFT,

2)для каждого значения k=1,…,N решают системы (11.107) и находят vk,m,

3)используя FFT, проводят обратное преобразование (11.104) сеточной функции vk,m для получения решения un,m.

Полное число операций необходимое для выполнения этой процедуры составляет Ns(2.5log2(Ns)+8), что значительно меньше, чем, например, в методе SOR.

Мы рассмотрели построение алгоритма, когда в граничных точках при x=0 и x=lx решение равно нулю. Можно легко обобщить данный метод и для граничных условий (11.4). Запишем уравнение (11.17) при n=1:

c1u1,m−1 + (c2u0,m )−c3u1,m + c2u2,m + c1u1,m+1 = f1,m .

Согласно граничному условию (11.101) слагаемое в скобках равно нулю. Если задано граничное условие u0,m=g1(ym), то мы можем перенести слагаемое c2u0,m в правую часть и при этом левая часть уравнения будет такой же как и для условия (11.18). Таким образом, применяя рассмотренный выше алгоритм с модифицированными правыми частями

f1,m → f1,m −c2g1(ym ),

fN −1,m → fN −1,m −c2g3(ym ).,

85

мы получим решение уравнения с граничными условиями (11.4). Естественно, что в конце вычислений мы должны положить u0,m=g1(ym) и uN,m=g3(ym).

Построение метода для других типов специальных граничных условий производится аналогичным образом. Дискретное преобразование Фурье в общем виде можно представить

в следующей форме

un,m = 1 ∑N vk,mϕk,n , m=0, …, M,

N k=0

где

ϕk,n = sin(πkn N ) - для граничных условий (11.87),

		1/2, k	= 0 , k = N
ϕk,n = cos(πkn
				- для граничных условий (11.88) и
	N	)	< k < N
		1 , 0

ϕk,n = exp(−2πikn N ) - для граничных условий (11.86).

Все эти преобразования вычисляются с использованием алгоритма FFT.

Для граничных условий Неймана (11.5) правая часть должна быть модифицирована следующим образом:

f0,m → f0,m −2c2hxg1(ym ), fN ,m → fN ,m −2c2hxg3(ym )..

Таким образом, мы полностью обсудили методы решения уравнения (11.84) с постоянными коэффициентами, граничными условиями (11.6)-(11.8) по x и граничными условиями общего вида по y. Решение уравнения с граничными условиями (11.89)-(11.91) по y и граничными условиями общего вида по x осуществляется совершенно аналогичным образом.

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

________________________________________________________________________________

Пример 11.9 (метод FFT)

Применим метод FFT к задаче из примера 11.8. Тогда система (11.107) примет вид

vk,0 = 0 ,

vk,m−1 −2(2 −cos(πk N ))vk,m + vk,m+1 = h2 dk,m ,

a

m=1, …, M-1,

86

vk,M = dk,M .

На сетке N=M=16 ошибка полученного решения

e

=

u −u(e)

2

≈ 5 10−5 ,

a

u(e)

2

а время вычислений примерно в 10 раз меньше по сравнению с методом SOR.

________________________________________________________________________________

11.4.4. Метод установления.

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

модельной задачи
	2			2
	∂ u		+	∂ u			= f (x,y),
a		2			2
	∂x			∂y
	{0≤x≤lx				, 0≤y≤ly},			(11.108)

u(0,y) = g1(y) , u(lx ,y) = g3(y) , u(x, 0) = g2(x) , u(x,ly ) = g4(x) .

Вначале мы должны сформулировать соответствующую нестационарную задачу. В

нашем случае она может быть записана в следующей форме:
∂v		2			2
		∂ v			∂ v			− f (x,y),
	= a		2	+		2
∂t		∂x			∂y
	{0≤x≤lx				, 0≤y≤ly},				(11.109)

v(x,y, 0) = g(x,y) ,

v(0,y) = g1(y) , v(lx ,y) = g3(y) , v(x, 0) = g2(x) , v(x,ly ) = g4(x),

где функция g(x, y) задаёт произвольные начальные условия. Так как граничные условия не зависят от времени, то естественно ожидать, что решение v(x,y,t) с течением времени будет

87

меняться всё медленнее и медленнее. Поэтому в пределе t→∞, это решение будет

стремиться к решению задачи (11.108), то есть

lim v(x,y,t) = u(x,y) .

t→∞

Следовательно, вместо стационарной задачи (11.108) можно решать нестационарную задачу (11.109) до некоторого момента времени t, когда решение v(x,y,t) перестанет изменяться в пределах заданной точности. Это и есть основная идея метода установления.

В соответствии с нашим рассмотрением нам теперь нужно построить разностную схему для задачи (11.109). Решение такого типа задач мы обсуждали в параграфе 11.2. Вначале рассмотрим явную схему (11.43), которая применительно к нашей задаче имеет вид

	vnk,+m1 −vnk,m	= a (					xvnk,m + yvnk,m )− f (xn ,ym ) ,	(11.110)
	τ
n=1, …, N-1 ; m=1, …, M-1 ; k=0, 1, …,
				vn0,m = g(xn ,ym ),
	v0,k+m1 = g1(ym ) , vNk+,m1 = g3(ym ) ,
				m=0, …, M,
	vnk,0+1 = g2(xn ) , vnk,+M1 = g4(xn ) ,
				n=0, …, N.
Здесь и далее используются следующие обозначения:
	xvnk,m =			1		(vnk+1,m −2vnk,m + vnk−1,m ),
				h2
				x
	yvnk,m =		1		(vnk,m+1 −2vnk,m + vnk,m−1 ).
			h2
				y
Тогда разностную схему для уравнения (11.108) можно записать в виде
	a ( xun,m +						yun,m ) = f (xn ,ym ),	(11.111)
n=1, …, N-1 ; m=1, …, M-1 ; k=0, 1, …,

u0,k+m1 = g1(ym ) , uNk+,m1 = g3(ym ) ,

m=0, …, M,

unk,0+1 = g2(xn ) , unk,+M1 = g4(xn ) ,

88

				n=0, …, N.
Решая	задачу (11.110), мы будем последовательно вычислять сеточную функцию
{	n,m }	, k=0,1,… . Когда эта функция уже не изменяется в пределах заданной точности,
vk = vk
то есть
				vk+1 − vk						≤ ε ,
					vk+1

процесс установления завершается и vk+1 есть приближённое решение стационарной задачи. Естественно возникает вопрос: всегда ли сеточная функция vk+1 будет приближаться к решению стационарной задачи? Для того чтобы ответить на этот вопрос нам необходимо провести анализ поведения разности

ek	= vk	−u	n,m	,	(11.112)
n,m	n,m

при k→∞. Задачу (11.111) можно записать в следующей эквивалентной форме

un,m −un,m = a ( xun,m + yun,m )− f (xn ,ym ) ,

τ

un,m = un,m ,

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

+ граничные условия. Вычитая эти уравнения из (11.110), получим

ek+1		−ek	= a (						yenk,m ),
	n,m	n,m			xenk,m +						(11.113)
		τ
e0		= v0	−u	n,m	= g(x	n	,y	m	) −u	n,m	,
	n,m	n,m

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

e0,k+m1 = eNk+,m1 = 0 ,

m=0, …, M,

enk,0+1 = enk,+M1 = 0 ,

n=0, …, N.

Для дальнейших оценок будем полагать, что N=M. Предположим, что функция g(x,y) удовлетворяет граничным условиям задачи (11.108). Тогда начальную ошибку можно представить в виде конечного ряда Фурье:

89

N −1 N −1
en0,m = ∑∑qr,s sin(πrn N )sin(πsm N ).
s=1 r=1
Подставляя это выражение в уравнение (11.113), можно получить, что
N −1 N −1
enk,m = ∑∑(λr,s )k qr,s sin(πrn N )sin(πsm N ),
s=1 r=1
где
λr,s = 1 −	4aτ			sin2 (πr 2N )− 4a2τ sin2 (πs 2N ).							(11.34)
	2
	hx								hy
Тогда относительную величину ошибки можно оценить как
			ek						≤(maxr,s | λr,s |)k	= αk .	(11.35)
			e0

Поэтому, для того чтобы ||ek||/||e0|| стремилось к нулю, когда k→∞, необходимо чтобы α<1. Так как α=α(τ), то при некотором значении τ величина α будет наименьшей и значит, будет наблюдаться самое быстрое уменьшение ошибки. Из соотношения (11.34) видно, что

λleft≤λr,s≤λright, где

λleft = 1 −2(β1 + β2 )cos2 (π2N )≤1,

λright = 1 −2(β1 + β2 )sin2 (π2N )≤1,

β =	2aτ	, β =	2aτ	.
	h2
1		2	h2
	x		y

Увеличение τ, начинаяя с τ=0, приводит к тому, что эти крайние значения сдвигаются влево на числовой оси. При некотором значении τ=τ* возникает ситуация, когда λleft и λright располагаются симметрично относительно нуля, то есть

-λleft=λright.

(11.116)

При этом α(τ*)=|λleft|=λright. Это и есть наименьшее из всех возможных значений α, так как при дальнейшем увеличении τ, λright будет приближаться к нулю, а λleft – удаляться от нуля и α(τ)=|λleft|>α(τ*). Таким образом, оптимальное значение τ определяется из условия (11.116) и равно

τ =	(hxhy )2
		.
	2a (hx2 + hy2 )

Тогда

α = 1 −2 sin2 (π2N ).

90