CC BY
12
УДК 517.956.32
И. В. Матюшкин1'2, Г. Я. Красников1'2, Н.В. Черняев1, Е. С. Горнев1'2,
Н. В. Евстратов1'2
ХАО «НИИ молекулярной электроники» 2Московский физико-технический институт (государственный университет)
Численное исследование параметрических свойств решений некоторого уравнения типа Штурма—Лиувилля
Объектом исследования было частное уравнение вида Штурма-Лиувилля, полученное ранее в результате решения уравнений Максвелла для плоской электромагнитной волны, распространяющейся в нестационарной среде. Для случаев линеаризованной, линейной, экспоненциальной и гармонической нестационарности исследовано влияние параметров нестационарности на гармонические и концевые свойства уравнения. Предложенная методика выделения амплитуды и фазы на правом конце отрезка позволила обнаружить явления квазипериодичности и детерминированного хаоса, названные нами «эффектом арки».
Ключевые слова: уравнения Максвелла, уравнение Штурма-Лиувилля, диэлектрическая проницаемость.
I. V. Matyushkin 1'2, G. Ya. Krasnikov1'2, N. V. Chernyaev 1, E.S. Gornev l'2,
N. V. Evstratov 1'2
1 Molecular Electronics Research Institute JSC 2Moscow Institute of Physics and Technology (State University)
Numerical Investigation of Parametric Properties of Particular Sturm-Liouville Equation Solution
A particular kind of the Sturm-Liouville equation, which is earlier derived from the solvation of Maxwell's equations for electromagnetic waves propagating in nonstationary medium is studied. The influence of nonstationary parameters on harmonic and end equation properties is studied in linearized linear exponential and harmonic cases. The proposed technique of magnitude and phase extraction at the right end of segment allows us to find the quasiperiodicity and deterministic chaos phenomenon, which we call «arch effect».
Key words: Maxwell's equations, Sturm-Liouville equation, permittivity.
1. Введение
Задача описания распространения электромагнитной волны (ЭМВ) в среде с изменяющимися материальными параметрами (диэлектрической и магнитной проницаемостями) является актуальной в связи с развитием устройств современной фотоники. Например, в работе [1] рассматривался случай константной магнитной проницаемости и линейной (и экспоненциальной) диэлектрической проницаемости; в формализме напряженностей поля < Е, В > найденное решение выражалось через функции Бесселя. При этом авторы даже учли протекание в среде токов.
Ранее в [2] мы вывели для плоской волны, распространяющейся в нестационарной среде, систему дифференциальных уравнений (в формализме потенциалов), в общем виде связывающую электромагнитные потенциалы с управлениями, т.е. парой динамически меняющихся проницаемостей среды е = е(Ь),^ = Системы уравнений Максвелла мы
свели к нескольким однотипным телеграфным уравнениям, а затем методом разделения переменных получили обыкновенные дифференциальные уравнения, которые формально классифицируются как частный случай уравнения Штурма-Лиувилля:
((2 X
-+ = ш2х, г е [0;т], X е [0; Ь].
Однако вопросы ставятся не для краевой задачи Штурма-Лиувилля, а для задачи Коши. В большей степени нас интересуют не глобальные свойства решения, такие как количество нулей на отрезке, а локальные свойства на правом конце отрезка. Такая постановка задачи актуальна для приборов фотоники [3], использующих интерференцию волн в качестве механизма обработки информации.
Там же в [2] мы представили методику нахождения локальных амплитуды, фазы и частоты ЭМВ. Таким образом, каждой дифференциальной задаче Штурма-Лиувилля, выделенной конкретными значениями параметров нестационарности, можно поставить в соответствие тройку введенных нами величин и исследовать уже поведение данной функции.
Хотя и на качественном уровне, нами будет проведен анализ зависимости фурье-образа решения дифференциальной задачи от способа и величины возмущений е = е(Ь),^ = Для ряда численных результатов указана непосредственная связь с теоретическими результатами, полученными в [2]. Преимущественный метод исследования - численный, но в конце статьи мы сделаем попытку аналитически объяснить ряд эффектов.
2. Общая постановка задачи
Целью работы является параметрическое исследование задачи Коши для уравнений вида [4] Штурма-Лиувилля:
T'j« + (1? - ^ + f) T"w = °- " í + f) T0<" = °- (1>
Индекс 0 относится к скалярному потенциалу, а1 — к векторному (j = 1, 2). Считаем, что ш = kc/^/e(0)^(0), c - скорость света, которую принимаем за 1, поскольку универсальные константы eo, Цо тоже приняты за 1.
Будем рассматривать три модели, которые наиболее естественно бы моделировали нестационарности среды - линейную (а), экспоненциальную (b) и гармоническую (с):
(а) : e(t) = e(0)(1 + at), n(t) = ц(0)(1 + bt), (b) : e(t) = e(0)exp(-at),^(t) = ц(0) exp(-bt),
(c) : e(t) = e(0) + а sinp^ , ¡i(t) = ц(0) + - sinp^ ,...; p ^ ш.
Проведем процедуру обезразмеривания. Нормируем время на период невозмущенной
волны t = —. Примем без ограничений общности с учетом вида (1) e(0)^(0) за единицу, 2п
т.е. нормируем проницаемости на их начальные значения. Соответственно коэффициенты a,b,p будут иметь размерность обратных секунд. Тогда уравнения (1) примут следующий вид, не содержащий ш:
e2T"+Шde)2 - ¿Seтцй=0
dt2 \4e2 \ dt) 2e dt2 ец.
d2To + ( * (дЛ\2--1+ Ш Tn(t) = 0, (2)
dt2 + - ~2ц~д¥ ^^)To(t) = 0,
e(0) = 1, ц(0) = 1,í e [0; N].
Без ограничения общности примем (^ = 1 в (2)) граничные условия на левом конце равными (3):
ТИ(0)=0,Т11(0) = 1,
To(0) = 0,T0 (0) = 1. ()
Таким образом, получаем математически рафинированную формулировку задачи (2) - (3). В дальнейшем проницаемости интересуют нас как безразмерные коэффициенты, зависящие от некого параметра. Для удобства и дальше будем именовать их проницаемо-стями, которые зависят от времени.
Достаточно решить лишь одно уравнение, для определенности выберем первое из (2), содержащее производные е по t и, соответственно, граничное условие из первого соотношения в (3). В отсутствие нестационарности его решением будет T(t) = — sin(2^t).
2п
3. Используемые численные методы
Для дифференциальных уравнений используем решатель ode45 в MATLAB. Он является одношаговым явным методом Рунге-Кутты 4-го порядка [5]. С помощью параметра odeset устанавливаем допуски ошибок (10-5) [6]. На промежутке времени шаг сетки составляет 0.001, интервал изменения независимой переменной - от 0 до N (Приложение, line 1-3).
Сравним в случае невозмущенной волны решение и производную, выдаваемые ode45, в начальный и конечный момент времени: Tend — To = 0.0002; T'end — TO = -0.0013.
Итак, за 100 периодов не произошло существенного накопления погрешности, что говорит о пригодности данного численного метода для нашей задачи. Общее число точек по оси времени - Mt = 105.
Для разложения решения уравнения T(t) в ряд Фурье используем функцию fft (Fast Fourier Transform, FFT) в MATLAB. Посредством этой функции производится преобразование (в нашем случае прямое) Фурье для одномерного массива длиной M с применением (4):
M
F(k) = J2 T(j)wM-1)(fc-1), k = 1: M, um = exp(-2ni)/M. (4)
В результате преобразования Фурье получают [7] комплексные значения, дающие результаты для амплитуды и фазы. На оси частот это выражается тем, что вместо одного пика для тестовой гармоники появляются два симметричных, отвечающих за отрицательные и положительные значения частот, однако нам достаточно рассмотрения одного из них. В MATLAB для FFT используется алгоритм Кули-Тьюки, согласно которому преобразование Фурье от M = M1M2 рекурсивно упрощается до преобразований более малых размерностей M1 и M2. В простейшем случае происходит деление надвое, поэтому M должно являться степенью двойки, а преобразование будет происходить особенно быстро [5,8]. В нашем случае это M = 210 (количество точек по оси частот и одновременно по оси времени). Это значение должно коррелировать с утверждением теоремы Котельникова о том. что частота сигнала должна быть в два или больше раз меньше частоты дискретизации
Fs = Ml = юз.
N
Так как после решения задачи Коши (2) - (3) имеем Mt точек, то интерполируем решение уравнения на M точек (Приложение, line 4). Для удобства и наглядности смещаем нулевую частоту в центр спектра с помощью функции fftshift [6] (Приложение, line 5).
Калибруем ось k так, чтобы максимальное значение было равно половине частоты дискре-
M M M тизации, для этого вычитаем — и получаем значения на оси от--до — (далее рассматриваем только положительные значения на оси), после чего растягиваем ось и получаем откалиброванную ось частот П (Приложение, line 6) с тем, чтобы для невозмущенной волны ~ sin(2^Qi) спектральная плотность была сосредоточена на частоте П = 1. В нашем случае П £ [0; Qmax = 5.12], а точки расположены на этом интервале однородно.
В дальнейшем будем графически отображать зависимость спектральной плотности W(П). Также будем анализировать зависимости амплитуды P и фазы Q в конечной точке от коэффициентов, отвечающих за нестационарность (a,b,p). Поиск проводился по следующей схеме, описанной в [2].
Вводились локальные амплитуды и частоты (P, u)t - у аргумента тильду писать не будем. Ограничивающим предположением методики является многоэкстремальность кривой T(t). Набег фазы считается от точки ближайшего (слева) k-го экстремума и равен произведению локальной частоты на разницу времен, текущего и локального экстремума: Q(t) = Uk(t — tk) £ [0; п], k : (tk < t) /\(t — tk < п). Аргументы экстремумов при численном решении ищем приближенно из условия y(i)y(i + 1) < 0, где y = T'(t) - предоставляемый ode45 вектор, а i - номер узла. В конце каждого периода введем значения амплитуды
п
Pk = P(tk) = \T(tk)| и частоты Uk = -. Величины Uk близки по смыслу к частоте
tk — tk-i
следования нулей в задаче Штурма-Лиувилля. Для простоты не будем маркировать последнюю точку индексами, относящимися к паре (P,Q)t=tk:k=max(k). Такой анализ назовем «концевым».
Опишем также методику анализа результата преобразования Фурье, т.е. спектральной плотности W(П), которая во всех рассмотренных нами случаях имела форму волнового пакета (ВП). Предположим, что при разложении Фурье П £ [0; Птах] и рассмотрим первообразную V(П) = JQw<Wmax W(n)dn. Введем для некоторого 0 ^ z ^ 1 определение Пх = П : (V(П)/V(Птах) = z). Тогда введем три характеристики ВП: к1 = П1/2 — 1 - цен-
Q Wmax Wmin I TIr ■ ттт/гл\ 1
, , „ - 9/10 — П1/10 - ширину, КЗ = —-——— Wmax,min = max, min W(П)
"max + Wmin ^ П1/10^П^П9/10 )
- видность (см. «видность» в оптике [5]). Поскольку исходными данными являются конечномерные вектора, то при численной реализации в MATLAB интегрирование заменяем простым суммированием, а равенство - двумя условиями <,> с последующей линейной интерполяцией.
Сводные результаты анализа фурье-разложения по описанной методике, исключая гармонический случай, приведены в табл. 1.
Примечание. На отрезок [0;1] на частоту П приходится 100 точек, т. е. числовая погрешность W(П) не менее 0.01. В строках П приведена номинальная частота в конце отрезка, т.е.
1 А2
— I T" + П(^2Т = 0, П2 = П^ = N). В строках к2 приведены соответственно ширина 2п
ВП, граничные частоты П1/10, П9/10 и число точек, вошедших в пределы ВП (указываем, как точно нами охарактеризован ВП).
4. Исследование линеаризованной модели
В первую очередь целесообразно исследовать линеаризованную модель нестационарности. Она является общей, поскольку при разложении в ряд Тейлора проницаемостей на малых временах все частные случаи (1) сводятся к (5):
( Ш 2 Ж + (1 + ^« = ° (5)
I т(0) = 0, т'(0) = 1, И < N.
Т а б л и ц а 1
Зависимость характеристик ВП от безразмерных параметров
2па 2па
нестационарности а =-,р =-
ш ш
Параметр а
0.0001 0.0063 0.0125 0.0376
Линеаризованный а 1 1.28 1.5 2.19
0.0023 0.1404 0.2534 0.5937
случай 0.0182 0.2135 0.3896 0.9447
К2 (1.0060, (1.0360, (1.0660, (1.1261,
1.0160), 1.2562), 1.4564), 2.0770),
2 точки 22 точки 39 точек 95 точек
КЗ 0.5490 0.1990 0.2728 0.2492
Параметр а 0.0001 0.0063 0.0125 0.0376
Линейный а 1 0.613 0.44 0.21
К1 -0.0012 -0.2578 -0.4230 -0.7088
а = в 0.0203 0.2890 0.3720 0.3391
К2 (0.9960, (0.6456, (0.4755, (0.2252,
1.0160), 0.9359), 0.8458), 0.5656),
2 точки 29 точек 38 точек 34 точки
КЗ 0.8080 0.6074 0.7760 0.8842
Параметр а 0.0001 0.0032 0.0063 0.0099
Линейный а 1 1.056 1.288 7.088
К1 0.0010 0.0157 0.0494 0.0951
а = —в 0.0083 0.0415 0.1989 0.6800
К2 (1.0060, (1.0060, (1.0060, (1.0160,
1.0060), 1.0460), 1.2062), 1.6967),
1 точка 5 точек 20 точек 68 точек
КЗ 0 0.5303 0.8041 0.9616
Параметр а 0.0001 0.0063 0.0125 0.02
Экспоненциальный а 1 1.88 3.49 7.4
К1 0.0087 0.3070 0.5574 0.7622
а=в 0.0203 0.6700 1.7127 2.9292
К2 (1.0060, (1.0660, (1.0961, (1.1161,
1.0260), 1.7367), 2.8077), 4.0490),
2 точки 67 точек 171 точка 292 точки
КЗ 0.8067 0.6065 0.8039 0.9829
Параметр а 0.0001 0.0063 0.0125 0.02
Экспоненциальный а 1 0.53 0.286 0.135
К1 0.0009 -0.1538 -0.3029 -0.4605
в = —2а 0.0088 0.2029 0.3531 0.4487
К2 (1.0060, (0.7657, (0.5756, (0.4054,
1.0060), 0.9659), 0.9259), 0.8558),
1 точка 20 точек 35 точек 45 точек
КЗ 0 0.4847 0.5745 0.7313
Мы будем исследовать (рис. 1) зависимости амплитуды и фазы конечной точки от коэффициента нестационарности а, а также фурье-разложение решения (5). Известно аналитическое решение (7) через функции Эйри [2].
Анализ табл. 1 при а = 0.0001 указывает нам на погрешность численного метода по оси частот - ~ 0.01 (критерий к2). ВП локализуется справа от центральной частоты невозму-
щенной волны. Наблюдается симметрия ВП относительно центроида. При малых значениях коэффициента нестационарности амплитуда максимального пика, т.е. тахШ(а), резко падает, затем ее изменения незначительные (рис. 1). Номинальная частота, положение центроида и ширина пика практически линейно пропорциональны а (см. табл. 1). Видность ВП слабо зависит от а и равна примерно 1/4. На абрисе спектральной плотности появляются биения, размах и число которых зависят от а. Экспериментальных данных недостаточно по точности и количеству, чтобы экстрагировать аналитическую зависимость огибающей абриса ВП; интерполяция по параболической гипотезе оказалась неэффективной.
Рис. 1. Зависимость спектральной плотности W от частоты Q для уравнения (5). Линеаризованный случай
Отметим превосходное совпадение частотных кривых, полученных по предложенной методике и непосредственно из уравнения (5). С помощью MATLAB Curve Fitting Toolbox интерполировали зависимость амплитуды от коэффициента нестационарности при следующей гипотезе: P(а) = m(1 + aN)-n, где m = ^ = 0.1591 (0.1591^0.1592), n = 0.25 (0.2499^0.2500). В скобках указаны доверительные интервалы для стандартного уровня ошибки 95. Такой результат является косвенным аргументом в пользу корректности предложенной нами методики экстракции параметров (см. п. 1), поскольку аналитическое решение задачи известно. В [2] была получена обратнопропорциональная зависимость квадрата амплитуды P ЭМВ от частоты Q при различных видах нестационарности, что и подтверждают рис. 2 и приведенная выше степенная зависимость. Действительно, если предположить, что номинальная частота в (5) приблизительно равна производной фазы по времени, то
\ 2 d
p (t)2 - (tm)
P (t)2Q(t) = const.
(6)
откуда искомое соотношение: Р =
const
const
(1 + aN)-1.
Ц (2п)2(1 + аЫ) лДЛ Соотношение (6) выполняется точно для мгновенных величин и приближенно для дискретных отсчетов, получаемых при «концевом» анализе.
Рис. 2. Зависимость экстрагированных из решения (5) набега фазы Q, амплитуды Р и частоты (все в конечной точке) от параметра нестационарности. Построена кривая для номинальной
частоты = а/1 + аМ в конечной точке (М = 100). Амплитуда волны нормирована на —. Набег
2п
фазы нормировался на 2п
На абрисе набега фазы наблюдаем феномен арки, когда набег фазы не достигает своих крайних значений, а огибающая имеет форму арки. Такое встречается дважды: при совсем малых а < 2 • 10_4 и на отрезке а = 0.02 ± 0.015. Инженерный смысл явления «арки», по-видимому, состоит в том, что именно в этих диапазонах нестационарности следует проводить измерения по деградации среды или налагать информационную составляющую сигнала в фотонных приборах.
5. Исследование линейной модели нестационарности
В случае линейной нестационарности имеем для Т\] (для То аналогично с перестановкой а и Ь, индексы для простоты писать не будем, Ь размерно):
e(t) = е(0)(1 + at), ¡i(t) = ^(0)(1 + bt),
T"(t) +
ш
+
T (t) = 0.
(7)
(1 + аЬ)(1 + Ы) 4(1 + аг)2/
При а = Ь дифференциальное уравнение имеет аналитическое решение вида ~ cos(ln Ь). Хотя вывод для общего случая приводился в первой части статьи, уместно привести более простой вывод для данной ситуации.
Будем искать решение уравнения Т"(Ь) +
ш + Т
—-т(t) = 0 как (1 + at)А, что даст нам
(1 I at)
1
квадратичное условие на Л, причем корни будут мнимые с ИеЛ = ^:
ш
4
гш
X - X + ^7 + Т =0, Ai.2 = - ± —.
2
(8
Суперпозиция двух комплекснозначных решений даст нам общее вещественное реше-
ние:
T(t) = CV1 + at cos (a ln(1 + at) + d) .
(9)
2
a
a
a
Здесь С, О - константы, которые можно найти из граничных условий при Ь = 0. В пределе а ^ 0 (9) естественным образом переходит в закон косинуса. Следует отметить, что даже в предположении ш ^ а,аЬ ^ 1 набег фазы в (9) может быть существенным на высоких частотах, что вполне соответствует выводу в [2]. Этот набег, если разложить логарифм в ряд, приблизительно равен (¿уШа)2. Здесь набег фазы отсчитывается от момента Ь = 0.
При а = Ь точного аналитического решения (7) найти не удается, но можно заменить уравнение приближенным:
ш2
+ (ТТОнТ+Ь)Т (() = (10)
При Ь = —а решение (7) вида (9) возможно, как показывает анализ (10), лишь при ш2 = 2а2, что выводит за пределы исходного предположения ш ^ а, и формально записывается как Т(¿) = 1 — а2Ь2.
Для численного исследования перепишем уравнение (7) в безразмерных величинах:
д2Т ( 2п2 а2 \ 2па 2пЬ
+ „ , , ^ + ^ , ЪЛ Т(Ъ) = 0, а =-,в =-. (11)
дЪ2 \(1 + аЪ)(1 + в- 4(1 + аЪ)2/ ' ш
Предположим, что а = в.
В отличие от линеаризованного случая ВП локализуется слева от единичной частоты, что и следует из уравнения (11), амплитуда наиболее выраженной гармоники при малых а падает медленнее, чем в линеаризованном случае, и растет при больших а. Появляющиеся на абрисе биения асимметричны (сдвинутые влево гармоники ВП вносят больший вклад, чем сдвинутые вправо), при больших значениях а сглаживаются. При возрастании а ширина пика уменьшается, однако появляется хвост ВП. При росте а огибающая пика приближается к оси абсцисс, появляется хвост. При экстремально больших значениях а ^ 1 мы получаем единичный пик при частоте, близкой к нулю, что можно трактовать как трансформацию ВП в одну инфракрасную моду. Видность растет (рис. 2) с ростом а, а вот номинальная ширина ВП меняется незначительно и противонаправленно.
Кривые двух частот а1 и а2 совпадают. Кривые Р(а), а2(а) коррелируют друг с другом, а именно соблюдается обратная зависимость квадрата амплитуды и частоты. Так же, как и в линеаризованном случае, при небольших значениях а набег колебания набега фазы происходят чаще. Арка в данном случае появляется при меньших значениях а, также можно увидеть еще одну слабовыраженную арку. Вогнуты арки в противоположную сторону (по сравнению с линеаризованным случаем). Обращает на себя внимание эффект обращения (зубцов) на кривой ^(а) в области больших а, когда происходит быстрый набор фазы и более медленное его падение. При этом плотность зубцов падает. Также отметим эффект дрожания амплитуды, видимый на рис. 4 на правом конце кривой Р(а), - его мы связываем особенностями метода экстракции, не исключая фактора приближенного выполнения (6).
1 1( а2
На рис. 3 построена кривая для номинальной частоты а2 = --тт—+ 1
(1 + (Ш) у \4(2п)2
в конечной точке N = 100. Амплитуда волны нормирована на —. Набег фазы в начальной
2п
точке, как и ожидается, равен п/2 и нормировался для удобства на 2п.
Обсудим линейный случай при в = —а.
Рассматривались значения а < 0.01, так как при остальных исходное уравнение не имеет смысла (проницаемость принимает отрицательные значения). Амплитуда наибольшей гармоники монотонно падает с ростом а (в отличие от синфазного изменения проницае-мостей). ВП локализуется справа от единицы. По сравнению со случаем а = в количество/плотность биений заметно меньше, в остальном абрис тот же - хвост по-прежнему сохраняется. В отличие от предыдущего случая ширина пика и видность растут при увеличении а, причем видность при максимальном значении близка к 1.
х10л(-3)
х10л(-4) В
7 6 5
з г 1
о.
- 1 х10л(-4)
16
14
- 12 -
- «=0.0125 10 (1=0.0376
- 3 8 6 ■ 1
\ 4 ] \
- V 1 V- 2 \
1 | 0.5 1 1.3 I I 1
06 п 1
15
2.5
35
Рис. 3. Зависимость спектральной плотности Ш от частоты О для уравнения (11). Линейный случай при а = в > 0
Что касается «концевого» анализа, то амплитуда падает, вначале медленно и затем более круто (при а > 0.01). Так же, как и при синфазном управлении, наблюдается начальная полуарка на графике набега фаз. Вторая арка довольно узкая и относится к форме абриса, не являясь огибающей, как в синфазном случае. Почти сразу за ней идет последняя, третья и весьма глубокая, т.к. падает почти до нуля, арка.
6. Экспоненциальная модель нестационарности
е(г) = е(0) ехр(-аг), /(г) = /(0) ехр(—Ы),
( а2\ (12) т"(г) + ш2 ехр(а + ъ)г - ^ т(г) =0.
4 ,
Считаем а = Ъ, и тогда в прежних обозначениях имеем:
г)2т ( а2 \
д-р + (^(2п)2 ехр2аг - ^ Т© = 0. (13)
Амплитуда максимального пика довольно быстро падает вследствие расширения ВП, локализованного справа от единицы (что подтверждается и центроидом, см. табл. 1). До некоторого значения частоты гармоники амплитуда биений падает (при увеличении их плотности), а после, наоборот, растет при уменьшении их плотности вплоть до последнего «росчерка». Видность и ширина ВП умеренно растут с увеличением а.
Рис. 4. Зависимость экстрагированных из решения (11) набега фазы Q, амплитуды Р и частоты а1 (в конечной точке) от параметра нестационарности
Рис. 5. Зависимость спектральной плотности Ш от частоты а для уравнения (11). Линейный случай при в = —а
С учетом сделанного выше замечания формально при коэффициенте нестационарности, близком к а > 0.02, начинаются биения спектральной плотности с элементами хаоса. Это отличает случай от других, так как некоторые высокочастотные гармоники могут быть «выключены» из имеющегося спектра.
При «концевом» анализе (рис. 7), аналогично предыдущим случаям, кривые частот совпадают. Арки появляются чаще (их 3, не считая начальной), выражены наиболее четко, т.к. в отличие от предшествующего, что принципиально, уже не являются огибающими линиями. Наиболее ясно выражена арка с минимумом а = 1/N. Поведение набега фазы обнаруживает черты самоподобия (фрактальности) в «межарочных» интервалах; нахождение здесь теоретических зависимостей представляет собой интересную задачу для дальнейших исследований. А с точки зрения практики наличие особых диапазонов значения коэффициента нестационарности, где набег фазы изменяется медленно, интересен возможностью кодирования информации.
х10Л(-3)
а
Рис. 6. Зависимость спектральной плотности Ш от частоты а для уравнения (13). Экспоненциальный случай при а = в
а2
Построена кривая для номинальной частоты а2 = у exp(2аN) — 4(2 )2. Амплитуда
волны нормирована на —. Набег фазы нормировался на 2п.
2п
Также был исследован случай Ь = —2а. Центроид при этом сдвигался влево от единицы более медленно, ср. с рис. 6. Ширина ВП была меньше, чем при синфазном управлении. Существует такая же точка перегиба, но она менее выражена, правый край так же заканчивается резким «росчерком».
Обратим внимание, что существует а = а* ~ 0.09723 (по-прежнему в = —2а ), см. (12) и (13), при котором номинальная частота обращается в ноль на конце интервала. «Концевой» анализ показывает небольшое расхождение частот, т.е. а1 = а2 ^ а1 ~ а2, более заметное при больших а; визуально частота падает линейно, хотя реализуется приближенно экспоненциальный закон (квадратичным слагаемым в (13) можно пренебречь) в области малых
аргументов. Феномен арки возвращается к абрису огибающей; первая полуарка приходится на начальный участок а < 0.002, при этом набег фазы (нормированный) меняется от стартовых 0.25 до почти максимального уровня в 0.5, а вторая наблюдается в промежутке а ~ 0.011^0.014. Эти диапазоны не слишком сильно отличаются от предыдущего случая, то уже третья арка не наблюдается. Плотность зубцов падает при увеличении а.
□ 0.002 0 004 0.008 0.03 0.01 0.012 0.014 0.016 0.018 0.02
а
Рис. 7. Зависимость экстрагированных из решения (13) набега фазы Q, амплитуды Р и частоты (в конечной точке) от параметра нестационарности
«Концевой» анализ проводился также для следующего участка 0.02 < а ^ 0.064 < а ^ 0.097. Амплитуда зубцов на кривой набега фазы становится все большей, временами превышая 0.5; одновременно с этим усиливается эффект дрожания кривой Р(а). Однако эти особенности получены в той области, где методика уже не вполне «работает», где кривая Т(¿) лишается «синусоидального» характера. Однако в этой области действителен фурье-анализ; типичный абрис состоит из резкого пика при О ~ 0 (+0.1) и небольшого хвоста справа, где практически незаметны биения.
Аналогично линеаризованному случаю интерполировалась зависимость амплитуды волны от коэффициента нестационарности. Для случая а = Ь получили Р(а) = ехр^а), где ¡1 = —50, (—50.01, —50), то есть результаты имеют высокую точность. Вместе с тем для а = —2Ь Р(а) = ехр(12а), где ¡2 = 24.82(24.79, 24.84). Результаты также соответствуют полученной обратнопропорциональной зависимости квадрата амплитуды волны от
I " О2"
частоты. Если в формуле для частоты ^ /ехр((а + в)^) — Т7—пренебречь вторым слагае-
у 4(2п)2
мым по сравнению с экспоненциальным, то получим следующую формулу для амплитуды: (а + в \
Р(а, в) = ехр I —4— N I, что удовлетворяет (6).
7. Гармоническая модель нестационарности
Существует два пути внесения гармонических нестационарностей при разных соотношениях уже пары, а не одного, коэффициентов нестационарности а, Ь, что делает анализ многовариантным. Рассматриваем сразу безразмерные уравнения (14) - (18), (20) с
2па п 2пЬ 2пр ^ .
а = -, в = -, 7 = -, где р << ш, и берем за основу пару < 7,а/^ >, где 7 -
ш ш ш а
частота модуляции проницаемостей, и — < 1 - амплитудный коэффициент, который пока-
7
зывает, насколько меняется диэлектрическая проницаемость (в относительных единицах). Для наглядности в уравнениях вынесем (2п)2 - так мы удостоверимся, что при нулевых коэффициентах нестационарности действительно получается уравнение для невозмущен-
д2Т
ной волны. Таким образом, решаемое уравнение имеет вид —+ ((2п)2 + q(a, !)) Т(-) = 0.
д-2
Напомним условие обезразмеривания е(Ь = 0) = = 0) = 1. 7.1. Синусоидальное управление, а = в
а \ ( Ь \
е(г) = е(0) ( 1 + - в'т(рг)) , /л(г) = /л(0) I 1 + - ят(рг)
q(а,i) =
О- — ( (2п)2 — (25 + 82)
(1 +
а
8 = — 8Ш(7- ).
1
(14)
Рис. 8. Зависимость спектральной плотности Ш от частоты а для уравнения (14). Гармонический случай (5.1) при а = в
При повышении частоты возбуждения 7 абрис ВП приобретает игольчатый вид, а в дальнейшем спектр ВП из непрерывного становится дискретным, что коренным образом отличает гармонический случай от предыдущих и затрудняет применение методики анализа фурье-спектра. На рис. 8б рядом с основной гармоникой появляется малоамплитудная, причем они практически симметричные. Например, при 7 = 1 наблюдается 7 гармоник и несколько малоамплитудных. Причем наиболее выраженная гармоника находится вблизи 1. Узость пиков сравнима с величиной численной ошибки по оси частот.
Физического смысла это не имеет, однако рассматриваем значение 7 = 4п, так как при
а
этом в (14) уйдут несколько членов, а при значениях — ^ 0 получится случай невозму-
7
щенной волны, что показывает рис. 8. С учетом (14) верна формула (15), показывающая
существование пар коэффициентов нестационарности, дающих нам гармоническое решение:
д2Т + ( y2 + а2 - y2 + 4(2п)2 \ 0
W + U + 4(1 + s)2 ) Т(t) = 0 (15)
А именно, при а2 — y2 + 4(2-/г)2 = 0 гармонический сигнал имеет частоту y/2, что соответствует выводам в [2]. В частности, сигнал сохраняется невозмущенным тогда и только тогда, когда гармоническая нестационарность имеет нулевую амплитуду. Мы не можем снизить частоту исходной ЭМВ, но можем её увеличить, например, в 1.5 раза с помощью высокоамплитудного (a/Y = л/б/3) возмущения; подобную ситуацию демонстрирует рис. 8в. Вместе с тем мы можем получить (рис. 8б) суперпозицию гармоник: на основную моду, имеющую частоту, большую частоты исходной ЭМВ, накладывается малоамплитудная гармоника, сдвинутая к меньшим частотам.
«Концевой» анализ проведем по двум параметрам только для амплитуды P(y, a/y). Общая картина для прямоугольника {(Y,a/Y)} = (0;2п) х (0; 1) неоднозначна, сеточная функция задана на слишком крупном шаге - 0.1 по первому и 0.01 по второму аргументу (а меньший шаг естественно требует ресурсоемкости). Тем не менее можно сделать ряд замечаний:
• нормализованная на начальную амплитуду, т.е. P/P(t = 0),P(t = 0) = 1/(2п), изменяется в пределах ~ 0.24^1.38:
• вплоть до y ~ 4.5 вертикальная структура контурного графика повторяется, при этом P сильно зависит от y, но слабо от значения a/Y; при y > 6 линейчатая структура разрушается, появляются островки, причем крупные характерны для нижней половины и выгнуты вправо;
• период по аргументу y составляет примерно 0.4.
Начальный фрагмент показан на рис. 11, где видна периодичность.
Рис. 9. Фурье-образ решения (15) при некоторых парах (a/Y, y)
Можно фиксировать один из параметров, например y, и проводить «концевой» анализ по обычной схеме. Для y = 0.01 примечательно, что амплитуда монотонно растет (и частота соответственно уменьшается), хотя гармоничность нестационарности среды не предполагает изначально направления изменения амплитуды; на кривой набега фазы феномен арки отсутствует. Однако более интересен случай y = 0.5 ^ 2п, см. рис. 10. Во-первых, изменилось направление изменения амплитуды; во-вторых, на кривых амплитуды и частоты становятся заметны «дрожания», более сильные в конце, т.е. при 0.5 < a/Y < 0.95, плот-
ность которых возрастает; в-третьих, на заключительном отрезке мы наблюдаем феномен арки (здесь наиболее выделены три арки, одна - по огибающей, две других - по самой кривой), но этот эффект отсутствует на кривой набега фазы. Теперь пусть 7 = 0.4. Что изменилось? Феномен арки стал наблюдаться и на кривой набега фазы, причем одновременно с двумя другими кривыми (амплитуды и частоты). «Дрожания» стали менее выраженными. Концевая амплитуда вновь стала возрастать, а не падать; её поведение сложное - так, в ситуации аналогичной рис. 10, мы наблюдаем рост амплитуды при 7 = 0.01,1.1 и её падение при 7 = 0.1, 0.6,1.5. При 7 ~ п и выше методика перестает работать, т.к. нормированный набег фазы превышает 0.5.
Наличие арок на последних кривых (рис. 10) заставляет подробнее рассмотреть точки (7,а/7) = (0.4, & 0.75) и (0.4, & 0.83), приходящиеся примерно на центры арок. Обратим внимание, что ранее совпадавшие кривые номинальной (О1) и экстрагированной (^2) частот начали расходиться, что объясняется «включением» слагаемого ¿О' в приближенной формуле (6).
Как показывает рис. 11, понятие «амплитудная модуляция» является хорошим представлением для анализа.
Рис. 10. Контурный график амплитуды в области малых значений аргументов. Амплитуда нормализована на 1/(2п). Пределы и шаг (в скобках) изменения по аргументам: 7 ^ 0.0001 : (0.002) : 0.2, a/j ^ 0.0001 : (0.002) : 0.05
Рис. 11. «Концевой» анализ для синфазного гармонического управления. Случай Y = 0.5, a/Y = 0.05 : 0.001 : 0.95. Обозначения рис. 7
t
Рис. 12. Графики решения (14) в точках (y, a/Y) = (0.4, « 0.75) u (0.4, « 0.83). Для удобства: знак тильды у аргумента не ставим, обе кривые нормированы на 1/2п , вторая кривая сдвинута вверх на +3
7.2. Синусоидальное управление, в = —а
(d \ ( b \ 1 + - sin(pt) I , ¡(t) = ¡(0)1 1 + - sin(pt M ^
p (16) д2Т + / y! + (4(2n)2 + a2 - Y2) + (4(2n)2 - a2 + y2)л T(t) = 0. lJ
дЪ2 V 4 4(1 + в)2(1 — в)
В этом случае не существуют значений коэффициентов нестационарности, дающих гармоническое решение, так как на всем отрезке числителя (16) нулю равняться не может. Тем не менее, как мы увидим дальше, решение может быть весьма похожим на гармоническое.
Рис. 13. Зависимость спектральной плотности Ш от частоты а для уравнения (16). Синусоидальная нестационарность при а = —в, т.е. противофазном управлении
В отличие от синфазного управления, исключая ситуации больших значений параметров (рис. 14) дискретности спектра, ВП не наблюдается, сам ВП компактен и имеет небольшую ширину. Даже при относительно значимых по величине параметрах нестационарности имеется один узкий пик, ширина которого близка к вычислительной ошибке (рис. 8).
Так, в предположениях рис. 9, но уже для противофазного управления, изменения амплитуды в пределах 4-го знака после запятой (причем всегда в отрицательную сторону). Требуется y ^ 1,a/Y & 0.1 , чтобы \P/P(t = 0) —1| & 0.0035. «Концевой» анализ показывает незначительность отклонений амплитуд от частот (по сравнению с 5.1), а феномен арки начинает проявляться при поздних значениях a/Y > 0.9. Кроме того, если в синфазном случае амплитуда изменяется практически линейно, то в противофазном случае амплитуда падает по выпуклой кривой. Небольшие значения y гарантируют незначительность «дрожаний» на кривых амплитуды и частоты, но не гарантируют малого разброса значений P, Q(y = const,a/Y): так, при y = 0.01 на отрезке a/Y Е [0.05; 0.95] нормированное значение P падает примерно на 0.22, а при y = 0.5 такое изменение составляет 0.03.
На рис. 14. показана ситуация достаточно сильной нестационарности, когда его частота составляет половину исходной частоты и амплитуда равна половине от исходной проницаемости среды.
(б)
Рис. 14. Сравнение результатов синфазного (а) и противофазного (б) управления в точке (7 = п, о/^ = 1/2). На левой стороне показан график решения Т(¿) (для удобства ограничен 50, а не 100), на правой - его фурье-разложение
7.3. Косинусоидальное управление
В рамках косинусоидального управления также можно выделить синфазный (17) - (18) и противофазный (20) случаи:
О \ / Ь \
е(г) = е(0) ( 1 + -(1 - еов(р«)П , ¡(г) = ¡(0) ( 1 + ~(1 - еов(р«))
- Т - (^п)2* (2 + Ч(2п)2 - Т) (1 + * - Ч3 + Р + *) а -
(а, - ) = --^-и-4 ,а -1-"-—--, ° = 1 - ^ сов(7- )•
(17)
Основное уравнение (2) с учетом подстановки (17) можно переписать в более удобной форме (18):
#т + ( £ + ) т (г> = о. (18)
д- V 4 4(а + а/7)2 Рассуждая аналогично (15), получаем условие гармоничности решения:
~ = 1 - ^ ° 7 = ^а2 + (4п)2 - а= 4п о ш- (19)
Из (18) и (19) следует, что точке минимума соответствует граница области, т.е. а/^ = т ~ 1. Если отбросить требование строгой гармоничности и принять слабую зависимость от времени знаменателя (18), то становится понятным такой численный результат - при (7,а/7) = (6п, 3/4) спектр содержит единственный узкий пик на частоте 0-35. Противофазное косинусоидальное управление приводит к (20):
О \ ( Ь \
е(г) = е(0) ( 1 + -(1 - еоё(рг))) , ¡(г) = ¡(0)1 1 + - (1 - еов(р«)) '
'± (т - т - т+(2п)2 (' + а)) - (т + £ + а
4 + / )
д2т д- +
у (1 + а (1 - ео«(7^) (1 - а (1 - еов(7-))] у
т (-) = 0-
(20)
Аналитическое исследование затруднено громоздкостью выражений, а численного исследования мы не проводили.
7.4. Сравнение синфазного и противофазного управления: концевая амплитуда
Остановимся на частном вопросе объяснения, исходя из вида (15) и (16) и представления уравнения Штурма-Лиувилля, почему в синфазном случае концевая амплитуда Р(а) = Р(7 = еоп8\,,а/^)1=и может как убывать, так и возрастать (дР/да >< 0), а в противофазном - только убывать?
Приведем несколько соображений. Заметим вначале, что, считая вполне верным (6), для относительного изменения амплитуды 5 можно записать линейное разложение (21):
Р М = N) / 2п д(а,М)
Р(•,- =0) у 2п + д(а^) ~ 4п - ()
Далее, если индексом 1 отметить синфазный, а индексом 2 - противофазный случаи, получим (22):
2(2п)2в а . АТ,
- = - (ТгЫг-о' 8 = а 5,п(7^ (22)
Сделав подстановку (22) в (21), получим независимость от а знака разности концевых амплитуд. С другой стороны, этот знак зависит от выбора концевой точки N и частоты возбуждения, что как будто снижает ценность этих замечаний. Однако надо учесть, что ¿(а) относится к дискретным отсчетам амплитуды Р, которая экстрагируется по определенной схеме, а не амплитуды Р, с которыми связано (22) и которые берутся из непрерывного соотношения (6). Поэтому если рассуждать для концевых значений, когда в течение полупериода ~ п/^ они не изменяются, то целесообразно усреднить (22). Тогда получим эллиптические интегралы (23), (24):
1 ('2п/т 2 С1 \/1 - х2йх а
<91 - 92 >=Ы*) - Ч2(Ъ))(И = 16пд —--> 0, д = -. (23)
2п/7 7о Уо (1 - д2 + д2*2)2 7
В тех же обозначениях среднее значение 91 зависит от обоих параметров:
т, а2 + 4Г Г1 (1+ д2г2)(1г . л,
< «м >= -г+-¿г- /о (1 -д+Ф)^. (24)
Мы видим, что знак <91 > (24) не определен, а знак <91 — 92 > положителен. Если абсолютное значение первой величины меньше второй, тогда становится ясным поведение концевой амплитуды в противофазном случае.
Можно посмотреть на уравнение Штурма-Лиувилля как на возмущенное поправкой 9(а, N) уравнение гармонического осциллятора, где частота 2п играет роль собственного значения оператора. Тогда попытаемся применить методы теории возмущений [9], ограни-
а
а
чившись членами второго порядка по параметру —
\7
диктует искать асимптотическое решение в виде (25):
7
< 1. Формальная техника
т(*) = То(*) + + (^2Т2(*), ЭД = ^^^М, Т1(0)= Т1(0) = Т2(0) = Т2(0) = 0.
7 \7/ 2п
(25)
Также предположим второй порядок в разложении возмущения. Введя два коэффициента, выпишем (26):
'(а) = 9(1,(^) а+а)2. (26)
Тогда, подставляя (25) в (14) и учитывая (26), приравняем к нулю коэффициенты разложения при первой и второй степени и получим условия для нахождения Т^© из дифференциальных уравнений (27):
Т1 + (2п)2Т1 = —9(1)То, Т2 + (2п)Т = -9(2,То - 9(1,Т1. (27)
Поскольку в правых частях неоднородных уравнений (24) стоят произведения синусов, то решением (27) будут сумма гармоник с некоторыми весовыми коэффициентами, а именно такая картина и наблюдается на рис. 12, 14. Эти веса могут быть выведены аналитически. Наибольший из них вес (и соответственно выделяющий основную моду), по-видимому, можно интерпретировать как концевую амплитуду.
Для синфазного случая, учтя (14), выпишем коэффициенты:
9(1) = -2Г8т(7*),9(2) = 2 (3Г + - 3ТС08(2Т*), Г ^ (2п)2 - ^. (28)
Для противофазного случая, учтя (15), также выпишем коэффициенты:
9(1) = ^ вш^*),?^ = 2г - 2 (г - ео8(27?).
,2
(29)
В синфазном случае из (25), (27), (28) следует выражение для Т1: со8(2п£) Г / еов(2п — 7)? еов(2п + 7)? ^
Т1
+
+
(30)
7 2^7 \ 4п — 7 4п + 7
Отметим, что, несмотря на видимое наличие в знаменателе (31) 7, Иш7^о ТД7) = 0, и это отвечает физическому смыслу. В правой части неоднородного дифференциального уравнения для Т2 находятся 5 мод со следующими весовыми коэффициентами:
В общем виде результат согласуется с экспериментом (рис. 7, 8, 13, 14). К сожалению, решение (31) будет всегда содержать, т.к. А0 = 0, растущее во времени резонансное слагаемое. Это говорит о том, данный асимптотический метод, даже если повысить порядок разложения, не совсем применим в нашем случае. Возможно, для лучшего анализа целесообразнее было бы применять метод усреднения Крылова-Боголюбова (формула (23) идейно с ним связана), а не разложение по малому параметру (25).
8. Заключение
Мы проанализировали несколько случаев управления ЭМВ через внесение нестационарности в среду распространения. Полученные нами результаты, к сожалению, не имеют прозрачного физико-технического смысла. Эмпирические результаты численного эксперимента на основании гармонического «концевого» анализа уравнения Штурма-Лиувилля, если суммировать наиболее важные из них, состоят в следующем:
• при гармоническом управлении спектр ВП обычно дискретный, что отличает этот случай от других;
• наблюдаются от 1 до 3-4 дискретных диапазонов изменения параметра нестационарности, когда набег фазы для концевой точки изменяется медленно; в этих диапазонах либо огибающая абриса, либо сама кривая набега фазы имеет форму арки;
• косинусоидальное управление, на наш взгляд, более перспективно с точки зрения практики, т.к. позволяет не только увеличивать исходную частоту ЭМВ, но и уменьшать её (но не более чем в у/2т + 1 ~ \/3 раз, считая, что технически проницаемость изменяема не больше чем в 2т + 1 ~ 3 раза); уменьшения частоты можно достичь и на большую величину, но лишь допустив небольшое уширение моды;
• направление перекачки энергии ЭМВ при гармоническом управлении: в синфазном случае определяется только частотой управления, а в противофазном случае - амплитуда волны падает (а частота возрастает).
Эффект арки, найденный нами впервые, представляется полезным для использования на практике, например, для информационной модуляции ЭМВ в волноводе. Происхождение эффекта арки и черты фрактальной геометрии на некоторых кривых набега фазы требуют дополнительного теоретического исследования.
9. Приложение
1; М1=1е+5; N=100; М=1024; Fs=Mt/N; ё=1е-5; Ь=0.001; 2; ор^оп8 = оёе8е^'Ке1То1',ё,'АЬ8То1',[ё ё]); 3; ^,Т] = оёе45(@шах№ге11,0:Ь^,[0;1],ор^оп8); 4; а=1т8расе(0^,М); а=Й'; х=Т(:,1); х1=1^егр1(^х,и,'8рНпе');
г=2
г=—2
5; Y=fft(x1,M)/M; Y=Y.*conj(Y); W=fftshift(Y); 6; Omega=((t1-N/2)*Fs/N); plot(Omega,W);
Литература
1. Pedrosa I.A., Petrov A.Y., Rosas A. On the electrodynamics in time-dependent linear media // Eur. Phys. J. D. 2012. V. 66.
2. Красников Г.Я., Матюшкин И.В., Черняев Н.В., Горнев Е.С., Евстратов Н.В. Электромагнитная волна в среде с нестационарными проницаемостями: часть 1 // Математическое моделирование. 2015. Т. 27, № 12. C. 48-64.
3. К'расников Г.Я. Конструктивно-технологические особенности субмикронных МОП-транзисторов. М.: Техносфера, 2002. 414 с.
4. Буфетов А.И., Гончарук Н.Б., Ильяшенко Ю.С. Введение в теорию Штурма-Лиувилля: Конспект курса «Обыкновенные дифференциальные уравнения». М.: Мех.-мат. фак. МГУ им. М. В. Ломоносова, 2015. 153 с.
5. Черняев Н.В., Матюшкин И.В., Евстратов Н.В. Численное и аналитическое исследование уравнений телеграфного типа, описывающих распространение электромагнитной волны в среде с деградацией // Материалы конференции «Математическая физика и ее приложения». 2014. С. 158-159.
6. Лазарев Ю.Ф. Начала программирования в среде MATLAB: учебное пособие Киев: НТУУ «КПИ», 2003. С. 94-95.
7. MatLab M. The language of technical computing // The MathWorks, Inc. 2012.
8. Родионов С.А., Вознесенский Н.Б., Домненко В.М. Численные методы в оптике: электронный учебник. СПб.: СПБГУ ИТМО, кафедра прикладной и компьютерной оптики.
Гл. 5, раздел 5.
9. Кельберт М. Что такое преобразование Фурье? // Математическое просвещение. 2000. Сер. 3, вып. 4. С. 188-202.
10. Федорюк М.В. Обыкновенные дифференциальные уравнения. 2-е изд., перераб. и доп. М.: Наука. Главная редакция физико-математической литературы, 1985. 448 с.
References
1. Pedrosa I.A., Petrov A.Y., Rosas A. On the electrodynamics in time-dependent linear media. Eur. Phys. J. D. 2012. V. 66.
2. Krasnikov G.Ia, Matiushkin I.V., Cherniaev N.V., Gornev E.S., Evstratov N.V. Electromagnetic wave's propagation in medium with nonstationary permittivity and permeability: Part 1. Matematicheskoe Modelirovanie. 2015. V. 27, N 12. P. 54-57. (in Russian).
3. Krasnikov G.Ia. Konstruktivno-tehnologicheskie osobennosti submikronnyh MOP-transistorov. M.: Tehnosfera, 2002. 414 p. (in Russian).
4. Bufetov A.I., Goncharuk N.B., IUashenko U.S. Vvedenie v teoriiu Shturma-Liuvillia. Konspekt kursa «Obyknovennye differentsialnye uravneniya». M.: Meh.-mat. Fac. MGU im. M. V. Lomonosova, 2015. 153 p. (In Russian).
5. Cherniaev N.V., Matiushkin I.V., Evstratov N.V. Chislennoe I analiticheskoe issledovanie uravnenii telegrafnogo tipa, opisyvayushchikh rasprostranenie electromagnitnoi volny v srede s degradatsiei. Materialy konferentsii «Matematicheskaya fizika i eyo prilozheniya». 2014. P. 158-159. (in Russian).
6. Lazarev Iu.F. Nachala programmirovaniya v srede MATLAB. Uchebnoe posobie. Kiev: NTUU «KPI», 2003. P. 94-95. (in Russian).
7. MatLab M. The language of technical computing. The MathWorks, Inc. 2012.
8. Rodionov S.A., Voznesenskii N.B., Domnenko V.M. Chislennye metody v optike. Electronnyi uchebnik SPb, SPBGU ITMO, kafedra prikladnoi i kompyuternoi optiki. Ch. 5, S. 5. (in Russian).
9. Kelbert M. Chto takoe preobrazovanije Furje? Matematicheskoe prosveshchenie. 2000. V. 3, N 4. P. 188-202. (in Russian).
10. Fedoruk M.V. Obyknovennye differentsialnye uravneniya. 2-e izd., pererab. I dop. M.: Nauka. Glavnaya redaktsiya fiziko-matematicheskoi literatury, 1985. 448 p. (in Russian).
Поступила в редакцию 20.10.2016
