Марчук Г.И., Лебедев В.И. Численные методы в теории переноса нейтронов

Подождите немного. Документ загружается.

Доказательство сформулированных выше утверждений везде ис-

пользует тот факт, что в силу инвариантности пространства U (UA)

относительно матрицы НА (АН) все векторы г|з* = и

— и* и p

будут

принадлежать пространству (/, а векторы

и Арь — пространству

UА.

Это обеспечивается специальным выбором начального приближе-

ния и

, для которого 1° = Аи° — / = (/ — АН) (Лср — /) £ UA-

В заключение рассмотрим частный случай обобщенного метода

сопряженных градиентов. Предполагая* что А и Н — симметричные

и положительно определенные матрицы, и выбирая D = Л"

, получаем

=(Щ*-*,Рк-дл/1Рн-г№

Э*

=(6*-*. Pfc-i)/|Pik|Я, 6=1,2,...

(10.11.16)

§ 10.12. МЕТОД РАСЩЕПЛЕНИЯ ОПЕРАТОРА

Начиная с работ Д. Писсмана, Дж. Дугласа, X. Рэкфорда [271,

320] большое внимание было привлечено к итерационным методам

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

физики, основанным на расщеплении сложных операторов на прос-

тейшие (см. также [272,341]). В Советском

Союзе это

направление пред-

ставлено исследованиями Н. Н. Яненко [247—250], Е. Г. Дьяконова

[76—78],

А. А. Самарского [202, 203], В. К. Саульева [209]

др.

Эти ра-

боты привели к созданию ряда эффективных алгоритмов решения ста-

ционарных и нестационарных задач.

Пусть в уравнении (10.2.16) оператор А положительно определен и

представим в виде

Л= 2 4*. (10.12.1)

а=1

Пусть A

, где {/

} — набор индексов (t

*.., /

), 1 < i

< я, являют-

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

вам к оператору A

. Образуем оператор

С = С (t,

..., xt

) = (/ + т/Д/J (1 + т,Д0 ... (/ + т/Д),

(10.12.2)

где T

(a = 1, 2, ..., s) — некоторые числовые параметры. Пусть мно-

жество допустимых операторов имеет вид ЭД = {С"

(т/

1Э

..., т/

В методах расщепления сложился следующий конструктивный спо-

соб построения итерационных операторов:

П (/ + ^1

)(а*

-м*)=-т(Ла*-/), (10.12.3)

где т — числовой параметр.

Уравнение типа (10.12.3) при Л

|у

= A

(s = п) могло появиться

как результат исключения вспомогательных функций в различных

предложенных методах расщепления. Введенный класс операторов

(10.12.2) родствен классу операторов, предложенному в работе

280

Е.

Г. Дьяконова [76], с той разницей, что в этой работе дополнительно

предполагается, что операторы A

являются разностными одномер-

ными операторами. На примере кинетических задач увидим (см. § 11.

22,

11.23), что требование одномерности операторов является излиш-

ним.

Запишем схему реализации, соответствующую каждому шагу итера-

ционного процесса (10.12.3), в которой предполагается, что можно до-

статочно просто найти обратные операторы (/ + т^Л^,)""

(/ + ^Л

)^+^/»=-т(Л^-/);

(10.12.4)

{I + *i

==l

k+{s

X)/s

Итерационная схема содержит s + 1 произвольных параметров т

г/>

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

алгоритма.

Проведем исследования схемы расщепления при s = п, А

= А

+ г«А

) "*

* +

= f (10.12.5)

в том случае, когда операторы А

заданы в гильбертовом пространстве

й коммутируют; тогда они

имеют

общую систему собственных функций.

Рассмотрим спектральные задачи

аФ

Х«*>ф;

Л;ф* =

А,

(а)

ф*,

а=1,2,..., п. (10.12.6)

Здесь Аа — оператор, сопряженный к Л

. Предположим далее, что

собственные функции задач (10.12.6) образуют нормированный биор-

тогональный базис {ф*}и{ф/}. Решение уравнения (10.12.5) будем

искать в виде

а*= § **Ч

«- (ЮЛ2.7)

Предполагается,

что

функция / также разложима

ряд

по

собствен-

ным функциям задачи:

/=£

/|ф|-;

/| =

(/.Ф/*).

(Ю.12.8)

/ = 1

Ряды (10.12.7), (10.12.8) подставим в уравнение (10.12.5) и резуль-

тат скалярно умножим на ф*. Тогда для коэффициентов Фурье ш

Приходим к уравнению

П (1+*сД!

а)

)

' +M?-/i =

(Ю-12-9)

a-I

281

где

= ^ М

а=1

Уравнение (10.12.9) разрешим относительно неизвестной

ut^

•

Тогда

П {1+г

Ь\

а)

)-гк

=J*=I-

"*+— ~

(Ю.12.10)

П (1+т

»,|

+ т

Хр>)

а=1

Если спектр

{%1

а)

}

положителен,

то

нетрудно сформулировать

критерий сходимости итерационного метода (10.12.10).

этой целью

рассмотрим знаменатель прогрессии

Я1=(П

(l +

TaOT-TX,)/

(1+т

а)

). (10.12.11)

\a=*l

/ / a=I

Нетрудно видеть,

что

сходимость будет иметь место

при

выполнении

неравенства

(l +

)>TV2. (10.12.12)

Очевидно, достаточным условием выполнения этого неравенства будет

следующее:

2т

>т>0,

а=1, ..., п.

(10.12.13)

Условие (10.12.13) позволяет осуществить

широкой области

по-

иск оптимальных значений параметров

, т. Эта

задача решается вся-

кий

раз

самостоятельно исходя

из

конкретной структуры операторов

задач.

Перейдем теперь

рассмотрению важного частного случая, когда

= А

+ А

(10.12.14)

Этот случай

литературе рассмотрен

большой полнотой

общностью.

Сначала исследуем частную схему, когда

+ т

; A

= A

; а = 1,2.

Будем иметь

(/ + tHi) (/ +

) (u

k+1

)ix

+ Au

= 1.

(10.12.15)

Полезно отметить, что это уравнение

помощью тождественных преоб-

разований переходит

следующее:

(/ + Mi) (/ +

)

и^

= (I -

i4i)

(/ -

) u

+ т/.

(10.12.16)

Возможны различные схемы реализации уравнения (10.12.16).

На-

пример,

a*+W*

(/-T

)a*;

и*+*'*

(1-т

AJu*+4*+xf;

(Ю 12 17)

*'*--=U

+ *'

;

(/ +

г,Л,)а*

=и*

/*.

J ' '

262

Анализ системы уравнений показываетг, чга первые Дба урайненйй

оказываются явно разрешенными относительно вспомогательных неиз-

вестных величин, а два последних уравнения являются уравнениями

с операторами, реализация которых предполагается достаточно прос-

той. Другая широко распространенная итерационная схема реализа-

ции, принятая в методе переменных направлений, следующая:

(/ + т

)^+1/2 = (/-т

)^ + т

/; 1

(/+т

)^+^(/-т

1/2

+ т

/. ( I ' • '

Представляя u

в виде (10.12.7), приходим к системе уравнений

вида (10.12.9), из которой получаем

и?

=<7|И? + Л. (10.12.19)

где

Л=(1-^

)(1-^>т

)/((1+^

>т

)(1+Я}»Ч,));

Л =

/((l

Х^>

т,)

)).

Таким образом, для ошибок г\ = щ

—

ui имеем

e?=^

e?-

(10.12.21)

Выражение (10.12.21) показывает, что метод (10.12.16) сходится при

sup|ft|< 1. (10.12.22)

Если XI

0, %\

и параметры x

>0, то условие (10.12.22) будет

выполнено.

Уточним этот результат. Пусть параметры xj (х

= т

) выбираются

независимыми от k. Выберем их в методе (10.12.16) оптимальным об-

разом, зная лишь только, что спектр каждого оператора Aj лежит в

круге

с центром в точке r

exp (i(p

)

радиусом р < г

.Для этого рас-

смотрим функцию комплексного переменного z

= (1

— xjz)

(1 +

%,г)-\

j = 1,2. (10.12.23)

Учитывая, что линия уровня с константой

> с >

для модуля

функции wj есть окружность радиусом 2 |

т/1"

с (1 —

)""

центром в точке (1 +

)

—

т/"

, заключаем, что \ij =

= min max |

| достигается при

т.

г е й

= exp(-i9o)(rg-p

)-i/2; lij^Pjiro +

VTf^Y

(10.12.24)

Тогда

<7*

< |л

. Если ф

= 0, то

(тМ)-и*;

^=-(М

-т

)(М

+ тЧ*)-\ (10.12.25)

где М = г

+ р; /п = го

—

р.

Перейдем к рассмотрению проблемы оптимизации, когда параметры

зависят от k. Пусть х

= а*"

, т

Р*

S где а

, р

—

меняющиеся

с номером k параметры. Пусть 0 <

rrij

< Лр < Mj

где /п^ и Mj

—

точные нижние и верхние границы спектров соответственно самосоп-

283

(10.12.20)

ряженных операторов A

. Тогда после N итераций для г

= и —

— и

имеем

|| < Q

где

<2N-

max \F

\; F

(K v, a, P)== f\ ^^T '

(10Л2

Параметры a

, p

выберем так, чтобы

была функцией, наименее

отклоняющейся от нуля в прямоугольнике D = {т

<1 M

/я

<! v < М

}. Решению этой задачи посвящены работы Е. И. Золо-

тарева [84], Е. Вакспресса

[342],

Ю. В. Воробьева [47] и др.

Обозначим

(М

1У

г%

, а,

р)

= max

(К v, а,

Р) |;

пп

9 9

£W (Afx, т

1э

, m

) = inf YN.

а, э

Для нахождения оптимального набора параметров a, р требуется

для заданного N найти среди функций вида (10.12.26) ту, для которой

в области D

I|?N

(Mi, т\, M

, a, p) = EN (M

, m

). (10.12.28)

Эта задача решается в несколько этапов. Сначала найдем дробно-ли-

нейное преобразование, переводящее область D в квадрат и не меняю-

щее вид функций (10.12.26). Его ищем в виде

Я.

= (

№

+ р)/(1 +

x), v = (qy - p)l (1 - sy), (Ю.12.29)

где параметры

р,

s выберем такими,

чтобы:

1)область {ri'<i ji, у ^ 1},

где rj — подлежащее определению число, отображалась преобразова-

нием (10.12.29) в область D; 2) преобразование (10.12.29) не меняло

тип каждого множителя в (10.12.26), т. е. чтобы

Я—ft

У—а

[л

—ft'

у—а'

К+ р

v+a

у + а' ' |х+Р' '

При этих условиях получаем, что если

0 = 2 (М

)

(М

~т

)/[(М, + т

) (М

то

+ е + (9(2+9))

/2)-1

;

(10.12.30)

s =

-m

)/[M

+ M

+(m

+ m

)r\]-

; (10.12.31)

-M

+ {M

+ M

)s)/2

= (M

+ (M

-M

)s)/2.

(ЮЛ 2.32)

Пусть

Pi = (?P/+p)/(l +

sP/);

a, = foa/-p)/(l~sa/)

1.....JV.

(10.12.33)

Тогда после замены переменных[(10.12.29), (10.12.33) в FN получим

a,P) = M|i,Y,a',p')

(Af

, т

, а,

р)

(1,

т],

а', р').

284

Можно показать, что inf T

(1, т),

т),

а', р') достигается при

а/ = Р/. Следовательно, экстремальная функция вида (10.12.26) в

новых переменных имеет вид

(JA,

a') f

(у, р'), где

М*,а) = П^р-. (ЮЛ 2.34)

Следовательно, EN

(М

1У

, т

) = Ем, где

= inf max | f

(t, а) | , (10.12.35)

а 0 = (ах, ..., 0N). Функцию, реализующую равенство (10.12.35),

обозначим~

(t). Решение задачи (10.12.35) было найдено Е. И. Золо-

таревым еще в 1877 г. [84] с помощью аппарата эллиптических функ-

ций,

а затем вновь найдено Р. Жорданом

[342];

параметры at совпада-

ют с множеством чисел

s.^dnf

K,k\, i = 1,2,..., N. (10.12.36)

Здесь dn(w,

k)—эллиптическая

функция, определяемая равенством

dn(M) = (l—

)

1/2

где

sin\p;

k = (l—

)

1/2

;

J(l— /j

sin

0)-

de =

Я/2

/C= f (1—A^sitfe)-

^^. При малых г) имеем

s, « 2с i (1 + с

- 'i>/(l +c

'i), (10.12.37)

где

с = 4/4; *

= p(tf-0+l)/W

При получении формул для параметров в

позаботимся о такой

перенумерации параметров s

при которой происходило бы интен-

сивное уменьшение ошибки преимущественно на первых итерациях.

Итак, пусть N =

2P,

а параметры, реализующие (10.12.35)

(обозначим их временно s

...

SN),

пронумерованы в возрастающем по-

рядке: s

при i<k. Их порядок употребления

итерационном методе

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

хлг

= (Д, /

,...,

/N),

< j

^ N, в которой

число j

означает, что на k-й итерации следует взять o

= s

Перестановки к/v определим рекуррентно: х

(1);

если известна

перестановка

K^-I—(/ь/*

•••»/*--0» то

щд

= (2*-i + l-j

29-

..., 2*-Ч-1-/ь, 2*-

/*,...)•

Например, v

= (6,11,3,14,7,10,2,15,5,12,4,13,8,9,1,16). Тогда

формулы для параметров а

, реализующих (10.12.35) и упорядочен-

ных в соответствии с перестановкой Х

Р, имеют вид, если t

= t,

^^/jVO + W

/|=(/

f+l

+ t|«

i/'i

i)/0 +

4i+i),

—

!• °-

(10.12.38)

285

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

/TLV /;

[ (10,12.39)

^+1=-

Йг. /=l,2,...

)

2';i =

)

l,...,/-l

о^

—

о^—оптимальный набор параметров задачи (10.12.35). Не-

трудно показать, что при этом для любого 0<>^р

(/) = f]-±l£.

<l (10.12.40)

и что

^=(1-Т1

)/(1+^

Упорядоченный набор параметров обозначим

Г?

{а

}^.

Пусть

Г?=К:2«(/-1)<Л<2-/};

Ф/(0= П |=£;

?(0

= М0-

а£е=Г?

Тогда

*J(0 =

(*>--*Г)/

(

^+°Г): Vs< ',-*<

*•

Обозначим

max

(0;

р**

= min

ф*.

(ft;

;

= max |сй (0|.

Справедливы следующие леммы

[135].

Лемма

10.12.1,

Если

>аГ"

>а?~

>]/"гь^ "^ Лр-в<<$~

<a?~

<]Ab-

> то m'<mj.

Лемма 10. 12.2. Если r\

<а}-

<1,

</</?, г)<1,

то

<V!

< <Vi < yTJ7<

а<,

а<

< 1.

Лемма 10.12.3. Для

любого

l<s<p и 1</<2"-

функция ф^(/)

имеет

на

отрезке

[г\,

1] 2

1 одинаковых

поло-

жительных максимумов

и 2

"™

одинаковых отрицательных

мини-

мумов; при

нечётном

^/^|Р/|,

яри

чётном:

.\^P

Лемма 10.12.4. Для любы* 0<s<p—

l^/*^-*-

^/-I^

772

!/»

а n

нечетных

/: ^2/~i^

2/+p

^/^2^-

2/+1>

2/+з> 1</<2^«-

-2.

Значение

достигается функцией

Ф^^О

при /=1, a

/raj/

функцией <p*

(/) яри

г].

286

Следствие

10.12.1.

mf-

mg-

mf-

mg-

<m§-

mj-

Лемма 10.12.5.

Для

каждого 0<s<p

— 1

min tn

}

достигается

при

j =

2P~

— 3,2"-

—

a max

т)—при

2P~

—

1,2P-

Пометим чертой сверху (соответственно снизу)

те

функции

Ф?(0,

для которых достигается minm^ (соответственно max

тр.

Тогда

после

итераций

/1 \

—г?—2

"~р—

2~о— 3

~'р—3 ~~p—k "~p—k ~~Q ~"0 1

М*,0) =

ф?

Ф2 Ф5 Фб

-ф5*-зФ2*-2-Ф

Р_

-2

-1-

(10.12.41)

При этом каждая

из

функций

ф/ (£)

устроена аналогичным образом.

Формула (10.12.41) показывает,

что

убывание ошибки

на

промежуточ-

ных этапах итерации происходит

определенном смысле оптимальным

образом

достаточно равномерным подавлением всех компонент ошиб-

ки.

Алгоритм упорядочения параметров можно обобщить

на

случай

= Пр

где p

—

любые простые числа.

работе

[136]

построены

бесконечные последовательности параметров {crjf

целых чисел

оо при п->- оо,

такие,

что для

каждого

N = N

функция

с параметрами

{crjf

наименее отклоняется

от

нуля

на

отрезке

[т|, 1].

Основой

для

получения таких последовательностей послужила теория

главных преобразований эллиптических функций

результаты работы

Е.

И.

Золотарева

[84].

Опишем коротко конструкцию такого алгорит-

ма.

Пусть

= (р

/?2>---> Ps)

—

некоторый набор простых нечетных

чисел;

—

любое целое число

из

, выбираемое

на п-м

шаге.

Пусть

= 2", где п > 0 —

задано,

a N

n +1

работе

[136]

построены классы бесконечных последовательностей

{<r*}f,

обладаю-

щие следующими свойствами.

Для

каждой последовательности имеем

7».

.(О =7*

(0 "П

Qj(t,N

),n

0,l,...

Функции

Qj (t, N

)

имеют

вид

287

Каждая функция Qj (t, N

) N

раз меняет знак на

[т],

1] и имеет [N

/2]+

одинаковых максимумов и N

—[NJ2] одинаковых минимумов на

[т|,

1]. Это свойство назовем почти чебышевским альтернансом.

Функции Q

^i (t, N

)

(t, N

), j = 1, ..., (p

— l)/2, имеют

почти чебышевский альтернанс.

Пусть б — любое из чисел 2*, k = 1, 2, ..., п, N

, п = 1,

2,...,

а Г

/ (

—О,

1, ....)—отрезок последовательности {ffj]/{.~i"

. Тогда для

любых б, / функции

8у

обладают почти чебышевским альтернансом.

Употребление в итерационном методе последовательностей

{(7/}~,

обладающих свойствами 1—4, позволяет вести итерационный процесс,

не ограничивая его априори каким-то числом итераций; при этом для

коммутативного случая итерационный процесс для всех k = N

бу-

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

он обеспечивает достаточно хорошее убывание ошибки. Наиболее

«плотные» и просто устроенные последовательности получаются при

Пусть п = О, т. е. Л/о = 1 и

УС

= (1). Если известны перестановки

« -

= Ц

Г\ /Г"

, ..., /V- *)•

то x

3" =

И1>

Й» •••• /а») определяет-

ся так:

Л # 1 о Qrt—1

//Н-г-з"-

—

(

—

^ , 6=1,2,..., 3'

(10.12.42)

Например, х

= (5, 2, 8, 6, 1, 7, 4, 3, 9).

Если п>0, то, зная перестановку и

<7- i

n = (Д, ..., /

<7 -1

«),

образуем перестановку х

л по правилу к

л = (2?-

+ 1 — Д,

2^—

+ Д,...) и полагаем

= 1, 2, ..., /г.

Проведенные расчеты диффузионных задач в прямоугольнике с не-

коммутативными операторами показали высокую эффективность бес-

конечных последовательностей параметров, упорядоченных предложен-

ными способами [15, 16].

§ 10.13. О ПРИНЦИПАХ, ПОЛОЖЕННЫХ В ОСНОВУ ОБЗОРА

МЕТОДОВ ПОСТРОЕНИЯ ИТЕРАЦИОННЫХ ПРОЦЕССОВ

К настоящему времени известно очень много различных итерацион-

ных методов. Классифицировать их все, отметить и исследовать ори-

гинальные находки в каждом методе — задача невыполнимая. В этом

деле лучше следовать по некой колее, пусть довольно узкой, но позво-

ляющей хоть отчасти систематизировать в определенном смысле извест-

ные методы. Очевидно, что при этом придется поступиться подчас

288

очень многим, в частности и приоритетными изобретениями тех авто-

ров,

чьи методы мы будем излагать.

Во-первых, будем условно считать эквивалентными итерационные

методы, отличающиеся друг от друга вычислительной схемой (поряд-

ком реализации операций), во-вторых, вне нашего внимания останут-

ся методы, специфичные только для какого-либо конкретного вида

уравнения. В-третьих, предположим, что как оператор А в уравнении

(10.2.16), так и операторы Т и Н в итерационных методах (10.4.3)

принадлежат коммутативной С*-алгебре скалярных операторов (см.

§ 10.2) 33 (X) с единицей.

При предположении, что 33 (X) есть коммутативная С*-алгебра с

единицей, можно интерпретировать и итерационные методы, и различ-

ные идеи, относящиеся к ускорению сходимости итераций, на уравне-

нии, в котором оператор А является просто скалярной функцией А (X),

где Я£Л. Мы уже воспользовались ранее этим приемом, проводя

доказательства сходимости итераций в § 10.12 [см. также формулы

(10.2.23) — (10.2.36), (10.3.16) — (10.3.19)]. Желающие ознакомиться

с обзорами итерационных методов, где в качестве модельной задачи ис-

пользуются матричные уравнения линейной алгебры, могут обратиться

к монографиям Д. К. Фаддеева и В. Н. Фаддеевой

[233],

Р. Варги

[341],

Е. Вакспресса [342] и др.

Рассмотрим теперь уравнение, заданное в 33 (X): Аи = /, где Л,

/ — заданные элементы, а и — искомый элемент. Согласно изложен-

ному в § 10.2, этому уравнению в пространстве

С(Л)

соответствует урав-

нение А

(Х)и(Х)

= /

(А,),

а итерационным формулам решения уравне-

ния в 33 (X) соответствуют формулы в пространстве С (Л).

Пусть для простоты Л — замкнутое множество, лежащее на дейст-

вительной оси, а X — элемент этого множества. Пусть функция

(X)

> 0 на Л. Норму элементов в этом пространстве определим как

|/|д,

= sup

/(Я) М

(Х)|;

/ е С (Л), (10.13.1)

А,еЛ

где | / | — абсолютная величина элемента /.

Наряду с пространством С (Л) рассмотрим гильбертово простран-

ство Нм (Л), построенное следующим образом. На множестве функций

из С (Л) введем скалярное произведение и норму согласно формулам:

если ф, я|>£С (Л), то

[ф,г|)]

Л1

= (Мф,1|))==|М(Я)ф(Я)115(Х)^а(Я), (10.13.2)

где da

(X)

— некоторая мера на Л. Тогда || ф

\\м

= [ф,

ф]]/

Замкнем множество С (Л) в норме (10.13.2); получим гильбертово

пространство функций, определенных на Л с конечной нормой (10.13.2),

его обозначим Нм (Л).

Таким образом, исследуем принципы построения итерационных ме-

тодов решения уравнения

А (к) и

(X)

= f(X) (10.13.3)

289