Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
Определим счетную область Gh. Будем считать точку (к, т) счетной, если (хк, ут) Є G. В этих точках определена сеточная функция ик т. Тогда Gh = U <оЛ т. Здесь, как и в дальнейшем, без
специальных указаний предполагается, что индексы (к, т) пробегают значения, соответствующие счетным точкам и только им. Приближенное решение ищем в виде функции
Используем билинейность формы I и вынесем коэффициенты Uk т и суммирование за пределы операции I. Обозначая Zij — (/, <р; у),
А\' ”1 — 1(ч>к, it>> Ф;, у)> получаем систему линейных алгебраических уравнений, аппроксимирующую уравнение (1) (или (4)):
Vk,m(t У]) = 0
Tl + |Т1 при О Sg I, Tl Sg 1,
к, т
и(х, у) = E Uk т ук т (х, у).
(6)
Составим уравнение (5), заменив V v на V ук т:
(7)
492
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. II
Очевидно, число уравнений равно числу неизвестных. В такой форме можно считать <р* m любым базисом, а не только тем, который
был описан выше.
Мы получили систему уравнений (7) с четырехиндексной матрицей. Хранение ее в памяти ЭВМ, если число базисных функций не очень мало, — проблема, впрочем, для современных ЭВМ (серии ЕС, например) уже преодолимая. Ho вот вычисление матрицы остается проблемой: ведь каждый элемент А — интеграл по G (хуже того, интеграл от функции с достаточно сильной особенностью). Именно это заставляет сузить общую галеркинскую конструкцию, используя специальные конечные элементы.
Заметим, что ядро интегрального уравнения (1) зависит не от четырех аргументов (лс, у, х’, у’), а только от двух (х — х', у — у'), что дает основание ожидать соответствующего характера матрицы: А^’ Jn = АЪ °. j . Это весьма полезное свойство действительно уда-
»» J Л- I1 у IfL '
ется получить, но только за счет использования одинаковых (с точностью до сдвига аргументов) конечных элементов. Очевидно,
У) = %,0(*-**> у — Ут)'
С учетом этого имеем /(<(>*, m, Vij) = tPo, о)'
Теперь уравнение (7) можно записать в форме
2Ahk-i,m-jUKm = fi,r Vi, j. (8)
к, т
Здесь Ah означает, что матрица А вычислена для сетки h х h. Несложный анализ (хотя бы размерностей) показывает, что можно вычислить «универсальную» матрицу на сетке Ixl один раз и после
этого пользоваться формулой Ah = h~lAl.
В универсальной матрице A1 наиболее сложно вычисляются элементы Ai j с малыми значениями |г| + |/'| <3, так как именно в этих случаях сказывается сингулярность подынтегральной функции. Для вычисления таких элементов были разработаны специальные программы аккуратного интегрирования, с помощью которых вычислены элементы Ai j для малых i, j (их можно найти в соответствующих работах). Остальные элементы Ai . легко вычисляются по простым асимптотическим формулам. При этом используется очевидное свойство симметрии: Ai j = A_t j — Ai _j, позволяющее хранить в памяти только «четверть» матрицы (/, / = 0, 1, ..., К, где К — целое число, такое, что Kh превосходит линейный размер G).
Итак, составление системы уравнений (8) проблем не содержит, если известны значения Ai у. при малых | /1 + | j \. Следующий вопрос — решение этой системы. Число уравнений не так уж мало. В расчетах, которые обсуждаются ниже, число узлов было от 500 до
§30]
ПСЕВДОДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ
493
1000, причем матрица системы полностью заполнена. Применение стандартных программ решения линейных уравнений потребовало бы порядка IO9 операций, IO6 памяти, что уже представляет серьезную проблему для современных машин ЕС (типа 1040, 1045, 1050), не говоря уже о БЭСМ-6, для которой это просто непосильная задача (а именно эта ЭВМ применялась в расчетах).
Достаточно эффективным средством решения системы (8) является простейший итерационный метод, известный под названием «релаксационный». Он состоит в поочередном пересчете ик т по формуле
Uk, т ~ Л>,0 (/*, т — 2 ^i-k,j-m Ui,j\- ^
' І
Здесь в сумме по і, j пропускается слагаемое і = к, j = m. Используется также ускорение сходимости методом сверхрелаксации. По формуле (9) вычисляется «предварительное» значение йк т, окончательное значение находится по формуле ик т = ик т + ш(йк т — ик т). Параметр ш определялся экспериментально.
В табл. 22 приведена эволюция нормы невязки уравнения (8) в зависимости от номера итерации v при разных значениях ш. Видно, что оптимальное значение 1.4, при этом скорость сходимости достаточно высока. Начиная с тривиального приближения Ui ^ 0, за 15 итераций можно получить невязку, существенно меньшую /. Такое приближенное решение системы (8) можно трактовать как точное, соответствующее незначительно измененной правой части.
Таблица 22
(О V 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7
1 4.6е1 5.Oel 5.4е1 5.8е1 6.4е1 7.Oel 7.8е1 8.8е1
4 1.5е1 1.8е1 1.5е1 І.Зеї I .Oel 7.3еО 5.6е0 6.5е0
7 8.8еО 6.8е0 4.8еО З.ОеО 1.4е0 4.7е— 1 8.3е— 1 2.OeO
10 4. IeO 2.6е0 1.4е0 6. Ie — 1 1.4е — 1 5.3е — 2 1.9е — 1 7.6е— 1
13 1.9е0 9.9е — 1 4.2е— 1 1.2е— 1 1.2е —2 8.4е —3 5.4е — 2 3. Ie — I
16 8.5е — 1 3.7е — 1 1. 2е — 1 2.3е —2 1.2е —3 1.5е-3 1.5е —2 1.4е — 1
19 3.9е — 1 1.4е — 1 3.5е —2 4.5е —3 I. Ie — 4 4.8е —4 4.4е — 3 5.7е — 2
