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

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

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

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

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

движения

d 2 x

m—7- = - q dt2

d 2 y

m = - q dt2

д U ( x , y , z ,

дx д U (x, y, z,

,

d 2 z

m   = - q dt2

д y д U ( x , y , z ,

д z

,

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

,

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

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

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

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

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

К счастью, в реальности ситуация не столь плоха. Более детальное рассмотрение показывает, что, во-первых, аномально большая локальная ошибка O ( h ) ~ C • h реализуется только для тех интервалов времени, которые попали на границу скачка функции, стоящей в правой части уравнения, в то время как для остальных интервалов локальная ошибка интегрирования по прежнему равна регламентному значению O ( h^{k + 1} ) , много меньшему, чем O ( h ) ~ C • h . Поэтому в реальности накопившаяся в процессе интегрирования ошибка равна не ~ CT , а ~ CT • ( h(H ) , где H — характерная временнáя длина одиночного импульса. (Здесь использована оценка ~ T H для числа скачков импульсов на интервале времени T ; эта величина при условии T/h > T]Н (т. е. когда h <  Н ) совпадает с числом тех критических интервалов длины h , на которых образуется аномально большая локальная ошибка интегрирования). Это означает, что вопреки первоначальному пессимистическому прогнозу ошибка интегрирования траекторий все-таки стремится к нулю с уменьшением шага h , хотя и очень медленно, и к тому же это происходит лишь начиная с достаточно мелких шагов интегрирования, удовлетворяющих условию h <  Н .

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

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

dxd

— = v ; ^— — v ;

dt    x dt dz dt  v

z ;

dvx _ q dU (x, У, z, t) dvy _ q dU (x, y, z, t)

----—---; ----—;

dt    m    dx      dt   m d vz _ q 5 U (x, y, z,r)

• dt     md

локальная ошибка O ( h ) , обусловленная неточностями алгоритма, не рассчитанного на импульсную функцию в правой части уравнения, формируется лишь в уравнениях для скоростей v x ( t ) , v y ( t ) , v _z ( t ) . В уравнениях же для координат x ( t ) , y ( t ) , z ( t ) накапливаемая локальная ошибка будет равна O ( h ² ) (результат интегрирования по интервалу времени длины h скоростей v x ( t ) , v y ( t ) , v _z ( t ) , вычисленных с ошибкой O ( h ) ). Поэтому в итоге для скоростей накопившаяся ошибка численного интегрирования, действительно, равна ~ CT • ( h/ H ) , но для координат накопленная ошибка ~ C' T • ( h²/H ) ведет себя гораздо лучше (т. е. существенно более обнадеживающим образом).

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

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

• 1.    Перед началом трассировки для всех импульсных сигналов, подаваемых на электроды рассматриваемой системы, выстраивается упорядоченная последовательность временных точек 0 < 4 k ) <  1 2 ^k ) <  t ^{( k} ) < ... <  T , в которых происходит переключение импульсного сигнала (см. рис. 3).

• 2.    Также перед началом трассировки все критические точки отдельных импульсных сигналов t 1 ^{( k} ) <  1 2 ^k ) <  t ⁽ k ) < ... <  T объединяются в одну объединенную глобальную упорядоченную последовательность 0 <  t * <  1 2 <  t * < ... <  T критических временных точек переключения импульсных сигна-лов.²⁾

• 3.    Процесс интегрирования траектории организован таким образом, что каждый интервал t * <  t <  t ** _{+ 1} обрабатывается независимо, а шаг интегрирования h j ® h корректируется так, чтобы вдоль текущего интервала t * <  t <  t * _{+ 1} уложилось целое число шагов длины h _j (но при этом не менее одного).

• 4.    При трассировке траектории на интервале t * <  t <  t * _{+ 1} подпрограмма, управляющая возвратом значения напряженности электрического поля в заданной точке пространства в заданный момент времени, должна возвращать правильные значения напряженности электрического поля для моментов времени t = t * и t = t * _{+ 1}. А именно, если запраши-

• вается значение напряженности электрического поля в точке t = tn, а tn входит в список критических точек переключений импульсов, то подпрограмма должна иметь информацию, используется ли в качестве текущего интервала интегрирования интервал t**_1 < t < t* либо интервал t* < t < t‘n+1. В первом случае при запросе текущего значения импульсного электрического поля в момент t = t** подпрограмма должна вернуть значение для t = t* - 0 (левый фронт импульса), во втором случае — для t = t* + 0 (правый фронт импульса).

Очевидно, что в таком модифицированном виде алгоритм трассировки будет обеспечивать для импульсных функций ту же декларированную глобальную точность O ( h^k ) , которая столь привлекательна для гладких функций. При этом существенной особенностью алгоритма являются переключения правила вычисления напряженности электрического поля при переходе от интервала t * <  t <  t ** _{+ 1} к интервалу t * _{+ 1} <  t <  t * _{+ 2}.

Действительно, имеется ошибочное мнение, что если все точки t* расположены регулярным образом, а шаг интегрирования h выбран таким образом, чтобы на каждом отрезке t* < t < t*+1 укладывалось целое число шагов длины h , то описанных выше проблем можно избежать. К сожалению, это не так: при переходе через точку t = t** либо левый интервал интегрирования t* - h < t < t*, либо правый интервал интегрирования t* < t < t* + h, либо оба сразу получат неправиль- для момента t = t** (это зависит от того, какими соображениями руководствовался программист при написании соответствующего фрагмента кода). Следовательно, по крайней мере один из интервалов интегрирования будет обработан неправильным образом, а результатом этого будет неправильное числовое значение, которое в конечном счете будет порождать ошибку порядка O(h) в окончательном результате. Совершенно ясно, что предложенная выше модифицированная версия стандартного алгоритма трассировки, позволяющая в зависимости от контекста процесса интегрирования траектории правильно принимать решение, какое именно значение должно быть приписано импульсному электрическому полю при t = t*, от этого недостатка свободна.

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

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

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

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

x ( t ) =

/ 0 ^- ( qE 0 Im ) t ,        t < ^т ;

_ V ‘ + ( qE 0/ m )( t - T ) , t > т ;

x 0 ⁺ v 0 t ^- ( qE 0 m i ) t 7^2,

_ x0+ v ‘(t- т) + (qE0 /m)(t- т )2/2, ное значение импульсного электрического поля

t <  т ;

t >  т,

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

где введены обозначения v ‘_ v ₀ - ( qE ₀ ]т ) т , x 0 _ x ₀ + v₀т - ( qE ₀ ]m ) т ²/2 для координаты и скорости заряженной частицы на границе переключения импульсного электрического поля. Иными словами,

v (t )_

_ v о ^- ( qE о/ m ) t ,         ⁰ < t < ^т ;

v ₀ + ( qE ₀/ m )( t - 2 т ) , т <  t <  h ;

x ( t ) _                                                 (3)

x 0 + v 0 1 — ( qE 0 Im ) 1 ²/2,              0 <  t <  т ;

x ₀ + v ₀ 1 + ( qE ₀ /m ) ( 1 ²/2 - 2 tт + т ² ) , т <  t <  h .

Без ограничения общности можно положить начальные условия x ₀ _ 0 и v ₀ _ 0, после чего в решении (3) остается единственная по-настоящему значимая часть — скачок электрического поля. Это возможно, поскольку любой метод Рунге—Кутты, которые мы будет рассматривать далее, а) вычисляет решение, устроенное по принципу линейной функции x ₀ + v ₀ 1 , аналитически точно; б) при наличии в решении линейной добавки x ₀ + v ₀ 1 аналитически точно приплюсовывает ее к итоговому результату. На рис. 4 показаны нормированные решения (при qE ₀ ]m _ 1) на интервале 0 <  t <  1 для разных значений т е { 0,1/8,1/4, 3/8,1/2, 5/8, 3/4, 7/8,1 } .

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

y = hF (Я, t о ), y _ hF(yo + yj2,10 + hl2),

У з = hF ( y o + y ₂/ 2, t o + h /2 ) ,                (4)

yA = ^hF ( y O ⁺ y 3 , t 0 ^{+ h} ) ,

y(t + h)« yo +1 (y + 2y + 2Уз + yA ).

Для достаточно гладких функций F ( y , t ) один шаг итераций, выполненный согласно алгоритму (4), должен продолжать приближенное решение y ( t ) системы обыкновенных дифференциальных уравнений y ' ( t ) = F ( y ( t ) , t ) от точки t _ 1 ₀ до точки t = 1 ₀ + h с погрешностью, не хуже чем O ( h ⁵ ) . Однако это справедливо для достаточно гладких функций, тогда как наша функция, стоящая в правой части дифференциального уравнения, вполне разрывна. Разбирая отдельно случаи 0 <  т <  h /2 и h /2 <  т <  h , для рассматриваемого модельного случая получим с помощью (4) следующий ответ:

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

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

v (h) ~ <

v ₀

v ₀

2 h

-—( qE ₀ (m ) при0 <  т <  h /2;

2 h

+ — ( qE ₀/ m ) при h /2 <  т <  h ;

x (h) ~ <

x ₀

x ₀

h ²

⁺ v о ^{h -} у ( qE о / m )

h ²

⁺ v о ^h+~r ( qE о/ ^m )

при0 <  т <  h /2;

при h / 2 <  т <  h .

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

v (h)-v (h) =

^- 3 ( qE о/ ^m ) ^h ( 5 ^- 6 ^Т ')

3 (qEо /m)h (-1 + 6Т') x (h) - x (h) =

при0 <  Т <  1/2;

при 1(2 <  Т <  1;

- ( qE ₀ jrn ) h ² ( 1 - 2 т ' + т'² )    при0 <  Т< 1/2;

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

На рис. 5 показаны нормированные графики для функций (6) (при qE ₀ / m = 1, h = 1) в зависимости от того, в какую именно точку внутри интервала интегрирования попадает момент релейного переключения электрического поля. Из проведенного анализа и формулы (6) можно сделать вывод, что вместо локальной погрешности O ( h ⁵ ) , гарантированной формулами Рунге—Кутты для достаточно гладких функций, для разрывных (импульсных) электрических полей погрешность локального шага по схеме (4) составляет O ( h ² ) для координат и всего O ( h ) для скоростей!

Учитывая, что локальные погрешности имеют тенденцию накапливаться, а не компенсировать друг друга, и что полное число шагов интегрирования траектории на интервале t е [ 0, T ] можно оценить как « T/h , при наличии локальной погрешности O ( h ) , имеется опасность накопления систематической ошибки численного расчета, которая не стремится к нулю при уменьшении шага интегрирования.³⁾ Такой же эффект наблюдается

—( qE ₀/ m ) h ² ( - 1 + 6 т '- 3 т ² ) при^2 < т <  1,

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

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

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

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

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

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

+ E a ( t - т ) при t <  т ,

- E a ( t - т ) при t <  т .

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

v(t ) = <

x (t) = <

v ₀

v 0

+ 2 ( qEalm ) t ( t ^{- 2 т} ) , t < ^т ;

^- 2 ( qEalm )( t ^{- т} ) , t > ^т ;

x ₀

x‘

+ v 0 t + 1 ( qEalm ) t ² ( t ^- 3 ^т ) ,     t < ^т ;

⁺ v o ( t ^{- т} )--( qEalm )( t ^{- т} ) , t > ^т , 6

где введены обозначения

^x ‘= ^x 0 ⁺ v 0 ^т — 3 ( qEalm ) ^{т 3}

v 0 = v 0 ^- ;2 ( qE al m ) ^{т 2},

для координаты и ско-

рости заряженной частицы в момент переключе-

ния импульсного электрического поля.

Как и раньше, без ограничения общности можно положить начальные условия x 0 = 0 и v 0 = 0 ,

после чего в решении (7) остается единственная

по-настоящему значимая часть — скачок электрического поля (это возможно, поскольку любой метод Рунге—Кутты с локальным порядком точности не ниже четырех обрабатывает не более чем кубические многочлены аналитически точно). На рис. 6 показаны нормированные решения (7) (при qE ₀ / m = 1) на интервале 0 <  t <  1 для разных значений т е { 0,18,14, 3/8,12, 5/8, 3/4, 7/8,1 } .

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

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

v (h) - v (h) =

• ^- 3 ( qEa/m ) ^hT '( 1 - ^{3 т} ' )      при0 <  ^T < 12;

• 1    ( qE a/ m ) h ² ( 2 - 5 т ' + 3 т '²) при1/2 <  т* <  1; x ( h ) - x ( h ) =

  - 1 ( qEa/m ) h ^Зт ' (1 - 3 т '+ т '²) при 0 <  Т <  1/2;

  - ( qEa/m ) h ³ ( 1 - т ')            при^2 <  т <  1,

где введена безразмерная нормировка т ' = т/h . На рис. 7 показаны нормированные графики для функций (6) (при qE 0/ m = 1, h = 1). Из формулы (8) можно сделать вывод, что, как и следовало ожидать, вместо локальной погрешности O ( h ⁵ ) , гарантированной формулами Рунге—Кутты для достаточно гладких функций, для электрических полей с изломом погрешность локального шага по схеме (4) составляет O ( h ³ ) для координат и O ( h ² ) для скоростей.

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

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

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

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

б)

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

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

У = hF(й,t0), y2 = hF ^ yo + 4 yi, t0 + 4 h j,

-   ,Л- з -   9 -     3, y. = hF yn + — y. +— y., tn + -h 3 I 0 32 1 32 2 0 8

У 4 = hF У 0

1932 -

+--- У

2197 1

7200 -    7296 -      12 j

-----У? + У з , tn +— h , 2197     2197 ⁰ 13 J

У 5 = hF У 0

8341 -   32832 -    29440 -    845 -

+ y,--y 7 +-- Уз--y ₄, tn + h

4104 1   4104 2   4104 3 4104 4 0

6080 -   41040 -    28352 -    9295 -    5643 -      1 ,j

-------У 1 +------- У 2-- У з +------- У д-- У5 , t o + ^-h , 20520    20520    20520    20520    20520   ⁰ 2 J

^y ( t 0 ^{+ h} ) ~ Уп ⁺

902880 - y

7618050 1

+

3953664 -

--------У з +

7618050 3

3855735 - _ 1371249 -    277020 - j

⁺ 7618050 ^{У 4} 7618050 ^{У 5 +} 7618050 ^{У 6} J ,

_

2090 -

--------У 1 +

752400 1

22528 -

-------У з +

752400 3

21970 -    15048 -    27360 - j

+-- y ₄--y--y₆ .

752400    752400    752400 J

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

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

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

б)

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

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

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

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

^E 0 ^h

—

32 + 270 т ' ^

135 J

v ( h ) — v ( h ) = <

f— 8176 + 12875 т '

2 E h

0 I 12875

f— 95066 + 141075 т '

2 E h

⁰ I 141075

^E 0 ^h

f — 59 + 50 т ' ^

I 25 ^J

при 0 <  Т <  3/8;

при 3/8 <  Т <  1/2;

при 1/2 <  Т <  12/13;

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

f-253 + 2160т'- 1080т '2)               „ , E 0 h ----------------------- при 0 < т’< 3 8; (          1080          J                         ' f-93667 + 205200r'- 102600т'2)       _   ,                  , х (h)-x (h ) = • E0 h 2 ----------------------------- при 3 8 < 1 < 12 13;       (10a) (          102600          ) f-51 + 100т'- 50т '2)                                  , E0 h 2 -------50-------J             при 12/13 < т < 1, в то время как априорные оценочные погрешности, вычисленные по формулам (9), равны:

	—Eoh 180 0	при	0 <  1 < 3/8;
Л г ( h ) .1	929 - 92^- Eoh 17100 0 3461 Eh 188100 0	при при	3/8 <  т < 12; 12 <  I <  12/13;	(11)
	- —Eoh 25 0	при	12/13 <  т' <  1;
	1 360 E 0 h	при	0 <  т' <  3/8;
А х ( h ) ~	- -161- Eh ’ 34200 0	при	3/8 <  т' <  12/13;	(11a)
	—Eoh2 [ 50 0	при	12/13 <  т' <  1.

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

ВЫВОДЫ

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

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