Конспекты. Blackboard / Тема 9.2 Учебное пособие по конечно-разностным методам решения краевой задачи

Суммируя по n=0; j; приходим к оценке решения задачи (4)

j
X
kvj+1kC kv0kC + k'nkC :	(34)

n=0

Схема с весами является монотонной при условии (33). Чисто неявная схема ( = 1) монотонна при любых h и . Схема Кранка-Николсона ( = 0:5) монотонна при условии h2, а явная схема

-ïðè 0:5h2 .

8.Метод энергетических неравенств. Одним из весьма общих и весьма эффективных спо-

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

Пользуясь тождествами [11, 12]

1				1
v=		(^v + v)		vt; v^=		(^v + v) +		vt;
	2		2		2		2
v^ + (1 )v=( 0; 5) vt + 0; 5(^v + v);
перепишем (4) в виде
vt ( 0; 5) vt 0; 5 (^v + v)=';									(35)
		v(x; 0)=0;		v(0; t)=v(1; t)=0:

Умножим уравнение (35) на 2 vt=2(^v v) и просуммируем полученное равенство по внутренним узлам xm=mh сетки !h:

2 kvt k2 2( 0; 5) 2 ( vt; vt) (^v + v); v^ v =2 ('; vt ):

(36)

Пользуясь разностной формулой Грина (см. с.29)

( V; W )=(Vxx; W )= (Vx; Wx]; V0=W0=0; VM =WM =0;

ïðè V =vt; W =vt è V =^v + v, W =^v v соответственно, учитывая затем, что v0=vM =0, будем иметь

( vt; vt)= kvtx]j2;
( (^v + v); v^ v)= (^vx + vx; v^x vx]= (kv^x]j2 kvx]j2):
Подставив это выражение в (36), получим энергетическое тождество
2 kvtk2 + ( 0; 5) kvtx ]j2	+ kv^x]j2=kvx]j2 + 2 ('; vt ):						(37)
Оно справедливо при любых . Предположим, что 0. Рассмотрим выражение
J =kW k2 + ( 0; 5) kWx ]j2;				ãäå W = vt:
Тогда с учетом оценки (см. с. 32)	4
kvx]j2			kvk2				(38)
		h2
имеем					h2
J =kW k2 + ( 0) kWx]j2 + ( 0 0; 5) kWx ]j2 kW k2						kWx]j2 0
					4
в силу (38). Отсюда и из (37) следует энергетическое неравенство
kv^x]j2 kvx]j2 + 2 ('; vt )				ïðè 0:

41

Åñëè '=0; v решение задачи (13), то kvxn+1]j : : : kvx0 ]j; т. е. схема (12) при 0 устойчива по

o

начальным данным в норме kvx]j, являющейся сеточным аналогом нормы W21:

Далее выясняем условия устойчивости по правой части. Пусть

;

=

1

(1 ")h2

;

0<"

1:

(39)

"

"

2

4

Покажем, что

J "kW k

ïðè ":

(40)

В самом деле,

J =kW k2 + ( ") kWx]j2 + ( " 0; 5) kWx]j2

k

W

k

2

(1 ")h2

k

W

]

2

k

W

k

2

(1

")

k

W ]

2 ="

W

k

2

:

4

x

j

j

k

Подставляя (40) в (37), получим энергетическое неравенство

2 "kvt k2 + kv^x]j2 kvx]j2 + 2 ('; vt );

":

(41)

Воспользуемся теперь неравенство Коши-Буняковского и "-неравенством:

2 ('; vt ) 2 k'k kvtk 2 "kvt k2 +				k'k2 :	(42)
			2"
После подстановки (42) в (41) будем иметь
kvxn+1]j2 kvxn]j2 +		k'nk2:
	2"

Просуммируем по n=0; j и учтем, что v0=0:

j

kvxj+1]j2 21" X k'nk2:

n=0

Òàê êàê kvkC 0; 5kvx ]j , òî

v

k

vj+1

k

C

1

j

'n

2;

":

2p2" un=0

k

k

uX

t

Применим эту оценку к задаче для функции z:

	j+1		n		T
kz		kC c 0 n j k		k; c=	2p2"		; ":
		max

Отсюда следует равномерная сходимость схемы с весами

kvj uj kC	8	M (h2 + ) ïðè >0; 5; u 2 C24;
		M (h2 + 2)	ïðè =0; 5; u 2 C34;
	<	M (h4 + 2)	ïðè = ; u	2	C36:
	:

Для явной схемы ( =0) из (39) не следует равномерная сходимость при условии h2=2: Однако

для нее можно непосредственно получить оценку

j 1

X

kvj kC kv0kC + k'nkC ; h2=2;

n=0

42

j

òàê ÷òî kzj+1kC P k'nkC : Отсюда и следует равномерная сходимость явной схемы со скоростью

n=0

O( + h2):

9. Исследование устойчивости методом гармоник. Познакомимся с еще одним распространенным приемом исследования разностных схем с постоянными коэффициентами, называемый методом гармоник. Хотя данный метод не является достаточно обоснованным, в частности, не учи-

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

Метод гармоник рассмотрим на примере схемы с весами [10, 12]
vt= v( v^ + (1 )v); v(0; t)=v(1; t)=0; v(x; 0)=u0(x):	(43)
Частные решения уравнения (43) ищем в виде
vmn (')=qneimh';	(44)

ãäå i мнимая единица, ' любое действительное число и q число, подлежащее определению. Подставляя (44) в уравнение (43) и сокращая на eimh', получим

q 1

=[q + (1

)]

eih'

2 + e ih'

;

h2

откуда найдем

h'

1 4 (1 ) sin2

q=

2

; = =h2:

h'

1 + 4 sin2

2

Начальные условия vm0 (') = eimh', соответствующие решениям вида (44) (их называют гармониками), ограничены. Если для некоторого ' множитель q станет по модулю больше единицы, то решение вида (44) будет неограничено возрастать при n ! 1. В этом случае разностное уравнение (43) называется неустойчивым, поскольку нарушается непрерывная зависимость его решения от на- чальных данных. Если же jqj 1 для всех действительных ', то все решения вида (44) ограничены при любом n и разностное уравнение устойчиво. jqj 1 для всех действительных '; åñëè

	1		h2		(45)
				:
	2		4

Отсюда видно, что все схемы с 12 абсолютно устойчивы. Схема повышенного порядка ап-

проксимации ( = ) также абсолютно устойчива. Явная схема условно устойчива при <0; 5h2 .

10. Трехслойные схемы для уравнения теплопроводности. Одной из первых схем, применявшихся для численного решения уравнения теплопроводности была явная трехслойная схема Ричардсона [11, 12, 13]

	vn+1 vn 1		= vn	èëè vo = v;					(46)
		2			t
ãäå	v^ v
vo =		; v^=vn+1	; v=vn; v=vn 1; v=vxx:
t	2
Эта схема имеет второй порядок аппроксимации по				è h (	= u ut	=O(	2	2	)): Однако она явля-
								+h
					o

ется абсолютно неустойчивой (т. е. неустойчивой при любом законе стремления h è ê íóëþ).

Уравнение можно переписать в виде

	vmn+1 vmn 1	=	vmn 1 2vmn + vmn +1	:	(47)
	2		h2

43

Если в правой части уравнения (47) заменить 2vmn суммой vmn+1 + vmn 1, то получим трехслойную

схему "ромб"(схему Дюфорта и Франкеля)

	vn+1 vn 1	=	vmn 1 vmn+1 vmn 1 + vmn +1	;	(48)
			h2
	2

которая остается явной относительно vmn+1 и является абсолютно устойчивой. Она может быть за-

писана в виде

vo +	2	v = v;	v =	vmn+1 2vmn + vmn 1	:	(49)
t	h2	tt	tt	2

Схема "ромб"получается из схемы Ричарсона добавлением к левой части (48) члена 2=h2v , îáåñ-

tt

печивающего устойчивость. Погрешность аппроксимации схемы есть

2

@2u @u 2 @2u

2

2

2 @2u

2

2

= u ut

utt =

+O(

+ h

)=

+O(

+ h

):

h2

@x2

@t

h2

@t2

h2

@t2

o

Отсюда видно, что схема "ромб"обладает условной аппроксимацией

=O( 2 + h2 + 2=h2)=O(h2)

ïðè =O(h2):

Если взять = h(1 + O(h)); ãäå =const, то, очевидно, что схема (49) аппроксимирует уравнение

âèäà @u=@t + 2@2u=@t2=@2u=@x2:
Обычно для уравнения (3) используются неявные трехслойные схемы с весами:
а) симметричные схемы	= ( v^ + (1 2 )v + v) + ';
vt				(50)
o
б) несимметричные схемы
	v	t	+ v = v + ':	(51)
			tt

Уравнения (50) и (51) содержат три слоя (tn 1; tn; tn+1). Значение v(x; 0) известно, значение v(x; ) надо задавать дополнительно; например, можно положить

v(x; )=v(x; 0) + (u000 (x) + f (x; 0)):

Тогда v(x; ) u(x; )=O( 2):

Òàê êàê u^ + (1 2 )u + u= u + 2utt= u + O( 2); то симметричная схема (50) при любом имеет второй порядок аппроксимации по и h: Выражение погрешности аппроксимации для схемы

(51) имеет вид

= u^+' (ut + utt)=Lu^+'

@u^

@2u^

@2u^

+O( 2+h2)=

(52)

0; 5

+

@t

@t2

@t2

@u^

@2u^

= Lu^+f^

+(' f^) ( 0; 5)

+O( 2+h2):

@t

@t2

Отсюда видно, что и для схемы (51)

2

2

^

=O(h2

+ )

ïðè = 0; 5; '=f:^

=O(h

+ )

ïðè =0; 5; '=f ;

6

Выписывая в (52) члены O(h2) и учитывая уравнение @u=@t= @2u=@x2+f; можно убедиться в том, что схема (51) имеет аппроксимацию O(h4+ 2) ïðè

			h	2		2 ^			^
2		^				@ f			@f
=0; 5 + h	=(12 );	'=f +	12		(	@x2		+	@t	):

44

Для определения v^ из (50) и (51) получаем трехточечные уравнения
Av^m 1 Bv^m + Cv^m+1= Fm	(53)

с правой частью Fm; зависящей от v; v; ' с обычными краевыми условиями при m=0; m=M: Ýòà

задача решается методом прогонки.

Достаточные условия устойчивости имеют вид

>1=4 для симметричной схемы (50);

0 для несимметричной схемы (51):

4.2Уравнение колебаний струны

1. Постановка задачи. Вычисление погрешности аппроксимации. Рассмотрим уравнение колебаний струны в виде [12]

	@u2	=	@2u		+ f (x; t);		0<x<1;			0<t T;	(1)
	@t2		@x2
удовлетворяющее начальным условиям
u(x; 0) = u0(x);						@u(x; 0)		= u0	;	0 x 1	(2)
						@t

(начальное отклонение u0(x) и начальная скорость u0(x)): Концы струны движутся по заданным

законам

u(0; t)= 1(t); u(1; t)= 2(t); 0 t T:

(3)

Введем в области G=(0 x 1; 0 t T ) прямоугольную сетку

!h =!h ! =f(mh; n ); m=0; M ; n=0; N g

с шагами h=1=M è =T =N: Так как уравнение (1) содержит вторую производную по t, число слоев

не может быть меньше трех. Воспользуемся следующие обозначения

v=vn; v^=vn+1; v=vn 1; v

=

v^ v

; v =

v v

; v=v

xx

;

t

t

v =

vt vt

=

v^ 2v + v

; vo =

vt + vt

=

v^ v

:

tt

2

t

2

2

Рассмотрим семейство схем с весами

vtt= ( v^ + (1 2 )v + v) + ';

'=f (x; tn);

(4)

v0= 1(t); vM = 2(t);

v(x; 0)=u0(x);

vt(x; 0)=~u0(x);

ãäå u~0(x) определим ниже.

@2u @2u

Дифференциальный оператор Lu=

был заменен разностным оператором Lh , исполь-

@t2

@x2

зующий значения сеточной функции в три момента времени t ; t; t + : Минимальным является

	2x	x	x3		x	x	x
пятиточечный шаблон		x		ïðè = 0, максимальным - девятиточечный шаблон	2x x x3
ïðè = 0.	4	x	5		4x	x	x5
6

45

Вычислим погрешность аппроксимации схемы (4) при '=f (x; tn): Пусть v - решение задачи (4), u = u(x; t) - решение задачи (1). Вводим погрешность решения z=v u: Подставляя v=z + u â (4),

получим
ztt= ( z^ + (1 2 )z + z) + ;		(5)
z0=zn=0;	z(x; 0)=0 zt(x; 0)= (x);
ãäå = ( u^ + (1 2 )u + u) + ' utt	погрешность аппроксимации для схемы (4) на решении

u=u(x; t), =~u0(x) ut(x; 0) погрешность аппроксимации для второго начального условия vt=~u0(x): Краевые условия и первое начальное условие u(x; 0)=u0(x) на сетке !h удовлетворяются точно.

Выберем u~0(x) так, чтобы = 0( 2). Из формулы
ut(x; 0)=	@u(x; 0)	+	@2u(x; 0)			+O( 2) = u0	(x) + 0; 5	u000(x) + f (x; 0) + O( 2)
	@t		2		@t2
видно, что u~(x) ut(x; 0)=O( 2); если положить
					u~(x)=u0(x) + 0; 5 (u000(x) + f (x; 0)):				(6)
Определим порядок аппроксимации .
Пользуясь формулами
	u^=u + ut; u=u ut; u^ + (1 2 )u + u=u + 2utt;
перепишем в виде
	= ( u^ + (1 2 )u + u) + ' utt = u+ 2 utt+' utt:								(7)
Подставляя сюда выражения

@u 2 @2u 3 @3u

+ O( 4);

u^ = u +

+

+

2 @t2

3 @t3

@t

@u 2 @2u 3 @3u

+ O( 4);

u = u

+

@t

2 @t2

3 @t3

@2u

h2

@2u

u =

+ 0( 2 ); u=Lu +

L2u + O(h4); Lu=

; L2u =

tt

@t2

12

@x2

@4u @x4 ;

получим

@2u

@2u

h2

= (Lu

+ f ) + (' f ) + 2L

+

L2u + O( 2 + h4):

@t2

@t2

12

@2u

h2

@2u

Отсюда видно, что

= 2L

+

L2u + O( 2 + h4) ïðè ' = f , òàê êàê

= Lu + f .

@t2

12

@t2

@2u

Учитывая, что L

= L2u + Lf , получим

@t2

h2

= (' f )+( 2 +

)L2u+ 2Lf + O( 2+h4):

(8)

12

Отсюда видно, что:

h2

1)

=O( 2 + h4)

ïðè

=

+ ; ' = f 2 f;

12 2

1)

=O( 2 + h2)

ïðè ' = f; 6= :

Здесь - постоянная, не зависящая от и h, которая должна выбираться так, чтобы схема была устойчивой (достаточно потребовать 0).

46

2. Устойчивость по начальным данным. Метод разделения переменных. Решение задачи

vtt= ( v^ + (1 2 )v + v)= v( ) ;

(9)

v0=vM =0; v(x; 0)=u0; ut(x; 0)=~u0(x)

будем искать методом разделения переменных.

Будем искать решение уравнения (9) в виде произведения функций, одна из которых зависит только от t а вторая только от x, полагая v(x; t)=X(x)T (t): Подставим это выражение в (9) и

учтем, что

v=vxx=T X; vtt = XTtt:

Тогда получим

X

T

=

tt

= :

(10)

X

T ( )

Отсюда и из краевых условий v0=vM =0 получаем для X(x) задачу на собственные значения

X + X=0;

x 2 !h;

X(0)=X(1)=0; X(x) 60:

Она имеет решения (см. с. 31)

4

2 kh

(k)

p

k =

h2

sin

2

; X

(x)= 2 sin kx; k = 1; M 1:

Собственные функции fX (k)g образуют ортонормированную систему.

Èç (10) äëÿ Tk (t) получим разностное уравнение второго порядка

(T )

+

T

( )

=0;

k tt

k

k

èëè

^

(1 +

2

2

k ]Tk + (1 +

2

k )Tk 2[1 + ( 0; 5)

k )Tk =0;

которое перепишем в виде

^

0; 5 2 k

(11)

Tk 2(1 k )Tk + Tk =0;

k =

1 + 2 k

:

Решение этого уравнения ищем в виде Tk =Tk (tn)=qn. Äëÿ q из (11) найдем квадратное уравнение

k p

q2 2(1 )q+1=0 (индекс временно опускаем). Его корни равны q1;2=1 2 2 . Åñëè 0< <2;

p

то корни q1;2=1 i (2 ) комплексные и jq1;2j=1. Введем новую переменную 'k ; полагая

p

cos 'k =1 k ; sin 'k = k (2 k ):

Тогда получим q1= cos 'k +i sin 'k =ei'k , q2= cos 'k i sin 'k =e i'k , (Tk )1=q1n=ein'k , (Tk )2=e in'k . Îá-

щее решение квадратного уравнения (11) имеет вид

Tk (tn)=Ck(Tk )1 + Dk (Tk )2=Ak cos n'k + Bk sin n'k ;

ãäå Ak è Bk произвольные постоянные.

Решение задачи (9) ищем в виде суммы частных решений

M 1
X
vn= (Ak cos n'k + Bk sin n'k )X(k)(x):	(12)

k=1

47

Пусть u0k ; u~0k коэффициенты разложений u0(x) è u~0(x)

M 1	M 1
X	X
u0(x)= u0k X(k)(x); u~0(x)=	u~0k X(k)(x):	(13)
k=1	k=1

Потребуем, чтобы сумма (12) удовлетворяла начальным условиям v0=u0, vt0 = (v1 v0)= =~u0(x): Тогда для определения Ak è Bk получим:

Ak =u0k ;

Ak

cos 'k 1

+ Bk

sin 'k

=~u0k :

Отсюда находим

1 cos 'k

Ak

=u0k ;

Bk=

u0k +

u~0k :

sin 'k

sin 'k

Подставив Ak è Bk в (12), после преобразований имеем

M 1

cos(n 0; 5)'k

sin n'k

vn=

u

0k

+

u~

X(k)(x):

k=1

cos 0; 5'k

sin 'k

0k

X

Потребуем, чтобы k удовлетворяло соотношению

k	2	; k=1; M 1:
	1 + "

Âэтом случае параметр удовлетворяет условию

1 + " h2

4 4 2 :

Тогда

cos 0; 5'k =p1 k =2

r

1 + "

:

Далее,

"

k =

k cos 0; 5'k

k r

1 + "

;

sin ' 2 sin 0; 5'

2 sin 0; 5'

"

è, òàê êàê

p

2 sin 0; 5'

2(1

cos '

)

k

2

=

p

k

=

k

= p k ;

то справедлива оценка

sin '

"

k

p k

:

1

+ "

Из (15) следует неравенство

(14)

(15)

(16)

(17)

(18)

vn	M 1	cos(n 0; 5)'k u			X(k)	+	M 1	sin n'k u~			X(k) ;
	X						X
k	k k	cos 0; 5'k		0k	k	k	k=1	sin 'k		0k	k
	k=1

подставляя в которое оценки (17) и (18), имеем

+ v

M 1

k

vn

k r

1 + "

k

u0

k

(~u0k )2

:

"

k

u

k=1

u X

t

48

Выражение

v

uM 1 2

u X (~u0k )

t

k=1 k

есть "негативная"норма (норма в HA 1 ):

ku~0kA 1 =(A 1u~0; u~0)1=2;

ãäå Av= v= vxx

Действительно,

и, следовательно,

в пространстве функций v, заданных на сетке !h и равных нулю при x=0; x=1:

M 1		M 1 u~0k
X		X
A 1u~0= u~0k A 1X(k)=				k	X(k);
k=1
		k=1
M 1	(~u0k )2
X
(A 1u~0; u~0)=				:
k=1		k

Итак, если выполнено условие (16), то для схемы (9) справедлива оценка

r

kvnk 1	+ "	ku0k + ku~0kA 1 ):	(19)
	"

4.3Задача Дирихле для уравнения Пуассона

1.Исходная задача. Рассмотрим уравнение Пуассона

	@2u			@2u		(1)
u =			+		= f (x; y); (x; y) 2 G:
	@x2			@y2
Будем искать его решение, непрерывное в прямоугольнике
G = G [ = f(x; y) : 0 x l1; 0 y l2g
и принимающее на границе заданные значения
				uj = (x; y):		(2)

	Покроем данную двумерную область
		G сеткой узлов (xi; yj ), образованных пересечениями пря-
	ìûõ x = xi è y = yj , ãäå		1
	xi = ih;	h =				;	i = 0; 1; :::; M;
				M
	yj = il;	l =	1		;		j = 0; 1; :::; N;
			N

(h è l - шаги сетки ïî îñÿì Ox è Oy соответственно). Будем называть узлы (xi; yj ): внутренними, когда i 2 f1; 2; :::; M 1g, j 2 f1; 2; :::; N 1g è граничными, когда i = 0 èëè i = M , à j 2 f1; 2; :::; N

1g, и когда j = 0 èëè j = N , à i 2 f1; 2; :::; M 1g. Множество внутренних узлов сетки обозначим !hl, граничных узлов hl, всех узлов (внутренних и граничных) - !hl. Заметим, что узлы (x0; y0),

(xM ; y0), (x0; yN ), (xM ; yN ), соответствующие вершинам прямоугольника						, при использовании
					G
шаблона типа "крест"и не относятся ни к внутренним, ни к граничным узлам [3].
Начнем с построения разностного аналога оператора Лапласа
	@2u		@2u
u = L1u + L2u; L1 =		; L2u =		:
	@x2		@y2

49

Для каждого внутреннего узла (xi; yj ) каждый из операторов L1u =

@2u

èëè L2u =

@2u

àï-

@x2

@y2

проксимируем трехточечным оператором 1 èëè 2: [9, 10]

1

(3)

L1u 1uij = uxx;ij

=

(ui+1;j

2uij + ui 1;j );

h2

1

(4)

L2u 2uij = uyy;ij

=

(ui;j+1 2uij + ui;j 1);

l2

где принято обозначение uij = u(xi; yj ).

Используя (3) и (4), заменим оператор Лапласа (2) разностным оператором

v = 1v + 2v = vxx

+ vyy ;

(5)

который определен на пятиточечном шаблоне "крест"2

3 ; состоящем из узлов (xi

h; yj ),

(xi; yj ), (xi; yj l).

4

5

Задаче Дирихле сопоставим следующую разностную схему

								(6)
vij = vxx;ij + vyy;ij	= fij ; i = 1; M 1; j = 1; N 1
								(7)
vi;0 = (xi; 0); vi;N = (xi; 1); i = 1; M 1;
								(8)
v0;j = (0; yj );	vM;j = (1; yj ); j = 1; N 1:

Не вызывает затруднений формальная аппроксимация разностной схемы со вторым порядком точности на той же сетке с тем же шаблоном "крест"более общего уравнения

	@2u		+ C@2u@y2		@u		@u
A				+ D		+ E		+ gu = F:
	@x2
					@x		@y

Аппроксимируем первые производные на шаблоне "крест"по симметричным формулам

@u				=	u(xi+1; yj ) u(xi 1; yj )		+ 0(h2 );
@x jxi ;yj						2h
	@u					u(xi; yj+1) u(xi; yj 1)		2
	@y		jxi ;yj	=		2h	+ 0(l		):

Тогда фиксируя в уравнении точку ((x; y) = (xi; yj ), после отбрасывания остаточных членов,

получаем соответствующее сеточное уравнение

	Aij	ui+1;j 2uij + ui 1;j	+Cij	ui;j+1 2uij + ui;j 1	+Dij	ui+1;j ui 1;j	+Eij	ui;j+1 ui;j 1	+gij vij = Fij :
		h2		l2		2h		2h

Легко видеть, что данное уравнение не имеет принципиальных отличий от уравнения (6) в том плане, что также представляет собой пятиточечное разностное уравнение и при всевозможных зна- чениях i 2 f1; 2; :::; M 1g, j 2 f1; 2; :::; N 1g - это система линейных алгебраических уравнений с

числом уравнений, совпадающих с числом неизвестных, если дополнить ее сеточными граничными условиями (7), (8), и с аналогичной структурой матрицы коэффициентов.

2.Аппроксимация. Введем погрешность решения zij = vij u(xi; yj ). Подставляя vij = zij + uij

âуравнения (1) (2), получим, что погрешность удовлетворяет уравнению

zxx;ij + zyy;ij = ij ; (xi; yj ) 2 !hl;	(9)
zij = 0; (xi; yj ) 2 hl;	(10)

50