Программирование, 2023, № 4, стр. 33-38
СРЕДСТВА КОМПЬЮТЕРНОЙ АЛГЕБРЫ ДЛЯ ГЕОМЕТРИЗАЦИИ УРАВНЕНИЙ МАКСВЕЛЛА
А. В. Королькова a, *, М. Н. Геворкян a, **, Д. С. Кулябов a, b, ***, Л. А. Севастьянов a, b, ****
a Российский университет дружбы народов
117198 Москва, ул. Миклухо-Маклая, д. 6, Россия
b Объединенный институт ядерных исследований
141980 Московская область, Дубна, ул. Жолио-Кюри 6, Россия
* E-mail: korolkova-av@rudn.ru
** E-mail: gevorkyan-mn@rudn.ru
*** E-mail: kulyabov-ds@rudn.ru
**** E-mail: sevastianov-la@rudn.ru
Поступила в редакцию 05.08.2022
После доработки 04.09.2022
Принята к публикации 21.10.2022
- EDN: MGAIHN
- DOI: 10.31857/S0132347423020127
Аннотация
При расчете оптических приборов в рамках геометризованной теории Максвелла используются широко известные формализмы общей теории относительности и дифференциальной геометрии. В частности, для подобных вычислений требуется знать аналитический вид уравнений геодезических. Что приводит к необходимости вычислять большое количество однообразных математических выражений. Одним из предназначений средств компьютерной алгебры является облегчение работы исследователя путем автоматизации громоздких символьных расчетов. Таким образом, использование систем компьютерной алгебры представляется вполне очевидным действием. В работе рассмотрено несколько свободных реализаций символьных вычислений для аппарата общей теории относительности. В конце статьи приводится практический пример символьных расчетов для геометризованной теории Максвелла.
1. ВВЕДЕНИЕ
Геометрический подход к уравнениям Максвелла прошел несколько этапов в своем развитии. Первоначальный интерес был вызван общей теорией относительности. Это работы Л.И. Мандельштама, И.Е. Тамма [1–3], В. Гордона [4]. В отсутствии практических приложений интерес к данной теме угас. Новая вспышка интереса возникла во время золотого века теории относительности (1960–1975 гг.). Это работы Е. Плебаньского [5], Ф. де Феличе [6]. Впрочем, следует заметить, что и в этот период не удалось определиться с приложениями данной теории.
Новая вспышка интереса к геометрическому подходу проявилась в середине 2000-х годов на фоне интереса к метаматериалам [7]. Это работы Д.Б. Пендри [8, 9] и У. Леонгардта [10, 11]. Эти работы породили целое направление – трансформационную оптику [12].
Геометризация материальных уравнений Максвелла позволяет изменить взгляд на прямую и обратную задачу оптики. Нахождение траектории лучей по параметрам среды можно назвать прямой задачей оптики, а нахождение параметров среды по заданным траекториям лучей – обратной. Естественно, что обратная задача сложнее прямой. В геометризованной максвелловской оптике эти задачи меняются местами. Прямая задача – нахождение диэлектрической и магнитной проницаемости по заданной эффективной геометрии (по траекториям лучей), обратная – нахождение эффективной геометрии по диэлектрической и магнитной проницаемости. Сложность обеих задач сопоставимая.
Геометризация заключается в построении индуцированной (электромагнитной) метрики. При этом индуцированная метрика должна быть связанна с метрикой лабораторного пространства. В прямой задаче геометризованной оптики (соответствующей обратной задаче при обычном подходе) задается траектория в лабораторной системе координат, далее по этим траекториям задаются геодезические в индуцированной (электромагнитной) системе координат. По уравнениям геодезических строится метрика в индуцированной системе координат. Такого типа задачи характеризуются большим количеством несложных операций. Ручные вычисления не только выматывающие, но и подвержены накапливающимся ошибкам. Поэтому естественно выглядит желание переложить эти вычисления на плечи системы компьютерной алгебры.
1.1. Структура статьи
В разделе 3 описываются основные положения методики геометризации уравнений Максвелла. Выбор конкретной библиотеки для символьных преобразований описано в разделе 3. Процесс применения средств компьютерной алгебры для исследования в рамках геометризованной теории Максвелла представлен в разделе 4.
2. МЕТОДИКА ГЕОМЕТРИЗАЦИИ УРАВНЕНИЙ МАКСВЕЛЛА
Под геометризацией уравнений Максвелла мы понимаем эффективную геометризацию. То есть мы устанавливаем соответствие между уравнением Максвелла в среде (в плоском пространстве Минковского) и вакуумным уравнением Максвелла в эффективном римановом пространстве. Мы не считаем, что электромагнитное поле искривляет пространство, но устанавливаем соответствие между средой и метрикой [13–15].
Это позволяет рассматривать траектории лучей как геодезические в римановом пространстве. Задавая же вначале пути распространения лучей и восстанавливая по геодезическим метрику, мы можем решать обратную задачу оптики.
Входными данными для задачи перевода материальных уравнений Максвелла в вакуумные уравнения Максвелла в искривленном пространстве являются параметры среды: диэлектрическая $\varepsilon _{i}^{j}$ и магнитная $\mu _{i}^{j}$ проницаемости, показатель преломления среды ${{n}_{{ij}}}$. В связи со спецификой геометризации на основе квадратичной метрики [1, 2] реально используется только показатель преломления.
Вид эффективного метрического тензора ${{g}_{{\alpha \beta }}}$ подбирается из соображений симметрии задачи [2, 3]. Сам метрический тензор зависит от параметров среды.
Имея эффективный метрический тензор, можно решить геометризованные уравнения Максвелла и получить искомые траектории распространения электромагнитных волн. В рамках статьи упростим задачу и будем использовать приближение геометрической оптики. В этом случае лучи будут распространяться по геодезическим:
Для случая изотропной среды метрический тензор имеет следующий вид [1, 5]:
(1)
${{g}_{{\alpha \beta }}} = \left( {\begin{array}{*{20}{c}} {\frac{1}{{\varepsilon \sqrt \mu }}}&0&0&0 \\ 0&{ - \sqrt \mu }&0&0 \\ 0&0&{ - \sqrt \mu }&0 \\ 0&0&0&{ - \sqrt \mu } \\ {}&{}&{}&{} \end{array}} \right)$Поскольку в методе геометризации на основе квадратичной метрики диэлектрическая и магнитная проницаемость равны [5], то дополнительно запишем:
Из уравнения (2) видно, что стандартная геометризация уравнений Максвелла подходит для исследования оптических задач в приближениях геометрической оптики и приближении эйконала.
Таким образом можно отметить, что исследование в рамках геометризованной теории Максвелла подразумевает операции с тензорами в римановом пространстве. Каждое такое действие достаточно простое, но количество таких манипуляций просто огромно. И это может привести исследователя в состояние уныния. И только возможность использования средств компьютерной алгебры сможет примирить его с таким подходом.
3. БИБЛИОТЕКИ СИСТЕМЫ КОМПЬЮТЕРНОЙ АЛГЕБРЫ SymPy ДЛЯ ИССЛЕДОВАНИЯ ИСКРИВЛЕННЫХ ПРОСТРАНСТВ
Нами было рассмотрено несколько вариантов библиотек для работы с формализмом общей теории относительности, который мы применяем для геометризации уравнений Максвелла. Ранее для таких вычислений мы применяли систему компьютерной алгебры Cadabra [16, 17]. Однако в данном исследовании нам было интересно рассмотреть библиотеки для системы компьютерной алгебры SymPy [18]. Заметим, что мы рассматриваем преимущественно свободнораспространяемые программные продукты, и предпочитаем не использовать проприетарные системы компьютерной алгебры, как многие другие исследователи [19].
Основной сложностью при выводе системы уравнений, задающих геодезическую кривую, является вычисление символов Кристоффеля $\Gamma _{{\beta \gamma }}^{\alpha }$. Задача несколько облегчается тем, что используется диагональный метрический тензор, однако все равно остается сложной в вычислительном плане. Поэтому разумно вычислить символы Кристоффеля с помощью системы компьютерной алгебры.
Все вычисления будут проводиться с помощью библиотеки символьной алгебры SymPy [20] и оболочки Jupyter Notebook или Jupyter Lab [21].
Мы выделили две библиотеки SymPy: EinsteinPy и gravipy.
3.1. EinsteinPy
Модуль [22] представлялся нам крайне перспективным. Однако при попытке задать метрику с помощью стандартной символьной матрицы SymPy мы столкнулись с ошибкой, которая заключалась в невозможности задать дробь $\frac{1}{{\varepsilon \sqrt \mu }}$. Так, вместо дроби в метрике (1) появляется артефакт, имеющий вид выражения $1\varepsilon \sqrt \mu $. Кроме того, библиотека EinsteinPy позиционируется как модуль для численных расчетов и символьные вычисления являются не главной ее задачей.
То есть при всей привлекательности данная библиотека оказалась просто непригодной для символьных вычислений (численные расчеты на базе этой библиотеки мы не рассматривали).
3.2. Gravipy
Так как использовать встроенные средства SymPy оказалось невозможным из-за указанной выше ошибки, авторы попытались найти альтернативу, которая работала бы в рамках SymPy. Такой альтернативой оказался небольшой модуль gravipy [23], в зависимостях которого указана только библиотека SymPy и, опционально, оболочка JupyterLab. Сам модуль фактически состоит из одного файла. Авторы работы использовали версию 0.2.0. Последнее изменение датируется сентябрем 2019 года.
При отсутствии альтернативы мы решили использовать именно эту библиотеку.
Из-за минимальных зависимостей и компактности, модуль можно удобно использовать без установки. Для этого следует в директории с Jupyter-блокнотом, в котором проходят вычисления, создать каталог gravipy и в него поместить файлы __init__.py и tensorial.py из каталога gravipy с официального репозитория [23]. После чего можно пользоваться модулем импортируя его в Jupyter-блокноте обычным способом.
4. ПРАКТИЧЕСКОЕ ПРИМЕНЕНИЕ СИСТЕМЫ КОМПЬЮТЕРНОЙ АЛГЕБРЫ ДЛЯ ИССЛЕДОВАНИЯ УРАВНЕНИЙ МАКСВЕЛЛА
Проведем необходимые вычисления.
Вначале импортируем SymPy и gravipy:
import sympy as sp
from gravipy.tensorial import *
Настроим отображение формул, включив кодировку юникод, отображение формул в окружении equation и их рендеринг с помощью mathjax.
sp.init_printing(use_unicode=True,
↪ latex_mode='equation' ,
↪ use_latex='mathjax' )
Далее определим необходимые символы:
# Координаты пространство-времени
t, x, y, z = symbols(' t, x, y, z ' )
# Характеристики среды, как функции
epsilon = sp.Function(' epsilon ' )(x, y,
↪ z)
mu = sp.Function(' mu ' )(x, y, z)
n = sp.Function(' n ' )(x, y, z)
Для работы модуля gravipy также необходимо отдельно задать координатный 4-вектор
# Создаем координатный 4-вектор
X = Coordinates('X' , [t, x, y, z])
Для создания метрики необходимо вначале определить символьную матрицу, а затем преобразовать ее в метрический тензор. Матрицу задаем с помощью функции diag, которая импортируется из SymPy автоматически при импорте модуля gravipy
# Metric = diag(1 / epsilon /
↪ sp.sqrt(mu), -sp.sqrt(mu),
↪ sp.sqrt(mu), -sp.sqrt(mu))
Metric = diag(1 / n / sp.sqrt(n),
↪ -sp.sqrt(n), -sp.sqrt(n),
↪ -sp.sqrt(n))
После задания вспомогательной матрицы можно задать уже метрический тензор, для чего помимо матрицы Metric нужно передать еще координатный объект X.
g = MetricTensor('g' , X, Metric)
Для отображения метрики следует выполнить ячейку с кодом g(All,All) в результате чего получим следующую матрицу11:
После создания метрического тензора, сразу же можно вычислить символы Кристоффеля. Делается это с помощью функции Christoffel, которой требуется передать метрический тензор:
Gamma = Christoffel(’Gamma’, g)
для распечатки результата также посчитаем в отдельной ячейке выражение Gamma(All, All, All). Получим список из четырех матриц (срезу компонентного представления символа Кристоффеля):
Функция Geodesic позволяет вычислить уравнения геодезической линии на основе метрики. Символы Кристоффеля вычисляются ею автоматически и отдельно передавать их в качестве аргументов не требуется.
w = Geodesic(’w’, g, t)
Результат выполнения функции выглядит крайне громоздко. Нет необходимости приводить здесь этот вывод.
При вычислениях функция Geodesic автоматически использует параметризацию по параметру t, однако не отождествляет этот параметр с координатой времени и поэтому не производит замены $\frac{{{\text{d}}t}}{{{\text{d}}t}} = 1$ и $\frac{{{{{\text{d}}}^{2}}t}}{{{\text{d}}{{t}^{2}}}} = 0$. Параметризацию можно отключить вызовом следующей функции:
Parametrization.deactivate(t)
# выводит сообщение о том включена ли
↪ автоматизация или нет
Parametrization.info()
Однако на работе функции Geodesic эта настройка не сказывается, из-за чего требуется произвести некоторые дополнительные преобразования уже вручную.
Прежде всего подставим $\frac{{{\text{d}}t}}{{{\text{d}}t}} = 1$ и $\frac{{{{{\text{d}}}^{2}}t}}{{{\text{d}}t}} = 0$ во все уравнения. Первое уравнение принимает наиболее простой вид:
(3)
$ - \frac{3}{2}\left( {\frac{{\partial n}}{{\partial x}}\frac{{{\text{d}}x}}{{{\text{d}}t}} + \frac{{\partial n}}{{\partial y}}\frac{{{\text{d}}y}}{{{\text{d}}t}} + \frac{{\partial n}}{{\partial z}}\frac{{{\text{d}}z}}{{{\text{d}}t}}} \right) = 0,$В каждом из этих трех уравнений в качестве левой части выступает выражение (3), что дает возможность заменить его на ноль. В результате получаем систему из трех уравнений следующего вида:
(4)
$\begin{gathered} \frac{1}{4}\frac{{\partial n}}{{\partial x}}\left[ {{{{\left( {\frac{{{\text{d}}x}}{{{\text{d}}t}}} \right)}}^{2}} + {{{\left( {\frac{{{\text{d}}y}}{{{\text{d}}t}}} \right)}}^{2}} + {{{\left( {\frac{{{\text{d}}z}}{{{\text{d}}t}}} \right)}}^{2}} + \frac{3}{{{{n}^{2}}}}} \right] = n\frac{{{{{\text{d}}}^{2}}y}}{{{\text{d}}{{t}^{2}}}}, \\ \frac{1}{4}\frac{{\partial n}}{{\partial y}}\left[ {{{{\left( {\frac{{{\text{d}}x}}{{{\text{d}}t}}} \right)}}^{2}} + {{{\left( {\frac{{{\text{d}}y}}{{{\text{d}}t}}} \right)}}^{2}} + {{{\left( {\frac{{{\text{d}}z}}{{{\text{d}}t}}} \right)}}^{2}} + \frac{3}{{{{n}^{2}}}}} \right] = n\frac{{{{{\text{d}}}^{2}}z}}{{{\text{d}}{{t}^{2}}}}, \\ \frac{1}{4}\frac{{\partial n}}{{\partial z}}\left[ {{{{\left( {\frac{{{\text{d}}x}}{{{\text{d}}t}}} \right)}}^{2}} + {{{\left( {\frac{{{\text{d}}y}}{{{\text{d}}t}}} \right)}}^{2}} + {{{\left( {\frac{{{\text{d}}z}}{{{\text{d}}t}}} \right)}}^{2}} + \frac{3}{{{{n}^{2}}}}} \right] = n\frac{{{{{\text{d}}}^{2}}x}}{{{\text{d}}{{t}^{2}}}}. \\ \end{gathered} $В свою очередь данную систему уравнений второго порядка можно преобразовать в систему ОДУ первого порядка из шести уравнений с помощью следующей замены переменных:
Тогда уравнения (4) приводятся к виду:
Следует отметить, что аналогичная замена получается при решении уравнения эйконала [24] методом характеристик [25–28].
5. ЗАКЛЮЧЕНИЕ
В статье приведён алгоритм расчёта характеристик оптических приборов в рамках геометризованной теории Максвелла. Одним из основных этапов расчётов является получение уравнения геодезической и преобразование его в систему обыкновенных дифференциальных уравнений, пригодную для дальнейших вычислений. Поскольку данная задача требует от исследователя крайне кропотливых вычислений, то для её решения используется система компьютерной алгебры. В статье даётся сравнительный обзор нескольких систем компьютерной алгебры, предназначенных для решения задач общей теории относительности. На примере одного такого пакета демонстрируются основные операции при практическом исследовании в рамках геометризованной теории Максвелла.
Список литературы
Тамм И.Е. Электродинамика анизотропной среды в специальной теории относительности // Журнал Русского физико-химического общества. Часть физическая. 1924. Т. 56. № 2–3. С. 248–262.
Тамм И.Е. Кристаллооптика теории относительности в связи с геометрией биквадратичной формы // Журнал Русского физикохимического общества. Часть физическая. 1925. Т. 57. № 3–4. С. 209–240.
Mandelstam L.I., Tamm I.Y. Elektrodynamik der anisotropen medien in der speziellen relativittstheorie // Mathematische Annalen. 1925. Bd. 95. H. 1. S. 154–160.
Gordon W. Zur Lichtfortpflanzung nach der Relativita‥tstheorie // Annalen der Physik. 1923. Bd. 72. S. 421–456.
Plebanski J. Electromagnetic waves in gravitational fields // Physical Review. 1960. V. 118. № 5. P. 1396–1408.
Felice F. On the Gravitational Field Acting as an Optical Medium // General Relativity and Gravitation. 1971. V. 2. № 4. P. 347–357.
Smolyaninov I.I. Metamaterial ‘Multiverse’ // Journal of Optics. 2011. V. 13. № 2. P. 024004.
Pendry J.B., Schurig D., Smith D.R. Controlling Electromagnetic Fields // Science. 2006. V. 312. № 5781. P. 1780–1782.
Schurig D., Pendry J.B., Smith D.R. Calculation of Material Properties and Ray Tracing in Transformation Media // Optics express. 2006. V. 14. № 21. P. 9794–9804.
Leonhardt U. Optical Conformal Mapping // Science. 2006. V. 312. № June. P. 1777–1780.
Leonhardt U., Philbin T.G. Transformation Optics and the Geometry of Light // Progress in Optics. 2009. V. 53. P. 69–152.
Foster R., Grant P., Hao Y. et al. Spatial Transformations: from Tundamentals to Applications // Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2015. 8. V. 373. № 2049. P. 20140365.
Kulyabov D.S., Korolkova A.V., Sevastianov L.A. A naive geometrization of maxwell’s equations // The 15th small triangle meeting of Theoretical Physics. Star Lesn, 2013. P. 104–111.
Кулябов Д.С., Королькова А.В., Севастьянов Л.А. Простейшая геометризация уравнений Максвелла // Вестник РУДН. Серия. Математика. Информатика. Физика. 2014. № 2. С. 115–125.
Kulyabov D.S., Korolkova A.V., Sevastianov L.A. et al. Algorithm for lens calculations in the geometrized maxwell theory // Saratov Fall Meeting 2017. V. 10717 of Proceedings of SPIE. Saratov : SPIE, 2018. 4. P. 107170Y.1–6.
Королькова А.В., Кулябов Д.С., Севастьянов Л.А. Тензорные расчеты в системах компьютерной алгебры // Программирование. 2013. № 3. С. 47–57.
Кулябов Д.С., Королькова А.В., Севастьянов Л.А. Новые возможности второй версии пакета компьютерной алгебры cadabra // Программирование. 2019. № 2. С. 41–48.
Sandon D. Symbolic Computation with Python and SymPy. Independently published, 2021. ISBN: 979-8489815208.
Диваков Д.В., Тютюнник А.А. Символьное исследование спектральных характеристик направляемых мод плавно-нерегулярных волноводов // Программирование. 2022. № 2. С. 23–32.
Sympy. 2022. URL: http://www.sympy.org/ru/index.html.
Project jupyter. 2022. URL: https://jupyter.org/.
Einsteinpy–making einstein possible in python. 2022. URL:https://einsteinpy-einsteinpy.readthedocs.io/en /latest/index.html.
Gravipy tensor calculus package for general relativity based on sympy. 2022. URL: https://github.com/wojciechczaja/GraviPy.
Bruns H. Das Eikonal. Leipzig: S. Hirzel, 1895. Bd. 35.
Borovskikh A.V. The two-dimensional eikonal equation // Siberian Math. J. 2006. V. 47. P. 813–834.
Moskalensky E.D. Finding exact solutions to the two-dimensional eikonal equation // Num. Anal. Appl. 2009. V. 2. P. 201–209.
Kabanikhin S.I., Krivorotko O.I. Numerical solution eikonal equation // Sib. Elektron. Mat. Izv. 2013. V. 10. P. 28–34.
Kulyabov D.S., Korolkova A.V., Velieva T.R., Gevorkyan M.N. Numerical analysis of eikonal equation // Saratov Fall Meeting 2018. Vol. 11066 of Proceedings of SPIE. Saratov: SPIE, 2019. 6. P. 56.
Дополнительные материалы отсутствуют.
Инструменты
Программирование
Пн.-пт.: 10.00-19.00