2.4. ИНТЕГРИРОВАНИЕ УРАВНЕНИЙ ДВИЖЕНИЯ

При моделировании физических процессов в плазме может потребоваться 10000 частиц, взаимодействующих в течение 1000 шагов по времени. Это означает, что уравнения движения должны интегрироваться 10 млн. раз. Поэтому необходимо воспользоваться самыми простыми и быстрыми методами вычислений, какие только возможны, сохраняя при этом приемлемую точность расчетов.

Кроме того, наш выбор метода вычислений должен учитывать возможности памяти компьютера, которые будут использоваться для хранения ряда величин для каждой частицы. Если бы мы наблюдали за траекторией одной частицы, движущейся в заданном поле (так называемые траекторные вычисления), тогда целесообразно было бы использовать метод интегрирования по времени высокого порядка точности с сохранением информации о скорости и координате за несколько предшествующих временных шагов. Однако минимальная информация, необходимая для расчетов, — это скорость и положение частицы, два машинных слова на частицу или 20000 слов для нашей задачи. Сохранение значений для нескольких предыдущих шагов существенно увеличит это число. Использование же метода высокого порядка, например Рунге — Кутта, увеличивает число операций, необходимых для каждой частицы. Следовательно, почти всегда мы будем выбирать для использования минимальную информацию и наибыстрейший метод.

Один из наиболее часто используемых методов интегрирования называется методом с перешагиванием. Два

Рис. 2.4. Схема метода интегрирования с перешагиванием, показывающая центрирование по времени силы при вычислении и при вычислении

Рис. 2.5. Плоскость показывающая силу нормальную к и приводящую к вращению без изменения модуля скорости с при

дифференциальных уравнения первого порядка, интегрируемые по отдельности для каждой частицы, имеют вид

где - сила. Эти уравнения заменяются конечно-разностными уравнениями

Схема интегрирования по времени и обозначения показаны на рис. 2.4, который поясняет идею временного центрирования. Компьютер будет вычислять используя хотя заданы в разные моменты времени. Для использования этого метода необходимо учесть два момента. Во-первых, начальные условия для скоростей и координат частиц, заданные в момент должны быть изменены ради возможности применения указанного способа интегрирования; в дальнейшем значение будет преобразовываться в с использованием вычисленной при силы Во-вторых, необходимо согласование способов вычислений кинетической и потенциальной энергии для проверки закона сохранения энергии в любой момент времени.

Метод с перешагиванием имеет погрешность, которая исчезает при стремлении Мы будем использовать этот метод интегрирования почти во всех наших программах, потому что он прост (легок для понимания и требует минимума запоминаемой информации) и удивительно точен, как будет показано позднее. Применяя этот метод к интегрированию простого гармонического осциллятора с круговой частотой (в следующем параграфе), получим, что при выполнении условия амплитудная ошибка равна нулю, а фазовый сдвиг на один шаг дается выражением

Из вида основного члена, определяющего фазовую ошибку, следует, что при условии возможно моделирование колебаний и волн в плазме с приемлемой точностью в течение нескольких десятков периодов.

Силу можно разбить на две части:

Электрическое поле и магнитное поле В должны вычисляться в точке нахождения частиц. Следовательно, используя пространственную сетку, необходимо интерполировать из узлов сетки в место нахождения частицы, точно так же, как это делается для определения координат частиц. Такая интерполяция рассматривается ниже, в § 2.6, где будет также показано, что электрическое поле, действующее на частицу, зависит не только от расстояния между частицами (физический эффект), но и от положения частицы внутри ячейки (нефизический эффект).

Рассмотрим одномерное расположение частиц вдоль оси х. Под действием одномерного магнитного поля направленного по оси каждая частица приобретает две составляющие скорости: Сила как показано на рис. 2.5, вызывает простое вращение, т. е. не изменяется по величине. Однако сила изменяет величину (т. е. пусть Следовательно, физически обоснованная схема с центрированием по времени будет следующей (вводятся вспомогательные переменные

Вычисляются скорости в момент (половина ускорения):

Затем вычисляются скорости в момент (вращение):

Конечные величины даются соотношениями (половина ускорения)

Угол вращения определяется формулой

т. е. , а циклотронная частота дается формулой

с сохранением знака Эта схема [Boris, 1970] тщательно разрабатывается ниже, в гл. 4.

Возникает одно усложнение в момент когда начальные

условия заданы в один и тот же момент времени. Следовательно, в начальный момент необходимо преобразовать значение используя величину Такой переход рассматривается ниже, в § 2.7, рис. 2.10.