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

УДК 517.3

Ю. Г. Смирнов, Д. И. Васюнин

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

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

Ключевые слова: диэлектрическая проницаемость материала, обратная краевая задача, итерационный метод.

Abstract. The article investigates a problem of permittivity determination of dielectric body located in rectangular waveguide. Nonhomogeneous body has an arbitrary shape. The authors offer an iteration method for numerical problem solving and proove its convergence. The article also introduces numerical results of determining dielectric body permittivity.

Key words: permittivity of dielectric body, inverse boundary value problem, iteration method.

Введение

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

В статье исследуется задача определения диэлектрической проницаемости неоднородных образцов материалов произвольной геометрической формы, помещенных в прямоугольный волновод с идеально проводящими стенками. В работе [4] задача сведена к решению нелинейного объемного сингулярного интегрального уравнения. Интегральное уравнение изучали, опираясь на результаты исследования соответствующей краевой задачи и теорему эквивалентности краевой задачи и интегрального уравнения [5]. Была доказана теорема о существовании и единственности решений нелинейного объемного сингулярного интегрального уравнения и обратной краевой задачи для определения эффективной диэлектрической проницаемости наноматериалов [6-8]. Численные результаты для случая однородного тела были получены в [9]. Некоторые особенности реализации численного алгоритма представлены в [10].

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

Пусть в декартовой системе координат Р = {х: 0 < xi < a, 0 < Х2 < b, —го < X3 < го} - волновод с идеально проводящей поверхностью дР . В волно-

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

ограниченной функцией в области Q , ее Ьх(0), а также е_1 е Ьх(0). Граница дQ области 0 кусочно-гладкая.

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

Как показано в [4, 5], рассматриваемая обратная задача может быть сведена к следующей задаче для нелинейного объемного сингулярного уравнения.

Введем ток:

е(х) _ 1 е0

(1)

где Е (х) - электрическое поле.

Тогда электрическое поле выражается через ток по формуле

Е (х ) =

е(х) _ 1

е0

_1

3 (х).

(2)

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

;(х) _ 1

-1

Г(х) = Е0 (х) + к01Ое (х,у)3(у)ёу -

б

^гаааіу |Ое (х,у3(у)<яу, хє б.

б

(3)

Здесь 0^ (х, у) - (известный) диагональный тензор Грина [4, 5, 10] с компонентами:

^ ^ —У хз —уз

е пт 3 3 пп . пт пп . пт

, , ----------СОЭ хі 8ІП-х2 008 — уі 8ІП------------------у2 ;

аЬП=0 т=1 У пт(1 + 80п) а Ь а Ь

^ ^ —У хз —уз

е чіт\ з -^і пп пт . пп пт

, , -------------8ІП хі 008---х2 8ІП у 008---------------------у2 ;

аЬ А т=0 У пт(1 +80т) а Ь а Ь

3 2 е Упт1хз у3 , пп .пт . пп .пт

Ое = — > > ----------------8ІП-хі8ІП — х2 8ІП-----------------------Л8ІП — у2 .

аЬ Упт а Ь а Ь

п=1 т=1 пт

В этих выражениях

у nm

пп

a

пт

2 — k0

при этом ветвь квадратного корня выбирается так, чтобы 1т упт > 0 и Яе упт > 0, если 1т упт = 0 . В формуле (3) к0 - волновое число свободного пространства, к^ = ю2Є0^0; ю - круговая частота. Параметры волновода выбраны так, чтобы п/а < к0 <п/Ь . В этом случае в волноводе может распро-

страняться только одна мода. E0 (x) = Є2A(+)iro^^—sin

nxi —iY(2)x3

——c *; 3

- извест-

п --8

а а

ное падающее поле (мода в волноводе); А+ - (известная) амплитуда пада-

ющей волны, у; нат.

(2) =

к2 п2 й й

к0-2 ; е2 - второй орт в декартовой системе коорди-

а

Дополнительное асимптотическое уравнение запишется в форме [4, 5]

Qi = A(

+ k0

i

Ьуі0ІПЮЦ0

Isi

Q

ny1 іу(2)у3 sin—- e 1 3

E(y)•e2dy, (4)

где уі0 =

П2 — k2

2

0

a

Коэффициент 01 ' считается известным из измерений. Требуется определить диэлектрическую проницаемость е(х),хе 0, посредством серии измерений.

Поскольку количество измерений должно быть конечным, то и неизвестных параметров также должно быть конечное число. Поэтому будем предполагать, что тело 0 состоит из N подобластей Qj таких, что

0=и , Qi П = 0, / Ф } . Мы предполагаем, что е(х) = е(}) при хе ,

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

При измерениях изменяются частоты ю(1),ю(2),...,ю(N) (происходит

сканирование по частоте); при этом волновое число изменяется по формуле

4° =®°УеоЦо .

2. Формулировка итерационного метода

Будем предполагать, что тело имеет форму параллелепипеда 0 = {х: а1 < х1 < ^2, 6 < х2 < &2, С1 < хз < ^2) . Выберем равномерную прямоугольную сетку в 0 , образованную элементарными параллелепипедами:

ПЫш {х • х1,к ^ х1 ^ х1,к+1, х2,1 ^ х1 ^ х2,/+1, х3,ш ^ х3 ^ х3,ш+1} ;

°2 _ а1 у т. Ь2 - Ь1 с2 - С

XI к = а + —---------- к, х2 1 = 01 + —-1, х3 ш = с + —------ ш,

1,к 1 Щ 2,1 1 N 3,ш 1 Щ

где к = 0,...,N1 -1, I = 0,...,N2 -1, ш = 0,...,N3 -1. Перенумеруем эти элементарные параллелепипеды с помощью одноиндексной нумерации П5, 5 = 0,...,N0 -1, N0 = N1N2Щ.

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

«М _1

е0 _

1п (х )J п (х)_ ^01 & (х, у р п (у )Ф-

(5)

е

Р — А + к)2

^гаааіу |&(х,у)п (у)dy — Е0 (х), хє Є ; (6)

е

Еп (х ) — ^п (х )J п (х); (7)

— 1вш Г^1 егу( ^уз Лп+1(У)Еп (у)• е2Ф , (8)

Уі0 ^ I а )

аЪЧ10 е

где

Р — е1(+), А — =0а(+) , Лп (х) —

п (х) _ 1

е0

^п(х)— Лп 1(х)•(9)

По этим формулам вычисление производится следующим образом. Сначала выбираем начальное приближение £0 (х ) = £е (п = 0), где £е = £е£-,

ее# - эффективная диэлектрическая проницаемость тела, вычисленная как

решение обратной краевой задачи с постоянной диэлектрической проницаемостью [4, 5]. Нельзя взять в качестве £е = £0, так как по формуле (7) нельзя определить электрическое поле. По формуле (5) вычисляется значение ^(х) .

Далее по формуле (6) определяется ток Зп (х) как решение интегро-дифференциального уравнения методом коллокации. Затем по формуле (7) по току определяем электрическое поле Еп (х) на сетке. Данную процедуру

проводим N раз при различных значениях к0 = к0(1), к0 = к0(2),..., к0 = к0(N). Таким образом, получаем N значений полей Е^\ еП2),..., еП^) при различных к0(1), к0(2),..., к0(N). На этом заканчивается вычисление на первом «слое».

На втором «слое» по известным значениям полей еП^ (х) ( = 1,...,N) из формулы (8) определяем новое значение ^п+1 (х). Для этого потребуется

произвести решение системы линейных алгебраических уравнений (СЛАУ), составленной из уравнения (8), относительно неизвестных параметров. При

этом «коэффициенты прохождения» ¥^ = ¥ (к^) находятся с помощью измерений. Считаем, что Л(+^ = 1.

Мы предположили выше, что Лп (х ) = Ли^') при X е Qj . Кроме того, пусть подобласти Qj состоят из объединения элементарных параллелепипедов сетки Qj = и П/ . Мы будем считать также, что ЕПг) (х) = ЕПг,/) при /

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

Формула (8) приводит к конечномерной СЛАУ. Образуем матрицы

г (// ^ N, М0 , ,м0,М

Е = |в2 • Е^ ^ 1/ 1 размера N X N и Н = \Н/.)/ 1 . 1 размера N х N, где

Ну = 0 при / таких, что П/ £ Qj. Тогда будем иметь СЛАУ с матрицей

ЛN = ЕН размера N х N:

^ Лп+1 = В, (10)

Т

которая решается относительно неизвестных лп+1 = (лП+1,. .,лП+1) .

Здесь и ниже мы будем отождествлять (кусочно-постоянную) функцию Лп+1 = Лп+1(х) и вектор (лП1_+1,...,лП+1) , так как они однозначно определяют друг друга.

Выпишем коэффициенты матрицы ЛN и правой части В в случае изотропного неоднородного тела. Они вычисляются по следующим формулам:

«у= Е е2 • Е^ )Нц ; (11)

/:П/ cQ;

in I m)е*4|0°1 V-,

H ! I ”У| ) I ,

Иц= I sin I—-I e ’ a ау^ау2 ауз;

П ^ a '

b,= 1 ' <l2)

Значение интеграла по параллелепипеду П/ можно вычислить аналитически:

xl 0 +hl/2 У10 +h2/2 zl0 +h3/2 / л izt I(k0,) I2 —

HU s I I I sin1 a\ V|0 1 a2 dxdydz =

xl0 -V2 У/0 -h2/2 zl0 -h3/2 a

4И2а 0

е

а

\

Далее проверяется выполнение неравенств ^+1 _ Ли)

<8 (/ = 1,..., Ы)

с заданной точностью 8(> 0) . Если требуемая точность достигнута для каж-

ются с формулы (6).

В качестве искомого выбирается значение для относительной диэлек-

Ключевым моментом в данном двухслойном итерационном процессе

то необходимо иметь по крайней мере результаты N различных измерений.

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

Решение уравнения (5) подробно описано в [10].

3. Теорема о сходимости итерационного метода

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

Запишем итерационный метод в операторной форме:

дого ^+1 ( = 1,...,N), то вычисления прекращаются. Если требуемая точность не достигнута, то ^п+1(х):=^~+1 (х), п := п + 1, и вычисления повторя-

трической проницаемости

является возможность определения £п+1 (х) по известному полю Еп (х) из формулы (7). Если искомая функция е(х) имеет N неизвестных параметров,

( + 5(і) ) ) = Е(0,і) (/ = 1,..., N);

Е£) = ^и Л*') (і = 1,..., N);

(13)

(14)

(15)

АЫЛп+1 = В ; ^п+1 = Ли+1,

где I - единичный оператор,

(16)

£1 = к01О (х, у ) (у )сіу + graddiv | О (х, у ) (у )сіу

(17)

Є

Є

линейный ограниченный оператор ^: ^((З) ^ -^(б) [4]; ^ ^(г) для

ко = 4°. Уравнения (13), (14) рассматриваются в пространстве ^(З.

Объединив эти уравнения, будем иметь

Е« = ^ (I + ^(0 ) Е(0,г) (/ = 1,...,N); (18)

Лп+1 = АкВ , предполагая, что существует АN = А^(лп). Тогда

Лп+1 = А^ ((I + Лп^Г1 Е0)В, (19)

или

Лп+1 = Р(Лп)В, (20)

где функция

Р(Лп) - АN ((I + Лп^)1 Е0) (21)

зависит от Лп .

Функцию (21) можно рассматривать как матрицу-функцию N пере-

Т

менных лп>),. .,Ли^). Пусть Лп ^Ли"*,-.,Лп^)) е ) - вектор из простран-

^ ~ II II

ства с нормой • ^ .

Пусть существует Р(ле) для некоторого Ле . Будем предполагать, что

Ле = Р(Ле)Ве . (22)

Тогда в силу непрерывности функции (I + л<5' ) 1Е0 по Л в некоторой

окрестности точки ле существует функция Р(л) = АN ((I + Л^) 1Е0). Выбе-

рем в этой окрестности замкнутый шар Вг (ле) := {л : ||л - Ле || — г} (уменьшая г , если необходимо). Рассмотрим оценки:

||р(Л)В - Ле|| = IР(Л)В - Р(Ле)Ве\\ = 11(Р(Л) - Р(Ле))В + Р(Ле)(В - Ве)|| —

— ||Р(Л) - Р (Ле )||||В1и + I |Р (Ле )||||В - Ве\1 . (22)

Отображение (20) действует в из шара в шар Р(•)В : Вг (ле) ^ Вг (ле) при условии

тахЛеВг(ле)11Р(Л) - Р(Ле|+||Р(Ле)\\\\В - Ве 1— г . (23)

Очевидно, что максимум в первом слагаемом в (23) существует и достигается на некотором элементе из шара, так как пространство конечномерное, а шар замкнутый. Оценка (23) будет выполняться при достаточно малом тахлев ( ) ||Р(Л) - Р(Ле )|| (чего можно добиться уменьшением г ) и

достаточно малом ||В - Ве 1.

Нетрудно видеть, что функция Р(л) = Ад/ ((I + л^) 1Е0) будет также

сильно дифференцируема (по Фреше) как суперпозиция сильно дифференцируемых функций, поэтому справедлива оценка [11]

||Р (л) - Р (Т)|| — Щ |л-^и (24)

для некоторого М(> 0) при л, те Вг (ле). Тогда отображение Р(•)В:

Вг (Ле) ^ Вг (ле) будет сжатием, если дополнительно выполнено условие

И»<М-1. (25)

Итак, доказано, что при перечисленных выше условиях отображение будет сжимающим, поэтому (см. [12]) верна

Теорема 1. Пусть существует Р(ле) для некоторого Ле и верно (22). Тогда найдется такое г > 0, что при выполнении условий (23) и (25) отображение Р(•)В: Вг (ле )^ Вг (ле) является сжимающим, итерационный

процесс (20) сходится к точному (единственному) решению Ле Вг (ле) уравнения

Л = Р (Л) В (26)

со скоростью геометрической прогрессии с показателем д := И» М (< •)

при любом начальном приближении Л0 е Вг (ле).

Эта теорема теоретически обосновывает двухслойный итерационный метод определения функции е(х). Наиболее сложным является обеспечение

условия (22). Оно означает, что должно быть известно решение обратной задачи с «близкими» параметрами. В качестве такой задачи можно выбрать, например, задачу определения эффективной диэлектрической проницаемости, когда Л(1) = .. = Л(N) =Ле. Эта задача подробно исследована в [4, 6, 7] (где также указаны условия существования ее решения).

4. Численные результаты

Описание решения интегродифференциального уравнения методом коллокации имеется в [10]. В качестве точек коллокации выбираем центры элементарных параллелепипедов. Параметры задачи: а = 2, Ь = 1, с = 2,

к0 = 2,5, N = 1, N0 = 8 . Коэффициент Р вычисляется с помощью аналитического решения прямой задачи дифракции [13, 14].

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

Точное значение равнялось е = 1,5 . В левом столбце табл. 1 указаны начальные значения £е; в правом столбце указаны вычисленные (приближенные) значения двухслойным итерационным методом.

В табл. 2 представлены результаты расчетов для случая неоднородного тела, состоящего из двух секций:

^1 = {х: 0 < Х1 < а, 0 < Х2 < Ь, 0 < Х3 < С1} ,

Q2 = {х: 0 < Х1 < а, 0 < Х2 < Ь, С1 < Х3 < с} .

Таблица 1

Начальное Вычисленное Начальное Вычисленное

значение є значение є значение є значение є

1,100000 1,507570 2,000000 1,507060

1,200000 1,507500 2,100000 1,507020

1,300000 1,507430 2,200000 1,506980

1,400000 1,507360 2,300000 1,506950

1,500000 1,507300 2,400000 1,506920

1,600000 1,507250 2,500000 1,506920

1,700000 1,507200 2,600000 1,507090

1,800000 1,507150 2,700000 1,507710

1,900000 1,507110

Параметры задачи: а = 2, Ь = 1, с = 2, С1 = 1, к0^ = 1,7, к|(2) = 1,6, N = 2, N0 = 10. Начальное значение относительной диэлектрической проницаемости равнялось ее = 1,15 в каждой секции. Точные значения равнялись е(1) = 1,1 в первой секции и е(2) = 1,2 во второй секции.

Таблица 2

Количество итераций Значение є

в первой секции во второй секции

1 1,1442211 - і0,0018885104 1,1547399 - і0,00087996185

2 1,1440666 + і0,00031775959 1,1574058 - і0,004368522

3 1,140544 + і0,0077741992 1,1624472 - і0,0096281959

4 1,1316639 + і0,011875111 1,1695469 - і0,012195842

5 1,1253693 + і0,010341899 1,1745528 - і0,012229604

6 1,1238371 + і0,010211105 1,1776618 - і0,013340704

7 1,1209914 + і0,013582285 1,181974 - і0,015370148

8 1,114805 + і0,015177354 1,1873796 - і0,015540706

9 1,1098722 + і0,013085815 1,1913597 - і0,013888457

10 1,1080044 + і0,011108461 1,1937901 - і0,01250952

11 1,1062838 + і0,010900297 1,1962147 - і0,011620387

12 1,1033408 + і0,010241321 1,1987555 - і0,010081171

13 1,1009931 + і0,0080135619 1,2004106 - і0,0078635895

14 1,1002929 + і0,0057195404 1,2010552 - і0,005890436

15 1,1001548 + і0,0044050605 1,2013838 - і0,0044480784

16 1,0996759 + і0,0034048383 1,2016867 - і0,0031334373

17 1,0992887 + і0,0020517684 1,2017066 - і0,0017845466

18 1,099517 + і0,00071302208 1,2013461 - і0,00066606828

19 1,1000676 - і0,00010505475 1,2008618 + і6,8213645е-005

Расчеты показывают высокую точность (порядка 0,4 %) при определении относительной диэлектрической проницаемости образца. Метод быстро сходится даже при выборе сильно отличающегося начального приближения от точного значения.

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

1. Solymar, L. Waves in Metamaterials / L. Solymar, E. Shamonina. - New York : Oxford University Press, 2009.

2. Zharova, N. A. Nonlinear Transmission and Spatiotemporal Solutions in Metamaterials with Negative Refraction / N. A. Zharova, I. V. Shadrivov, A. A. Zharov, Yu. S. Kivshar // Optics Express. - 2005. - V. 13, № 4. - P. 1291-1298.

3. Shestopalov, Yu. V. Development of Mathematical Methods for Reconstructing Complex Permittivity of a Scatterer in a Waveguide / Yu. V. Shestopalov, Yu. G. Smirnov, V. V. Yakovlev // Proceedings of 5th International Workshop on Electromagnetic Wave Scattering (Antalya, Turkey, October 22-25). - Antalya, 2008.

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

5. Kobayashi, K. Investigation of Electromagnetic Diffraction by a Dielectric Body in a Waveguide Using the Method of Volume Singular Integral Equation / K. Kobayashi, Yu. V. Shestopalov, Yu. G. Smirnov // SIAM Journal of Applied Mathematics. - 2009. -V. 70, № 3. - Р. 969-983.

6. Смирнов, Ю. Г. О существовании и единственности решений обратной краевой задачи для определения эффективной диэлектрической проницаемости наноматериалов / Ю. Г. Смирнов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2009. - № 1. - С. 11-24.

7. Смирнов, Ю. Г. О существовании и единственности решений обратной краевой задачи определения диэлектрической проницаемости материалов / Ю. Г. Смирнов, Д. А. Миронов // Журнал вычислительной математики и математической физики. -2010. - Т. 50, № 9. - С. 1587-1597.

8. Smirnov, Yu. G. Existence and uniqueness of a solution to the inverse problem of the complex permittivity reconstruction of a dielectric body in a waveguide / Yu. G. Smirnov, Yu. V. Shestopalov // Inverse Problems. - 2010. - V. 26, № 105002. -P. 1-14.

9. Smirnov, Y. Analysis of Inverse Scattering in a Waveguide using the Method of Volume Singular Integral Equation / Y. Smirnov, Y. Shestopalov, D. Mironov // URSI Internetional Symposium on Electromagnetic Theory (EMTS-2010) (Berlin, Germany, August 16-19). - Berlin, 2010. - P. 532-534.

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

11. Колмогоров, А. Н. Элементы теории функций и функционального анализа / А. Н. Колмогоров, С. В. Фомин. - М. : Наука, 1976.

12. Канторович, Л. В. Функциональный анализ / Л. В. Канторович, Г. П. Акилов. -М. : Наука, 1984.

13. Гурина, Е. Е. Численное и аналитическое решение задачи дифракции электромагнитного поля на диэлектрическом параллелепипеде, расположенном в прямоугольном волноводе / Е. Е. Гурина, М. Ю. Медведик, Ю. Г. Смирнов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2010. - № 2. - С. 32-43.

14. Гришина, Е. Е. Численное и аналитическое решение задачи дифракции электромагнитного поля на двух секциях с разной диэлектрической проницаемостью, расположенных в прямоугольном волноводе / Е. Е. Гришина, Е. Д. Деревянчук,

М. Ю. Медведик, Ю. Г. Смирнов // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2010. - № 4. - С. 73-81.

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

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

E-mail: [email protected]

Васюнин Денис Игоревич аспирант, Пензенский государственный университет

E-mail: [email protected]

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

Vasyunin Denis Igorevich Postgraduate student, Penza State University

УДК 517.3 Смирнов, Ю. Г.

Итерационный метод определения диэлектрической проницаемости неоднородного образца материала / Ю. Г. Смирнов, Д. И. Васюнин // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2011. - № 1 (17). - С. 20-30.