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

Всероссийская открытая научная конференция «Современные проблемы дистанционного зондирования, радиолокации, распространения и дифракции волн» - «Муром 2021»

Особенности численного решения бихарактеристической системы уравнений методами компьютерной математики.

Ю.И. Бова, А.С. Крюковский

АНО ВО «Российский Новый Университет»,

Москва, ул. Радио, д.22, kryukovsky56@yandex.ru, julia_bova@mail.ru

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

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

Features of the numerical solution of bi-characteristic equations by methods of computer mathematics

Yu.I. Bova, A.S. Kryukovsky

Russian New University, Moscow, st. Radio, 22,

kryukovsky56@yandex.ru, julia_bova@mail.ru

The features of applying the methods of computer mathematics to the numerical solution of the bi-characteristic system of equations describing the propagation of radio waves in the ionospheric plasma are considered. Based on the comparison of the numerical and exact solutions obtained for the parabolic plasma layer, the accuracy of the numerical solution was investigated. Key words: ionosphere, bicharacteristic system, propagation, radio waves, testing, plasma, computer mathematics.

Для определения лучевых траекторий, характеризующих величину и направление распространения энергии в задачах радиосвязи и радиозондирования, необходимо построение численных методов и алгоритмов решения системы бихарактеристических уравнений

dr _ ЗГ dk _ ЗГ dt _ ЗГ dm _ ЗГ

dr З к dr З r' dr dm' dr З t'

в которой т - параметр вдоль траектории, г - координаты луча, t - групповое время, к - волновой вектор, m - круговая частота, Г - гамильтониан. Такая задача решалась во многих исследованиях различными способами [1-4]. Главной особенностью настоящей работы является привлечения языков программирования четвёртого уровня и объединение численных и аналитических методов с помощью применения символьных вычислений.

Рассмотрим основные этапы такого подхода. Остановимся на структуре бихарактеристической системы (1). Система состоит из восьми обыкновенных дифференциальных уравнений первого порядка. Система нелинейная. Первые три уравнения определяют производные координаты луча r = (x, y, z) относительно

параметра вдоль траектории r. Вторые три уравнения определяют производные

волнового вектора к = (кх, ку, к2 ) относительно параметра вдоль траектории, седьмое

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

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

В данной работе рассматривается бихарактеристическая система с Гамильтонианом вида:

СО1

Г = кХ + к2у + к] -Б —. (2)

с 2

Определяющим для вычислений является вид эффективной диэлектрической проницаемости Б , которая зависит от е- заряда электрона, те - массы электрона и значения электронной концентрации. Даже в простейшем случае:

Г -- V - 2

Б = 1 - V , V =

— У

е N

-5-, (3)

ш„—

Б зависит как от координат (так как Б зависит от электронной концентрации, которая в свою очередь зависит от координат), так и от частоты —. В случае магнитоактивной среды эффективная диэлектрическая проницаемость зависит и от волнового вектора к . Такой средой является ионосферная плазма, находящаяся в магнитном поле Земли [5,6]. Кроме того, диэлектрическая проницаемость б может зависеть от группового времени I, поскольку в ионосфере могут распространяться гравитационные волны. Если же диэлектрическая проницаемость не зависит от группового времени, то, так как частота ш остается постоянной, можно исключить параметр т и преобразовать систему к «стационарному» виду, известному как бихарактеристическая система Гамильтона-Лукина:

с1г дГ /дГ ¿к дГ /дГ

— =--- -, — = — -. (4)

Ж д к / дш & д г I дш

Тогда параметром интегрирования становится групповое время

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

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

непригодно, когда плотность электронной концентрации задана интерполяционной функцией.

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

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

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

г\ т=0 = (х0. Уо> * о) (5)

- начальные условия для первой группы уравнений.

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

ш I- ~ , ш I— . „ 7 ш

kx(0) = —^/^cos^cos^, ky(0) = -^sm Lcos^, kz(0) = — д/^^т^ (6)

>0 cus^ cos^, ky (О) = --WCoSin L, cos^, kz (О) = -Л/й0

c ' c c

Вообще-то начальное значение эффективной диэлектрической проницаемости S0 само

зависит от начального значения волнового вектора, но поскольку эффективная

2

диэлектрическая проницаемостьзависит от волнового вектора через cos а , а

2 (Н0xkx + Н0yky + Н0zkz )

cos а =--, (7)

-*• 2

k

н 02

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

Если бихарактеристическая система стационарная (см. (4)), то шести начальных условий (5) и (6) оказывается достаточно. Они определяются при t=0. Для нестационарных задач необходимо задать ещё два условия:

ш(0) = ш0 (8)

- начальное значение круговой частоты, и

х (0) =

(9)

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

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

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

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

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

' 2 2

X * ,2*0 + 2и] 0= * - - Ч (10)

^ 0, * €[ ,2 + ]

В формуле (1) 0)0 это круговая плазменная частота в максимуме слоя, -полутолщина слоя, а - высота от поверхности земли до нижней границы слоя (см. рис. 2). Здесь и ниже предполагалось, что =100 км, а =200 км.

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

_X_

Определение основных констант

п, е, с, е, те, Я2 и др.

_X_

Определение управляющих констант кг, ЯК, Г, (Хо, Ус, п, С и др.

_X_

Определение параметров модели электронной концентрации ионосферы и магнитного поля Земли

Задание модели электронной концентрации ионосферы и магнитного поля Земли

Задание модели магнитного поля Земли

_X_

Задание модели эффективной диэлектрической проницаемости и определение гамильтониана

_X_

Вычисление гамильтониана бихарактеристической системы

V

0Г ог ог ог

Вычисление , д г д к дю' 1д7

Формирование списка начальных данных бихарактеристической системы

Формирование уравнений бихарактеристической системы

I

Вход в цикл]к

Вычисление значения угла выхода луча (или начальной частоты для частотно-модулированного излучения)

_X_

Вычисление лучевой траектории

_4_

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

_х_

Конец

Рис. 1. Базовый алгоритм построения решения бихарактеристической системы.

Можно показать [13], что уравнение лучей как функция г(х) представимо в виде:

* = ^

с

х-,

¿0 + 2Ь + (С бИ Ж - а еИ Ж),

а

2 **- Аа ь

а-с

а + с

с

- х —

0 < х < х

I

хI < х < х р

х р < х < х р ^ хI

(11)

где

Ж=а

2а

х

~8

ъ л с)

, c=sin(ц), 5=cos(ц), а = ш0 / ш.

ш0 =

4л е 2 N

0

шг

(12)

Л"

700 Р 600 -

0.0 0.5 1.0 1.5 2.0

И, 10б см"3

Рис. 2. Зависимость электронной концентрации от высоты.

В формулax (12) ш=2л/- круговая рабочая частота,/- рабочая частота, N0- значение электронной концентрации в максимуме слоя, ц- угол выхода луча с горизонтальной

осью х, то есть угол между вектором к в источнике излучения и положительным направлением оси х. При вычислениях предполагалось, что ^=2х106 см-3, /0=&0 / (2 л)= 12.6979 МГц,¡= 14 МГц, а= 0.906993.

На рис. 3 красным цветом приведена лучевая структура в параболическом слое в плоскости (х,г).

В области, которую можно условно ограничить по оси х интервалом (300 км, 500 км), а по оси г интервалом (80 км, 180 км), наблюдается сгущение лучей [14], что соответствует каустической фокусировке типа А5 («бабочка») [15,16].

Для тестирования выбрано четыре луча: пологий луч (угол выхода 30°), стандартный луч, отражающийся от ионосферы (угол выхода 60°), стандартный луч, проходящий ионосферу (угол выхода 70°) и луч, отражающийся от ионосферы, но близкий к лучам, проходящим ионосферный слой (угол выхода 64°) (см. рис. 4).

О 200 400 600 800 1000 1200 1400

X, ЕМ

Рис. 3. Лучевая структура в параболическом слое.

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

Рис. 4. Лучевая траектории с углами выхода: 30°- зелёный цвет, 60°- синий, 64°- фиолетовый, 70°- тёмно-жёлтый.

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

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

точного решения и тогда находится горизонтальное отклонение Ах у = ху — ху , либо по

х у определяется значение координаты г( 2у ) для точного решения и тогда находится

вертикальное отклонение А2 у = 2у — 2у .

На рис. 5 приведено отклонение А2 в метрах для луча с углом выхода 30°.

™ I, км

Рис. 5 а Рис. 5 б

Рис. 5 Отклонение А2 , угол выхода луча 30°, стационарная система - а, полная система - б

Здесь и ниже рассмотрены два варианта бихарактеристической системы: стационарная система из шести уравнений (а) и полная система из восьми уравнений (б). В данном случае отклонение положительное. На рис. 5,а можно выделить три участка: начальный участок (до входа луча в слой), на котором отклонение практически равно нулю, средний участок (луч в слое), на котором отклонение резко возрастает (ошибка накапливается) и конечный участок (луч выходит из слоя), на котором отклонение убывает. В случае полной бихарактеристической системы (рис. 5, б) зависимость ошибки иная: на первом участке ошибка практически постоянна и близка к нулю, на втором участке она осциллирует, а на последнем участке она линейно нарастает. Причём максимальное значение ошибки на порядок меньше. Таким образом, для пологих лучей полная бихарактеристическая система предпочтительнее.

На рис. 6 показано относительное отклонение А2 / 2, показывающее какую реальную точность (выбирающуюся «автоматически») вместо 16 разрядов обеспечивает встроенная функция, интегрирующая бихарактеристическую систему.

Рис. 6 а Рис. 6 б

Рис. 6. Относительное отклонение А* / 1, угол выхода луча 30°, стационарная система - а, полная система - б

На рис. 7 приведено отклонение А в метрах для луча с углом выхода 60°. Как и предыдущем случае, на рис. 7 а можно выделить три участка: начальный участок (до входа луча в слой), на котором отклонение практически равно нулю, средний участок (луч в слое), на котором отклонение резко возрастает (ошибка накапливается) и конечный участок (луч выходит из слоя), на котором отклонение возрастает значительно медленнее.

Г, KM X, ЕМ

Рис. 7 а Рис. 7 б

Рис. 7. Отклонение Az, угол выхода луча 60°, стационарная система - а, полная система - б

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

На рис. 8 приведено отклонение Az в метрах для луча с углом выхода 64°. Этот угол близок к критическому углу, при котором луч не выходит из слоя, но асимптотически приближается к линии z=zo. При углах, больших критического, лучи проходят ионосферный слой насквозь. В этом случае для обеих бихарактеристических систем отклонения меньше нуля, а абсолютные значения несколько больше и достигают 21 м и 43 м соответственно.

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

Рис. 8 а Рис. 8 б

Рис. 8. Отклонение А2, угол выхода луча 64°, стационарная система - а, полная система - б

Максимальное отклонение стало меньше, причем отклонение опять положительное. На рис. 9 б показано отклонение по оси хтакже для стационарной бихарактеристической системы.

Рис. 9 а Рис. 9 б

Рис. 9. Номер угол выхода луча 70°, стационарная система, отклонение А2 — а, отклонение Ах — б.

Заключение

Возникает вопрос: достаточна ли точность вычислений для задач, посвященных моделированию распространения радиоволн в ионосфере Земли [17]. На этот вопрос следует ответить утвердительно, поскольку модель ионосферной плазмы обычно известна с точностью не выше нескольких километров [5, 18], а максимальное отклонение порядка 40 м, то есть две длины волны в пустоте ( Л = 21.4138 м). Однако, если возникает необходимость увеличить точность вычислений, то в первую очередь следует обратить внимание на процедуру, обеспечивающую решение бихарактеристической системы, поскольку потеря заявленной точности возникает именно при численной решении системы обыкновенных дифференциальных уравнений. На рис. 10 приведены значения отклонений при угле выхода луча 64° с повышенной точностью. Параметр «рабочая точность» увеличен до 32. Видим, что точность вычислений резко возросла и теперь порядка 1 мм.

Рис. 10 а Рис. 10 б

Рис. 10. Отклонение А , угол выхода луча 64°, стационарная система - а, полная система - б, (рабочая) точность 32

Литература

1. Казанцев А.Н., Лукин Д.С., Спиридонов Ю.Г. Метод исследования распространения радиоволн в неоднородной магнитоактивной ионосфере. // Космические исследования, 1967. Т. 5. Вып. 4. С. 593-600.

2. Лукин Д.С., Спиридонов Ю.Г. Применение метода характеристик для численного решения задач распространения радиоволн в неоднородной и нелинейной среде. // Радиотехника и электроника, 1969. Т. 14. № 9. С. 1673-1677.

3. Крюковский А.С., Куркин В.И., Ларюнин О.А., Лукин Д.С., Подлесный А.В., Растягаев Д.В., Черняк Я.М. Численное моделирование амплитудных карт для скорректированной модели IRI-2012 с плавными возмущениями ионосферы // Радиотехника и электроника. 2016. Т. 61. № 8. С. 794-799.

4. Ипатов Е.Б., Крюковский А.С., Лукин Д.С., Палкин Е.А., Растягаев Д.В. Методы моделирования распространения электромагнитных волн в ионосфере с учетом распределений электронной концентрации и магнитного поля Земли // Радиотехника и электроника. 2014. Т. 59. № 12. С. 1180-1187.

5. Дэвис К. Радиоволны в ионосфере М.: Мир.1973. 502.

6. Крюковский А.С., Скворцова Ю.И. Математическое моделирование распространения радиоволн в нестационарной плазме с учетом кривизны поверхности Земли и ионосферных слоев // Вестник Российского нового университета. Серия: Сложные системы: модели, анализ и управление. 2016. № 1-2. С. 34-40.

7. Крюковский А. С., Лукин Д. С., Палкин Е. А., Растягаев Д.В. Теория катастроф в проблемах стационарной и нестационарной дифракции // Труды X школы - семинара по дифракции и распространению волн .7 - 15.02.1993./М.: МФТИ. 1993. С. 36-111.

8. Крюковский А.С., Лукин Д.С. Краевые и угловые катастрофы в равномерной геометрической теории дифракции. Учебное пособие. М.: МФТИ, 1999. 134 с.

9. Кравцов Ю.А., Островский Л.А., Степанов Н.С. Геометрическая оптика неоднородных и нестационарных движущихся сред.//ТИИЭР. 1974. Т.62. № 11. C.91 -112.

10. Felsen L.B. Transients in dispersive media, part 1: theory // IEEE Trans. on Ant. and Prop. 1969. AP-17.№ 2. P.191 - 200.

11. Lewis R.M. Asymptotic theory of transients //In: Electromagnetic Wave Theory. Part 2. Ed. by J. Brown / N.Y.: Pergamon Press. 1967. P.845-869.

12. Анютин А.П. Асимптотическая теория распространения радиосигналов в неоднородной плазме. //Распространение радиоволн в ионосфере. М.: ИЗМИР АН СССР. 1978. C.29-36.

13. Крюковский А.С., Лукин Д.С. К вопросу о поле в окрестности каустического острия в ионосферном плазменном слое. // РЭ, 1981. Т. 26. № 6. С. 1121 - 1126.

14. Орлов Ю.И. Особенности лучевых и каустических картин в неоднородном параболическом слое. // Изв. вузов: Радиофизика, 1977. Т. 20. № 11. С. 1669-1675.

15. Карепов С.Л., Крюковский А.С. Расчет волнового поля методом интерполяционной локальной асимптотики. //Радиотехника и электроника. 2001. Т. 46. № 1. С. 40 - 46.

16. Крюковский А.С. Локальные равномерные асимптотики волновых полей в окрестности основных и краевых каспоидных каустик. // Радиотехника и электроника. 1996. T.41. № 1. C. 59-65.

17. Бова Ю.И., Крюковский А.С., Лукин Д.С. Распространение частотно-модулированного излучения электромагнитных волн в ионосфере Земли с учетом поглощения и внешнего магнитного поля // Радиотехника и электроника. 2019. Т. 64. № 1. С. 3-14.

18. Куницын В.Е., Терещенко Е.Д., Андреева Е.С.Радиотомография ионосферы. М.: Физматлит, 2007. 345 с.