Особенности численного расчета траекторий заряженных частиц в импульсных электрических полях Текст научной статьи по специальности «Физика»

ISSN 0868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2014, том 24, № 3, c. 62-74

- МАТЕМАТИЧЕСКИЕ АНАЛИЗ И МОДЕЛИРОВАНИЕ В ПРИБОРОСТРОЕНИИ

УДК 537.534.7:621.319.7 © А. С. Бердников, Н. Р. Галль

ОСОБЕННОСТИ ЧИСЛЕННОГО РАСЧЕТА ТРАЕКТОРИЙ ЗАРЯЖЕННЫХ ЧАСТИЦ В ИМПУЛЬСНЫХ ЭЛЕКТРИЧЕСКИХ ПОЛЯХ

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

Кл. сл.: численное решение дифференциальных уравнений, трассировка заряженных частиц в электрических полях, артефакты численных алгоритмов

КРАТКОЕ ОПИСАНИЕ ПРОБЛЕМЫ

При трассировке заряженных частиц в электрических полях [1, 2] численно решаются уравнения движения

U(t)

1.0

m

m

m

d2 x_ dU ( x, y, z, t)

dt2

d2y dt2

d!z dt2

= ~ci

dx

dU ( x, y, z,

t)

dy

dU ( x, y, z, t) dz

III

(1)

Рис. 1. Непрерывно меняющееся высокочастотное напряжение, приложенное к электродам ионно-оптической системы (пример взят из [10])

где х, у, z — декартовы координаты частицы, х (, у ({), г (t) — траектория частицы, т — масса частицы, q — заряд частицы, и (х, у, г, t) — потенциал электрического поля, t — время. Для решения уравнений (1) обычно используются типовые рецепты численного решения обыкновенных дифференциальных уравнений [3-9]. В случае, когда напряжения на электродах представляют собой достаточно гладкую функцию времени (см. рис. 1), их применение не представляет особых технических проблем.

U(t)

1.0

0.5

t

•1.0

Рис. 2. Импульсное высокочастотное напряжение, являющееся цифровым аналогом напряжения с рис. 1

t

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

Целью настоящей работы является исследование работоспособности типовых алгоритмов численного интегрирования обыкновенных дифференциальных уравнений применительно к трассировке заряженных частиц в электрических полях, характеризуемых импульсными функциями времени со скачками (см. рис. 2). Для достаточно гладких правых частей дифференциальных уравнений используемые численные алгоритмы должны обеспечивать оговоренную численным методом локальную точность О (кк+1) (где h — шаг

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

Однако будет показано (см. следующий раздел), что для правых частей со скачками локальная точность интегрирования равна всего лишь О (h)« С • h , независимо от порядка использованного численного метода (где С — константа, пропорциональная скачку функции в правой части

1) Порядок метода в целом на единицу меньше, чем порядок локальной ошибки интегрирования одиночного шага: если 0(кк+1) — локальная ошибка на одиночном шаге, а N = Т/к = О(1/И) — полное число шагов интегрирования вдоль интервала 7 е [70,70 + Т], то накопившаяся в конце интервала интегрирования ошибка будет равна о(кк+1 )• О(\/к ) = о(кк ).

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

К счастью, в реальности ситуация не столь плоха. Более детальное рассмотрение показывает, что, во-первых, аномально большая локальная ошибка О (к)« С • к реализуется только для тех

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

равна регламентному значению О (кк+1), много

меньшему, чем О (к)« С • к . Поэтому в реальности накопившаяся в процессе интегрирования ошибка равна не ~ СТ, а ~ СТ • (к/Н), где Н — характерная временная длина одиночного импульса. (Здесь использована оценка ~ Т/Н для числа скачков импульсов на интервале времени Т ; эта величина при условии Т/к > Т/Н (т. е. когда к < Н) совпадает с числом тех критических интервалов длины к, на которых образуется аномально большая локальная ошибка интегрирования). Это означает, что вопреки первоначальному пессимистическому прогнозу ошибка интегрирования траекторий все-таки стремится к нулю с уменьшением шага к, хотя и очень медленно, и к тому же это происходит лишь начиная с достаточно мелких шагов интегрирования, удовлетворяющих условию к < Н .

Во-вторых, в силу того что компьютерный код для решения уравнения (1) на самом деле использует эквивалентную запись (1) как систему уравнений первого порядка:

U(t) 1.0-

0.5

-0.5

-1.0

тг

Рис. 3. Определение и упорядочивание критических точек времени переключения импульсного сигнала (импульсного напряжения, приложенного к электроду)

t

/1 12 h

t

■ = v.

. dy

dt = VJ;

dz dt

=v

q dU (x, y, z,t) dvy q dU(x, y, z, t) (2)

m

dx

dt

m

cy

dx

¥

dt

¿V, _ q ди (х, у, z, t) dt т дz

локальная ошибка О(к), обусловленная неточностями алгоритма, не рассчитанного на импульсную функцию в правой части уравнения, формируется лишь в уравнениях для скоростей Ух ),

Уу ^), V, ^). В уравнениях же для координат х^), ), , (t) накапливаемая локальная ошибка будет равна О (к2) (результат интегрирования по интервалу времени длины И скоростей vx ^), vy ^),

V, (^), вычисленных с ошибкой О(к)). Поэтому в итоге для скоростей накопившаяся ошибка численного интегрирования, действительно, равна ~ СТ • (И/И), но для координат накопленная ошибка ~ С 'Т •( И2 /И) ведет себя гораздо лучше (т. е.

существенно более обнадеживающим образом).

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

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

1. Перед началом трассировки для всех импульсных сигналов, подаваемых на электроды рассматриваемой системы, выстраивается упорядоченная последовательность временных точек

0 < ) < t2k) < ^) <... < Т, в которых происходит переключение импульсного сигнала (см. рис. 3).

2. Также перед началом трассировки все критические точки отдельных импульсных сигналов

^к) < t2k) < t3k) <... < Т объединяются в одну объединенную глобальную упорядоченную последова-

„ * * * ^

тельность 0 < ^ < t2 < ^ <... < Т критических временных точек переключения импульсных сигна-лов.2)

2) Эта процедура, впрочем, может выполняться динамически в процессе трассировки траекторий. Перед тем как трассировать траекторию на протяжении очередного интервала t е [t*, t*+1J, для текущей стартовой точки

t* программа "опрашивает" все импульсные сигналы,

в какой именно момент времени t* + At, ожидается

j k

очередное переключение соответствующего импульса. После этого устанавливается значение t*+1 =

= t* + min Atk, уравнения движения благополучно ин-

тегрируются на интервале t е

[t; f *+iJ,

и процедура по-

вторяется для новой стартовой точки t .

3. Процесс интегрирования траектории организован таким образом, что каждый интервал

< t < ¿*+1 обрабатывается независимо, а шаг интегрирования hj «h корректируется так, чтобы

* *

вдоль текущего интервала tj < t < tj+1 уложилось целое число шагов длины hj (но при этом не менее одного).

4. При трассировке траектории на интервале

* *

tj < t < tj+1 подпрограмма, управляющая возвратом

значения напряженности электрического поля в заданной точке пространства в заданный момент времени, должна возвращать правильные значения

напряженности электрического поля для моментов

* *

времени t = t.. и t = tj.+1. А именно, если запрашивается значение напряженности электрического

* *

поля в точке t = tn, а tn входит в список критических точек переключений импульсов, то подпрограмма должна иметь информацию, используется ли в качестве текущего интервала интегрирования интервал < t < либо интервал t¡* < t < tn+1.

В первом случае при запросе текущего значения

*

импульсного электрического поля в момент t = tn подпрограмма должна вернуть значение для t = - 0 (левый фронт импульса), во втором случае — для t = tn + 0 (правый фронт импульса).

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

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

электрического поля при переходе от интервала

* * * *

tj < t < tj+1 к интервалу tj+1 < t < tj+2.

Действительно, имеется ошибочное мнение, что если все точки tJ■ расположены регулярным образом, а шаг интегрирования h выбран таким образом, чтобы на каждом отрезке t* < t < tj+1 укладывалось целое число шагов длины h, то описанных выше проблем можно избежать. К сожале-

*

нию, это не так: при переходе через точку t = t.

**

либо левый интервал интегрирования t. - h < t < tj, либо правый интервал интегрирования

< t < + h, либо оба сразу получат неправильное значение импульсного электрического поля

для момента t = ^ (это зависит от того, какими

соображениями руководствовался программист при написании соответствующего фрагмента кода). Следовательно, по крайней мере один из интервалов интегрирования будет обработан неправильным образом, а результатом этого будет неправильное числовое значение, которое в конечном счете будет порождать ошибку порядка О(h) в окончательном результате. Совершенно ясно, что предложенная выше модифицированная версия стандартного алгоритма трассировки, позволяющая в зависимости от контекста процесса интегрирования траектории правильно принимать решение, какое именно значение должно быть приписано импульсному электрическому полю при t = tj, от этого недостатка свободна.

Обоснованию высказанных здесь тезисных положений посвящена оставшаяся часть работы.

АНАЛИЗ ОШИБКИ ИНТЕГРИРОВАНИЯ ТРАЕКТОРИИ, КОГДА СКАЧОК ИМПУЛЬСА НАХОДИТСЯ ВНУТРИ ШАГА ИНТЕГРИРОВАНИЯ

Рассмотрим модельный случай. Пусть выполняется одномерное движение заряженной частицы вдоль оси OX в однородном электрическом поле E(х) = const, и пусть в момент t = т происходит импульсное переключение электрического поля от значения —E0 в моменты t < т к значению + E0 в моменты t > т . Исследуем, как такой скачок поля влияет на результат, посчитанный типичным алгоритмом Рунге—Кутты [3-9, 11] при условии, что момент переключения электрического поля приходится в точку внутри текущего интервала интегрирования траектории заряженной частицы.

Аналитически точное решение уравнений движения для координаты заряженной частицы х (t)

и скорости заряженной частицы v (t) для описываемого случая имеет вид

/ч k — (яЕо/m)t, t <т;

v (t Н.,

:(t ) = •

+ (qE„/m)(t — т), t > т;

+ v0t — ( qE0/ m ) t2/2, х0+ v0(t — т) + (qE0/m)(t — т)2/2, t >

t < т;

т

Рис. 4. Нормированные решения системы уравнений движения (с разными моментами т релейного переключения электрического поля) вдоль нормированного интервала времени 0 < t < 1 для модельной задачи с импульсным переключением электрического поля: а — для скорости зараженной частицы, б — для координаты заряженной частицы

где введены обозначения VI = v0 -(qE0/т)т, Х0 = Х0 + V0т

-(qEo/т )т2/2 для координаты и скорости заряженной частицы на границе переключения импульсного электрического поля. Иными словами,

* ) =

0 < t < т;

V,-( qEo/т ) t, v0 +(qE0/т)(t - 2т), т < t < h;

(3)

0 < t < т;

Для системы обыкновенных дифференциальных уравнений вида у' = F (у, t) рассмотрим классический метод Рунге—Кутты [11], вычисляющий приближенное решение от точки t = t0 до точки

t = ^ + h :

* ) =

| х0 + v0t - (qE0/т)12/2, х0 + v0t + (qЕ0 /т)^2/2 - 2т + т2), т < t < h.

Без ограничения общности можно положить начальные условия х0 = 0 и v0 = 0 , после чего в решении (3) остается единственная по-настоящему значимая часть — скачок электрического поля. Это возможно, поскольку любой метод Рунге—Кутты, которые мы будет рассматривать далее, а) вычисляет решение, устроенное по принципу линейной функции х0 + v0t, аналитически точно; б) при наличии в решении линейной добавки х0 + v0t аналитически точно приплюсовывает ее к итоговому результату. На рис. 4 показаны нормированные решения (при qE0/ т = 1) на интервале 0 < t < 1 для разных значений

у1 = ^ ( y0, Ч ),

у2 = № (у0 + У/2, t0 + Ч2) ,

.Уз = № (у0 + Л/2, ^ + V2),

у4 = ^ (у0 + ^ t0 + h),

(4)

т е

{0,1/8, 1/4, 3/8,1/2, 5/8, 3/4, 7/8,1} .

у ^ + h )« у0 + 1 (у1 + 2 у2 + 2 у3 + у4 ). 6

Для достаточно гладких функций F (у, t) один шаг итераций, выполненный согласно алгоритму (4), должен продолжать приближенное решение у ^) системы обыкновенных дифференциальных уравнений у' ^) = F (у ^), t) от точки t = t0 до точки

t = t0 + h с погрешностью, не хуже чем О (h5) .

Однако это справедливо для достаточно гладких функций, тогда как наша функция, стоящая в правой части дифференциального уравнения, вполне разрывна. Разбирая отдельно случаи 0 < т < h|2 и h|2 < т < h , для рассматриваемого модельного случая получим с помощью (4) следующий ответ:

Ах

Рис. 5. Погрешность численного интегрирования методом Рунге—Кутты (4) одномерной системы уравнений движения для модельной задачи с импульсным переключением электрического поля. Показана нормированная разность между аналитическим и численным решением в конечной точке шага интегрирования t = t0 + И в зависимости от нормированного момента г' = г/И переключения электрического поля: а — для скорости заряженной частицы, б — для координаты заряженной частицы

б

т

К И у

:( И )>

2И

v0 qEo/ т ) при0 < г < И/2; 2И

v0 + qE0|m) при И/2 < г < И; И2

х0 + v0И--(qE0|m) при0<г <И/2;

(5)

И2

х0 + v0И +--(qE0|т) при И/2 < г < И.

6

Сравнивая между собой точное решение (3) и приближенное решение (5), получаем, что разность этих двух решений имеет вид:

V (И) - V (И) =

-1 (qE0/т)И(5 - 6г') при0 <г'< 12; 1

-( qE0/ т ) И (-1 + 6г') при1/2 <г'< 1;

х (И) - х (И) =

-(qEJт)И2 (1 - 2г'+ г'2) при 0 < г' < 12;

(qEJт)И2 (-1 + 6г'- 3г'2) при1/2 < г' < 1,

(6)

решения действительно несущественны, ничтожны и не имеют отношения к делу).

На рис. 5 показаны нормированные графики для функций (6) (при qE0/т = 1, И = 1) в зависимости от того, в какую именно точку внутри интервала интегрирования попадает момент релейного переключения электрического поля. Из проведенного анализа и формулы (6) можно сделать вывод, что вместо локальной погрешности

О (И5), гарантированной формулами Рунге—Кутты

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

локального шага по схеме (4) составляет О(И2) для

координат и всего О (И) для скоростей!

Учитывая, что локальные погрешности имеют тенденцию накапливаться, а не компенсировать друг друга, и что полное число шагов интегрирования траектории на интервале t е[0,Т ] можно оценить как «Т/И, при наличии локальной погрешности О (И) , имеется опасность накопления

систематической ошибки численного расчета, которая не стремится к нулю при уменьшении шага интегрирования.3"1 Такой же эффект наблюдается

3

где введена безразмерная нормировка г' = г/И . (В частности, из (6) следует, что значения констант х0 и v0, определяющих начальные значения, для рассматриваемого нами аспекта численного

3) То, что на самом деле ситуация не столь плоха, подробно разбирается в первом разделе статьи. Но этот факт следует отнести скорее к редкостной удаче, свя-

занной с использованием специфической формы уравне-

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

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

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

АНАЛИЗ ОШИБКИ ИНТЕГРИРОВАНИЯ ТРАЕКТОРИИ, КОГДА ВНУТРИ ШАГА ИНТЕГРИРОВАНИЯ ИМЕЕТСЯ ИЗЛОМ ЗАВИСИМОСТИ ПОЛЯ ОТ ВРЕМЕНИ

Для сравнения рассмотрим, что происходит, когда внутрь шага интегрирования по времени попадает не скачок, а излом напряженности электрического поля. В качестве модельного примера снова возьмем одномерное движение в релейно переключающемся в момент времени t = т электрическом поле, однако закон изменения электрического поля во времени теперь будет другим:

Е ( х, t ) =

[+Еа (t -т) при С < т, \-Еа (с -т) при с < т.

Аналитически точное решение уравнений движения для координаты заряженной частицы х (С)

и скорости заряженной частицы V (с) для описываемого случая имеет вид

* ) =

:(С ) =

V + 2(qEa|m)С(С - 2т), С < т;

1 2

V;-2^Еа/т)(С-т) , С > т;

хс + v0t + 6(qEa|m)С2 (с - 3т), С < т;

1 3

х;+ ^(С-т)-6 ( qEa|m )(С-т) , С > т,

(7)

где введены обозначения V = v0 - —(qEa|m)т2,

1

_3

х0 = х0 + ^т--(qEa|m)т для координаты и скорости заряженной частицы в момент переключения импульсного электрического поля.

Как и раньше, без ограничения общности можно положить начальные условия х0 = 0 и у0 = 0 , после чего в решении (7) остается единственная по-настоящему значимая часть — скачок электрического поля (это возможно, поскольку любой метод Рунге—Кутты с локальным порядком точности не ниже четырех обрабатывает не более чем кубические многочлены аналитически точно). На рис. 6 показаны нормированные решения (7) (при qEс/т = 1) на интервале 0 < С < 1 для разных зна-

чений т е

{0,18,14, 3/8, 12, 5/8, 3/4, 7/8,1} .

-0.4

а

х(С)

-0.05

-0.10

т = 1/2

-0.15

Т= 1/4

-0.20

т = 1/8

-0.25

= 0

-0.30

б

Рис. 6. Нормированные решения системы уравнений движения (с разными моментами т излома закона изменения электрического поля) вдоль нормированного интервала времени 0 < С < 1 для модельной задачи с негладким переключением электрического поля: а — для скорости зараженной частицы, б — для координаты заряженной частицы

С

После применения итерационного процесса (4) получаем следующий ответ для расхождения между истинным аналитическим решением (7) и приближенным решением в точке t = И , полученным с помощью алгоритма (4):

V (И) - V (И) =

--(qEa|m)И2г'(1 - 3г') при0 < г' < 12; -(qEa|m)И2 (2 - 5г'+ 3г '2) при1/2 < г' < 1;

х (И) - х (И) =

1 (qEa|m)И3г'(1 - 3г '+г' 2) при 0 < г' < 12;

(8)

3

1 ( qEa/m) И3 (1 -г')3

при 12 <г' < 1,

где введена безразмерная нормировка г' = г/И . На рис. 7 показаны нормированные графики для функций (6) (при qE0/т = 1, И = 1). Из формулы (8) можно сделать вывод, что, как и следовало ожидать, вместо локальной погрешности О (И5),

гарантированной формулами Рунге—Кутты для достаточно гладких функций, для электрических полей с изломом погрешность локального шага по схеме (4) составляет О (И3) для координат

ИСПОЛЬЗОВАНИЕ ЧИСЛЕННЫХ

АЛГОРИТМОВ С ПРЕДСКАЗАНИЕМ ШАГА

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

В качестве иллюстрации рассмотрим метод Рунге—Кутты—Фельдберга 5-го порядка с автоматическим контролем шага интегрирования [12]:

и О (И2)

для скоростей.

б)

0.08

0.04

0,02

Рис. 7. Погрешность численного интегрирования методом Рунге—Кутты (4) одномерной системы уравнений движения для модельной задачи с негладким переключением электрического поля.

Показана нормированная разность между аналитическим и численным решением в конечной точке шага интегрирования t = t0 + И в зависимости от нормированного момента г' = г/И переключения электрического поля: а — для скорости заряженной частицы, б — для координаты заряженной частицы

т

а

у1 = ИР(^ to),

у2 = ИР ( у0 + 1 ^ to + 4 И ],

у = ИР ^ у0 + 32 у1 + 32 У2, t0 + 3 И ],

у , ГУ 1932 У 7200 У 7296 у 12, ^

у4 = ИРI у0 +-у--у2 +-у3, t0 +— И I,

4 ^ 0 219/ 1 219/ 2 2197 3 0 13 )

у У 8341 у 32832 у 29440 у 845 у

у5 = ИРI у0 +-у1--у2 +-у3--у4, t0 + И

410^ 4104 2 4104 3 410^ 0

у У 6080 У 41040 У 28352 У 9295 У 5643 У 1,

у, = Ир I у„--у +--у.--у +--у.--у, ^ + — И

У 20520 20520 20520 20520 20520 0 2

У, ,ч У ( 902880 У 3953664 у

у+ И) « у0 +1-у1 +-у3 +

У0 ' ^ 7618050 761805^

3855735 у 1371249 у 277020 у

+--у л--у +--У л

7618050 7618050 7618050

Ау, ,ч ( 2090 у 22528 у

Ау (t0 + И)«!--у1 +-у3 +

У0 ; ^ 752400 752400

21970 у 15048 у 27360 _

+--у л--у--У I.

752400 752400 7524006 1

(9)

Здесь предпоследняя строчка вычисляет приближенное значение у ^) для момента времени t = t0 + И, а последняя строчка используется для контроля точности полученного значения у (t0 + И) . Впоследствии, зная оценку Ау (t0 + И) , можно обоснованно принять решение об уменьшении шага интегрирования, если оценка погрешности слишком велика, либо об увеличении шага интегрирования, если оценка погрешности чрезмерно мала.

Посмотрим, однако, как будут работать эти формулы применительно к импульсным электрическим полям. На рис. 8 показано нормированное расхождение между аналитическим решением и решением, посчитанным численно в соответствии с алгоритмом (9) для рассмотренной ранее модели движения заряженной частицы в постоянном поле со скачком в момент t = г . Как и раньше, абсолютная (не нормированная) локальная ошибка интегрирования для координаты х ^) имеет порядок

О(И2) ~ (qE0/т)• И2, а для скорости V(t) ошибка

интегрирования имеет порядок О (И) ~ (qE0/т) • И .

Графики на рис. 8 показывают зависимость истинной нормированной ошибки счета от момента переключения электрического поля г =г' И (0 <г' < 1). Одновременно на рис. 9 показана априорная оценка ошибки интегрирования, рассчитанная в соответствии с формулами (9) (последняя строчка). Из графиков видно, что априорная и апостериорная ошибки очень сильно расходятся и что автоматически посчитанная оценка погрешности на самом деле на порядок отличается от настоящей погрешности (при этом содержательны даже не столько сами графики, сколько масштабы шкал графиков). Детальный анализ интервалов, на который отрезок г е [0,И] разбивается реперными точками формул Рунге—Кутты— Фельдберга (9), дает для апостериорных погрешностей рассматриваемой модельной задачи следующие формулы:

б)

Ду

-1.0

т

б

Рис. 8. Погрешность численного интегрирования методом Рунге—Кутты—Фельдберга (9) с автоматическим выбором шага системы уравнений движения для модельной задачи с импульсным переключением электрического поля в зависимости от положения момента переключения электрического поля. Показана нормированная разность между аналитическим и численным решением в конечной точке шага интегрирования С = С0 + h в зависимости от нормированного момента т' = т/к переключения электрического поля: а — для скорости заряженной частицы, б — для координаты заряженной частицы

Ду

0.04 0.02

-0.02 -0.04 -0.06 -0.08 -0.10

0.2

0.4

|.»г

Ах

0.020

0.015 0.010 0.005

0.2

-0.0051

0.4

0.6

0.8

1.0 т

а

Рис. 9. Оценки погрешности численного интегрирования при использовании метода Рунге—Кутты— Фельдберга (9) с автоматическим выбором шага системы уравнений движения для модельной задачи с импульсным переключением электрического поля в зависимости от положения момента переключения электрического поля.

Показана нормированная априорная оценка разности между аналитическим и численным решением в конечной точке шага интегрирования С = С0 + к в зависимости от нормированного момента т' = т/к переключения электрического поля: а — для скорости заряженной частицы, б — для координаты заряженной частицы

V (к ) - V (к ) =

Е0к\ -

32+270т

2Е0к

2Е0к

Е0 к

135

-8176+12875т

12875 -95066+141075т

141075

-59+50т 25

при 0 <т'< 3/8;

при 3/8 < т' < 12;

при 12 <т' < 12/13;

при 12/13 < т' < 1;

х (h) - х (h) =

Е0 Ъ

Ео Ъ

Ео Ъ

-253 + 2160т'- 1080т 1080

'2 Л

-93667 + 205200т'- 102600т 102600

12 Л

-51 + 100т'- 50т 50

2

при 0 < т' < 38; при 3/8 < т' < 12/13; при 12/13 < т' < 1,

(10а)

в то время как априорные оценочные погрешности, вычисленные по формулам (9), равны:

ДУ (h):

Ах(h)'

1_ 80 - 929 17100 3461 188100

Е0h при 12/13 < т'< 1;

-Е0Ъ2 при 0 < т' < 3/8;

-Е0Ъ при 0 < т' < 3/8;

180 0

Е0Ъ при 3/8 <т' < 12; Е0Ъ при 12 < т ' < 12/13;

(11)

360

161 34200

—Е0 Ъ2 50 0

Е0Ъ2 при 3/8 <т' < 12/13; при 12/13 < т' < 1.

(11а)

Этот результат достаточно неожиданный и вполне неприятен, поскольку обычно алгоритмы с автоматическим контролем шага интегрирования (в частности, рассмотренный здесь алгоритм Рунге—Кутты—Фельдберга [12]) вполне обоснованно считаются надежным инструментом, позволяющим автоматически измельчать шаг интегрирования для интервалов с нерегулярным поведением правой части уравнения и/или с нерегулярным поведением решения уравнения. Однако, к сожалению, если нарушаются априорные предположения о непрерывности (гладкости) функций, на основе которых эти алгоритмы сконструированы, работоспособность алгоритмов оказывается под вопросом: алгоритм послушно отслеживает деликатные нерегулярности в поведении правой части и/или решения, но перед грубыми разрывами в правой части уравнения определенно пасует.

ВЫВОДЫ

При интегрировании траекторий заряженных частиц в импульсных электрических полях, к со-

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

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

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Хокс П., Каспер Э. Основы электронной оптики. Пер. с англ. В 2 т. М.: Мир, 1993.

2. Yavor M.I. Optics of charged particle analyzers // Ser. Advances of Imaging and Electron Physics. Amsterdam, Elsevier, 2009. Vol. 157.

3. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. 4-е изд. М.: БИНОМ. Лаборатория знаний, 2006. 636 c.

4. Калиткин Н.Н. Численные методы. М.: Наука, 1978. 512 c.

5. Демидович Б.П., Марон И.А., Шувалова Э.З. Численные методы анализа, 3-е изд. М.: Наука, 1967. 368 с.

6. Бабенко К.И. Основы численного анализа. М.: Наука, 1986. 452 с.

7. Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М.: Мир, 1990. 512 с.

8. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999. 685 с.

9. Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical recipes: the art of scientific computing. Third edition. Cambridge University Press, 2007. 1235 P. URL: (http://www.nr.com).

10. Андреева АД., Бердников А.С. Масс-спектро-метрические устройства на основе радиочастотных электрических полей с архимедовыми свойствами // Масс-спектрометрия. 2011. Т. 8, № 4. С. 293-296.

11. Kutta M. W. Beitrag zur näherungsweisen integration totaler differentialgleichungen // Zeitschrift für Mathematik und Physik. 1901. Vol. 46. P. 435-453.

12. Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. М.: Мир, 1980. 280 c.

Институт аналитического приборостроения РАН, г. Санкт-Петербург (Бердников А.С., Галль Н.Р.)

Физико-технический институт РАН

им. А.Ф. Иоффе, г. Санкт-Петербург (Галль Н.Р.)

Санкт-Петербургский политехнический университет (Галль Н.Р.)

Контакты: Бердников Александр Сергеевич, asberd@yandex.ru

Mатериал поступил в редакцию: 28.07.20(4

UDK 537.534.7:621.319.7

SPECIFICS OF NUMERICAL SIMULATIONS OF THE TRAJECTORIES OF CHARGED PARTICLES IN PULSED ELECTRIC FIELDS

A. S. Berdnikov1, N. R. Gall1'2'3

1 Institute for Analytical Instrumentation of RAS, Saint-Petersburg A.F. Ioffe Physico-Technical Institute of RAS, Saint-Petersburg Saint-Petersburg State Polytechnical University

The problem of numerical simulation of the trajectories of charged particles as applied to the pulsed electric fields, is considered. Typical numerical recipes imply that the right hand side function of the differential equations is smooth and can be differentiated up to necessary order; however, when applied to the functions with sharp jumps these recipes can result to serious inaccuracies of the calculations. The paper considers the true inaccuracy artifacts which exists for the pulsed functions (pulsed electric fields) when some "smooth" numerical recipes are used for trajectory integration. It also considers some simple modifications of the standard tracing algorithms which enable to eliminate these artifact effects.

Keywords: numerical solutions of differential equations, trajectory tracing in pulsed electric fields, artifacts of numerical algorithms

REFERENСES

1. Hawkes P.W., Kasper E. Principles of electron optics. London, 1989.

2. Yavor M.I. Optics of charged particle analyzers. Ser. Advances of Imaging and Electron Physics, Amsterdam, Elsevier, 2009, vol. 157.

3. Kutta M.W. Beitrag zur näherungsweisen integration totaler differentialgleichungen. Zeitschrift für Mathe-

matik und Physik, 1901, vol. 46, pp. 435-453.

4. Press W.H., Teukolsky S.A., Vetterling W.T., Flan-nery B.P. Numerical recipes: the art of scientific computing. Third edition. Cambridge University Press, 2007, 1235 p. URL: (http://www.nr.com).

Contacts: Berdnikov Alexander Sergeevich, asberd@yandex.ru

Article arrived in edition: 28.07.2014

HAYHHQE nPHEQPQCTPQEHHE, 2014, tom 24, № 3