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

Просуммировав это неравенство по переменной t = , 2 , : : :, n и учитывая, что kV ( )k = 0,

получают оценку (10).

В качестве примера рассмотрим несимметричную трехслойную схема для решения уравнения

теплопроводности

v

+ v

= v:^

t

tt

Пользуясь формулами

vt = vt

+

; vt = vt

; A = ;

2 vtt

2 vtt; v^ = v + vt +

2 vtt

o

o

o

приводим ее к каноническому виду

2

E

A

(E + A)vo +

(

+ +

)v

+ Av = 0;

t

2

tt

ò.å. B = E + A; R = E= + + A=2. Отсюда видно, что условия R > A=4; B E выполняется при

любом A > 0 è 0.

5. О регуляризации разностных cхем. Теорию устойчивости разностных схем можно использовать для формулировки общего принципа (принципа регуляризации) получения схем заданного качества, т. е. устойчивых, обладающих аппроксимацией и удовлетворяющих дополнительному требованию экономичности (минимума арифметических действий, достаточных для решения на ЭВМ получающихся разностных уравнений).

Требование экономичности применительно к нестационарным задачам математической физики обычно означает, что число арифметических действий, затрачиваемых для решения разностных уравнений при переходе со слоя на слой, пропорционально числу узлов сетки !h.

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

(E + R)vt + Av = ';		заданоv0;	(11)
Bv _ + 2Rv + Av = ';		заданы v0; v1	(12)
t	tt

было обнаружено, что за устойчивость отвечает оператор R (регуляризатор). Достаточные условия

устойчивости имеют простой вид:

R oA; o =			1			1			для двухслойных схем;	(13)
			2			kAk
R	1	A èëè R	1 + "				A	для трехслойных схем:		(14)
	4				4

Устойчивость или неустойчивость схемы (из семейства) зависит только от выбора оператора R. С точки зрения теории устойчивости произвол в выборе R ограничен лишь двумя требованиями:

1) cхема должна принадлежать иcходному cемейcтву, т. е.

B = E + R > 0 äëÿ (11) R = R > 0 äëÿ (12);

2) должны быть выполнены уcловия (13), (14).

Для получения устойчивой схемы заданного качества необходимо добиться, чтобы она имела аппроксимацию заданного порядка и была экономной, т. е. для решения уравнений (E + R)^v = F (äëÿ (11)) èëè (B + 2 R)^v = F (для (12)) требовалось бы минимальное, в некотором смысле, число

действий.

Заметим прежде всего, что если схема (11) или (12) с некоторым оператором R устойчива, то и

~

схема с оператором R R также устойчива. Обычно при построении разностных схем поступают

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

61

Оcновная идея регуляризации разностных cхем заключаетcя в cледующем: cхемы заданного ка- чеcтва надо иcкать в клаccе уcтойчивых cхем, отправляяcь от некоторой иcходной cхемы и заменяя ее, путем изменения оператора R другой cхемой нужного качеcтва, принадлежащей клаccу уcтой-

чивых cхем.

Запиcь cхемы в каноничеcкой форме удобна и для оценки порядка аппрокcимации. В (11) при R còîèò , à â (12) - 2. Поэтому, еcли при изменении R в cлучае двухcлойных разноcтных cхем оcтаетcя выполненным уcловие kRvtk = 0(1); то погрешноcть аппрокcимации меняетcя на величину, а в трехcлойных cхемах kRvttk = O(1) гарантирует получение cхемы, отличной от прежней на

0( 2):

Основным является вопрос о выборе регуляризатора R. Так как условия устойчивости имеют вид операторных неравенств, то в качестве R естественно выбирать операторы возможно более простой структуры энергетически эквивалентные оператору A. Пуcть например, A è A0 - энергетичеcки эквивалентные операторы c поcтоянными 1 è 2, òàê ÷òî:

1A0 A 2A0; 1 > 0; 2 > 0:

Полагая затем R = A0; получим уcтойчивые cхемы при 0 2 (èëè 0; 5 2) в cлучае (11), при

> 2=4 - в cлучае (12).

Простейшим видом R является оператор R = E (A0 = E). Условия устойчивости выполнены, если 0kAk äëÿ (11), kAk=4 äëÿ (12).

Укажем еще один способ выбора R. Пусть A0 = A1 + A2, A2 = A1 > 0. Выберем R так, чтобы

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

B = (E + A1)(E + A2) = E + ( A0 + 2A1A2);

òàê ÷òî

R = A0 + 2A1A2:

Òàê êàê (A1A2v; v) = (A2v; A2v) = kA2vk2 > 0; то эта cхема уcтойчива при 0 2:

Схемы с факторизованным оператором

B = (E + !R1)(E + !R2); R2 = R1

применяются в качестве итерационных схем для решения уравнений Av = '.

VI. Итерационные методы решения сеточных уравнений

6.1Введение

Применение разностных методов для решения дифференциальных уравнений приводит к системе линейных алгебраических уравнений специального вида к разностным уравнениям. Эта система обладает следующими специфическими свойствами: 1) она имеет высокий порядок, равный числу узлов сетки; 2) система плохо обусловлена (отношение максимального собственного числа матрицы системы к минимальному велико; так, для разностного оператора Лапласа это отношение обратно пропорционально квадрату шага сетки; 3) матрица системы является разреженной в каждой ее строке отлично от нуля несколько элементов, число которых не зависит от числа узлов; 4) ненулевые элементы матрицы расположены специальным образом матрица является ленточной.

Разностные эллиптически задачи в случае оператора общего вида или областей сложной фор-

мы решаются в основном при помощи итерационных методов.

Сеточные уравнения можно трактовать как операторные уравнения первого рода

Au = f

(1)

62

с операторами, заданными на пространствах H сеточных функций. В пространстве H вводятся скалярное произведение (; ) и энергетические нормы kukD = p(Du; u), D = D > 0, D : H ! H, ãäå D - некоторый линейный оператор в H.

Итерационные методы решения операторного уравнения Au = f можно трактовать как опера-

торно - разностное (разностное по фиктивному времени или по индекс-номеру итерации) уравнения с операторами в гильбертовом пространстве H. Если новая итерация vk+1 вычисляется через m предыдущих итераций vk, vk 1, ..., vk m+1, то итерационный метод (схема) называется m + 1- слойными (m-шаговым). Отсюда видна аналогия итерационных схем с разностными схемами для

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

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

Любой двухслойный (одношаговый) итерационный метод записывается в следующей канониче-

ской форме	vk+1 vk
B		+ Av	k	= f;	k = 0; 1; :::; v	0	2	H;	(2)
	k+1

ãäå B : H ! H - линейный оператор, имеющий обратный B 1, 1, 2, ::: итерационные параметры, k - номер итерации, vk итерационное приближение номера k.

Параметры f k g и оператор B произвольны, и их следует выбирать из условия минимума числа итераций n, при котором решение vn уравнения (1) приближает в HD точное решение u уравнения Au = f с относительной точностью " > 0:

kvn ukD "kv0 ukD :	(3)
На операторы накладываются следующие ограничения:
A = A > 0; B = B > 0; 1B A 2B; 1 > 0:	(4)

Операторные неравенства означают, что заданы постоянные 1, 2 энергетически эквивалентности операторов A è B или границы спектра оператора A в пространстве HB ( 1 минимальное, 2 максимальное собственные значения обобщенной задачи на собственные значения: Av = Bv).

В методах Зейделя и верхней релаксации оператор B несамосопряжен, и потому нельзя пользоваться методом (2). Оператор B можно сделать самосопряженным, если положить его равным

произведению сопряженных друг другу операторов:
B = (E + !R1)(E + !R2); R2 = R1;	(5)

ãäå ! > 0 параметр. В качестве R1 è R2 можно взять операторы, имеющие треугольные (нижнюю R1 и верхнюю R2) матрицы, так что R1 + R 2 = R : H ! H, R = R > 0. В частности, можно

положить

R1 + R2 = A;				R2 = R1:	(6)
Типичным является предположение
R E; R1R2			A;	> 0; > 0:	(7)
	4

p

Выбирая затем ! = 2= из условия min n0("), находим параметры 1, 2 и вычисляем параметры f k g. Определение vk+1 через vk è f сводится к последовательному решению двух систем уравнений

с нижней и верхней треугольными матрицами.

Построенный итерационный метод (2) с факторизованным оператором вида (6) называют попеременно - треугольным методом. Данный метод является универсальным, так как представление A в виде суммы R1 + R2 = A; R2 = R1 возможно всегда.

63

(x; y) 2 G;

= f (x; y);

При использовании общей теории итерационных методов формулируют общие правила решения разностных уравнений. Сначала изучаются общие свойства оператора A и устанавливается, например, его самосопряженность и положительность, A = A > 0, затем строится оператор B = B > 0 и вычисляются постоянные 1, 2 и, наконец, находят n = n0(") и параметры f k g.

6.2Модельная задача

Рассмотрим задачу Дирихле для уравнения Пуассона [13]

@2u @2u @x2 + @y2

u(x; y) = 0; (x; y) 2 ;

в единичном квадрате G(0 < x < 1; 0 < y < 1) с границей . Введем в G квадратную сетку с шагом h, т.е. множество точек

! = f(xi; yj ); xi = ih; yj = jh; i; j = 1; 2; :::; N; hN = 1g:

Пусть, как обычно, !h множество внутренних точек, h множество граничных точек сетки

Заменим исходную дифференциальную задачу разностной задачей

(8)

!h.

vxx;ij + vyy;ij = fij ; (xi; yj ) 2 !h;	(9)
vij = 0; (xi; yj ) 2 h;

которую будем рассматривать как модельную при изучении методов решения сеточных уравнений. Задачу (9) можно записать в виде системы

	vi 1;j 2vij + vi+1;j	+	vi;j 1 2vij + vi;j+1		=	f		;	(10)
	h2
				h2			ij
vi0 = viN = 0; v0j = vN j = 0;				i; j = 1; 2; :::; N 1:

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

Действительно, порядок системы (9) совпадает с числом точек сетки !h и равен (N 1)2. Äàæå ïðè h = 0; 1 имеем (N 1)2 = 81, т. е. матрица системы (9) является квадратной матрицей 81 порядка. На более типичной сетке, когда шаг h = 0; 01, порядок системы равен примерно 10000.

Сильная разреженность видна из того, что каждое уравнение системы (10) содержит не более пяти отличных от нуля коэффициентов. Тем самым отношение числа ненулевых элементов данной матрицы к общему числу ее элементов не превосходит 5=(N 1)2 = O(h2).

Отношение наименьшего собственного числа 1 = 8=h2 sin2 h=2 к наибольшему собственному числу 2 = 8=h2 cos2 h=2 равно

= 1 = tg2 h :2 2

Отношение является величиной второго порядка малости при h ! 0, а именно

= 2h2 + O(h4): 4

Следствием малости величины является плохая обусловленность системы (9). По этой причине

явные итерационные методы для системы (9) сходятся медленно.

64

6.3Методы простой итерации и Зейделя

Запишем разностное уравнение (9) в операторной форме Av = f , где оператор A определяется

следующим образом:

(Av)ij = vxx;ij vyy;ij ; (xi; yj ) 2 !h:

(11)

В дальнейшем будем рассматривать для этого уравнения одношаговые итерационные методы,

записанные в каноническом виде			vn+1 vn
		B			+ Avn = f:					(12)
Начнем с наиболее простых методов метод простой итерации и Зейделя.
Метод простой итерации для системы (10) записывается в виде
1
vijn+1 =		(vin 1;j + vin+1 + vi;jn		+1 + vi;jn		1 + h2fij );		(xi; yj ) 2 !h;		(13)
	4
		vijn = 0; (xi; yj ) 2 h:
Здесь vijn значение решения в точке (xi; yj ) 2 !h íà n-й итерации, = h2=4.
Скорость сходимости метода (13) определяется числом =							1		, = tg2	h	. Метод требует
							1	+		2

O(h 2) итераций для достижения заданной точности. Это очень медленная сходимость.

Реализация метода Зейделя для системы (10) приводит к следующему итерационному методу:

vin+11;j 2vijn+1 + vin+1	+	vi;jn+1+1 2vijn+1 + vi;jn	1	= fij ; (xi; yj ) 2 !h;	(14)
h2		h2

vijn = 0; (xi; yj ) 2 h:

Вычисления начинаются с точки i = 1; j = 1 и продолжаются либо по строкам, либо по столбцам

сетки !h.

Метод Зейделя сходится несколько быстрее, чем метод простой итерации, однако число итераций также является величиной порядка h 2.

6.4Метод верхней релаксации

Система Av = f представляется в виде

(L + U + D)v = f;

ãäå L; U - нижняя треугольная и верхняя треугольная матрицы с нулевыми диагоналями, D - äèà-

гональная матрица. Для ускорения сходимости метода Зейделя его модифицируют, вводя в итерационную схему итерационный параметр !, òàê ÷òî

(D + !L)

vn+1 vn

+ Avn = f:

!

Метод верхней релаксации определяется уравнениями

n+1

n

n+1

n

4

1 n+1

1 n

vi 1;j + vi+1

vi;j+1

+ vi;j 1

(15)

+

[

vij

+ (1

)vij

] = fij ; (xi; yj ) 2 !h;

h2

h2

h2

!

!

Способ нахождения значений vn+1

vijn = 0; (xi; yj ) 2 h:

на новой итерации тот же, что и в методе Зейделя.

ij

Метод верхней релаксации сходится при ! 2 (0; 2), а число итерации является величиной O(h 1).

65

6.5Попеременно-треугольный метод

Оператор A = A > 0 представляют в виде суммы A = R1 + R2; ãäå

vx

vy

vx

vy

R1 =

+

; R2 =

:

h

l

h

l

Для определения vk+1 íà k + 1-й итерации получаем уравнение

Bvk+1 = (E + !A1)(E + !A2)vk+1 = Fk ;

o

o

ïðè (x; y) 2 h):

Fk = Bvk + k+1( vk + '); (vk = ; vk = 0

Значения vk+1 находятся последовательно из уравнения

o

(E + !A1)W = Fk ; (E + !A2)vk+1 = W:

Отсюда получаем формулы

1Wi 1;j + 2Wi;j 1

+ Fk (i; j)

!

!

Wi;j =

; 1 =

; 2

=

;

1 + 1 + 2

h2

l2

o

o

o

1vk+1(i + 1; j) + 2vk+1(i; j + 1) + Wi;j

(16)

vk+1(i; j) =

:

1 + 1 + 2

Чтобы определить W1;2, выбираем узел i = 1; j = 1 в левом углу прямоугольника. Зная W1;1, последовательно находим W ïðè i = M1 1 è j = 1 (на первой строке). Далее полагаем j = 2 и находим последовательно W на второй строке при i = 1; M 1; т. е. проводим вычисления по столбцам сверху вниз. Счет vk+1 начинаем с правого верхнего угла и можно вести счет по строкам

справа налево или по столбцам сверху вниз.

Вычисления ведутся по рекуррентным формулам (16). Данный алгоритм называется алгоритмом бегущего счета.

Границы и оператора A; т. е. числа, для которых выполнено неравенство E A E,

равны

4

h

4

l

4

4

2

2

= min (A) =

sin

+

sin

; =

+

;

h2

h2

h2

l2

2l1

2l2

! = 2=p

p

; =

1

; 1 =

; =

; 2 =

:

1

p

4

+ 2

2(1 +

)

66

Список литературы

[1]Березин И.С., Жидков Н.П. Методы вычислений, том II/ - М.: Физматгиз, 1962. 640 с.

[2]Бахвалов Н.C.Жидков Н.П., Кобельков Г.М. Чиcленные методы. -М.: Лаборатория базовых знаний, 2000. 624 с..

[3]Вержбицкий .

[4]Волков Е.А. Чиcленные методы. -М.: Наука,1987. 248 с.

[5]Годунов C.К.,Рябенький В.C. Разноcтные cхемы. -М.:Наука, 1977.

[6]Калиткин Н.Н. Чиcленные методы. -М.:Наука, 1978. 512 с.

[7]Козлов В. Ф. О применении монотонных разностных схем при диагностических расчетах морских течений. -Изв. АН СССР, Физика атмосферы и океана, 1977, т.13, • 7, с.728-737.

[8]Охлопков Н.М. Методологические вопросы теории и практики разностных схем. -Иркутск: Издво Иркут. ун-та, 1989.

[9]Рябенький В.С. Введение в вычислительную математику. -М.: Физматоит, 1994.

[10]Cамарcкий А.А. Введение в чиcленные методы. -М.: Наука, 1987. 286 с.

[11]Cамарcкий А.А. Введение в теорию разноcтных cхем. -M: Наука, 1971.552 с.

[12]Cамарcкий А.А. Теория разноcтных cхем. -M: Наука,1983. 616 с.

[13]Cамарcкий А.А., Гулин А.В. Чиcленные методы.-М.: Наука,1989. 432 с.

[14]Турчак Л.И. Основы численных методов. -М.: Наука, 1987. 320с.

67

Учебное издание

Лилия Александровна Молчанова

РАЗНОСТНЫЕ МЕТОДЫ РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Учебное пособие для cтудентов математических cпециальноcтей

В авторской редакции

Технический редактор Л.М. Гурова Компьютерный набор и верстка автора

Подписано в печать ..08

Формат 60 84 1/16. Усл. печ. л. , . Уч.-изд. л. , . Тираж экз.

Издательство Дальневосточного университета 690950, Владивосток, ул. Октябрьская, 27.

Отпечатано в лаборатории

кафедры компьютерных технологий ИМКН ДВГУ 690950, Владивосток, ул. Октябрьская, 27, к. 132.