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

УДК 517.9, 519.6

Д. В. Валовик, Ю. Г. Смирнов

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

Аннотация. Предложен метод коллокации как альтернатива методу Галерки-на, для решения псевдодифференциального уравнения электрического поля. Ключевые слова: прямое и обратное преобразование Фурье, псевдодифферен-циальный оператор, псевдодифференциальное уравнение, метод коллокаций.

Abstract. Collocation method (alternative to Galerkin method) for solving pseudodifferential equation of electric field is suggested.

Keywords: direct and inverse Fourier transform, pseudodifferential operator, pseudodifferential equation, collocation method.

1. Постановка задачи

Рассмотрим следующую задачу дифракции. Пусть в свободном пространстве R расположено объемное тело (область) Q с границей dQ класса C ™, хаpактеpизующееся постоянной магнитной проницаемостью Цд и 3 X 3-матрицей-функцией (тензором) диэлектрической проницаемости s (x ). Компонентами тензора s(x) являются Zÿ (x) - бесконечно гладкие функции

в Q , т.е. Zÿ (x)е C™(Q), причем Zÿ (x) = eÿ (х)/eo, где £q - диэлектрическая проницаемость свободного пространства.

Из условия конечности энергии необходимо [1], чтобы Eе L2 (Q) =

= L (Q )X L2 (Q )X L2 (Q ).

Требуется определить электромагнитное поле E, Hе L2 (Q), возбуж-

Е0 ТТ0 ~ —ifôt

, H с временной зависимостью вида e . Будем искать электромагнитное поле E, H, удовлетворяющее уравнениям Максвелла, условиям непрерывности касательных компонент поля при переходе через границу тела и условиям излучения на бесконечности [1].

Задача отыскания E, H сводится к решению интегродифференциаль-ного уравнения [1]

0 (x )J (x ) = E0 (x ) + ko J G E (x, y )J (y )dy +

Q

+ graddiv JGe (x,y)J(y)dy, xе Q, (1)

Q

где J(y) = ( (y), J2 (y) J3 (y))T.

Считаем, что 0 (x ):=( (y ) — I ) 1 существует при всех x е Q и J (y ):=

:=((y)—1 )E(y) тогда E(y) =0(y)J(y); E(y) = ((y)e2(y)E3(y))T -

(комплекснозначный) вектор электрического поля и у = (, у2, Уз) - точка в пространстве Л3; I - единичная 3х3-матрица; 6Е (х,у) =

о ОЕ о

0 0 о\

V

1 еК\х-у| _

тензорная функция Грина, где О™ (х, у ) =---------:-----¡- + gm (х, у), х, у е 2.

4л х - у

gm е С“(х2) - гладкая функция, ( = 1,2,3); ко - волновое число свободного пространства.

Уравнение (1) как псевдодифференциальное запишется в виде [2]1

где

AJ = E0,

AJ = —1T JJei(x-y)^(0(x)- dt(П))(y)dydI, (2я)

(2)

(3)

(r- 2

и dt (£)=-1 (П)

П2 - ^2 П1П3

^1^2 ^2 - k0 ^2^3

^1^3 ^2^3 ^2 - k0

' (П)=• pi = і-

"=1

in

В работе [2] относительно оператора A доказаны следующие теоремы. Теорема 1.1. Если выполняются условия:

1) матрица 0(х):=((y)-1) 1 существует при всех xє Q ;

3 _

2) Д(х )= І cos щ cos ay £j (x )^ 0, x є Q , и cos2 «1 + cos2 «2 +

i, J =1

+ cos2 «3 = 1, то оператор A: #^omp (Q)^ ^loc (Q), определенный по формуле (3), является эллиптическим псевдодифференциальным оператором

(пдо).

Теорема 1.2. Пусть выполнены условия теоремы 1.1. Тогда, если дополнительно выполнено одно из двух условий:

1) Res(x)v • v > (C2 + 1)|v|2, при xє Q и C2 > 0 ;

2 _________

2) Ims(x)v • v > C3 v , при xє Q ,

то оператор A: L2 (Q)^ L2 (Q), определенный формулой (3), является

фредгольмовым с нулевым индексом.

Нас будет интересовать часть главного символа, которая определяется выражением

1 Далее во всех интегралах, где пределы интегрирования не указаны явно, считаем, что интегрирование ведется по всему пространству.

1

Й?

5? 5l5? 5і5з ^ -.2

5і5? 52 5?5з 2

(4)

2. Метод коллокации

Сначала кратко опишем общую схему метода коллокации, а затем применим ее к уравнению (1) (или (2)).

Для уравнения Аф = f (, f е X) в гильбертовом пространстве X

рассмотрим метод коллокации, который формулируется следующим образом. Приближенное решение фи е Xn определяется из уравнения PnApn = Pnf . Здесь pn е Xn (Xn есть n -мерное подпространство пространства X), Pn : X — Xn - оператор проектирования на конечномерное подпространство, который определяется ниже.

Разобьем область Q на элементарные подобласти Q- с кусочногладкими границами dQ- так, чтобы выполнялись условия Q- n Qj = 0 при

i Ф j и Q = U Qi . Выберем в каждой подобласти Q- точку (узел) коллокации i

i [1,х е Qi

x . Рассмотрим базисные функции v- = < . Пусть подпространства Xn

[0,x g Qi

являются линейными оболочками базисных функций: Xn = span{v/,...,vn} .

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

Vxе X lim inf ||x - x|| = 0.

n——^ xfeXn

Проектор Pn : X — Xn определим так: (Pnp)(x) = p(x- ), xе Qt. Заметим, что при таком определении проектора не определены значения функций (Pnp)(x) при xе dQi, но это не будет важно, так как в нашем случае X = L2 . Уравнение PnApn = Pnf эквивалентно следующему:

(APn )(xj ) = f (( ) j = 1,^,n-

Представим приближенное решение в виде линейной комбинации ба-

n

зисных функций: pn = ^ c^v^ . Подставив это представление в схему метода

к=1

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

k=1

Основная трудность применения метода коллокации в данной работе связана с тем, что в качестве пространства X рассматривается пространство = ¿2 X ¿2 X ¿2, в котором значения функции в точке, вообще говоря, не определены. Таким образом, оператор проектирования Рп : X ^ Хп определен не на всем пространстве X и, вообще говоря, не ограничен. Это приводит к тому, что нельзя применить стандартные утверждения о сходимости проекционных методов. Однако в нашем случае правая часть / является

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

Определение 2.1. Метод коллокации будем называть сходящимся для оператора Л и / е 1т Л , если существует число N такое, что приближенные

уравнения (Лфп)(ху) = /(ху), у = 1,..,п, имеют единственное решение

Фп е Хп для всех п > N , и если эти решения сходятся фп ^ ф при п

к единственному решению ф уравнения Лф = / .

Рассмотрим вопрос о построении схемы для метода коллокации для уравнения (1).

Представим уравнение (1) в виде системы трех скалярных уравнений:

3

2(х)- к0 |°1Е (х,УУ1 (у)Ф -м <2

-Xdivx JGE (x,y)J(y)dy = E0l (x), l

dxi

= 1,2,3.

Q

(1 2 3 \

Jn, Jn, Jn ) следующим образом:

n n n

Jn = Zakfk (x) J« = Zbkfk (x) Jl = Zckfk (x) k=1 k=1 k =1

где fk - базисные функции-«ступеньки».

Ниже проводится построение функций fk . Будем считать, что Q - параллелепипед: Q = {x: a1 < x1 < a2, b1 < x2 < b2, C1 < x3 < C2 } .

Разобьем Q элементарными параллелепипедами Qj = n^m :

Пklm ={x : x1,k < x1 < x1,k+1, x2,l < xl < x2,l+1, x3,m < x3 < x3,m+1} ;

x1,k = a1 + a2 a1 (k -1), x2,l = b1 + Ь2 Ь1 (l - 1), x3,m = C1 + C2 C1 (m - 1),

n n n

где k, l,m = 1,...,n .

Получим формулы для fklm :

f = J1, x ^ ^klm,

Jklm = ] р. _ гг

l0, x ^nklm.

Построенное множество базисных функций удовлетворяет необходимому условию аппроксимации в Ь2 = ¿2

X ¿2 X ¿2.

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

Л11 сч Л13 Бл

Л Л22 Л23 Б2

Л Л32 Л33 Б3 і

где элементы столбцов и матриц определяются из соотношений

4 = 4 (х );

ЛЫ = 8у - 8к/ко |°Е (,У) (у)Ф |-Х°Е (,У) (у)dУ, (5)

Ч д х1

б

а координаты точки коллокации имеют вид

Х = {х11,х12,х13 ) хп = (71 + 1/2)^1, х72 = (72 + хг3 = (з +

И = °2 а1, ^2 = Ъ^-Ъ-, Из = —2 —, к,I = 1,2,3; 7 ,7 = 1,...,#; N = п3 .

п п п

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

дифференциального уравнения.

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

этого достаточно рассмотреть интегралы вида —— Г—-С1Е (х^,у) (у)<^у,

дхк^ дх1 К '

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

д г д

1

Эхк ■> Эхі

к б I

х] - У

(6)

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

М

(х У-11

1

Эхк б дхі Iх - у|

(7)

в виде действия ПДО на базисную функцию и вычислить его аналитически. Учитывая формулу (4), можно показать, что интеграл (7) будет иметь

вид

Ikl íy) = _L fff+^2 +y3^3) sin^l^Lsinsin ^ lkhd12 ,

JJJ 2 2 2 ^ |É|2

, s , s , s , , ^/2 hh/2 h/2

ГДе -^Sin ^l^Sin ^sin ^ =1 f f f ^^+^2+^3)) есть

^2 ^3 2 2 2 8 J J J

123 -h^/2-h2/2 -h,/2

преобразование Фурье элементарного параллелепипеда П, центр которого

расположен в начале координат, а ребра параллельны координатным осям.

Интегралы Ik в точке коллокации y = У2 = У3 = 0 вычислены в п. 3.

Приведем здесь значения этих интегралов:

Т11 2 h2 h3 2 h2 h3

I = — arctg----, = = — arcsin- z J •

UlVlg I---------------- — UlVOlll I---------- I--------- ,

n h^yj h2 + h| + h2 я д/hj2 + hf д/ hj2 + h2

T22 2 . hjh 2 . hjh

1 = — arctg-------------------------------------------------------, — = — arcsin^ ;

я h2^¡ hj2 + h| + h2 я д/hj2 + hf д/h| + h2

T33 2 hh2 2 hh2

I = — arctg-----------------------------------------------------------. 1 2 = = — arcsin^ 1 2 ;

л h3yj hj2 + h| + h2 л ^/hj2 + h2^ hf + h32

I !2 = 121 = i 23 = i 32 = i 13 = i 31 = g

> g, i = 1,2,3.

что совпадает со значениями из [1, с. 121]. Отметим, что |

Таким образом, интегралы (6) (значения интегралов (7) в точках коллокации) ограничены константой, не зависящей от шагов Ъ\_, ^>, И3. Остальные части в формуле (5) представляются непрерывными функциями и также могут быть ограничены константой, не зависящей от шагов Ъ\_, ^>, И3. Следовательно, мы получаем следующее

Утверждение 2.1. Существует константа М такая, что для коэф-

фициентов (5) верно неравенство

Akl

< M, причем M не зависит от

\, Й2, Нъ и /, ], к, I.

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

Введем «-мерные пространства Rj, Rf и R« с нормами, соответст-

венно,

I, - ZI ь

І=1

1

ZЫ ■ И~- ."HM■

. ! 1<i<j

І=1

Т

где Ь = (*!,...,ьп ) .

Будем рассматривать конечномерные (матричные) операторы Ап : Д™ ^ ДД. Соответствующая операторная (матричная) норма имеет вид

14

n I1_Z

= max

1<i, j <n

где a\j - коэффициенты матрицы An = j a;y

Anb = c , b є R1, c є RZ, то

і,j=1

. Действительно, если

n ( ) f { \ \

тах (Anb - = max Z aA" bi j =1 < max max (n ) aij

1<г<п 1<i<n v 1<i<n1< j <n J

III •

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

Рассмотрим матричное уравнение

(An + Bn ) b c , где An : Rn _ RZ, Bn : Rn _ RZ, b є Rn , c є RZ.

(8)

Лемма 2.1. Если существует обратная матрица Ап и для всех п вер-1 , то уравнение (8) имеет единственное реше-

на оценка ||В,

,-1

°_1

ние при всех п .

Доказательство. При выполнении условий леммы уравнение (8) мож-

-1

но переписать в виде Ап + Ап^Бп | Ь = с. Так как

An Bn

<

1_1

°_1

В

n I11_Z

< 1, то решение уравнения (8) существует и единственно,

и имеет вид b = (l + Ап 1В;

-1 A-1c

'ni n ь ■

Если Ап : ДД ^ Д™ и все а— ' > 0 (или все а\- ' < 0) или матрица явля-

„ , (п) й (п )ч

ется диагональной (а—' = о-а^- '), то для соответствующей операторной

(матричной) нормы верна формула

п

IlAn IL_n Z

h м

aj

(9)

где ajj”- - коэффициенты матрицы An = \аА‘

i, j =1

Действительно, если Anb = c , b e , c e , то

ICI l = ï |(Л,ь ) = ï

n n f n \

Z X"1 (n )a Z aij bj < Z (n ) aij

і=1 j=1 1i,j=1 J

n

n

Если все компоненты вектора b равны между собой, то указанное вы-

(п) ^ п (п) ^ п

ше неравенство перейдет в равенство с учетом условия а); ' > 0 или а)■ ' < 0.

У У

Если матрица является диагональной, то можно выбрать вектор b с компонентами bi = sign | a(n)). Тогда снова неравенство перейдет в равенство, что и

(10)

доказывает формулу для нормы.

Рассмотрим конечномерные уравнения метода коллокации:

Л-хы = Ь,

где

A12 A13 ^ f B11 f J11

AN - A21 A22 A23 , b - B2 , u - J2

1A31 A32 A33 J VB3 J V J3 J

Теорема 2.1. Пусть тензор г(х)е С(2) диагональный, вещественнозначный и £ц (х)> 1, хе 2 (I = 1,2,3). Тогда существует N такое, что при N > Ы0 решения уравнений (10) существуют и единственны.

Доказательство. Из условий теоремы сразу следует, что матрица

((х)-1)е С(2) обратима в 2 , ((х)-1) 1 е С(2).

Представим матричные коэффициенты (5) в виде

ЛЫ = Ск1 + Бк1, где Ск1 = % (+ ^ (^ )), д г д

Dh - _5k/ko j ge (xJ, y f (у )dy j xX~GE (xJ , у f (у - (xJ f

Q 4 Q xi

Запишем уравнение (10) в виде

(СЫ + )ы = ь

с матрицами Сы и Бы , имеющими коэффициенты СЦ и соответственно. Из утверждения 2.1 следует, что для коэффициентов матрицы Бы выпол-

няются неравенства

DJ

< М. Матрица СN является диагональной и, очевидно, обратима (здесь мы учли, что Iі |х^' |> 0, I = 1,2,3.). Неотрицательны

будут и все элементы диагональной матрицы СN . Тогда для нормы обратной

= О XN) при N —— ^

матрицы верна формула (9). Ясно, что

С-1

CN

(здесь символ Е = О (Ы) означает, что существуют не зависящие от N константы С1 > 0 и С2 > 0 такие, что верно неравенство С^ < |е| < С2N). По-

скольку интегралы, входящие в Б1^ , берутся либо от ограниченных функций, либо от функций, имеющих особенность О (г-1), для нормы будет верна

оценка ||Бы = О (-2) при N ^^. Действительно, интеграл

I

-dy легко оценить, заменив интегрирование по параллелепипеду Qj

- , J - - - - }

Qj I - y

интегралом по описанному около параллелепипеда Qj шару (от той же

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

Для доказательства теоремы остается заметить, что левая часть неравенства в оценке в лемме 2.1 стремится к 0 быстрее, чем правая часть, поэтому начиная с некоторого No эта оценка будет выполняться и, следовательно, уравнения (10) будут однозначно разрешимы.

Доказанная теорема 2.1 устанавливает разрешимость конечномерных уравнений в методе коллокации при некоторых ограничениях на тензор s(x). Заметим, что эти ограничения выделяют широкий класс диэлектриков, используемых на практике.

3. Вычисление интегралов

Рассмотрим интеграл

IiJ =_L fff /Ml +y2^2 +y3S3 )sm Äsm ^2 sin ^ Si Sjd S 2 . (11)

JJJ 2 2 2 S1S2S3 |S|2

Мы считаем, что Il = Il (, y2, y3). Достаточно вычислить интегралы

только двух типов Ikl при к Ф l и Ikk для любых конкретных к, l = 1,2,3, значения остальных интегралов можно получить из соображений симметрии. Далее мы будем работать с интегралами 11 = IlJ (0,0,0), т.е. со значениями рассматриваемого интеграла в точке коллокации. Легко видеть, что замена hi := h /2, i = 1,2,3 не изменяет Il (0,0,0), будем этим пользоваться для сокращения записи.

12

Рассмотрим I , поскольку подынтегральная функция нечетна по S1 (и по S2), то

I12 = “7fff sin(S1h1)sin(S2h2 )sin(S3h3 )-dS2 = 0 .

n3JJJ S3 |S2

Отсюда получаем, что 112 = 121 = 123 = 132 = 113 = 131 = 0 .

Заметим, что интеграл I12 (, y2, y3) можно вычислить точно даже в произвольной точке y .

Еще один интеграл типа Iг/ , который необходимо вычислить, это

122 =

-ГШSln(1h )Sln(2h2 )sin(^3h3) ^2 • (12)

Л Sb3 |s|

Сначала вычислим через вычеты интеграл по ^2, получим

I = | sin (i^ (iff2 = »-«', где 5'^?2 + Й • i tf

Теперь в оставшемся двойном интеграле используем формулу Эйлера

eix _ e_ix

sln x =-------- и переходим к полярным координатам: ^ =р cos ф,

2i

^2 = Р sln ф • Здесь, учитывая формулу

I e_U* _e-vx v

J----------dx = ln-, (13)

J x u

0 r

где Reu >0 и Rev > 0 [3, с. 348] и проводя простые преобразования, получаем

122 = _1_42ln h22 ( + tg2 ф) + (1 + h3tgф)2 dф _

2я2 0 h2 ^1 + tg2 ф) + (( _ h3 tgф)2 tgфcos2 ф 1 Ч2 h2 ( + ctg2 ф) + (^1 _ h3ctg ф)2 dф

_^2 J l^

2л 0

h2 (1 + ctg2 ф) + ( + h3 ctgф)2 ct^sin2 ф

Теперь в первом интеграле делаем замену t = tg ф, а во втором -t = ctg ф, после некоторых преобразований получаем, что

122 = J_ Jln h2 ^ +12 ^ + ^ + ^^ dt

^2 0 h2 (1 + t1 ) + (h1 _ h3{)2

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

7 1п xdx (1п Р)2-(1п у)2 Г3 5471

-----—------ = ----,---^— 13, с. 5471 и помня о том, что

0 (х + Р)(х + у) 2(р-у) ,

1п (-1) = л/ + 2гак , к е Я, получаем

J 22 = — = Я2

12+h2

2га(1 + к +1)ln 1—2-2 + 2%(l _к)arctg

h1h3

h2 + h3 h2^j hj2 + h2 + h2

22

Легко видеть из (12), что Im I = 0, отсюда получаем, что 1 + k +1 = 0 .

22

Можно показать, что I

=1 (также см. [1]). Учитывая последнее,

1\ =Н- =к3 =1 3

получаем, что I - к = 1. Отсюда находим, что к = -1, I = 0. Окончательно имеем

122 = -агав ^

Замечание. Интеграл 122 не отражен в известных справочниках [3, 4]. Теперь можно выписать значения остальных интегралов в точках кол-локации, просто циклически переставляя индексы (см. п. 2).

В качестве модельного примера можно рассматривать куб со стороной h = 1, т.е. h = ^2 = Нз = 1. Из предыдущих формул получаем

111 = 122 = 133 = 1

3 .

Поясним кратко, как вычисляется интеграл

1 = -1 f sin ?! sin ?2 Sin 53 ^2^2 = 3 >

^3 |5|2 3

j22

который есть I

. Схема такова: сначала берется через вычеты ин-

h —h^2 —h?3 —1

теграл по , затем вводятся полярные координаты и используется формула

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

—/ 2 2

Г ln(1 + psinx)—— = ——1 (arccosp)2;

J v ’ sin x 8 2 ’

0

—2 2

f 1 \ dx — 1 / ч2 2 ^ i

I ln(1 + pcosx)-------=-(arccosp) , p < 1

cos x 8 2

0

2

(последние интегралы при p < 1 приведены в [3]).

Список литературы

1. Самохин, А. Б. Итерационные методы в электромагнитном рассеянии / А. Б. Самохин. - М. : Радио и связь, 1998.

2. Валовик, Д. В. Метод псевдодифференциальных операторов для исследования объемного сингулярного интегрального уравнения электрического поля / Д. В. Валовик, Ю. Г. Смирнов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2009. - № 4. - С. 70-84.

3. Градштейн, И. С. Таблицы интегралов, сумм, рядов и произведений / И. С. Градштейн, И. М. Рыжик. - М. : Физматгиз, 1962.

4. Прудников, А. П. Интегралы и ряды. Элементарные функции / А. П. Прудников, Ю. А. Брычков, О. И. Маричев. - М. : Наука, 1981. - Т. 1.

Валовик Дмитрий Викторович

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

E-mail: dvalovik@mail.ru

Смирнов Юрий Геннадьевич

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

E-mail: mmm@pnzgu.ru

Valovik Dmitry Viktorovich Candidate of physical and mathematical sciences, associate professor, sub-department of mathematics and supercomputer modeling,

Penza State University

Smirnov Yury Gennadyevich Doctor of physical and mathematical sciences, professor, head of sub-department of mathematics and supercomputer modeling, Penza State University

УДК 517.9, 519.6 Валовик, Д. В.

Метод коллокации для решения уравнения электрического поля /

Д. В. Валовик, Ю. Г. Смирнов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2010. - № 4 (16). -С.89-100.