Применение ГРИД-технологий для решения объемного сингулярного интегрального уравнения для задачи дифракции на диэлектрическом теле субиерархическим методом Текст научной статьи по специальности «Математика»

МАТЕМАТИКА

УДК 518.1

М. Ю. Медведик, Ю. Г. Смирнов

ПРИМЕНЕНИЕ ГРИД-ТЕХНОЛОГИЙ ДЛЯ РЕШЕНИЯ ОБЪЕМНОГО СИНГУЛЯРНОГО ИНТЕГРАЛЬНОГО УРАВНЕНИЯ ДЛЯ ЗАДАЧИ ДИФРАКЦИИ НА ДИЭЛЕКТРИЧЕСКОМ ТЕЛЕ СУБИЕРАРХИЧЕСКИМ МЕТОДОМ1

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

Введение

В работе исследуется задача дифракции стороннего электромагнитного поля на локально неоднородном теле, помещенном в свободном пространстве. Актуальность работы определяется применением результатов исследования, например, при решении задач дифракции в СВЧ-диапазоне. Для численного решения задачи возможно использование метода конечных элементов. Однако прямое применение метода конечных элементов встречает ряд трудностей. Во-первых, краевая задача для системы уравнений Максвелла не является эллиптической, поэтому «не работают» стандартные схемы доказательства сходимости проекционных методов. Во-вторых, для получения приемлемой точности расчета поля в теле с большой диэлектрической проницаемостью необходимо выбирать достаточно мелкую сетку, что влечет также выбор мелкой сетки и в объеме вне тела (выбор же сетки разного масштаба внутри и вне тела ведет к неверным результатам). А это, в свою очередь, с учетом трехмерного векторного характера задачи приводит к разреженным матрицам очень больших порядков в методе конечных элементов.

От этих недостатков свободен метод объемных сингулярных инте-гральнь1х уравнений. Здесь оператор получается эллиптическим, а интегральное уравнение решается только внутри тела (в области неоднородности).

Мы изучаем интегральное уравнение, опираясь в основном на теорему эквивалентности краевой задачи и интегрального уравнения, а также на результаты исследования свойств объемного сингулярного интегрального уравнения [1]. На этом пути удается доказать теорему о существовании и единственности решений в ¿2 интегрального уравнения, сходимость численного метода Галеркина, получить некоторые результаты о гладкости решений.

1 Работа выполнена при поддержке РФФИ, грант № 06-07-89063а.

2

Постановка задачи для системы уравнений Максвелла

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

магнитной проницаемостью и положительной (3 х 3 )-матрицей-функцией (тензором) диэлектрической проницаемости ё(х). Компоненты ё(х) являются ограниченными функциями в области Q , ё (@), а также е-1 е (Q). Граница дQ области Q кусочно-гладкая. Точнее, предположим, что для

каждой точки границы хо е дQ существует окрестность 0 (в Я ) и

2 3

С -диффеоморфизм этой окрестности на Я , при котором точка хо переходит в точку 0, а образом множества 0nQ является множество одного из

2

следующих типов (ниже (Х1, Х2, Х3) - декартовы; (г, 0), г > 0, 0е S - сфери-

3

ческие координаты в Я ). Либо Х1 > 0 (Х0 - точка гладкости границы); либо

3

Х1 > 0, Х2 > 0 (Х0 - точка на «выходящем» ребре); либо Я \{Х1 > 0, Х2 > 0}

2

(Х0 - точка на «входящем» ребре); либо г > 0, 0 е Q', где Q' с S - односвязная область с кусочно-гладкой границей дQ' (Х0 - вершина «конуса с ребрами»). В частности, если дQ' - гладкая, то Х0 - коническая точка; если дQ' образована дугами больших окружностей, то Х0 - вершина многогранного угла. Пусть Q - ограниченная область и каждая точка Х е дQ принадлежит одному из этих типов. Тогда будем говорить, что Q - область с кусочно-

гладкой границей.

Требуется определить электромагнитное поле Е, Н е ¿2^), возбуждаемое сторонним полем с временной зависимостью вида в~гШ. Источник стороннего поля - электрический ток ] 0 или падающая плоская волна.

Будем искать электромагнитное поле Е, Н , удовлетворяющее уравнениям Максвелла, условиям непрерывности касательных компонент поля при переходе через границу тела и условиям излучения на бесконечности

гоН = -г'юеЕ + ]<° ; гОЕ = г'юц0Н; (1)

[ Е ]1 ьа =[ Н ]1 = 0; (2)

д ( Е ^ ., ( Е \ _-К ( Е \ „-1

V Н

- їк

V Н

= о( я-1),

V Н

= 0(я~1), я := х (3)

дя

Здесь к0 - волновое число свободного пространства (вне Q ).

Объемное сингулярное интегральное уравнение

Сведем краевую задачу к объемному сингулярному интегральному уравнению. Перепишем (1) в эквивалентной форме:

гоШ = -їюєо Е + Іе ; го^Е = їюЦо Н, (4)

3

где

^ ^ 0 , ^ р

]е = ]е + ]е ■

(5)

В последнем равенстве ]Р = -гт(£(х) - £01)Е - электрический ток поляризации. Нетрудно проверить, что решение последней краевой задачи име-

ет вид

где

Е = г'юц0 Ае--------1—Бгаё Шу Ае ;

г'юе0

Н = го^е ,

Ае = |0(г)]е (У)йу, г = \х - у\

(6)

(7)

векторный потенциал электрического тока. Потенциал АЕ удовлетворяет уравнению

ААе + к0 Ае =- ]е .

(8)

Таким образом, потенциал АЕ представляет собой свертку с функцией Грина свободного пространства для уравнения Гельмгольца, обеспечивающей выполнение требуемых условий для полей. Функция Грина имеет вид

1 -'ко\х-у\

0( х, у) =— ------ . (9)

4п \ х - у \

Определим диагональный тензор Грина О (г) = Ша§(0,О,О).

Формулы (6) не дают явного решения задачи дифракции, т.к. ток }е

зависит от Е . Из полученных соотношений для поля Е следует интегро-дифференциальное уравнение

Е(х) = Е0(х) + ко 10(г)

Q

£( У)

-1

Е(у)йу+

+§гаё Шу |О (г) Q

ё( у) £0

£0

-1 Е(у)ёу, хе Q.

(10)

Кроме того,

Е(х) = Е0(х) + ко 10(г)

+§гаёШу 10(г) Q

Q

£( у) £0

ё( у) £о

-1

Е(у)йу +

-1

Е(у)йу, хе Я \ Q.

р

Последняя формула дает представление решения Е(х) в области

3 ^ ^

Я \ Q , если Е(у), у є Q - решение уравнения (10). Поле Н выражается через решение (10) в виде

H(X) = H0(X) - z'ro£0rot IG(r) Q

g( У) e0

-1

E(y)dy, xє Q.

Сведем полученное выше интегро-дифференциальное уравнение к объемному векторному сингулярному интегральному уравнению.

Представим функцию Грина в виде

G(r) = Go(r) + G:(r),r =| x- y |; Go(r) =

eik0r -1 4nr

,Gi(r) =

1

4nr

Пусть xi,xj,X3 и yi,У2,Уз - координаты точек x,y . Рассмотрим второй интеграл в уравнении (10) для электрического поля и исследуем вопрос о возможности внесения операции grad div под интегралы

JG(r)U(y)dy (p = 0,1,2).

Q

В декартовой системе координат

ґ \

д

grad div JG(r)U(y)dy

Q h

dxl

X iG(r )Un (y)dy

n=1 °Xn Q

, l = 1,2,3.

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

Пусть а = (у -х)/ г. Сформулируем важное утверждение о дифференцировании интегральных операторов, ядро которых имеет особенность порядка 1/ г2.

Лемма 1 [1]. Пусть функция ^(х, а) имеет в Q непрерывные первые производные по декартовым координатам точек х и а, а и є Ь} (@). Тогда интеграл

т( х) = |¥( Х2 а) и( у )ёу Q Г

имеет обобщенные производные дт/дхк є Ь>^), к = 1,2,3. Эти производные определяются из формул

дю г д

------= v.p. ---------

dxk I дхк

Q

у ( x, а)

u (y )dy - u( x) Jу (x, а) cos(^k )dS,

где £,к — угол между вектором г , направленным от точки х к у, и ортом, который соответствует координате Хк .

Применяя данную лемму к сингулярному интегралу с ядром

_Э______1_

dxn 4nr ’

придем к известному представлению

д 1 рГ~ I~Un (y)dy = v'p'i‘

j ПГ 4nr J

1 1

д Д—Un (y)dy ~~^nlUn (x) ,

dxi dXn 4nr 3

дхі дхп 4пг

1 Q п Q

Ьпі - символ Кронекера.

Используя полученные соотношения, переходим от интегро-дифференциального уравнения (10) к векторному сингулярному интегральному уравнению:

E (х) + — 3

£( х) е0

-1

E (х) - v.p. J Г—( х, y) Q

g( y) є0

-1

E(y)dy -

J Г( x, y)

Q

e( y) e0

-1

E( y)dy = E 0( x).

(11)

Здесь функции Г, Г— имеют вид

Г( х, y) = k0G (r) + (• ,grad)grad G0( r), Г—(х, y) = (• ,grad)gradG—(r).

Вопрос о разрешимости уравнения (11) и об эквивалентности краевой задачи дифракции и сингулярного интегрального уравнения подробно изучен в [2].

Метод Галеркина

Одним из наиболее эффективных методов численного решения интегральных уравнений является метод Галеркина. Для уравнения Аф = f, (ф, f е X) в пространстве X метод формулируется следующим образом. Приближенное решение фп е Xn определяется из уравнения РпАфп = Pnf . Здесь фп е Xn (Xn есть n -мерное подпространство пространства X ), Pn : X ^ Xn - оператор ортогонального проектирования на конечномерное подпространство. Пусть подпространства Xn являются линейными оболочками базисных функций: Xn = span{vi,..., vn}. Потребуем, чтобы для выбранных базисных функций выполнялось условие аппроксимации

Vxе X lim inf ||x - x|| = 0. (12)

n^<™ xeXn

Уравнение РпАфп = Pnf эквивалентно следующему:

(Афп , vl )X = (f, vl) X , l = 1 •••, n. где ( • , • )x - скалярное произведение в X .

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

n

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

k=1

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

Xck(Avk,vl)X - (/,vl)X ,l -1 •••, к-1

и.

Определение 1. Метод Галеркина будем называть сходящимся для оператора A, если существует число N такое, что для каждого f е Im A

приближенное уравнение PnAфn = Pnf имеет единственное решение

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

единственному решению ф уравнения Аф = f . В этом случае имеет место

квазиоптимальная оценка скорости сходимости:

-ф||< С inf ||¥-ф||.

VeXn

Рассмотрим вопрос о сходимости метода Галеркина для векторного интегрального уравнения электрического поля (11). Сформулируем важное утверждение о сходимости проекционного метода.

Лемма 2 [Кресс]. Предположим, что А: X ^ X есть ограниченный оператор, имеющий ограниченный обратный, и что проекционный метод сходится для А. Пусть B — линейный ограниченный оператор, А + B инъ-ективен. Оператор B удовлетворяет любому из двух условий:

а) sup A-lpn ||B|| = q <1;

neN

б) B компактен.

Тогда проекционный метод также сходится для оператора A + B .

Перепишем интегральное уравнение (11) для электрического поля в виде

(I + S) E = E0, где оператор S определяется в соответствии с (11):

(SE)(х) -

£( х) е0

-1

E (х) - v.p. J Гі( х, у) Q

g( У) е0

-1

E (y)dy.

(13)

(14)

Аналогично [1] получаем

Предложение 1. Пусть тензор диэлектрической проницаемости таков, что

ess sup

xeQ

3

I

l ,n-1

1/2

eln( x)

-5

є0

nl

<

1

и выполнено условие аппроксимации. Тогда уравнение (11) однозначно разрешимо для любой правой части Е0 е Ь^(О) и метод Галеркина сходится для уравнения (11).

Рассмотрим теперь вопрос о построении схемы Галеркина для рассматриваемой задачи дифракции. Будем формулировать метод не для сингуляр-

П

ного интегрального уравнения (11), а для интегро-дифференциального уравнения (10). Этот подход оказывается эффективным в силу более удобного

представления интегралов. Будем предполагать, что матрица

£( х) £0

-1

об-

ратима в Q ,

£( х) £0

-1

є Lx (Q), I - единичная матрица.

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

£ =

£ (х) £0

-I

, J :-

£ (х) £0

-I

E,

перейдем от (10) к следующему уравнению:

AJ -J(х)-к2 J0(х,у)J(y)dy-grad div J0(х,у)J(y)dy-E0(x). (15)

Q Q

Представим это уравнение в виде системы трех скалярных уравнений:

3 д -

X\liJl (х) - ко JG(х, y) J (y)dy - — divх JG(х, y) J(y)dy - E0l (х), l -1, 2, 3.(16)

i-1 Q х Q

Определим компоненты приближенного решения J следующим образом:

N N N

J1 - Xак/1(x), J2 - Xbkfk(x), J3 - Xck/k(x),

к-1

к -1

к-1

где /к - базисные функции-«крышки», существенно зависящие лишь от переменной . Ниже проводится построение функций /1. Будем считать, что Q - параллелепипед: Q = {х: а < Х1 < «2 Ъ < Х2 < Ъ2, с < Х3 < . Разобьем

2 параллелепипедами

Пк/т = {х: х1,к-1 < Х1 < х1,к+1, х2,1 < Х1 < х2,/+1, х3,т < х3 < Х3,т+1};

«2 - «1 , „ Ъ2 - Ъ| , „ с2 - с

х1к = а + —-----1 к, х2/ = Ъ + 2—----1 /, х3^ = С1 + 2—--1 т,

п п п

где к = 1,...,п-1; /,т = 1,...,п/2-1.

Обозначив й1 :=| х1 к - х1 к-11, получим формулы для /1/т :

/Ыт

1 -“ТІ х1 - х1,к І, хє ПЫт, h1

0, х і П

1

Ыт.

2 3

Функции /к/т , /кт , зависящие от переменных х2 и х3 соответственно, определяются аналогичными соотношениями. Из определения базисных

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

Перенумеруем базисные функции:

Л1, /к, /к, к = 1,..., N,

где N = —(п3 - п2).

4

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

Г Аі А12 А13 В1'

А22 А23 В2

-2Р А32 А33 В3 і

элементы колонок Вк и матриц Ли определяются из соотношений:

Вк = (4,),

| О (X, у)/■ (у )ёу, $

2

ч— к— О( X У) /)(УЖ ії дхк •> ох/

к 2 1

(17)

к, I = 1,2,3; г, ; = 1,..., N.

Преобразуем второе скалярное произведение в (17). Применяя к внешнему интегралу формулу интегрирования по частям, получим:

г л с \

IОЕ(х, У)/1](У) ^, /к

дхк •> дх/

к д і

2

I дг~ ОЕ(х, у )/1/ ( у ^, -т" /і

~ дХ/ —х

дхк ■"

. (18)

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

Ак/ = (^к///, ) - §к/к0

|О (х, у)/-(у )dy, /к

2

[О( х, у) /• (у )dy, /і

* д-к —х

Предложенный метод Галеркина реализован для решения ряда задач дифракции.

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

H (x) = H 0(x) + JG(г)

Q

Д( у) ^0

-1

H(y)dy +

■ grad div JG(г)

Q

Д( y) ^o

-1

H(y)dy, xe Q,

E(x) = E0(x) + гюцоrot JG(r) Q

Д( y) ^0

-1

H(y)dy, xe P.

(19)

Формирование матрицы

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

Аы = I £и// (х)/к (хК- 8ико I I °(x,У)/] (у)/к (х)^х +

пк П';

1 J

+11 ° (х у')дг/1] (у)д^" /гк (х)дУёх.

П1; пП k

J 1

Пк П j

1 J

(20)

Будем предполагать сетку равномерной с шагом h по всем трем направлениям. Будем также предполагать, что ai = bi = ¿i = 0, Ü2 = ¿2 = ¿2 = 1. Тогда носитель базисной функции состоит из двух примыкающих друг к другу кубиков Пк . Индекс к определяет направление, в котором располагаются примыкающие кубики. При к = 1 направление вдоль оси 0x ; при к = 2 - направление вдоль оси 0y ; при к = 3 - направление вдоль оси 0z . Индекс i

является тройным мультииндексом i = (1, i'2, Í3 ), он определяет положение первого из двух кубиков ¿1, ¿2, ¿з - номера кубиков по направлению осей 0x ,0y ,0z соответственно (рис. 1). Индексы l и j вводятся аналогично.

Число к0 является волновым числом, символ £¿i определяет диэлектрическую проницаемость среды, а 8И - символ Кронекера, равный нулю при к Ф l. Если вместо параллелепипедов рассматривать кубики, то построенная сетка примет следующий вид. Вся матрица заполняется путем вычисления

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

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

Рассмотрим три случая (рис. 2): пусть пересечение происходит по паре

носителей, изображенных в левой части рисунка, тогда интеграл равен

2Г

~з

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

нъ нъ З 1

первого интеграла равны — и — соответственно. Здесь п - длина грани

6 4

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

П(1), П(2) от функции /о (х, у) вычисляется по формуле

п,п+1

| | /о(х, у)ёхёу = X /(х(к\ у(т))А!3д2 , (21)

п(1) п(2)

к ,т

здесь X = (хьХ2,Хз), у = (У1,У2,Уз), Д1 = -, д2 =—Ц-,

п п + 1

х(к) = „(с) _ 1 + . + А1 у(т) = у(с) _ п + т +А1

х1 = „1 2 + к' А1 + ~, у1 = у1 2 + т' +~,

х(к) = х(с) _ п ' А + А1 (т) = л (с) _ п ^ + А2

Л^'~у — Л-'

2 ~ Л2

_2 + к'А1 + -^, у2т) = у2с) _2 + т'А2 + -2-,

х(к) = х(с) _ п + к 'А +А1 у (т) = у(с) _ п + т 'А I А 2

х3 = х3 2 + к 1 + ~, у3 = у3 2 + т ' А2 + ~,

здесь (х|с), „2с), „Зс)), (у|с), у2с), уЗс)) - центры кубиков П(1), П(2). В случае

совпадения носителей при интегрировании второго интеграла в (20) возникает слабая особенность, поэтому при численном интегрировании в одном из носителей выбирается п точек интегрирования, а во втором п +1 по каждому из направлений. Это приводит к избавлению от особенности за счет несовпадения узлов интегрирования.

Рис. 2

Для заполнения всей матрицы нас будут интересовать элементы, вычисляемые по двум комбинациям носителей: первый - оба носителя ориентированны по оси 0х; второй - один ориентирован по оси 0х, а второй - по оси 0у . Остальные элементы матрицы вычисляются на основе этих за счет внутренней симметрии матрицы.

Субиерархический, параллельный алгоритм

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

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

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

Применение ГРИД-технологий

Как отмечалось раннее, решение поставленной задачи может оказаться очень емким с точки зрения вычислительного процесса. Использование су-биерархического подхода позволяет решать задачи самой разной вычислительной сложности от очень простых, когда выбрано несколько носителей, до очень сложных - порядка ста тысяч. Для быстрого решения простых задач достаточно ресурсов небольшого кластера, в то время как решение сложных задач требует использования ресурсов самых современных кластеров. В связи с этим естественно использовать ГРИД-технологию к решению столь сложной вычислительной задачи. Используя web-интерфейс, позволяем пользователю задать его параметры счета и геометрию задачи, и перед решением задачи на кластере производим анализ сложности запускаемой задачи. В случае если задаче требуются небольшие вычислительные ресурсы, она решается на мини-кластере, в случае использования больших вычислительных ресурсов решение производится на более мощных кластерах. В качестве кластера для небольших вычислительных задач предполагается использовать мини-кластер ПГУ, для задач со средней сложностью кластер ЮУрГУ, для самых мощных задач предполагается использование кластера НИВЦ МГУ. Функциональная схема использования ресурсов при решении задачи отображена на рис. 3.

Использование данных кластеров и ресурсов СКИФ-ГРИД полигона

будет производиться в соответствии с проектом союзного государства

СКИФ-ГРИД, государственный контракт от 16 июля 2007 г. № СГ-2/07.

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

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

2. Смирнов, Ю. Г. Существование и единственность решения объемного сингулярного интегрального уравнения в задаче дифракции / Ю. Г. Смирнов, А. А. Цу-пак // Дифференциальные уравнения. - 2005. - Т. 41. - № 9. - С. 1190-1197.

3. Колтон, Д. Методы интегральных уравнений в теории рассеяния / Д. Колтон, Р. Кресс. - М. : Мир, 1987.

4. Медведик, М. Ю. Параллельный алгоритм расчета поверхностных токов в электромагнитной задаче дифракции на экране / М. Ю. Медведик, Ю. Г. Смирнов, С. И. Соболев // Вычислительные методы и программирование. - 2005. - Т. 6. -С. 99-108.

5. Медведик, М. Ю. Субиерархический параллельный вычислительный алгоритм и сходимость метода Галеркина в задачах дифракции электромагнитного поля на плоском экране / М. Ю. Медведик, Ю. Г. Смирнов // Известия высших учебных заведений. Поволжский регион. - 2004. - № 5. - С. 3-19. - (Естественные науки).

6. Медведик, М. Ю. Субиерархический параллельный вычислительный алгоритм для решения задач дифракции электромагнитных волн на плоских экранах / М. Ю. Медведик, Ю. Г. Смирнов // Радиотехника и электроника. - 2007. - Т. 6. -№ 4. - С. 1-6.