Файл: Вычислительные методы в физике плазмы..pdf

ВУЗ: Не указан

Категория: Не указан

Дисциплина: Не указана

Добавлен: 11.04.2024

Просмотров: 298

Скачиваний: 6

ВНИМАНИЕ! Если данный файл нарушает Ваши авторские права, то обязательно сообщите нам.

§ 2. Прямые методы

167

Можно показать, что соотношения (44) — (47) алгебраически иден­ тичны более простым соотношениям (41) и (42). Дополнительная

привлекательность новых

рекуррентных

соотношений состоит

в том,

что и прямой, и

обратный ход рекурсии записываются

в одной и той же форме (а именно В~ХР +

Q), и многие вычисли­

тельные операции можно использовать совместно.

На

каждом шаге прямой и обратной

редукции приходится

решать

21 трехдиагональных систем уравнений. Каждая такая

трехдиагональная система решается методом рекурсивной цикли­ ческой редукции, описанной в п. 3, в. В методе DCR, как и в ме­ тоде FAGR, необходима только одна сетка, которая используется как для заданного заряда, так и для неизвестного потенциала.

б. Подсчет числа операций и быстродействие

Рассмотрим разностную сеть из п рядов по т узлов в каждом. Общее число операций на 1-м уровне прямой и обратной редук­ ции определяется их числом в квадратных скобках выражений

(44)

и (46),

умноженным на п/21. Это составит І8тп/21 операций.

■Суммируя

по

I от

1

до бесконечности, получаем (с точностью

не

хуже

1%

при

п >

128) 18?шг операций.

 

На 1-м уровне нужно дважды для каждого ряда (во время

прямой

и

обратной редукции) решить 2(1~ѵ> трехдиагональных

вистем. Решение одной трехдиагональной системы уравнений требует примерно 6т операций. Решение всех систем на 1-м уровне требует 2 -6т •2<г-1>-п/21 = Qmn операций. Число необходимых уровней зависит от типа граничных условий в ^-направлении. При заданном потенциале на границе (ІВСХ — 1) нужно log2 п—1

уровней, а при периодических граничных условиях (ІВСХ =

3)

их нужно log2 п.

 

 

 

 

 

Полное число операций в методе DCR, таким образом, соста­

вит

 

 

 

 

 

6пт [loga п +

2],

ІВСХ =

1,

 

 

6ш?г [log2 п +

3],

ІВСХ =

3.

^

'

Соответствующий подсчет для метода FACR из и. 3, е дает

Ъпт [log2 п +

0,3],

ІВСХ =

1,

 

 

2,Ъпт [loga п +

2,4), ІВСХ =

3.

^49

Интересно также привести число операций в очевидном методе двумерного преобразования Фурье с последующим делением и обратным двумерным преобразованием Фурье [15]. Для такого метода нетрудно написать программу, располагая лишь библио­ течной подпрограммой преобразования Фурье. Назовем этот


168

Гл. 4. Методы расчета потенциала

 

метод DFA и, пренебрегая делением, запишем число операций:

 

 

10пт [log2 пт — 3,2],

ІВСХ = IBCY — 1,

 

 

Ъпт [log2 пт — 2,8],

ІВСХ = IBCY — 3.

'

Проведенное выше сравнение показывает, что в случае перио­ дических граничных условий метод FACR может в 2 раза пре­ взойти метод DCR по быстродействию, а при заданном значении потенциала на границе — лишь на 20%. В последнем варианте граничных условий метод DCR в 3 раза превосходит метод DFA, а для периодических условий на 50%.

Часть этих преимуществ теряется при учете вспомогательных операций. В этом отношении метод DCR гораздо предпочтитель­ нее, чем FACR, ввиду простоты программы. Программа Бунемана, применяющая метод DCR, состоит из единственной под­ программы и содержит 59 операторов языка Фортран (2464 байта), в блоке, решающем уравнение Пуассона. Этому следует сопо­ ставить 428 операторов Фортрана (11 120 байтов) и 14 подпро­ грамм в написанной Хокни по методу FAGR программе РОТ1. Общие выводы, основанные на подсчете числа операций, довольно хорошо подтверждаются значениями фактического машинного времени, сведенными в табл. 5.

Таблица 5

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

Время в секундах для IBM 360/67-1

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

компилятора '

(уровень H)

с Фортрана.

 

 

 

Сетка

DCR °

PACK 6

FACR 8

DFA 8

32X32

0,297

0,219

0,265

0,529

64x64

1,363

0,879

1,113

2,155

128x128

6,152

3,527

4,711

9,021

256x256

27,488

14,539

20,187

38,359

а

Подпрограмма Бунемана XYPOIS при ІВСХ = 1 , IBCY = l.

6

Подпрограмма POTI при ІВСХ =

3, IBCY =

I .

8

Подпрограмма POTI при ІВСХ =

I, IBCY =

I.

eПодпрограмма POT3 при ІВСХ = 3, IBCY — 3.

в. Выбор алгоритма

Для решения большой задачи на ЭВМ с большим объемом памя­ ти, по-видимому, предпочтителен алгоритм FACR. При этом могут быть удовлетворены увеличенные требования к объему


§ 2. Прямые методы

It,9

памяти и одновременно приобретаются преимущества в скорости от 20 до 90%, что при решении крупных задач может означать значительную экономию машинного времени. Алгоритм DCR становится оптимальным для малых ЭВМ с ограниченным объемом памяти (например, для IBM 1800 с памятью в 4096 слов или IBM ИЗО с памятью в 8000 байтов), так как, по-видимому, это единственный прямой метод, который может быть приспособлен к таким системам. В любом случае задачи, решаемые на малых ЭВМ, должны быть небольшими, так что более длительное время действия алгоритма DCR несущественно.

5. Оптимальный метод FACR (?)

Обобщим метод FACR, включив в него /-кратную нечетно­ четную редукцию перед выполнением разложения Фурье, и назо­ вем такой метод методом FACR (/). Во второй главе описан метод FACR (1) с однократной редукцией, а метод DCR можно рас­ сматривать как метод FACR (log2 п). Кроме того, метод FACR (0)г т. е. метод FACR без применения нечетно-четной редукции, был использован Веронисом [16] в связи с задачей о ветровых цирку­ ляциях океана.

Как мы видели, метод FACR (1) требует меньшего числа опе­ раций, чем DCR, и нетрудно убедиться, что он экономнее, чем FACR (0), поэтому встает вопрос о существовании оптимального уровня редукции.

Будет

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

/ = 1 или

I = 2 и что проблема переполнения не возникает при

использовании для редукции простых рекуррентных соотно­ шений (41) и (42). Как и раньше, рассмотрим сетку (п X т} по координатам (х , у) с разложением Фурье в ^-направлении.

Преобразование свободного члена по формуле (41) требует выполнения 1,5пті + пт (1 — 2~1) операций. В результате полу­ чим т!21 уравнений, которые решаются методом разложения Фурье с последующей циклической редукцией, для чего необхо­ димо соответственно 2пт [2,5 log2 п — 3,5п]/2г и ЬптІ21 операций. На стадии обратной редукции необходимо 2пт [1 — 2~1] опера­ ций для формирований свободных членов и 3пті операций для решения трехдиагональной системы.

Полное число операций на один узел сетки в методе состав­

ляет

 

3 + 4 ,5 /+ (51 ° у —4).

(51)

На фиг. 6 показана зависимость числа операций на один узел сетки от числа уровней редукции / для типичного случая п = 128. Неглубокий минимум достигается при двукратной


170 Гл. 4. Методы расчета потенциала

редукции, и оценки показывают, что введение следующего уровня редукции в подпрограмму РОТ1 может привести к повышению числа операций на 10%.

Очевидно, что оптимальный уровень редукции зависит от отно­ сительной эффективности используемых подпрограмм разложения Фурье и циклической редукции. Ключевой вопрос состоит в том, как быстрее решить систему Ві1)х = Ъ — разложением Фурье дли циклической редукцией. По мере увеличения I и заполнения

»5s

5S

«

I

У

£cd

и-

§

I

•5

а

<0

5S

to

•Ф и г. 6. Количество операций на один узел сетки 128 х 128 для алгорит­ мов FACR (I), DFA и DCR.

Дано также эквивалентное (по затратам времени) число итераций для метода SOR в пред­ положении, что число операций в одном узле на одну итерацию составляет наименьшее значение, равное 7. Во всех случаях граничные условия выбраны в форме, наиболее благоприятной для данного алгоритма. Заметим, что минимальное число операций дости­ гается при I = 2 .

матрицы B {h ненулевыми элементами разложение Фурье стано­ вится все более и более подходящим методом решения. Выбрав метод разложения Фурье на 1-м уровне, получаем алгоритм FACR (I). Во время разработки исходного алгоритма FACR (1) на IBM 7090 [6] Голуб и Хокни обсуждали желательность сле­ дующего шага нечетно-четной редукции, однако затем эта возмож­ ность была отвергнута, так как тщательный хронометраж имев­ шихся тогда программ разложения Фурье и циклической редук­ ции показал, что алгоритм FACR (2) будет более медленным.

Алгоритм типа FACR (I) можно применять для более общих уравнений, чем рассмотренные выше. Всесторонний математиче­ ский обзор таких обобщений был дан Базби и др. [14] при описа­

§ 2. Прямые методы

171

нии семейства алгоритмов GORF (методы циклической нечетно­ четной редукции и факторизации).

Установлено (см., например, [17]), что выгодно использовать девятиточечную аппроксимацию уравнения Пуассона (см. [18]). Последняя имеет более высокий порядок точности по сравнению с пятиточечной схемой, и ее применение может оказаться необхо­ димым, когда используется крупная сетка. Описанные выше прямые методы можно применять и для девятиточечных уравне­ ний. Объем вычислений и машинное время для них возрастают по сравнению со случаем пятиточечной аппроксимации.

6. Учет электродов

Может показаться, что описанные здесь прямые методы реше­ ния уравнения Пуассона имеют довольно ограниченную сферу

применения, в

частности при расчете электронной пушки, так

как

в рассмотренной задаче электроды не могут находиться

во

внутренней

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

ее границу. Один из возможных способов обхода этого ограни­ чения [19] состоит, например, в предварительном вычислении матрицы емкости С, связывающей потенциал и заряд в некото­ ром числе точек во внутренней области. Эти точки следует рас­ сматривать как электроды, на которых поддерживается заданный потенциал, и вычислить поверхностный заряд, индуцированный на этих электродах окружающим пространственным зарядом. После этого результат будет определен путем двукратного реше­ ния уравнения Пуассона следующим образом. Сначала уравнение Пуассона решается при нулевом заряде на электродах и запоми­ нается отклонение потенциала на электродах от требуемого значения. Умножив это отклонение на матрицу емкости и на —1, определяют искомый поверхностный заряд в каждой точке элек­ трода. Уравнение Пуассона решается снова с найденными значе­ ниями поверхностного заряда, и это решение будет правильным всюду, включая и точки на электроде.

Для определения матрицы емкости строят сначала обратную матрицу А = С-1. Каждый столбец матрицы А находят, помещая единичный заряд по очереди в каждую точку электрода (при нуле­ вом заряде в остальных точках электрода) и решая уравнение Пуассона. Значения потенциала в точках электрода образуют элементы столбца матрицы А. Матрица емкости вычисляется путем обращения матрицы А, и это ограничивает применимость метода 100—200 точками для CDC 6600. Матрица емкости сим­ метрична, так что для I точек электрода необходимо I (I + 1)/2 ячеек памяти. Так как матрица емкости зависит только от гео­ метрии задачи и не зависит от распределения зарядов, ее доста­ точно вычислить только один раз в начале расчета и затем хранить