Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
Таким образом, нужно обеспечить равенство
QtFGP = IFGQ. (!7)
Оно выполняется, если Q = Г. В'схеме с гарантированным сохранением импульса из операторов Q, I только один строится независи-
13*
388 ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ [Ч. II
мо, например оператор интерполяции I. Оператор раздачи заряда Q после выбора I определяется автоматически. Поясним смысл сказанного. Пусть оператор 1(х, у; i, j) реализован как кусочно-били-нейный. Это означает, что для вычисления
7(*. >0 = 27(*’ У! i, J) Іи ' (18)
U
нужно найти индексы I, т, при которых xG[lh,lh + l), у Є [mh, mh + h), и вычислить соответствующие коэффициенты интерполяции wa р (a, P = O, 1). Тогда
/(*> У) = "о, О A, m + wI1 О // + I, m + W0, . Л, m + l + Wl, 1 Л + 1, m + 1- (19)
где W0 0 = h~2(lh + h — x)(mh + h — у) и т.д. Элемент «матрицы»
/*(/, х, у), очевидно, вычисляется так: по точке (х, у) нужно найти индексы I, т и коэффициенты интерполяции, после чего определяются элементы матрицы Г:
Г(1 + а, т ¦+ х, у) = и>а р, а, 0 = 0, 1.
Все остальные элементы Г равны нулю.
Итак, операция раздачи заряда qk, находящегося в точке (Xk, Yk) состоит в следующем. По значениям Xk, Yk находится ячейка сетки (I, т), в которой эта точка находится, вычисляются соответствующие коэффициенты интерполяции Wa р, и доли заряда qkwa^ раздаются в узлы (7 + а, т + (3):
Р/ + а, m + p Р/ + а, m + P + Р‘
Если бы использовалась более точная интерполяция, мы столкнулись бы с таким «нефизичным» явлением, как раздача отрицательной доли заряда в какой-то узел. Видимо, поэтому на практике ограничиваются простейшим кусочно-билинейным оператором интерполяции, когда все Wa р > 0 и 2 wCt р ~ 1 •
Таким образом конструируются схемы интегрирования, обеспечивающие непрерывность плотности заряда р; т по положениям частиц (Xk, Yk) и сохранение импульса.
Быстрое дискретное преобразование Фурье. Имеется некоторая сеточная функция f(k) (к = 0, 1, ..., N— 1). Она может быть представлена разложением в сумму Фурье:
JV-I
f(k) Л(п) wkn, w=e2*ilN, к = 0, 1, ..., N- 1. (20)
п=о
§ 24] ПРИБЛИЖЕННОЕ ИНТЕГРИРОВАНИЕ УРАВНЕНИЯ ВЛАСОВА 389
Коэффициенты Фурье А(п) вычисляются по аналогичным формулам:
N-I
Л(«)=тя-2/(*) w~kn, я = 0, 1,...,TV-I. (21)
/t=0
Формулы (21), (20) представляют собой прямое и обратное преобразования. Выполнение каждого из них требует 0(N2) операций, если программировать непосредственно формулы (20), (21). Заметим, что нужно аккуратно вычислять степени w, чтобы избежать лишних операций. Ниже это будет сделано. Однако наиболее существенное сокращение объема вычислений связано с тем, что при вычислении сумм многие слагаемые и даже группы слагаемых вычисляются неоднократно, и этого повторения можно избежать. На этом основан алгоритм FFT (Fast Fourie Transform).
Рассмотрим реализацию (20), опуская, для простоты, множитель I/VN. Заметное сокращение операций происходит тогда, когда число N можно разложить на множители. Пусть N = MN1. Представим к и п в виде
к = XN1 + &j, п = V + ПуМ. (22)
Здесь х = 0, I, ..., M — I; = О, I, ..., TV1 — I; v = О, I, ..., M — 1;
M1 = 0, 1, ..., JV1 — 1. Показатель степени кп можно записать так:
кп = &(v + njM) = к V + к IilM = kv + л1М(хА^1 + ^1) =
= kv + Zt1XMTV1 + ^jrt1M.
Поскольку wMNі = Wjv= 1, слагаемое HiXMNl можно из кп убрать и формула (20) примет вид
JV-I JV-I
f(k) = ^ А(п) wkv WkЛм = ^ А(п) wkv WkSh W1 = wM. (20а)
п=0 п-0
Введем двухиндексную нумерацию, используя одно и то же имя для функций с разным числом целочисленных аргументов:
f(x,kl)zzf(xNl + kl), A(VjItl)=A(V^niM). (23) Сумму по п = 0, 1, ..., TV — 1 представим как двойную:
M-I Ni~l
/(и, &i) = ^ ^ А( v> rti) wkv (206)
v = 0 M1=O
где X = 0, 1, ..., M — I; = 0, 1, ..., TV1 — 1. Внутренняя сумма в (206) вычисляется при ПОСТОЯННЫХ V и t, поэтому множитель Wkv
390
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. II
можно вынести за скобки:
M-I Г nT1 1
/(*, *i) = 2 E лі) wI1"1
(20в)
v = 0 M1=O
Обозначим внутренние суммы, зависящие лишь от v и A1:
(24)
Тогда
/(X, ^) = 2A(V- *i) wkv’
(25)
х = 0, I, ..., М— I, =0, 1, ..., JV1 — I, Ic = XNl +Icl.
Итак, преобразование (20) распалось на два последовательных. Сначала исходные величины А преобразуются в A1 по формуле (24); это стоит CMN1N1 операций (постоянную С вычислим ниже: C = O(I)). Затем за CMMN1 операций выполняется преобразование (25). Итого мы имеем CN (М + N1) операций вместо CNMN1. При больших N и M«N1^VN это уже заметный выигрыш.
Если число N1 может быть разложено на множители, этот прием может быть использован снова. Будем понимать (24) как M преобразований типа (20). В самом деле, единственное свойство основания w в (20), которое было использовано, есть соотношение Wn = 1. Ho И WjrI = WwjvI = I. Пусть M = M1 — простое число, a TV1 =M2N2. Тогда при каждом v = 0, 1, ..., M1 — 1 нужно вычислить N1 коэффициентов ^j(v, ^1), а это можно сделать за CNl(M2 + N2) операций. Таким образом, мы имеем CTV1M1 (M2 + W2) операций. Вместе с CM1M1N1 операциями на вычисление (25) получаем CN(M{ + M2 + N2) операций. И т.д.
