ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 11.04.2024
Просмотров: 300
Скачиваний: 6
162 |
Гл. 4. Методы расчета потенциала |
рифмической зависимости. Взяв N X = N Y = 128 и 2 мкс на опе рацию, получаем 3,8-ІО5 операций и 0,77 с. Фактическое время счета на CDC 6600 составило 0,75 с.
Приведенные выше оценки относятся к случаю периодических по координате х граничных условий ІВСХ — 3, который наиболее удобен для применения разложения Фурье. Для условий с задан ными значениями на границе или с нулевым градиентом (ІВСХ = = 1 или ІВСХ = 2) необходимое для преобразования Фурье время примерно в 2 раза больше и число операций принимает значение
5N X -NY [log2 N X + 0,3]. |
(36) |
Однако, как уже упоминалось, путем усовершенствования программы, выполняющей преобразования Фурье, это число, по всей вероятности, можно приблизить к соответствующему числу в периодическом случае (35).
ж. Полный объем памяти
На всех этапах расчета а — д вновь вычисленные величины можно записывать на место тех, из которых они получены. Сле довательно, в основной памяти хранится сама сетка из N X -NY узлов. Кроме того, в программе преобразования Фурье исполь зуется 3N X вспомогательных ячеек (для хранения массивов У, Z, INDEX, SI) и еще N X ячеек для хранения коэффициента каждой гармоники (массив АКХ), что составляет всего
N X -NY + 4NX машинных слов. |
(37) |
з. Время расчета
Втабл. 3 дано машинное время для сетки 128 X 128 при всех девяти возможных типах граничных условий. В колонке ошибок дано максимальное отклонение известного точного решения разностных уравнений от решения, найденного подпрограммой РОТ1. Точное решение лежит в пределах от —Ѵ2 до Ѵ2. Найден ная ошибка примерно в 500 раз превышает точность представ ления чисел на CDC 6600. Это примерно согласуется с ожидае мым результатом случайного сложения 3,8 НО5 ошибок округле
ния х), поскольку Y 3,8 НО5 = 618.
Табл. 4 дает затраты машинного времени в зависимости от размера сетки при фиксированном типе граничных условий. Время счета растет несколько медленнее, чем линейно, с ростом числа узлов сетки. Хорошо оправдывается следующая грубая оценка: время счета составляет 50 мкс на один узел сетки.
В При N — 128 требуется 3,8 -ІО5 арифметических операций.
§ 2. Прямые методы |
163 |
Таблица 3
Затраты машинного времени и величина ошибки при использовании подпрограммы РОТ1 на сетке 128x128 для различных типов граничных условий
|
CDC |
6 6 0 0 а |
IB M |
3 6 0 /6 7 6 |
І В С Х |
I B C Y |
|
|
О ш и бк а 8 |
|
В р ем я , с |
О ш и бк а 8 |
В р ем я , с |
1 или 2 |
1 |
0,91 |
2,0-10-1'' |
4,71 |
|
2 |
0,94 |
2,2-10-12 |
4,85 |
|
3 |
0,93 |
1,4-10-12 |
4,79 |
3 |
1 |
0,75 |
1,6-10-12 |
3,53 |
|
3,64 |
|||
|
2 |
0,77 |
3,7-10-12 |
|
|
3,58 |
|||
|
3 |
0,76 |
1,4-10-12 |
|
|
|
а А в т о к о д C O M P A SS .
0,8-10-4
2,2-10-4
|
К* |
to |
н* О1 |
1,2-10-4
1,9-10-4
3,2-10-4
6 |
В т о р о й в а р и а н т (у р о в ен ь Н ) к ом п и л я то р а с |
Ф ор тр ан а XV |
д л я м о дел и 6 7 - 1 . |
в |
Р а з л и ч и е в в е л и ч и н е ош и б о к о т р а ж а е т р а з |
н и ц у в д л и н е |
сл ова обы ч н ой точ н ости |
у р а с с м а т р и в а е м ы х м аш и н .
Таблица 4
Машинное время, необходимое для решения двумерного уравнения Пуассона методом FACR на сетках различного размера
при использовании подпрограммы РОТ1 |
для І В С Х = 3, I B C Y = 1 |
||||||||||
|
|
|
|
СХЭС 6 6 0 0 а |
|
|
|
IB M |
3 6 0 /6 7 |
6 |
|
С етка |
В р е м я , с |
О ш и бк а |
В р ем я , с |
|
О ш и бк а |
||||||
|
|
||||||||||
3 2 x 3 2 |
0 |
, 0 |
5 |
6 |
1 , 7 - 1 0 - 1 3 |
0 |
, 2 |
1 |
9 |
4 , 5 - 1 0 - 5 |
|
6 4 x 6 4 |
0 |
, 1 |
9 |
6 |
4 , 4 - 1 0 - 1 3 |
0 |
, 8 |
7 |
9 |
6 , 5 - 1 0 - 5 |
|
1 2 8 x 1 2 8 |
0 |
, 7 |
4 |
6 |
6 , 7 - Ю - і з |
3 |
, 5 |
2 |
7 |
|
1 , 1 - 1 0 - 4 |
2 5 6 X 2 5 6 |
2 |
, 9 |
5 |
4 |
— |
1 4 , 5 3 9 |
|
5 , 4 - 1 0 - 4 |
|||
а А в т о к о д C O M P A SS . |
|
|
|
|
|
|
|
||||
6 В т о р о й |
в а р и а н т |
(у р о в ен ь |
Н ) к о м п и л я то р а с |
Ф о р т р а н а |
д л я м одел и |
6 7 -1 . |
|||||
4. Метод DCR Бунемана
а. Алгоритм
Бунеман [13] разработал интересное обобщение процесса нечетно-четной редукции — метод двойной циклической редук ции (DCR). Он заметил, что уравнения для четных рядов тожё
11*
164 Гл. 4. Методы расчета потенциала
могут быть сведены к уравнениям для каждого четвертого ряда, а последние в свою очередь можно свести к уравнениям для каж дого восьмого ряда и т, д. Наконец, остается единственное урав нение для одного ряда. Это уравнение решается, а затем на остав шихся рядах потенциал определяется из промежуточных урав нений, соответствующих каждому уровню редукции. Этот алгоритм исключает необходимость преобразования Фурье и приводит к очень короткой машинной программе, но, как мы увидим в даль нейшем, эта программа не является самой быстродействующей.
Переобозначим для удобства исходные разностные уравнения
на некотором ряде сетки: |
|
Фяі ~ В<0>(Рі + Фя-і = — 2рГ т |
(38) |
где ВЫ — трехдиагональная матрица с элементами {—1, 4, —1},
апропорционально известному свободному члену.
Первый шаг нечетно-четной редукции приводит к новой систе ме уравнений для четного ряда узлов сетки,
фу-2— £ ШФ/ + Фя-2= —2p)1\ |
(39) |
Где ВЫ = (5(°>)2 - 21 и № = В») p f + р<°Л + |
р<•%. |
Новое уравнение (39) имеет ту же форму, что и уравнение (38), и, следовательно, процесс можно повторять по рекурсивной формуле
ф і-г^-В ^ф і+ ф і+ г1— — 2р(^, |
(40) |
|||
1 Б<Ч-і>= (5<«))*_ 2/| |
||||
log2 n' |
|
|
|
|
рГ |
^ |
рР + р^ |
+ рЙ -,.. |
(41) |
где I — уровень редукции, |
а п — исходное количество |
рядов |
||
сетки. |
|
редукции |
сводится к модификации |
|
Алгоритм однократной |
||||
свободного члена р® по формулам (41), до тех пор пока не оста нется единственное уравнение. На фиг. 5 показаны ряды, свя занные на различных уровнях редукции, в случае 16 узлов.
Для получения решения уравнения (40) определим обратный шаг алгоритма рекуррентной формулой
Фі = (B(0)_1 [2pf + фі+2г + Фі_2*] Для Z= log2«, . . 0 . |
(42) |
||
Отметим, что значения |
потенциала в правой части известны, |
||
так как они были вычислены на предыдущем уровне. |
|
||
Матрицы 5® быстро |
заполняются ненулевыми |
элементами |
|
(#(°) — трехдиагональна, |
ВЫ — пятидиагональна, |
ВЫ — девя- |
|
тидиагональна и т. д.), и, следовательно, необходимое для реше ния уравнения количество операций растет вместе с уровнем
§ 2. Прямые методы |
165 |
рекурсии. Процедура, однако, существенно упрощается, если представить В<г> в виде произведения 2l трехдиагональных
I
о
!
г
з
Ф |
Х |
Х |
Х |
Х |
Х |
Х |
Х |
Ф |
ф |
|
|
|
|
X |
|
X |
|
|
|
|
|
|
X |
|
|
|
Ф и г . 5. Столбцы, связанные на различных уровнях редукции в алгоритме
DCR.
С толбцы н еи зв ест н ы х вел и ч и н п о к а за н ы к р ести к ам и , |
а р я ды |
и зв естн ы х гр ан и ч н ы х зн а |
||||||
|
|
чен и й — к р у ж к а м и . |
|
|
|
|
||
матриц: |
- 4 |
1 |
0 |
0 ..,. 0 |
|
0 |
0 |
|
|
|
|||||||
В ш = |
1 |
- 4 |
1 |
0 .,.. 0 |
|
0 |
0 |
|
0 |
1 - 4 |
1 . .. 0 |
0 |
0 |
||||
|
о ’ |
о • |
о . |
о • |
о ■ |
1 |
•‘äf . і ■ |
|
Сначала рассмотрим матрицу |
|
|
|
|
|
|
||
В а>= |
2 і = (jgcw + у 2 1 ) (Вш - У |
2І), |
||||||
которая является произведением двух трехдиагональных матриц. Затем рассмотрим
В ш = (Я11’)2 - 21 = (Ваі + У 21) (Ва) — У Щ = |
|
= ((Вт)2— (2 — У 2)1) ((Bw>)2— (2 + У 2) I). |
(43) |
Следовательно, матрица
Б (2> = (Вт - у У ‘2,— У 2 1) {Bw - V 2 — Ѵ 2 I) х
X (В<т-{-У2 У~2 I) { В ^ - Ѵ 2 + У 2 І )
будет произведением четырех трехдиагональных. Очевидно, что процесс повторяется, и нужно лишь вычислять и запоминать
166 |
Гл. 4. Методы расчета потенциала |
диагональные элементы матриц-сомножителей. Элементы верхней и нижней диагоналей всегда равны —1. Общая формула такого факторизованного представления была дана Базби и др. [14].
Диагональный элемент матрицы Z?« быстро возрастает и имеет
порядок 42*. Для сетки 128 X 128 I = 7 и диагональный эле мент будет иметь порядок ІО78, что приведет к переполнению для некоторых вычислительных машин (например, для серии IBM 360 с максимальным представимым числом ~ 1075). Во всяком случае, для весьма небольших п переполнение произойдет на любой из существующих ЭВМ, если в ней использована встроенная длина действительного числа. (В GDC 6600 с максимальным вещественным числом ~10307 переполнение случится для сетки 512 X 512.) Можно хранить показатель экспоненты в виде отдель ного целого числа и иметь специальную подпрограмму для пере множения чисел в форме с плавающей запятой, однако лучше использовать исправленную форму рекуррентного соотношения, которое содержит умножение на матрицу, обратную к В, вместо бамой матрицы В. В этом случае может произойти скорее обра щение величины в машинный нуль, чем переполнение, что не создает трудностей, так как в большинстве ЭВМ обратившаяся в нуль величина автоматически заменяется точным нулем и вычис ления продолжаются.
Примененные Бунеманом [13] исправленные рекуррентные соотношения должны на каждом уровне редукции строить пере
менные |
|
|
<*> |
|
„(г-i) |
<г-1) |
со I |
|
|
|
|
|
|
|
|||||
|
Р?+1) = :(Я<0)-Ч2рУ |
-Pj-ft |
— Pj+ft |
+Vi-2h + |
|
||||
|
|
+ P;+2ft— P j- 3 f t + |
Pj+3ft] + |
1 |
(44) |
||||
|
|
rn® |
— |
„ ( И |
I |
„Й |
|
||
при 1 = |
1, . |
. ., IP; — P j-ft |
PjJ-/( |
+ |
P j - 2 f t + |
Pj+2/iJ |
|
||
|
|
log2 n, где |
h |
= |
2<!-1>, |
|
|
||
|
P? = ( В ^ Г 1[2pn |
+ |
[pf- 1+ Ря-i] |
при l = 0. |
(45) |
||||
Если |
для |
заряда и потенциала |
имеются раздельные |
сетки, |
|||||
то можно использовать несколько более простую форму этих рекуррентных соотношений.
Неизвестные значения потенциала вычисляются во время
обратного хода |
рекурсии по соотношению |
|
— Pi+ft*] |
(46) |
(В№) |
|
\-h |
||
при I45'=—log2 п,1 [2. р.^-.,^-j—ф7-_2Н —Н ФЛ-2Л.] + [ р ^ — Р |
|
|||
|
1, где h = 2(1_15 |
|
|
|
Здесь' опять на всех стадиях значения потенциала, входящие |
||||
в правую часть, |
известны. При 1 = 0 используются исходные |
|||
уравнения |
|
|
|
|
45' = ( В <0>) 1 [2 р )0) 4- 45'-і -|- 45-и] • |
(47) |
