Тензоры: опыт создания пользовательского пакета программ
Создание объекта с помощью описания его свойств и правил. Посмотрим, как реализован в пакете объект tensor и подобные ему. Создание нового тензора (тензорного поля) и работа с уже имеющимися происходит всегда в некотором базисе (системе координат), поэтому первым всегда должен быть выполнен один из операторов: makeBasis или makeBasisDif. Первый оператор имеет, как уже говорилось выше, следующий формат: makeBasis[<имя базиса>,<размерность>,<откуда>,<как>], где первый аргумент - имя, присваиваемое элементам базиса, второй -размерность пространства. Два оставшихся необязательных аргумента связывают новый базис с уже имеющимися: третий задает номер базиса, из которого осуществляется переход (по умолчанию - равен номеру текущего базиса), а четвертый - матрицу перехода, столбцы которой суть компоненты векторов создаваемого базиса в базисе, указанном в третьем аргументе:
In[174]:=
makeBasis[e_Symbol, n_Integer, num_Integer: - 1,
mat_List: {} ] : =
Module[{nnum},
If[num == -1, nnum = curBas, nnum = num] ;
numBas++; curBas = numBas;
If[mat != {} && coor[nnum] == {},
matCC[{nnum, curBas}] = mat;
matCC[{curBas, nnum}] = Inverse[mat];
graphBas = Join[graphBas,
{nnum -> curBas, curBas -> nnum}]];
basis[curBas] = Join[Table[Overscript[e, i], {i, n}] ,
Table[Underscript[e, i], {i, n}]]; coor[curBas] = {};
initScan];Сам модуль присваивает базису очередной номер, устанавливает номер текущего базиса в номер вновь созданного, сохраняет, если заданы, матрицу перехода и обратную к ней в массив матриц matCC и заносит информацию о том, что переход задан в граф graphBas. Наконец, в массив базисов basis помещается список символов базисных векторов и ковекторов. Оператор makeBasisDif устроен аналогично.
После того как базис определен, можно выполнить операцию makeTensor или какую-либо ее разновидность.
Для написания команд такого типа мы воспользовались возможностью Mathematica задавать значения функции только для тех значений аргумента, для которых мы этого хотим. При выполнении команды makeTensor[T] происходит следующее:
-
Устанавливаются значения на аргументе T для функции tensor, возвращающей True на именах заданных тензоров, функции type, возвращающей тип, и функции inBasis, хранящей номер базиса, в котором будут сейчас заданы компоненты тензора:
tensor[T] = True;type[T] = {p, q}; inBasis[T] = curBas; -
Затем определяется
val[T] = <таблица значений тензора в текущей системе координат - просто символы с индексами в нужном месте и числе>;
- Также в некоторых случаях переопределяются значения функций symm и skew - признаков симметричности или кососимметричности тензора.
Кроме того, написаны аналоги функций tensor и type, работающие с произвольным количеством аргументов:
In[175] := tensorQ[T__] := MatchQ[And @@ tensor /@ {T}, True];
typeQ[T__] := SameQ @@ type /@ {T};Первая проверяет, что последовательность ( sequence ) аргументов T_ состоит из тензоров, а вторая - что тип у всех тензоров из T_ одинаков.
Теперь можно реализовать тензорные операции в виде рекурсивных правил пересчета этих функций. Например:
In[177]:=
(*сумма тензоров - тензор, при условии, что типы совпадают*)
tensor [Plus [T__]] := Head [T] = = = Plus && typeQ @@ T;
(* произведение тензора и чего угодно - тензор*)
tensor [ (а_) * (T_) ] /; tensor[T] := True;Чтобы правильно вычислялись значения компонент тензоров, нужно правильно написать функцию val:
In[179] := (*сложение тензоров одного типа - покомпонентное*)
val [Plus [T_]] /; Head [T] = = = Plus && typeQ @@ T : =
Plus @@ val /@ List @@ T;
(*умножение тензора на что угодно - покомпонентное*)
val [ (а_) * (Т_) ] /; tensor [T] := a* val [T];SetDelayed с условиями как средство избежать бесконечной рекурсии. Пересчет значения тензора в новых координатах также выполняется с помощью функции val:
In[181]:=
val[T_Symbol] /; inBasis[T] != curBas :=
Module[{inB, res}, inB = inBasis[T]; inBasis[T] = curBas;
res = changeVal[T, inB, curBas];
val[T] = res[[l]];
If[ ! res[[2]], inBasis[T] = inB]; res[[l]]];Казалось бы, мы вычисляем val[T] через val[T] , и программа должна зацикливаться. Чтобы понять, почему этого не происходит, обратим внимание на условие, участвующее в нашем определении. Напомним, что оператор f[x_] /; cond[x] := < тело функции > ; задает функцию, действующую только при выполнении условия cond[x] . При этом функции с одним и тем же именем, но разными условиями, воспринимаются Mathematica как разные функции. Поэтому приведенная выше функция val[T Symbol] /; inBasis[T] != curBas работает, только если текущий базис curBas отличен от базиса inBasis[T] , в котором были вычислены в предыдущий раз компоненты тензора T. При обращении же к val[T] внутри модуля это условие не выполнено, так как уже произошло присвоение inBasis[T] = curBas. Поэтому на самом деле вызывается другая функция. Отметим, что сам тензорный закон спрятан в более громоздкую фу нкцию changeVal, которая пользователю недоступна. Эта функция возвращает список, первый элемент которого - пересчитанные или непересчитанные значения компонент тензора (последнее происходит, если связь между координатами неизвестна), а второй - логическая переменная, равная True, если пересчет произошел, и False в противном случае. Это позволяет вернуть начальное значение inBasis[T] , если пересчет не произошел. Отметим, что если просто написать условный оператор присвоения If[<можно пересчитать>, inBasis[T] = curBas] , то в случае, когда условие <можно пересчитать> не выполняется, произойдет зацикливание программы.
Использование незанятых операторов и программирование их свойств. Для таких операций, как тензорное произведение, внешнее произведение, симметричное произведение мы использовали предусмотренные в Mathematica незанятые операторы. Такие операторы удобно программировать с помощью свойств и атрибутов. Разберем более простой пример внешнего произведения (незанятый оператор Wedge, верхняя галка ^, которая с клавиатуры вводится Esc ^ Esc ), действующего на разложенных по базису формах (базисные ковекторы считаются разложенными по базису). Напомним, что запись аргумента определяемой функции в виде <имя> ?<условие> используется для задания функции только на тех аргументах, которые удовлетворяют этому условию.
In[182] :=
(*Устанавливаем атрибуты: Flat дает ассоциативность,
Oneldentity дает возможность записать произведение
одного элемента: Wedge[х] теперь равно х*)
SetAttributes[Wedge, {Flat, Oneldentity}];
(*Делаем Wedge на формах дистрибутивным по сложению*)
( (T_) ?skewBasQ + (S_) ? skewBasQ) ^ (R_) ? skewBasQ : =
T ^ R + S ^ R;
(R_) ? skewBasQ ^ ( (T_) ? skewBasQ + (S_) ? skewBasQ) : =
R ^ T + R ^ S;
(* разрешаем выносить скалярные множители*)
(R_) ?skewBasQ ^ ( (a_) * (S_) ? skewBasQ) := а* R ^ S;
( (a_) * (SR_) ?skewBasQ) ^ (S_) ? skewBasQ := a*R ^ S;
(* умножение на нуль*)
(R_) ? skewBasQ ^ 0 := 0 ;
0 ^ (R_) ? skewBasQ := 0;
(*косая симметрия на базисных элементах, при повторах нуль *)
Wedge[(T__) ?basisQ] /; Signature [ {T} ] != 1 : =
If[Signature[{Г}] ==0, 0,
-(Wedge[Sequence @@ Sort[{T}]])];В результате внешнее произведение фактически индуктивно определено на формах, разложенных по базису форм. Теперь этим можно воспользоваться для вычисления компонент внешнего произведения кососимметричных тензоров, заданных обычным образом, т. е. с помощью makeSkew:
В результате внешнее произведение фактически индуктивно определено на формах, разложенных по базису форм. Теперь этим можно воспользоваться для вычисления компонент внешнего произведения кососимметричных тензоров, заданных обычным образом, т. е. с помощью makeSkew:
In[l90]:=val[Wedge[T__]] /; skewQ[T] : =
Module[{n, q, prod}, n = nn[basis[curBas]];
q = Last[type[Wedge[T]]];
prod = Wedge @@ toBasisSkew /@ {T};
If[prod = = = 0,
ConstantArray[0, ConstantArray[n, q]] ,
toListSkew[prod]]] ;А именно, мы сначала раскладываем сомножители по базису форм (оператор toBasisSkew применяется к последовательности аргументов функции Wedge ), там выполняем внешнее умножение, определенное приведенными выше правилами, а потом - или возвращаем массив нулей нужного размера, или переводим результат обратно в массив компонент с помощью функции toListSkew.
Операции тензорного произведения ( CircleTimes ) и симметричного произведения ( SmallCircle ) реализованы аналогично. Операция внешнего дифференцирования реализована так же, хотя оператор для нее не понадобился, она названа обычной буквой d. Это, впрочем, вызвало дополнительные трудности, возникшие при написании правил для основных функций типа tensor. Сравним следующие две строки кода:
In[191] := tensor [Wedge [T]] := skewQ[T];
tensor[h_[T_]] /; h === d := skewQ[T];Дело, видимо, в том, что функцию d в выражении типа tensor[d[T]] (в отличие от Wedge в аналогичном выражении) Mathematica стремится сразу выполнить, что приводит к ошибкам при обработке определений типа
In[193] := tensor [d [ T_] ] := skewQ[T];
Чтобы избежать их, был придуман трюк с переменным именем функции (h_) и условием равенства этого переменного имени конкретному значению.
Обратим также внимание на операцию присвоения, написанную нами для тензоров. Так как оператор Set - встроенный, то прежде чем его менять, приходится выполнять команду Unprotect. Благодаря условию на второй аргумент встроенная операция Set остается неизменной для всех случаев, кроме того, когда второй ее аргумент определен нами как тензор. Сама операция, по сути, представляет собой модификацию описанной выше операции makeTensor:
In[194]:= Unprotect [Set] ;
(A_ = (B_) ?tensorQ) : =
(tensor[A] = True; type[A] = type[B];
inBasis[A] = inBasis[B];
val[A] = val[B]; If[skewQ[B], skew[A] = True];
If[symmQ[B], syram[A] = True];
If[crist[B] === True, crist[A] = True]);
Protect[Set];Возможности дальнейшего развития
Еще раз оговоримся, что наш пакет носит учебный характер. Наверняка он не свободен от погрешностей, которые будут устраняться в процессе тестирования пакета пользователями. Кроме того, в нем реализованы далеко не все функции и понятия, нужные для работы с тензорными полями. Написание, на основе имеющихся, некоторых естественных функций, таких как поднятие/опускание индекса на римановом многообразии, операция "звездочка Ходжа" на дифференциальных формах, вывод и численное решение уравнений параллельного переноса и уравнений геодезических на римано-вом (аффинном) многообразии, вычисление компонент тензора кривизны и тензора Риччи, а также полной кривизны, на наш взгляд, не вызовет принципиальных трудностей и предлагается вам в качестве задач.
Следующее расширение пакета должно предоставить пользователю возможность работать с несколькими векторными пространствами или многообразиями (пока состоящими из одной карты). Для этого можно было бы, например, завести еще один массив данных, скажем Object, в котором хранить информацию о заводимых объектах, а каждой из реализованных в пакете функций приписать необязательный аргумент - номер или имя объекта, на котором функция собирается применяться.
Следующий шаг - задание отображений между пространствами или многообразиями, реализация дифференциала отображения, переноса форм, проверки регулярности отображения и т. д.
И, наконец, создание настоящих многообразий, состоящих из нескольких карт, с помощью стыковки однокартовых многообразий при помощи функций перехода.
Реализация этой программы, конечно, потребует труда и времени.