Метод сплайновой интерполяции для решения уравнения векторного потенциала в задаче транспортировки электронных пучков Текст научной статьи по специальности «Математика»

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

1. Вылегжанин О.Н., Шкатова Г.И. Учет ограничений равенств при решении оптимизационных задач с линейными ограничениями // Известия Томского политехнического университета. - 2008. - Т. 312. - № 5. - С. 76-78.

2. Васильев Ф.П., Иваницкий А.Ю. Линейное программирование. - М.: Факториал, 1998. - 323 с.

3. Вылегжанин О.Н., Шкатова Г.И. Сравнительная оценка двух методов выбора наилучших линейных регрессоров // Применение математических методов и ЭВМ в медико-биологиче-

ских исследованиях / Межвузовский научно-техн. сб. под ред. В.А. Кочегурова. - Томск: Изд-во Томск. политехн. ин-та, 1988. - С. 18-22.

4. Гантмахер Ф.Р. Теория матриц. - 5-е изд. - М.: Физматлит, 2004. - 559 с.

5. Cline R.E. Representation for the Generalised Inverse of a Partitioned matrix // J. Soc. Industr. and Appl. Mathem. - 1964. - V. 12. -№ 3. - P. 588-600.

Поступила 15.04.2009 г.

УДК 681.31:533.95

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

В.П. Григорьев, И.Л. Звигинцев, А.В. Козловских

Томский политехнический университет E-mail: [email protected]

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

Ключевые слова:

Сплайн-интерполяция, сплайновая аппроксимация, уравнения векторного потенциала, математическое моделирование.

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

В данной работе предлагается новый метод для решения электродинамических задач. Апробация и тестирование этого метода проводится на примере решения проблемы транспортировки пучков к мишени в газе низкого давления.

Обычно подобного рода задачи решаются в пакете МаИаЬ, если другие пакеты с готовыми встроенными алгоритмами (например, Сош8о1) не могут справиться с поставленной задачей из-за расходимости вычислений. МаНаЬ имеет достаточно большой набор солверов, основанных на различных численных методах для решения систем обыкновенных дифференциальных уравнений. Для решения жестких систем с высокой точностью подходит солвер, основанный на многошаговом методе Гира, который допускает изменение порядка точности [3]. При вычислении производных в МаНаЬ ис-

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

Исходная задача математически представлена системой нелинейных уравнений в частных производных:

dAz

и

ef

dnt ~dt

dt

1 A

r dr

d2A _ 4n

dr2

_— evbnb

(1)

_ G, ■ V. ■ И ■ П. + П ■ K ■ П

где А(г,1) - векторный потенциал поля, и;(г,0 -плотность ионов, ие(г,/) - плотность плазмы, г0 -классический радиус электрона, цДг,/) - частота соударений электронов и ионов плазмы, КС*,0 -коэффициент ионизации, - плотность газа, нь(г,1) - плотность пучка, уъ - скорость пучка, е -заряд электрона, ст; - сечение ионизации.

Эти нелинейные уравнения связаны через сложную зависимость приведенных выше параметров от Аг(г,/) и и;(г,0 [1].

Краевые условия:

дА

дг

= 0, 4 (г = Яс) = 0, п. (/„) = 0.

Система двух уравнений в частных производных сводится к одной системе нелинейных обыкновенных дифференциальных уравнений, размерность которой 2N, где N — число разбиений пространственной координаты (представляем пространство дискретным). Для каждой точки пространства получается дифференциальное уравнение первого порядка по времени. Решение такой системы уравнений представляет значительные вычислительные трудности, связанные как с жесткостью системы, так и с большой размерностью задачи. При этих условиях особое значение имеет вычислительная эффективность алгоритма. Содержание предлагаемого метода заключается в разделении задачи на две подзадачи: вычисление вектора правых частей и его интегрирование. Векторы правых частей вычисляются отдельно для векторного потенциала поля и для плотности ионов плазмы. Для этого был разработан метод аппроксимации производной по пространственной координате, основанный на сплайн-функциях. После вычисления эти вектора объединяются в один вектор размерностью 2N, который интегрируется с помощью солвера. Таким образом решаются сразу все уравнения системы. Ускорение вычислений обеспечивается за счет применения постоянной матрицы при вычислении значений производных от векторного потенциала в узлах дискретной сетки. Для нахождения этой матрицы рассмотрим первое уравнение системы (1). Определим первую и вторую производные по пространственной координате от искомой функции в каждом узле.

Будем использовать кубический сплайн для интерполяции зависимости /(х), заданной в узлах х0,хь...,х« значениями /0/1,.../. Кубический сплайн на интервале [х;,х;+1] записывается в виде [4]:

Я, (х) = а, + Ь1 (х - х,) + с, (х - % )2 + 4 (х - х , )3, (2)

где а, Ь, с1 и 4 - коэффициенты сплайна, определяемые из дополнительных условий; /=0,1,...,«— 1 -номер сплайна.

Коэффициенты сплайна определяются из условия непрерывности функции-интерполятора и двух первых ее производных в узлах исходной дискретной сетки:

(х,-о) = Я(х+оХ ЯXх-.о) = (х+ о)>

5 '(х,_о) = (х,.+о), (3)

а также равенства значений сплайна £(х) и интерполяции функции /(х) в узлах:

5,.(х) = /, Я(х,.+,) = /+г (4)

Для доопределения системы линейных алгебраических уравнений используются краевые условия (значения первой и второй производных сплайна в точках х0 и х„). Например, естественные краевые условия:

50(хо) = 0, ЯП(хп) = 0. (5)

Получим алгоритм определения коэффициентов кубических сплайнов. Условия (4) в узлах х1 и х;+1 после подстановки -го сплайна принимают вид:

а, = /,, а, + ЬА + сА2+4А3 = /, +

(6)

где к==хн1—х, 0< /<«—1.

Продифференцируем дважды сплайн (2) по переменной х:

Я'(х) = Ь + 2с, (х- х) + 34 (х- х )2,

Я'( х) = 2с, + 64 (х - х). (7)

Из условий непрерывности производных (3) при переходе в точке х, от ( /—1)-го к /-му сплайну с учетом выражений (7) получим следующие соотношения:

Ь,-1 + 2с,-Л-1 + 34-Л-1 = Ь,, С-1 + 34,- А-1= С • (8)

Соотношения (6) и (8) представляют собой полную систему алгебраических уравнений относительно коэффициентов сплайнов а, Ь, С и с/. Но, прежде чем решать эту систему, выгодно преобразовать ее так, чтобы неизвестными была только одна группа коэффициентов с. После несложных преобразований получим уравнение, в которое входят только неизвестные коэффициенты с,, пропорциональные второй производной функции интерполятора в узлах сетки:

к,-1с ,-1 + Щ-1 + к,)с, + кс+1 = = 3-[(/+ - /)/к - (/ - £-1)/к Л, (9)

где 1< /<«—1.

Из граничных условий (5) получим, что:

с0 = 0, сп = 0.

(10)

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

В каждое из уравнений типа (9), входит только три неизвестных с последовательными значениями индексов: с—1, с;, с/+1. Следовательно, матрица системы линейных алгебраических уравнений относительно с является трехдиагональной, т. е. имеет отличные от 0 элементы только на главной и двух примыкающих к ней диагоналях. Для решения систем с трехдиаго-нальной матрицей наиболее эффективно применять так называемый метод прогонки, являющийся частным случаем метода исключения Гаусса [5].

Будем считать, что шаг по х постоянен, т. е.

+1 =к для всех /. Тогда выражение (9) примет вид:

к2 (с,-1 + 4с, + с,+1) = 3 - (/+1 - 2 / + /-1). (11)

Воспользовавшись выражениями (7), получим значения производных в узлах:

(х,) = Ь, ЯДх) = 2с. (12)

Первую производную сплайна в -ом узле можно вычислить, воспользовавшись выражениями (8), по рекуррентной формуле (для всех 1< /<«):

Ь1 = (с, + с,-1) - к+Ь,-1, (13)

или, используя совместно выражения (6) и (8):

г=0

b = f+1 - f + (C+1 + 2c, )h ' h 3

b, = f - f-i + (2c, + С-i)h

(14)

2C = 6 C-lFF, h2

(15)

B =

-1 1 0 .. 0 0 0 -1 1 .. 0 0

1 0 -1 1

H =

1 1 0 0 1 1

0 0 0 0

0 0 0 0

dt

= 0. Следовательно:

52 A

dr2

cA Rc dr

= 0

(17)

Cn-1 + Cn (2 + 6n) = h^(an-1 - an ).

Окончательный вид матриц для вычисления второй производной сплайнами будет иметь вид:

(

к 3

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

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

C =

2 1 1 4 0 1

0 0 0 0 0 0

0 0 0 0 0 0

F =

а выражение (13) так:

В = кВ ~1ИС. (16)

Значение второй производной вычисляется через значения функции во всех точках сетки, а не в трех соседних, как для метода конечных разностей.

Теперь найдем дифференцирующие матрицы для рассмотренного алгоритма аппроксимации сплайнами. Каждая строка матриц соответствует узлу по пространственной координате. Приведенные выше выражения описывают внутренние узлы. Для построения матриц необходимо также знать ее строки, соответствующие граничным узлам. Определим матрицы выражения (16), отвечающие нашей задаче. Учитывая, что на левой границе первая производная принимает нулевое значение, получим:

(1 о о .. о о^ (о о о .. о о^

-1 1 1 -2 0 1

0 0 0 0

0 1 -2

0 0

0 0

0 0

-2 1

где а?=2+6и.

Для построения матрицы, аппроксимирующей дифференциальный оператор ЦЛ) необходимо знать матрицу, соответствующую 1/г - матрицу Я1: Ь¥ =2 С +Яг1Б. На левой границе матрицы имеет место неопределенность, которая разрешается по правилу Лопиталя. Чтобы существовала обратная матрица Яг1, сделаем не нулевой первую строку у матрицы Я. Это ни на что__не повлияет, т. к. в результирующем векторе ЯгБ придется заменить нулевой элемент на нулевой элемент вектора 2 С:

(1 о о .. о о ^

R = h •

0 1 0 0 0 2

0 0 0.. n-10

0 0 0

0

Для определения матриц выражения (15) возьмем выражение (14) и, учитывая равенство нулю первой производной на левой границе, получим:

3

С1 + 2со = Тг(а1 - «о). к

Так как на правой границе векторный потен -циал принимает постоянное нулевое значение, то дА

для всего временного интервала. Если записать это в дискретном виде, то -2 Д,-с„=Ь„. Взяв выражение (14) и учитывая предыдущую формулу, запишим:

Полученная матрица, аппроксимирующая дифференциальный оператор L(A), используется в ode-функции при интегрировании полученной системы уравнений. Эту матрицу достаточно вычислить один раз перед вызовом солвера Matlab (на каждом временном шаге она будет постоянной), т. к. используется сетка с постоянным шагом. Уравнения векторного потенциала поля и плотности ионов плазмы решаются одновременно на каждом временном шаге. Сначала внутри ocfe-функции вычисляется вектор правых частей (размерность N) векторного потенциала поля. Сразу же, используя полученный вектор, вычисляется вектор правых частей (размерность N для плотности ионов плазмы. Затем происходит их объединение в один вектор размерности 2N - вектор правых частей всей системы, который используется для интегрирования в солвере.

Предложенный алгоритм для подобного рода задач работает на треть быстрее встроенных в Mat-lab функций для вычисления сплайнов. При использовании стандартных функций Matlab на каждом временном шаге решается полная задача ин-

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

— lb=12 kA lb- 15 kA lb= 17 kA

* t j *

J 'f.' •

i

О 0.1 0.2 0 3 0.4 0.5 0.6 0.7 0.8 0.9 1

5 х Ю-8

Рис. 1. Зависимость векторного потенциала Az от времени t в центре пучка при различных значениях тока пучка Ь и давлении газа ~0,3 Па

t- 2.5-10"7 s t= 5.0-10"7 s t= 7.5 10"7 s

____

\ \

\ \

V 's

4 \. V N

\ \

^s. "V

i

i "

012345В789

г ст

Рис. 2. Зависимость векторного потенциала Azот радиуса г в различные моменты времени при токе пучка Ь=15 и давлении газа ~0,3 Па

В случае решения системы уравнений (1), при количестве интервалов разбиения пространственной координаты «=36 (система уравнений 72 порядка) эффективность алгоритма составила 33 %; при «=72 достигла 36 %, а при «=144 - 38,4 % при одинаковой погрешности. Под эффективностью понимается, насколько разработанный алгоритм вычисления производных в узлах сетки ускоряет решение задачи по сравнению с использованием стандартных функций для сплайнов пакета МаНаЬ. Результаты решения задачи при использовании предложенного метода и числе узлов «=72 представлены на рис. 1-3.

Рис. 3. Зависимость полного тока П и плотности электронов плазмы Ne в центре пучка от времени t при давлении газа ~0,3 Па

Выводы

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

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

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

1. Григорьев В.П., Поташев А.Г., Шулаев Н.С. Численное моделирование формирования плазменного канала при инжекции мощного электронного пучка в нейтральный газ // Физика плазмы - 1979. - Т. 5. - № 2. - С. 376-381.

2. Grigoriev V.P., Koval T.V., Potashev A.G. Plasma channel forming and an electron beam transporting in the air with the pressure gradient in the external magnetic field // Proc. 15th Intern. Conf. on High-Power Particle Beams. - Saint-Petersburg, 2005. -P. 151-154.

3. Ануфриев И.Е., Смирнов А.Б., Смирнова Е.Н. Matlab 7 в подлиннике (Наиболее полное руководство). - СПб.: БХВ-Петер-бург, 2005. - 1104 с.

4. Алберг Дж., Нильсон Э., Уолш Дж. Теория сплайнов и ее приложения. - М.: Мир, 1972. - 316 с.

5. Мудров А.Е. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. - Томск: МП «Раско», 1991. - 272 с.

Поступила 03.03.2009 г.