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

ISSN 0868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2007, том 17, № 1, c. 83-90

ИССЛЕДОВАНИЯ, ПРИБОРЫ, МОДЕЛИ-И МЕТОДЫ АНАЛИЗА

УДК 517.925.4: 62 - 408.64

© С. И. Шевченко

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

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

ВВЕДЕНИЕ

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

Решение уравнения Лапласа методом граничных интегральных уравнений означает переход от дифференциального уравнения Лапласа к эквивалентному ему интегральному уравнению [3]

ф(го) = J G(r, r0)dS .

(S)

(1)

где (7 (г)— функция плотности поверхностного заряда (ППЗ); С(г, г0)— ядро интегрального уравнения (ядро для задачи Дирихле для уравнения Лапласа); точка г0— точка, в которой ищется потенциал (точка наблюдения, ТН); интеграл в правой части берется по площади поверхности всех электродов (границы); коэффициент 1/4яе0, который должен быть перед интегралом в правой части, включен в функцию о.

Для решения уравнения (1) эффективным оказался метод коллокации, одним из возможных путей реализации которого является метод прямого интегрирования — непосредственное представление стоящего в правой части (1) интеграла в виде интегральной суммы (чаще всего по формуле Гаусса [4]). Для вклада в потенциал от одного отрезка границы можно выписать [4]

VP(го) = dp £АgЧг,) G(г,,Го),

(2)

где гi — точка интегрирования; — длина отрезка границы (электрода), по поверхности которого идет интегрирование (далее индекс р будем опускать); Л(п8)— п-й коэффициент квадратурной формулы Гаусса.

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

В результате выполнения первой стадии решения уравнения (1) вблизи оси получается резкий пик ППЗ, достигающий 1.5-2 % от теоретических значений. Ясно, что этот пик является следствием вычислительных, а не физических причин. Поэтому без внесения в алгоритм поправок, учитывающих влияние особенности оператора Лапласа (или соответствующего ему ядра интегрального уравнения G) вблизи оси, ни о какой приличной точности вычислений потенциала и компонент напряженности электрического поля не может идти речи.

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

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

Описание решения проблемы оси в литературе отсутствует. Некоторое касательство к проблеме оси (связь с проблемой оси) можно заметить в [6], где рассмотрена задача нахождения потенциала и его производных на оси вблизи некоторого касающегося оси электрода. При этом функция ППЗ а аппроксимируется полиномом четвертого порядка без всякого упоминания о получающемся пике ППЗ вблизи оси, что можно объяснить или тем, что проблема оси считается уже решенной, или, что авторы [6] считают необходимым учет особенности ядра только на втором этапе решения уравнения (1). Однако результаты, которые получались при вычислениях с помощью известного комплекса программ "Ро1880п2" [1], когда вблизи оси получается резкий пик ППЗ, достигающий 1.5 % от теоретических значений, нельзя признать удовлетворительными.

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

В представленной работе будут освещены особенности учета проблемы оси на первой стадии решения уравнения (1) применительно к алгоритму, основанному на методе коллокации и прямом интегрировании в методе граничных интегральных уравнений.

1. ОСНОВЫ АЛГОРИТМА УЧЕТА ЭФФЕКТА ОСИ

Ядро для задачи Дирихле для уравнения Лапласа в случае аксиальной геометрии имеет вид [1]

в(т, г0) = с • К(к):

(3)

где

4г

с =

у1(г + Го)2 + (г - 20)2 '

2 = 4ГГ0 .

— 2 2 ? (г + Го)2 + (г - го)2

К(к)— полный эллиптический интеграл первого рода, имеющий интерполяционное представление

[7]

К (к) = Х акЩк +|Х Ы

< .п

т1

т1(г, г) = 1 - к2 =

(г - Го)2 + (г - го)2 (г + Го)2 + (г - го)2

После несложных преобразований ядро (3) может быть представлено [2] в виде

б(т, Го) =

= Л (г,г) + /2(г,г) 1п[(г - Го)2 + (г - ^)2],

(4)

где

/Л г, Г) = -

4г

7(г + Го)2 + (г - го)2 х{& (г,г) + g2(г,г) 1п[(г + Го)2 + (г - го)2]};

/2( г, г) = -

4г

у1(г + Го)2 + (г - го)2

гг) = Е акт1;

к=о

4

g2( г, г) = Е ькт1к;

значения коэффициентов ак и Ък можно найти в справочнике [7].

Входящие в (4) функции /, /2 являются гладкими везде в области г Ф о. Ясно, что сингулярность дифференциального уравнения при переходе к интегральному уравнению трансформируется в некоторую сингулярность ядра G(r, го) этого интегрального уравнения.

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

с = -

4г Г + го

. = 4_х + 1о

Г - Г

х + 21о

У С

г + го

х + 2/о

где х = I - 1о, 1о — длина отрезка от начала электрода до ТН, I — длина отрезка от начала электрода до точки интегрирования (ТИ).

Графики последних выражений для функций с, т1 представлены на рис. 1, 2. Эти функции имеют особенности: функция с(х)— полюс первого порядка в точке х = -/о и функция т (х) — полюс второго порядка в точке х = -21о.

т1 =

C(x)

°'0. 0.1 0.2 0.3 0 .А 0.5

Рис. 1. Графики зависимости функции c(x) для точек наблюдения, расположенных в первых трех (считая от начала отрезка границы, находящегося на оси) точках коллокации. ТН расположена:- в первой точке коллокации, ---во второй, ............ в третьей точке

коллокации

mJ(x)

Рис. 2. Графики зависимости функции mj ( x) для точек наблюдения, расположенных в первых трех (считая от начала отрезка границы, лежащего на оси) точках коллокации. Типы линий соответствуют таковым на рис. j

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

Это особенно касается функции mj(x), для которой (для ТН, расположенной в первой точке коллокации) "отросток" (левее первой точки коллокации) вообще невозможно аппроксимировать некоторым полиномом. Поэтому этот "отросток" не может быть учтен при численном интегрировании по формуле Гаусса.

Чтобы преодолеть эту трудность, в данной работе используется предложенный в [8] и широко в настоящее время применяемый метод понижения степени особенности интегрируемой функции. Это метод был обобщен на случай интерполяции полиномом второго порядка в [9] и полиномами более высокого порядка в [10].

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

{f ( x)L = Х f ( x )i q x.

1=0 j=0

Эту формулу можно получить, записав интерполяционную формулу Лагранжа [11] и представив коэффициенты при значениях функции в виде полинома. В [j0] описан менее трудоемкий метод получения этой интерполяции.

Для реализации метода коллокации запишем уравнение (j) в каждой точке коллокации ( r0 = rn; n = J +M ; М — число точек коллокации). Считаем, что находящаяся под знаком интеграла функция о является гладкой вдоль длины отрезка интегрирования.

Поэтому, если l0 мало по сравнению с длиной отрезка интегрирования (т. е. l0 является первой или второй точкой коллокации при отсчете номеров от начала отрезка границы, опирающегося на ось), то на участке [0^ l0] (или [0^210]) функции c, mJ изменяются довольно резко (быстро). Поэтому численное интегрирование этих функций (и функций от этих функций), проводимое вдоль длины отрезка интегрирования, может не обеспечить удовлетворительную точность.

Из рис. J, 2 видно, что для ТН, совпадающей с первой точкой коллокации (считая от начала отрезка границы, опирающегося на ось), функции c(x) и mj(x) имеют вблизи этого конца границы довольно резкий характер. Более того, если смотреть на точки коллокации (отмеченные на рисунках жирными черными точками), то видно, что

2. ФОРМИРОВАНИЕ МАТРИЦЫ

При переходе к случаю непрямолинейного отрезка представим функции с, т1 в виде произведения "прямолинейной" части на некие гладкие поправочные функции, которые допускают интерполирование полиномами и которые близки по своим значениям к единице:

c = Pc ( x);

x + 2l0

(

mJ =

l + 2l0

Pm (x).

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

4 .,о

Фр(Го) = Хак I х)с(х)т1к(х)<3х-

к=о -1о

л-ю

С 1 >

т

dx.

(5)

+ ЕЬк / х)с(х)т1 (х)1п

к=о Ч I ,

Если входящий в выражение (5) логарифм представить в виде

1п

С 1 >

т

1

= 21п( х + 2/о) - 21п| х| - 1п Рт (х), то выражение (5) примет вид

4 л-о

Фр (Го) = Х I х)с(х)т1 (х) [ - Ьк 1пРт(х)]] +

к=о -¡о

4 лЧ

+ 2^Ьк I сг(х)с(х)тк (х)1п(х + 2/о^ -

к=о Ч

4 л - ¡о

2^Ьк I сг(х)с(х)т1к (х)1п| х| dx. (6)

к=о

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

л-¡о

Ф

.(го) = 4Е Ц-( х) хх++^ Рс (х)

к=о -

'о

х + 21о

С ^

х + 21о

х

хРтк (х) [ак - Ьк 1пРт(х)] [dx +

л-¡о

8Е Ьк 1>( х) Рс (х)

к=о -/о

х + 2¡0

С ^

х + 24

хРтк (х) 1п(х + 2/о) № -

х +1 С х ^

8Х Ьк I н х) Р< (х)

к=о -¡

'о

х

х + 24

хРк (х)1п| х| №.

(7)

Обратим внимание на то, что в первом члене присутствует только одна сингулярность (полюс

1/(х + 2^)2к 1), находящаяся вне пределов интегрирования и оказывающая влияние тем, что вблизи точки х = о (точнее говоря, на промежутке х = [-¡о, ¡о] при условии, что ¡о << 1) подынтегральная функция не может считаться достаточно

гладкой для проведения численного интегрирования по формуле Гаусса.

Во втором члене в точке х = -2^ расположен полюс и логарифмическая сингулярность.

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

В третьем члене в точке х = -2^ расположен полюс и логарифмическая сингулярность в точке х = о . Т. е. имеется совсем другая ситуация. Для первых двух точек коллокации, если отсчитывать от начала рассматриваемого отрезка границы, шаблоны для обеих сингулярностей совпадают. А для третьей точки коллокации шаблоны отличаются (сдвинуты на единицу). Численные эксперименты показали, что рассмотрение особенности оси оказывает заметное влияние только для первых трех точек коллокации (если для данного отрезка границы рассматривать не более одиннадцати точек коллокации). Чтобы обойти сложность с двумя разными шаблонами, были применены шаблоны по четырем точкам (эти шаблоны являются несимметричными относительно точки кол-локации), причем шаблон интерполяции, применяемый для понижения степени логарифмической сингулярности, сдвинут на один номер в сторону начала отрезка границы, в результате чего шаблоны совпали.

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

Процедура для 1-го члена из (7)

Первый член имеет вид

Ф =

4Е /М х)

к=о -

х + ¡о х + 24

С

Р (х)

2к

х + 24

хРкт (х) [ак - Ьк 1п Рт (х)]

х

dx.

Введем обозначение

Ро(к)(х) = Рс (х) • Рк (х).

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

{(х)^)(х)[ -Ък 1пРт(х)] =

= Ъ\°(Хт(i))Ро(к)(Хт(i)) X

к=0 -/

-{( Х) Р0( к)( Х) [ - Ък 1п Рт (х)] )■

(х +10) х2к

2к +1

(х + 2/0)

4^ / {(х)Р0(к)(х)[ -Ък 1пРт(х)]}

dх +

к=0 -/

X

(х + /0) х ■dх.

Х[к - Ък 1ПРт (хп)]

(хп +Ц х2к (хп + 2/0) 2к +1

хт(,)) Р0(к)(хт(,)) X

х[ ак - Ък 1п Рт (хт(,)]Х С" ((1к - < к)[,

где хп = ё • х(п8); х(п8) — п-й узел (абсцисса) квадратурной формулы Гаусса;

9 (1) = V Л( 8) хп (хп + /0) х«к .

= Ъ Лп х (хп + 2/0)2к+1;

п=0; пФт(})

ё-/0

Х[ак - Ък 1П Рт (хт(1>) ]Ъ Сi, х Г. (8)

3=0

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

Используя интерполяцию (8), представим первый член в виде

€ = 4Ъ I Их)Р0к)(х)[ -Ък 1пРт(х)]-

Б(1) = Г" хз (х + /0)х2к Бк ^ х (х + 2/0)2к+1

dх.

Интеграл Б(к можно полностью привести к сумме

аналитических (табличных) интегралов (см. Приложение).

Из выражения для €1 можно выделить вклад

первого члена в квазидиагональные элементы матрицы (выражения, которые стоят при ст(х), п = т(1)):

с(1) =

= -4ЪЪРок)(хп)[-Ък 1пРт(хп)]С® (^к-Б(,)

к=0 3=0

и во все остальные элементы матрицы, относящиеся к рассматриваемому отрезку интегрирования (п Ф m(i)):

С% = 4ёлп8) £ Р0(к)(хп) [ - Ък 1п Рт (хп ^^ .

к =0 (хп + 2/0)

В приведенных формулах остается еще определить полином Р0(к)(х) . Для этого рассмотрим ра-

венство

[к - Ък 1п Рт (хп )]С( хп )т\(хп ) =

(х + 2/(1У

Как обычно (см. [1]), в первом члене в правой части последнего выражения (точнее говоря, к стоящему там интегралу) применяем численное интегрирование методом Гаусса, а во втором члене проводим аналитическое интегрирование:

4 Г

€ = 4Ъё Ъ] лп8) а(хп)Р0(к)(хп) х

к =0 п=0; [ пФт(г)

= 4

хп + 10

х„ + 2/0

V

х„ + 2/0

Р0(к )(хп).

Отсюда видно, что искомый полином Р0к (х) можно представить в точках шаблона интерполяции (за исключением точки хп = 0) в виде

Р0к)( хп ) =

с(хп )т1к (хп ) [к - Ък 1п Рт (х)]

Г

х„ + /0

п_0_

х,. + 2/„

V

ч 2к

х„ + 2/„

В правой части этого выражения в точке хп = 0 существует неопределенность типа 02к/02к . Из интуитивных соображений ясно, что в точке хп = 0 полином Р0к (х) имеет значение Р0к (хп = 0) = [ - Ък 1п Рт (хп )] = ак, т. к. Рт (х = 0) = 1 и 1п Рт (х = 0) = 0.

3 =0

4

Процедура для 2-го члена из (7)

Второй член имеет вид

л-¡о

=

8ХЬк I Йх)х^Р(х)

к=о -¡

х + 24

С 4

х + 24

2к

Процедура для 3-го члена из (7)

Третий член имеет вид

л-¡о

=

8Х Ьк I Йх) хх7+2гРс (х)

С

2к

к=о -¡

х + 24

х + 24

х

х Рк (х)1п (х + 2^ ) dx.

Используем интерполяцию по четырем точкам {(х) Ро( к)(х)} =

= Х*( хт (,.)) Ро(к)(х„(. ))Х С

г=о 1=0

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

Фф2) = 8ХЬкл ] \ ^Й(хи)Ро(к)(хи)х

(хп + ¡о) х2

к=о п=о;

п Фт (г)

'(хп + 2^2 к+1

1п (хп + 2¡о) +

+ 8] Ьк хт(1)) Ро( к)(хт(1))

к=о 1=0 [

х]с® ((2) - ^2) ),

х

1 =о

где ^к = )хп ^^ 1п(хп + 2¡о);

п=о; ((хп + 2/о)

п Фт (г)

1 = Т ^ ((хХ ^ (хп + 2¡о ) .

Из последнего выражения можно выделить вклад второго члена в квазидиагональные (п = т(1)):

с^=8] Ьк^)(хп): с® (л • « 2)- в 2))

к=о 1=о

и во все остальные члены матрицы (п Ф т(1)), относящиеся к отрезку интегрирования:

с1я

4 (хп + ¡о) хп2к

X ЬкРо(к)(хп)

к=о

(хп + 2«

2 к+1

1п (хп + 2¡0 ).

х Ркт (х)1п|х| | dx.

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

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

Фф3) = 8ЕЬка X ^Й(хп)Ро(к)(хп);

к=0 п=0;

пФт(г)

, (хп + ¡о)хп2к

Ч( хп + 2¡о)2 к+1

1п х\ К

+ 8] Ьк хт(1)) Ро( к)(хт(1)):

к=0 1=0 [

хХ(3) - ^3) )),

1=0

где

«(3) = У )х1 (хп + ¡о)хп 1п |х I-

1 = Х А хп (хп + 2^+1^1^ ;

п=0; пФт(г)

л-¡о

(х + ¡0) х2к

в(з) = I х1 ^ +11п|хdx.

1 ,к ^ (х + 2^)2к+1 1 1

Из выражения для ф(р^) можно выделить вклад

третьего члена в квазидиагональные и во все остальные члены матрицы, относящиеся к отрезку интегрирования:

п = т(1):

с® = 8]ЬкР(к)(хп)ХС® (л• 3к -В<3),

к=0

1=0

п Ф т(1):

с® = 8^)XЬкР0(к)(хп) (хп + ¡о)х2" 1п|хп|.

г ,п п к о уп'/|Т 1 \ 2 к+1 I п|

к=0 (хп + 2/0)

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

Матричное уравнение решалось методом Гаусса с выделением ведущего элемента.

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

3. ТЕСТОВЫЙ ПРИМЕР И ОБСУЖДЕНИЕ

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

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

Основным результатом реализации приведенного алгоритма явилось заметное увеличение точности вычисления плотности поверхностного заряда. В среднем для всех тестовых примеров были получены примерно одинаковые результаты по точности функции ППЗ на оси системы — не хуже 0.1 %.

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

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

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

ЗАКЛЮЧЕНИЕ

Таким образом, в данной работе представлен алгоритм решения проблемы оси (учета особенности оператора Лапласа вблизи оси) применительно

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

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

Тестовые вычисления показали, что описанный алгоритм является работоспособным и позволяет понизить получаемую в результате решения уравнения Лапласа величину пика функции ППЗ на оси до 0.1 % от теоретического значения.

Приведенный в данной работе алгоритм был реализован в виде программ, включенных в очередную версию пакета прикладных программ "Shift" [12].

Приложение. Нахождение аналитических интегралов

В первом члене выражения (7) возникает ин-

d(1) f° j (Х + l0)x2k j

теграл Byk = J x --^ )2k+1 dx, который с помощью замены переменных y = x + 2l° и некоторого преобразования приводим к совокупности простых табличных интегралов вида

d

const-J ym dy,

0

где m = -9 ^ 3 .

Во втором члене выражения (7) возникает ин-

d -l° (x +1 ) x2k

теграл B(2k = J x 2° +1 ln (x + 2l°)dx, кото--l° (x + 2l0) рый с помощью замены переменных y = x + 2l° и некоторого преобразования приводим к совокупности табличных интегралов вида

d

const-J ym ln(y) dy .

0

В третьем члене выражения (7) возникает интеграл Bfk = J xJ ^^21 )2k+1 ln| x|dx , который после замены переменных y = x + 2l° и некоторого преобразования приводим к совокупности интегралов вида

d

const-J ym ln|y — 2Z0| dy .

0

Если учесть, что l0 < d, то точка y = 2l0 принадлежит отрезку интегрирования. Последний интеграл тоже является табличным и может быть выписан в конечном виде, за исключением случая m = —1, для которого имеется решение в виде ряда [13].

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

1. Иванов В.Я. Методы автоматизированного проектирования приборов электроники. Новосибирск: Изд-во СО АН СССР, 1986. 193 с.

2. Шевченко С.И. Алгоритм получения предельной точности в электростатических расчетах элементов электронно- и ионно-оптических приборов, имеющих плоскую симметрию // Научное приборостроение. 1997. Т. 7, № 1-2. С. 45-53.

3. Соболев С.Л. Уравнения математической физики. М., 1966. 400 с.

4. Крылов В.И., Шульгина Л.Т. Справочная книга по численному интегрированию. М.: Наука, 1976. 370 с.

5. Scherzer O. // Z. Phys. 1936. V. 101. P. 593-603.

6. Монастырский М.А., Иванов В.Я., Куликов Ю.В. Устранение особенностей в методе интегральных уравнений при расчете осевых распределений потенциала и его производных вблизи границы // Новые методы расчета элек-

тронно-оптических систем. М.: Наука, 1983. С.187-191.

7. Абрамовиц М., Стиган И. Справочник по специальным функциям. М.: Наука, 1979. 830 с.

8. Канторович Л. В. О приближенном вычислении некоторых типов определенных интегралов и других применениях метода выделения особенностей // Математический сборник. 1934. Т. 41, № 2. C. 235-244.

9. Фрейкман Б.Г. Вычисление электростатического поля вблизи заряженной поверхности // ЖТФ. 1979. № 11. С. 2464-2472.

10. Шевченко С.И. Алгоритм получения высокой точности в расчетах аксиально-симметричных электростатических полей // Научное приборостроение. 2002. Т. 12, № 1. С. 40-45.

11. Корн Г., Корн Т. Справочник по математике для научых работников и инженеров. М.: Наука, 1970. 720 с.

12. Шевченко С.И. Пакет прикладных программ "Shift" для решения уравнений Лапласа и Пуассона // Материалы Всес. семинара "Методы расчета электронно-оптических систем". Алма-Ата, 1992. С. 11.

13. Двайт Г.Б. Таблицы интегралов и другие математические формулы. М.: Наука, 1973. 228 с.

Институт аналитического приборостроения РАН,

Санкт-Петербург

Материал поступил в редакцию 6.09.2006.

ON SOME PECULIARITIES ARISING IN CALCULATING NEAR-AXIS AXIAL ELECTROSTATIC FIELDS. I. FORTH INTEGRATION METHOD

S. I. Shevchenko

Institute for Analytical Instrumentation RAS, Saint-Petersburg

The paper presents a method of solving the axis problem related to the boundary integral equation method where the integral equation is transformed into a matrix equation by collocation and forth integration. In forming the matrix, especial attention has been given to the near-axis peculiarities of the integral equation kernel.