ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 11.04.2024
Просмотров: 280
Скачиваний: 6
112 |
Гл. |
3. |
Модель «водяного мешка» |
|
Типичные |
значения |
в |
случае 48-разрядного |
кода х) следующие: |
а = 30, I = 6, р = 46. Целая константа р выбирается большой, |
||||
чтобы в памяти ЭВМ |
координата X всегда |
оставалась положи |
||
тельной, как бы далеко не сдвинулась точка. Это делается для того, чтобы арифметические операции не зависели от способа представления отрицательных чисел в любой конкретной ЭВМ. Имеется несколько преимуществ в использовании степеней 2. Например, становится легче использовать операции типа MASK в Фортране (если они имеются в применяемом варианте языка) или перейти на автокод ЭВМ, чтобы наиболее трудоемкие участки программы считать как можно быстрее. Восьмеричная запись легко интерпретируется. Показатель а характеризует точность, с которой можно определить координату. Имеется несколько свободных разрядов для согласования между тремя частями «а»,
«б», и «в» |
кода, но |
практически |
этот метод так же |
точен, как |
|||
и метод, использующий числа с плавающей запятой. |
|
||||||
Перед |
вычислениями |
все точки |
лежат в |
исходном периоде |
|||
(Р ^ X < |
Р + AL), |
а |
затем |
они |
изменяют |
свое |
положение |
в соответствии с |
разностными уравнениями |
|
|||||
Х г [(п + |
1) Д*]= |
|
1(п - |
1) At] + |
Vi |
[nAt], |
|
Ѵі [(тг + |
1) At] = |
Vi [(n - |
1) At] + |
Ri |
\nAt], |
(21) |
|
где V и R — целые константы, представляющие соответственно 2vAt и ArАР, а г — ускорение. Такое представление позволяет избежать всех ненужных умножений. По истечении некоторого времени большинство точек выйдет за пределы начального периода через одну из его границ х = 0, L и попадет в соседний период. Однако нет необходимости проверять этот уход точек (на что понадобилось бы машинное время), так как для расчета плотности заряда, электрического поля и ускорения частиц мы просто про ектируем координаты на основной период 0 ^ X < AL, затирая все значащие разряды, кроме а + I последних. Это можно сделать на Фортране с помощью оператора
IXFUND = IXTRUE. AND. MASK, |
(22) |
|
где MASK = 2a+l — 1, |
IXTRUE — истинная |
координата |
и IXFUND — координата |
внутри основного периода. Маскиру |
|
ющие операторы присваивания очень удобны при решении задач, подобных задаче Власова или других, где имеют дело с частицами, так же как и при манипуляциях с цифрами. К сожалению, этих операторов нет на всех ЭЕМ и нет в ASA FORTRAN.
*) Программа была написана для CDC 3600 в Калифорнийском универси тете Сан-Диего и переработана для ICL KDF9 в Калэмской лаборатории. Обе машины используют 48-разрядный код. Программа также использова лась на CDC 6600 в Лоуренсовской радиационной лаборатории (Ливермор). У этой ЭВМ 60-разрядный код. Все три варианта Фортрана очень схожи.
§ 3 р Ч и с л е н н ы е м е т о д ы |
ИЗ |
Для расчета эйлеровского интервала, внутри которого лежит точка, используется деление целых чисел в версии Фортрана. Так как делитель — степень двух, то такое деление эквивалентно логическому сдвигу вправо на а разрядов, при котором зати рается дробная часть координаты. Если используется автокод ЭВМ. то применяются стандартные операции MASK и SHIFT. Очевидно, что решение уравнений движения (21) занимает малую долю полного времени счета. Основное время тратится на расчет распределения заряда, что обсуждается в и. 3.
2. Образование цепочек
Постоянно приходится иметь дело с двумя типами кривых: с кривыми, образующими замкнутую петлю на фазовой плоскости (тип 1), и с кривыми, которые являются непрерывными функция
ми при изменении аргумента х |
от —оо до +оо (тип 0). |
Кривая |
каждого типа определяется цепочкой последовательных |
точек: |
|
{%Ві Ев)> (З-сі» Ѵ<х) ••• |
(*£р1 Ѵр), (Хеі Ѵе)і |
(23) |
где символы В и Е относятся соответственно к начальной и конеч
ной точкам цепочки. За точкой (хв , |
ѵЕ) следует точка (хв, ѵв) |
в случае цепочки типа 1 и точка (хв + |
L , ѵв) — в случае цепочки |
типа 0, где L — длина периода. Условимся считать, что замкну |
|
тые петли (тип 1) имеют направление по часовой стрелке, а беско нечные кривые (типа 0) направлены слева направо. Это согла шение инвариантно относительно топологических преобразований фазовой плоскости. Если же кривые надо разрезать и соединить заново, чтобы получить новые замкнутые петли, то направление
некоторых |
участков кривых приходится менять на обратное |
(см. и. 5). |
|
Точки хранятся в памяти ЭВМ в ячейках со случайными номе
рами: |
|
|
|
^В> |
• • |
• ) Ірі ІЕ- |
(24) |
Введем функцию, определяющую преемственность точек |
|||
|
Sin = Г. |
(25) |
|
Эта функция определяет номер |
ячейки I" для точки, |
следующей |
|
в цепочке за точкой Г, и три массива bj, е$, rij, которые определяют соответственно номера ячеек для первой и последней точек на кривой і и число точек. Все неиспользованные точки объединяются, и из них образуется кривая о, или список «свободных мест», который характеризуется величинами Ъ0, е0, п0.
С течением времени каждая кривая вытягивается, и если мы хотим, чтобы точность расчета оставалась неизменной, то необ-
8 -0 1 2 3 6
114 |
Г л . 3 . М о д е л ь «в о д я н о г о м е ш к а » |
ходимо вводить новые точки. Мы можем ввести точку у между точками а и ß следующим образом:
Z7: = b0, b0: — S (b0), n0: = n0 — 1, nf. = tij + 1 |
(26) |
(изъять первую ячейку из списка «свободных мест» и включить ее в число ячеек, в которых хранятся точки кривой Cj) и
S (Іа): = Zv, S (1У): =./р |
(27) |
(вставить новую ячейку между ячейками Іа и Zp в цепочке). Новая точка у размещается в середине сегмента (а, ß), т. е. в общем случае
Xy: = j ( X a + Xf>), F7: = y ( F a + Fß). |
(28) |
Нужна некоторая осторожность в случае, если сегмент совпадает с (Е , В), но мы этого касаться не будем.
Точку ß, которая лежит между а и у на кривой /, можно изъять, а ячейку, в которой она хранилась, включить в список «свобод ных» мест:
В (Za): — ly, |
. —■Tij |
1, В (e0). — Zp, |
. — Zß, ^o* ■— |
1 * |
|
|
|
|
(29) |
|
3. |
Расчет электрического поля |
|
|
Чтобы рассчитать плотность заряда, нужно сначала вычислить |
||||
интеграл |
|
|
|
|
|
|
j / (X, v) dv |
|
(30) |
как функцию X. Так как фазовая плоскость поделена на области постоянной плотности, то / (х 0, ѵ) — ступенчатая функция ѵ для каждого значения х 0. Пусть контур Cj разделяет две области Rj и Zij-i с фазовыми плотностями fj и /;-_4. Область Rj должна оставаться справа, если двигаться по контуру Сj в положительном направлении. Определим
Mi — ij fj-1- |
(31) |
Предположим, что Сj пересекает вертикальную линию х = х0 произвольное число раз (возможно, и ни разу) в точках Vj,q, и определим функцию Sj,q, которая равна + 1 , если контур пере секает эту линию в направлении увеличения х, и равна —1 в про тивном случае. Тогда
f (х0, v) = — |
vj,q) Sj, q, |
(32) |
3 |
9 |
|
§ |
|
3. Численные методы |
115 |
||
где Ѳ (V) — ступенчатая |
функция Хевисайда |
|
|||
Ѳ (ѵ) = 1, |
V > |
О, |
Ѳ (к) = 0, V < 0, |
(33) |
|
и интеграл принимает |
вид |
|
|
|
|
/ (аг0. |
ѵ) = |
2 |
2 A/jty. Ä . g- |
(34) |
|
39
Врасчетах мы будем предполагать, что каждая кривая — много
угольник, т. е. ломаная линия. Нам нужно вычислить интеграл
хт+і |
|
j dv j dxf(x, v) |
(35) |
xm
в каждом из интервалов
Xm -А Ж Xm+ii 0 YTT<d L* |
(36) |
Добавляя новые точки в случае, если сегмент пересекает одну из вертикальных линий на границах интервалов, мы можем пред ставить каждый член в (35) в виде суммы по всем подсегментам, лежащим внутри интервала (36):
2 |
M j { x " - x ' ) ^ ± ! L , |
(37) |
по всем |
|
|
подсегментам |
|
|
где ж" следует за точкой х' |
(принимаются в расчет и новые точки). |
|
Нетрудно проверить, что выражение (37) приводит к правильным результатам, какова бы ни была форма сегментов.
Следовательно, вычислительная процедура довольно проста: мы движемся вдоль каждой кривой в положительном направлении и с помощью операций MASK и деления целых чисел (см. п. 1) определяем ж-интервалы, в которых помещаются две концевые точки каждого сегмента. Затем каждому интервалу, внутри которого имеется сегмент, приписывается в соответствии с форму лой (37) некоторое добавочное количество заряда, которое вычис ляется в режиме с плавающей запятой.
Когда кривые ломаные, этот расчет приводит к точной величи не полного заряда внутри каждого эйлеровского интервала, если не считать ошибки округления. Электрическое поле на границе каждого интервала легко затем рассчитать по интегральной фор муле
(АЕ)т = 4я {Qm — g0), ' (38)
где Qm — заряд внутри т-то интервала, д0 — однородно распре деленный положительный заряд фона, который подбирается таким,
8*
116 |
Гл. 3. Модель «водяного мешка» |
чтобы суммарный заряд в точности равнялся нулю, и
(39)
Уравнение (39) точное, но, чтобы наложить условие периодичности
L
(40)
о
мы должны предположить, что внутри каждого эйлеровского интервала поле Е (х) изменяется линейно. Эта аппроксимация также нужна для расчета значений электрического поля в лагранжевских точках x t, которые подставляются в уравнения (18).
Следовательно, на каждом шаге по времени два небольших фактора приводят к ошибкам в расчете электрического поля. Во-первых, мы предположили, что кривые Cj — ломаные линии, но это свойство кривых нарушается с течением времени, так как сегменты искривляются. Во-вторых, мы приняли, что величину Е (X) можно рассчитать, пользуясь линейной интерполяцией значений в узлах эйлеровской сетки, в то время как в действи тельности электрическое поле имеет более сложную структуру: величина cPEldx2 в каждой точке x t изменяется от одного постоян ного значения к другому. Первый тип ошибок можно свести к минимуму, если увеличить число лагранжевских точек на кон турах, а второй — если увеличить число эйлеровских интервалов. Имеется и третий источник ошибок при расчете движения конту ров, который зависит от величины шага по времени At. Этот ис точник рассматривается в § 4.
4.Подправление контуров
Впроцессе счета два контура или две части одного и того же контура часто сближаются друг с другом, пока площадь между ними не становится исчезающе малой на отрезке конечной длины.
Вэтом случае две кривые превращаются в одну «двойную» кри вую, которую не имеет смысла проходить дважды. Если величина / на каждой стороне «двойной» кривой одинакова, то, исключив совсем «двойной» участок кривой, мы не изменим состояния
системы.
Фиг. 9 иллюстрирует три возможных случая образования «двойных» кривых. Проще разобрать диаграмму на фиг. 9, а; с нее мы и начнем наш анализ. В этом случае цепочка Сt переделывается следующим образом: точку С соединяют с точкой D и образуют новую замкнутую кривую, начинающуюся в точке А и оканчиваю
