CC BY
20
Доклады БГУИР
2012 № 5 (67)
УДК 519.62:621.385
ЭФФЕКТИВНЫЙ АЛГОРИТМ ЧИСЛЕННОГО ИНТЕГРИРОВАНИЯ УРАВНЕНИЙ ДВИЖЕНИЯ КРУПНЫХ ЧАСТИЦ В ПРИБОРАХ СВЧ
А О. РАК
Белорусский государственный университет информатики и радиоэлектроники П. Бровки, 6, Минск, 220013, Беларусь
Поступила в редакцию 18 июня 2012
Сформулирован алгоритм численного интегрирования уравнений движения крупных частиц. Алгоритм представлен в трех формах: явной, предсказания-коррекции и модифицированной, предназначенной для моделирования движения крупных релятивистских частиц в электромагнитном поле. Проведен анализ порядка точности предложенного алгоритма. Представлены результаты сравнения со стандартными алгоритмами. Показана эффективность сформулированного алгоритма.
Ключевые слова: численное интегрирование, уравнение движения, временной слой, четные и нечетные производные, предсказание, коррекция.
Введение
При моделировании электронно-вакуумных приборов необходимо численно решать уравнения движения крупных частиц. Особо важную роль скорость и точность этих вычислений играют при использовании Р1С-методов, когда требуется многократно вычислять траектории и скорости миллионов частиц. В работе представлен эффективный метод численного интегрирования уравнений движения, основанный на принципах, описанных в [1].
Движение крупной частицы описывается в соответствии со вторым законом Ньютона:
F
а = —. т
Основные параметры, описывающие состояние частицы: х(г) - координата частицы, ), а(г) - соответственно ее скорость и ускорение. Эти величины связаны следующим образом:
X = и;х = и = а.
Формулировка алгоритма
Значение функции, ее четных производных и ее нечетных производных определим на разных временных слоях [1], как это показано на рисунке. На основании этого:
xt +Дt xt
= » 1 , (1)
Дt г+2Д
2
■+дг = 1 Д , (2)
г+-д 2
и 1 -ч
t+1д/ 2
= а
1 , t+-дt
-Дt 4
2
1
и 1 = и + —а 1 Дt.
/+—дt 2 t+—дt
2 ^ 4
Конечно-разностная схема
Значение а 1 экстраполируем по значениям ускорения с двух предыдущих времен-
t +1 дt
4
ных шагов (а/-д/ и аt):
а 1 д = 1 (5at - at-Дt).
4
(4)
Подставив (4), (3) в (2) получим выражение для хи
1
х+д = х Дt+^ (5а- at-дt) Д/2.
Для уточнения значения хг+дt интерполируем а 1 по значениям аг и а(+м (xt+дt):
t +-д/
4
а 1 д = 1 (at+Дt + ^ ) .
t+—Дt 4
л 1
(5)
Подставив (5), (3) в (2) получим следующее выражение для коррекции координаты:
= х Д/+8 (а+д.+3at )Дt 2.
Теперь найдем выражение для с использованием at+дt (х{+д):
(6)
и/+дt 1 А, t+ -дt
_^ = а
1 3
1 , t+ -Дt
-Дt 4
2
(7)
t+- дt 4
4
=4 (3at+д/+а). Используя (1) и (8), приведем (7) к следующему виду:
(8)
4+ Д/Дt = xt+Д/ - X + Г (3at+ Д/ + а) Д/2.
(9)
а
Окончательно алгоритм в форме предсказания-коррекции формулируется следующим образом:
Предсказание:
х+м = х г + Д + 1 (5а - а-Д, ) Д2 . (10) Коррекция: 1
Ч+д,Д = Х+ы - х + 8(3а»+д» + а) Д2 • (12)
Количество итераций коррекции может контролироваться по разности между значениями х гна разных шагах коррекции.
Используя (10) для определения координаты и подставив (10) в (12), приведем алгоритм к явной форме:
= х, + ид + -( а+д. + ^ )д 2, (11)
= х+ид+-(5а - а-д, )д \ (13)
8
1
Ч+д, = + 8 (+ 6а, - )д • (14)
Алгоритмы (10)-(12) и(13), (14) могут быть использованы только в том случае, когда аг не является функцией от и,. При движении электронов в магнитном поле либо с релятивистскими скоростями это условие не выполняется. В таких случаях необходимо также предсказывать значение :
4+д = а 1 =1 (3а, - а(Д,) ,
и,+д, = и, + 2 (3а, - а,-д, )Д • (15)
Для коррекции значения (15) может быть использовано выражение (14). Полученный алгоритм весьма близок к алгоритму Бимана (Бимана-Шофилда). Для удобства сравнения оба алгоритма представлены в табл. 1.
Определение порядка точности алгоритма
Определим порядок точности предложенного алгоритма. Для этого выразим с помощью ряда Тейлора основные выражения. Для выражения (2) представим х,+д, их, в следующем виде:
х+д=х 1 д 1 д =х 1 д + 2и 1 д д,+8а 1 д д,г+48ь 1 д д,+о (д4), (1б)
,+-д/+—д/ ,+—д, 2 /+-д/ 8 ,+—Д, 48 ,+—Д, 4 '
2 2 2 2 2 2
х. = х , , = х , -1 и , Д, + 1 а , Д,2 - — Ь , Д, + 0 (д,4) . (17)
1 ,+_Д1-_Д1 ,+1дг 2 l+_Дl 8 l+_Дl 48 l+_Дl
2 2 2 2 2 2
Вычтем (17) из (16):
:+и 1 Д, + — Ь 1 Д,3 +0(Д,4) . (18)
_ Д , О Л l+_Д , V '
х+д, = х +и
д 1 1
,+—д, 24 ,+—Д,
2 2
Выполним аналогичные операции для выражения (3):
, + - а , At + — Ь , At2 + О (к,3) ,
— к, А ,+-к, 11 ,+—к, V ) '
и 1 = и 1 1 = и 1 +—а — —
,+—к, ,+—к,+—к, ,+-к, 4 ,+-к, 32 ,+—к,
2 44 4 4 4
, -—а , А, + — Ь , А,2 + О (А,3) .
и — Л, Л ,+—Аt 11 ,+—к, V '
и, = и — — = и — — а — —
,+—к,—к, ,+—к, 4 ,+—к, 32 ,+—к,
4 4 4 4 4
После вычитания получим:
+—а . к, + О (к,3) . П9)
1 ,+—к, V / у 7
и — = и, + — —
,+—к, 2 ,+—к,
2 ^ 4
Для (4) имеем:
— = а, + — Ь,к, + О (к,2) , (20)
а — Л л
,+—к, 4
4
= а, - Ь к, + О (к,2) . (2!)
Используя представления (20) и (2—), выразим а,+к двумя способами: первый - умножим (20) на 4 и прибавим (2—); второй - умножим (20) на 3 и прибавим (2—). В результате получим следующие выражения:
—к= 4 (5а, - а, к,) + О (к,2), (22)
4 ' 4
,+!к, = 3 (4а, - а,-к,) - —2 Ь,к + О (к,2) . (23)
Подставляя поочередно (22), (23) в (—9), а результат в (—8), получим выражения для предсказания координаты частицы в следующем виде:
х+к, = х+ик+—(5а, - а-к)к 2 + 2т ЬА, 3 + О (к 4)
8 24 (24)
= х, + и,к, + 8 (5а, - а,^ ) к,2 + О (к,3)
X= X + и,к, + 6 (4а, - а,-к) к,2 + О (к,4) . (25)
Выражение (25) соответствует алгоритму Бимана. Как видно из (24), предложенный алгоритм определяет позицию с третьим порядком точности, однако при малых значениях к, и
— з
Ь, (зависит от характера движения) слагаемое — к, может стремиться к нулю (например,
при равноускоренном движении). Также следует отметить, что в (24) используется коэффициент —8 = 0Д25, являющийся конечной десятичной дробью, в отличие от коэффициента —6 в алгоритме Бимана, приводящего к появлению ошибки округления из-за замены бесконечной десятичной дроби конечной.
Аналогичным образом можно определить порядок точности для оставшихся выражений.
Для (6) и (9):
х+к = х+ик,+—(а+к + 3а )к2 + ^ Ьк,ъ+О (к,А), (26)
8 24
и,+д,А = X+д, - X +1 (За,+д, + а ) д2 - ^ ьА3 + О (а/4) .
(27)
После подстановки (26) в выражение (27) (которое имеет второй порядок точности) получим выражение для (14), обеспечивающее третий порядок точности:
+8 (За+д + 6а, - а-д)А+О (д 3).
и,+д = и, + 81 За ■
Результат оценки порядка точности для всех форм алгоритма представлен в табл. 1.
Как видно из таблицы, алгоритм Бимана определяет координату с четвертым порядком точности, а скорость - с третьим, т.о. алгоритм имеет общую погрешность третьего порядка. Предложенный нами метод имеет второй порядок точности в форме предсказания-коррекции и третий порядок в явной форме и модифицированной форме предсказания-коррекции. При этом, как было показано ранее, предложенный алгоритм более устойчив к ошибкам округления за счет использования коэффициента 18 вместо 1/6 . Также отметим, что достоинством алгоритма Бимана является хорошее сохранение энергии и низкая чувствительность к ошибкам округления, благодаря тому, что в алгоритме вычисляется разность значительно отличающихся чисел (4а, - аг-д,). В предложенном алгоритме вычисляется разность между числами, которые сильнее отличаются по значению (5а, - аг-д,).
Таблица 1. Предложенный алгоритм и алгоритм Бимана
Предложенный алгоритм
Алгоритм Бимана
Форма предсказания-коррекции
Предсказание:
х,+а = X + и,д +1 (5а, - а,-а)а2 + О (д3) 8
Коррекция:
X,+а = X, + и,д +1 (а,+д, + За, )д 2 + О (д 3) 8
°'+а= Х±1-Х' + 8 (За+Д+ а )Д + О (А2)
Предсказание:
х,+а = X + и,д + 6 (4а, - а,-д)д 2 + О ( д4 ) Коррекция:
х,+а = х, +и,д + 6 (а,+д + 2 а, )д 2 + О (д 4 )
и,+А = (х1д_-^ + 6 (2а,+д+ а, )д, + о (а3) А 6
Явная форма
X,+а = X, +и,д +1 (5а, - а,-а)а2 + О (д3) 8
и,+А = и, +1 (За,+д, + 6а, - а,)д + О ( д3)
X,+а = X + и,д + 6 (4а, - а,-д)д 2 + О ( д4 ) и,+д = и, +1 (2а,+д, + 5а, - а,-д)д + О (д3)
Модифицированная форма предсказания-коррекции
Предсказание:
X,+а = X, + и,д +1 (5а, - а,-а)а2 + О (д3) 8
и,+д = и, + 2 (За, - а,-д)д + О ( д3) Коррекция:
X,+а = X +и,д +1 (а,+д, + За, )д 2 + О (д 3) 8
и,+д = и, +1 (За,+д, + 6а, - а,-д,)д + О ( д3)
Предсказание:
X,+а = X + и,д + 6 (4а, - а,-д)д 2 + О ( д4 )
и,+а = и, + 2 (За, - а,-а )А + О ( А3)
Коррекция:
X,+а = X +и,А + 6 (а,+а + 2 а, )А 2 + О (А 4 ) и,+А =и, +11 (2а,+а + 5а, - а,-д,)А+ О(А3)
Численное сравнение с другими алгоритмами
Для практической оценки характеристик алгоритма и сравнения с другими методами выполним численное моделирование гармонического осциллятора. Его движение в одномерном пространстве описывается следующим уравнением:
a = -ш x .
(28)
Аналитическое решение этого уравнения: x = Cj cos (Ш) + C2 sin (üt) .
Константы Cj и C2 зависят от начальных условий.
Вычисления выполнялись с использованием следующих методов:
- Верле (в скоростной форме);
- Бимана (в явной форме);
- предложенный алгоритм (в явной форме);
- Рунге-Кутты 4-го порядка (уравнение второго порядка (28) приводилось к системе уравнений первого порядка).
Скоростной алгоритм Верле - численный метод, предназначенный для решения уравнений движения. Алгоритм Бимана считается его вариацией, лучше сохраняющей энергию.
Записывается скоростной алгоритм Верле в следующем виде:
xt+At = xt+ч At+- a At
j
— i 2
j
= ч+2(a
+ at )At.
Алгоритмы Верле, Бимана и предложенный алгоритм являются специальными алгоритмами для интегрирования дифференциальных уравнений второго порядка. Основным их преимуществом относительно алгоритмов для уравнений первого порядка является скорость вычислений. В методе Рунге-Кутты вычисление силы (или ускорения) на каждом временном шаге выполняется четыре раза, в приведенных алгоритмах - только один раз. Поэтому когда основное время затрачивается на вычисление силы специальные алгоритмы оказываются в четыре раза быстрее. Если же для определения силы требуется меньше вычислений, чем для самого алгоритма, то специальные алгоритмы оказываются в десятки раз быстрее. Так для рассматриваемого примера специальные алгоритмы более чем в 60 раз быстрее метода Рунге-Кутты.
В рассмотренных методах по-разному проявляется погрешность вычисления координаты: в методе Рунге-Кутты ошибка накапливается и при больших значениях шага и расчетного интервала решение вырождается; в специальных методах погрешность вычисления координаты приводит к изменению частоты колебаний, но амплитуда и энергия остаются практически неизменными, позволяя выполнять расчеты с крупным шагом и на сколь угодно протяженных расчетных интервалах.
Важным параметром при моделировании движения частицы является сохранение энергии. Специальные методы имеют погрешность определения энергии, не выходящую за определённые границы на всем интервале вычислений. Метод Рунге-Кутты имеет накапливающуюся погрешность, из-за чего его на больших интервалах вычислений он хуже сохраняет энергию.
В табл. 2, 3 представлены погрешности определения энергии в зависимости от шага и интервала интегрирования соответственно. Из представленных данных видно, что предложенный алгоритм значительно лучше сохраняет энергию, по сравнению с другими специальными алгоритмами. Следует отметить, что метод Рунге-Кутты при значительном увеличении шага или интервала интегрирования уступает предложенному алгоритму.
Таблица 2. Максимальная погрешность определения энергии при t = 1000 с, %
dt, c Метод 0.001 0.01 0.02 0.08 0.1 0.3
Верле 2,5-10"5 2,5^10-3 0,01 0,15 0,25 2,3
Бимана 810"° 8,5-10-4 0,05 0,08 0,65
Предложенный 610-9 6,5-10-6 4,810-5 3,810-3 510-3 0,3
РК4 °-10-12 1.5-10-' 4,5-10-6 4,5^10-3 0,015 3,3
и
t, c Метод ^^^^^ 102 103 104 105
Верле 4-10"2 410-2 410-2 4-10"2
Бимана 1,410-2 1,410-2 1,4^10"2 1,410-2
Предложенный 4,5-10"4 4,5-10"4 4,5-10"4 4,5^10"4
РК4 1,4-10"5 1,410-4 1,4-10"3 1,410-2
Заключение
Представленные в работе результаты указывают на целесообразность использования предложенного алгоритма в качестве замены алгоритма Бимана, а при необходимости большого количества расчетов за ограниченное время и в качестве альтернативы другим распространенным методам.
EFFICIENT ALGORITHM FOR NUMERICAL INTEGRATION OF MOTION EQUATIONS OF LARGE PARTICLES IN MICROWAVE DEVICES
AO. RAK Abstract
The algorithm for the numerical integration of the motion equations of large particles is formulated. The algorithm is presented in three forms: explicit, prediction-correction and a modified, designed to simulate the motion of large relativistic particle in an electromagnetic field. The analysis of the order of accuracy of the proposed algorithm is presented. The results the comparison with the popular algorithms is shown. The efficiency of the algorithm is formulated.
Список литературы
1. Батура М.П., Кураев А.А., Попкова Т.Л. и др. Нерегулярные электродинамические структуры. Теория и методы расчета. Минск, 2011.
2. ГулдХ., ТобочникЯ. Компьютерное моделирование в физике. М., 1990.
