Диссертация (Компьютерная реализация геометрических методов в максвелловской оптике), страница 16

И она естественно разбилась на две части. Вначале191использовалась система Cadabra, поскольку эта система и была разработана для задач манипулирования тензорными объектами. Все манипуляции проводились исключительно с абстрактными индексами. Использование системы компьютерной алгебре на этом этапе позволилополностью исключить бумагу и ручку из теоретической работы. Навтором этапе возникла, собственно, необходимость создания программного комплекса для пакетного решения однотипных задач. Для этогобыла применена система FORM.

Данная система специализируется навекторных и тензорных вычислениях. Её спецификой является использование неинтерактивного (пакетного) подхода, а также возможностьиспользования распараллеливания вычислений и внешней памяти.Продемонстрируем вычисление выражения (6.29) в системе компьютерной алгебры Cadabra1 .Зададим список индексов :{\alpha,\beta,\gamma,\delta} : :Indices(vector).Введём частную и ковариантную производные2 :\partial_{#} : :PartialDerivative.\nabla_{#} : :Derivative.При выводе выражения (6.29) для нас будут иметь значения симметрийные свойства символов Кристоффеля, которые зададим с помощьюдиаграммы Юнга [182] :\Gamma^{\alpha}_{\beta \gamma} : :TableauSymmetry(shape={2}, indices={1,2}).1Точка в конце команды подавляет вывод результата, а точка с запятой заставляет напечататьрезультат.2Знак # является шаблоном для любого выражения.192Для тензора Fαβ нет необходимости приводить диаграмму Юнга,достаточно указать, что он антисимметричный.F_{\alpha \beta} : :AntiSymmetric.Запишем выражение для ковариантной производной в виде подстановки1 :nabla :=\nabla_{\gamma} A ?_{\alpha \beta} ->\partial_{\gamma}{A ?_{\alpha \beta}} A ?_{\alpha \delta}\Gamma^{\delta}_{\beta \gamma} A ?_{\delta \beta}\Gamma^{\delta}_{\alpha \gamma} ;nabla := ∇γ A?αβ → (∂γ A?αβ − A?αδ Γδβγ − A?δβ Γδαγ );Мы используем постфикс « ?», превращающий предыдущую букву вшаблон (в противном случае подстановка могла быть использованалишь с фиксированной переменной).

Заметим, что применение модификатора « ?» для индексов необязательно.Теперь запишем уравнение (2.6) :maxwell1 := \nabla_{\alpha}F_{\beta \gamma} +\nabla_{\beta} F_{\gamma \alpha} +\nabla_{\gamma} F_{\alpha \beta} ;maxwell1 := ∇α Fβγ + ∇β Fγα + ∇γ Fαβ ;Выполним подстановку nabla в выражение maxwell1 :@substitute !(maxwell1)(@(nabla)) ;maxwell1 := ∂α Fβγ − Fβδ Γδγα − Fδγ Γδβα ++ ∂β Fγα − Fγδ Γδαβ − Fδα Γδγβ +1Выражение := служит для задании метки.193+ ∂γ Fαβ − Fαδ Γδβγ − Fδβ Γδαγ ;В выражении maxwell1 приведём подобные :@canonicalise !(%) ;@collect_terms !(%) ;maxwell1 := ∂α Fβγ − ∂β Fαγ + ∂γ Fαβ ;В результате мы получили в точности уравнение (6.28).Теперь запишем уравнения Максвелла (2.6) и (2.7) в среде в декартовых координатах с метрическим тензором ηα β = diag(1 , −1 , −1 , −1) :∂α Fβγ + ∂β Fγα + ∂γ Fαβ = 0,∂α Hαβ4π β=j .c(C.1)Аналогично запишем вакуумные уравнения Максвелла в эффективном римановом пространстве с метрическим тензором gαβ :∂α fβγ + ∂β fγα + ∂γ fαβ = 0, 4π β√1√ ∂α −ghαβ =j̃ .−gc(C.2)Здесь fαβ и hαβ суть тензоры Максвелла и Минковского в эффективномримановом пространстве1 .В Cadabra для этого введём ковариантную дивергенцию :div :=\nabla_{\alpha}A ?^{\alpha \beta}->1/\sqrt{-g}\partial_{\alpha}(\sqrt{-g} A ?^{\alpha \beta}) ;1Заглавные буквы будем использовать для величин в пространстве Минковского, строчные — длявеличин в римановом пространстве.194√1div := ∇α A?αβ → √ ∂α ( −gA?αβ );−gЗдесь опять постфикс « ?» превращает предыдущий символ в шаблон.Уравнение (2.7) запишем для случая эффективного риманового пространства :riman :=\nabla_{\alpha} h^{\alpha \beta}=j^{\beta} 4 \pi / c ;πriman := ∇α hαβ = 4j β ;cПодстановка ковариантной дивергенции в уравнение с меткой rimanдаст нам уравнение (C.2).@substitute !(riman)(@(div)) ;√1πriman := √ ∂α ( −ghαβ ) = 4j β ;−gcПоскольку для вакуума выполняется соотношение :fαβ = hαβ ,(C.3)f αβ = g αγ g βδ hγδ .(C.4)то, подняв индексы, получим :В Cadabra для этого сделаем соответствующую подстановку :fh :=h^{\alpha \beta} ->g^{\alpha \gamma}g^{\beta \delta}f_{\gamma \delta} ;f h := hαβ → g αγ g βδ fγδ ;195@substitute !(riman)(@(fh)) ;√1πriman := √ ∂α ( −gg αδ g βγ fδγ ) = 4j β ;−gcСделав почленное сравнение, получим :√Fαβ = fαβ ,j α = −gj α ,√H αβ = −gg αγ g βδ Fγδ .(C.5)(C.6)Отсюда получаем соотношение :1Fαβ = √ gαγ gβδ H γδ .−g(C.7)Из соотношения (C.7) получаем явную запись для вектора электрической индукции :√Di = −−g ij1 ijkg Ej +ε gj0 Hk .g00g00(C.8)При этом геометризованная диэлектрическая проницаемость будетиметь вид :√ijε =−−g ijg .g00Проведём аналогичные выкладки и для вектора магнитной индукции.196Соответственно получим вектор магнитной индукции√Bi = −−g ij1 ijkg Hj −ε gj0 Ekg00g00(C.9)и геометризованную магнитную проницаемость√ijµ =−C.4.6.−g ijg .g00(C.10)Пример расчёта на FORMПроделанные с использованием системы Cadabra символьные манипуляции позволили получить результат, который можно непосредственно применить для компьютерных моделирования и проектирования элементов систем максвелловской (тензорно-векторной) оптики, допускающих простое представление в криволинейных координатах.

Дляданных расчётов мы выбираем другую систему. Предлагается использовать для этого систему FORM.Продемонстрируем вычисление определителя метрического тензораg для задачи из раздела 5.1.2 на примере системы FORM. Будем использовать преобразование (5.7). Для компактности обозначим b := R2 ,a := R1 .Отключим дополнительную информацию по процессу вычисления(ресурсы, время). Для того, чтобы вывод результатов уместился в статью, уменьшим ширину вывода до 40 символов (по умолчанию 80).Off statistic ;Format 40 ;197Укажем явно, что мы работаем в 4-мерном пространстве. Далее зададим основные элементы : индексы (i, j, k, l)1 , тензоры (gαβ ) и объекты, не имеющие дополнительной семантики (a, b) :Dimension 4 ;Indices i, j, k, l ;Tensors g ;Symbols a,b ;Зададим формулу для определителя (поскольку в FORM такаяфункция не задана).

Для этого используем символ Леви–Чевиты :det{gαβ } =: g = ε0123 εαβγδ g0α g1β g2γ g3δ .Local detG = e_(0,1,2,3) * e_(i,j,k,l)* g(0,i) * g(1,j) * g(2,k) * g(3,l) ;contract ;Print ;.sortРезультат будет представлен в следующем виде :detG =g(0,0)*g(1,1)*g(2,2)*g(3,3) - g(0,0)*g(1,1)*g(2,3)*g(3,2) - g(0,0)*g(1,2)*g(2,1)*g(3,3) + g(0,0)*g(1,2)*g(2,3)*g(3,1) + g(0,0)*g(1,3)*g(2,1)*g(3,2) - g(0,0)*g(1,3)*g(2,2)*g(3,1) - g(0,1)*g(1,0)*g(2,2)*g(3,3) + g(0,1)*g(1,0)*g(2,3)*g(3,2) + g(0,1)*g(1,2)*g(2,0)*g(3,3) - g(0,1)*g(1,2)*g(2,3)*g(3,0) - g(0,1)*g(1,3)*g(2,0)*g(3,2)+ g(0,1)*g(1,3)*g(2,2)*g(3,0) +g(0,2)*g(1,0)*g(2,1)*g(3,3) - g(0,2)*g(1,0)*g(2,3)*g(3,1) - g(0,2)1В качестве индексов в FORM будем использовать латинские буквы вместо греческих для большейкомпактности записи.198*g(1,1)*g(2,0)*g(3,3) + g(0,2)*g(1,1)*g(2,3)*g(3,0) + g(0,2)*g(1,3)*g(2,0)*g(3,1) - g(0,2)*g(1,3)*g(2,1)*g(3,0) - g(0,3)*g(1,0)*g(2,1)*g(3,2) + g(0,3)*g(1,0)*g(2,2)*g(3,1) + g(0,3)*g(1,1)*g(2,0)*g(3,2) - g(0,3)*g(1,1)*g(2,2)*g(3,0) - g(0,3)*g(1,2)*g(2,0)*g(3,1)+ g(0,3)*g(1,2)*g(2,1)*g(3,0) ;Конкретные значения компонентов в FORM задаются с помощьюподстановки по шаблону (в данном случае gαβ = diag[1 , −b2 /(b −a)2 , −1 , −1]).

Постфиксный модификатор « ?» превращает предыдущийсимвол в шаблон :id g(0,0) = 1 ;id g(1,1) = - b^2/(b-a)^2 ;id g(i ?,i ?) = - 1 ;id g(i ?,j ?) = 0 ;Print ;.sort.endВ результате, собственно, получаем значение определителя :b2det{gαβ } =: g = −.(b − a)2detG =- 1/(b^2 - 2*a*b + a^2)*b^2 ;C.4.7.Пример тензорных вычислений в MaximaПоскольку Cadabra на данный момент не поддерживает компонентные вычисления, продемонстрируем их в системе Maxima. В качестве199примера рассмотрим запись уравнений Максвелла в цилиндрическихкоординатах в голономном базисе (см. раздел 2.7.1).Загрузим вначале написанный нами небольшой пакет, содержащийопределения для дифференциальных операторов :(%i1) load("diffop.mac")$Зададим цилиндрическую систему координат :(%i2) ct_coordsys(polar cylindrical)$Посмотрим компоненты метрического тензора gij :lg ;1 0 0 2(%o3) 0 r 00 0 1Вычислим и посмотрим компоненты метрического тензора g ij :(%i3)(%i4)cmetric()$ug ;1 0 0 1(%o5) 0 r2 00 0 1Зададим необходимые векторы и определим их зависимость от ко(%i5)ординат.

Из-за ограничений Maxima обозначим j 1 через j1, а j1 черезj_1 :(%i6)(%i8)j :[j1,j2,j3]$depends(j,cons(t,ct_coords))$B :[B1,B2,B3]$depends(B,cons(t,ct_coords))$200(%i10) D :[D1,D2,D3]$depends(D,cons(t,ct_coords))$(%i12) H :[H_1,H_2,H_3]$depends(H,cons(t,ct_coords))$(%i14) E :[E_1,E_2,E_3]$depends(E,cons(t,ct_coords))$Теперь вычислим все части уравнений Максвелла в цилиндрическихкоординатах :(%i16) Div(B) ;dddB1(%o16)B3 +B2 +B1 +dzdθdrr(%i17) Div(D) -4*%pi*rho ;dddD1(%o17)D3 +D2 +D1 +− 4πρdzdθdrr(%i18) Rot(H)+ diff(transpose(matrix(D)),t)/c; - 4*%pi/c*transpose(matrix(j))ddθH_3− ddz H_2|r|ddtD1c− 4 πcj1  d H_3− d H_1d4 π j2 drdzd t D2(%o18) + c − c −|r| ddd4 π j3d r H_2− d θ H_1d t D3+ c − c|r|(%i19) Rot(E); + diff(transpose(matrix(B)),t)/cddθE_3− ddz E_2|r| d B2dt(%o19)  c −dddr++ddtB1c E_3− ddz E_1 |r|E_2− ddθ E_1|r|ddtB3c+Таким образом, результат совпал с полученными в разделе 2.7.1 анаdrлитическими выражениями.201Список иллюстраций5.15.25.35.45.55.65.7C.1C.2C.3C.4Исходная система координат .