Тензоры: опыт создания пользовательского пакета программ

Функции для работы с таблицами

Как мы уже выясняли в прошлом семестре, в Mathematica имеется ряд команд, манипулирующих списками. К наиболее популярным из них относятся команды, имена которых начинаются с Map: Map, MapThread, MapIndexed и др. Команда Map настолько широко используется, что для нее была придумана специальная форма, а именно, f /@ list (этой формой записи команды Map мы будем пользоваться многократно):

In[26]:=Map[f [#] , {{а, b}, {с, d}}] 
           f[#] /@{{a, b}, {с, d}} 
           MapThread[f[#l, #2, #3] ,
            {{{a}, {b}}, {{c}, {d}}, {{e}, {f}}}]
               (* #l соответствует элементам
               первого списка, #2 - второго и т.д.*) 
             MapIndexed[f [#1] [#2] , {{а, b} , {с, d} } ] 
               (* #2 соответствует номеру текущего
                 элемента списка взятого в фигурные скобки *)
Out[26]={f [{a, b}], f [{с, d}]} 
Out[27]={f [{a, b}], f [{с, d}]} 
Out[28]={f[{a}, {с}, {е}], f[{b}, {d}, {f}]} 
Out[29]={f [{a, b}] [{1}], f[{c, d}] [{2}]}

Эти команды допускают явное указание уровня списка, к элементам которого и нужно применять функцию:

In[30]:=Map[f [#] , {{а, b} , {с, d} } , {2}] 
           MapThread [ f [ #l, #2, #3] , 
                 {{{a}, {b}}, {{c}, {d}}, 
               {{e}, {f}}}, 2] 
             (* #l соответствует элементам 
               первого списка, #2 - второго и т.д.*) 
             MapIndexedff[#1][#2] ,
                 {{а, b}, {с, d}}, {2}] 
           (* #2 соответствует номеру 
                 текущего элемента списка, 
               взятого в фигурные скобки *)
Out[30]={{f [a], f [b]}, {f [с], f [d]}}
Out[3l]={{f [а, с, e]}, {f[b, d, f ]}}
Out[32] = 
         {{f[a] [{1, 1}], f[b] [{1, 2}] }, {f[c] [{2, 1}], f[d] [{2, 2}]}}

Заметьте, что теперь номер текущего элемента является списком длины 2. Такая общность интерпретации команды MapIndexed и объясняет, почему номер элемента берется в фигурные скобки.

Следующая команда Tuples строит или n -кратное декартово произведение списка на себя:

In[33] := Tuples [{а, b, с}, 3]
Out[33] =
       {{а, а, а},{а,а, b},{а, а, с},{a, b, a},{a,b, b}, {а, b, с},
        {а, с, а},{а,с, b},{а, с, с},{b, a, a},{b,a, b}, {b, а, с},
        {b, b, а},{b,b, b},b, b, с},{b, с, а},{b,с, b},
        {b, с, с},{с,а, а},{с, а, b},{с, а, с},{с,b, а},
        {с, b, b},{с,b, с},{с, с, а},{с, с, b},{с,с, с}}

или, более общо, декартово произведение списков:

In[34]:= Tuples [{{a1, a2 } , {b1, b2 , b3}, {c1, c2}}]
Out[34] = {{a1, b1, c1}, {a1, b1, c2 }, {a1, b2 , c1}, {a1, b2, c2 }, 
           {a1, b3, c1}, {a1, b3, c2}, {a2, b1, c1}, {a2, b1, c2}, 
           {a2, b2, c1}, {a2, b2, c2}, {a2, b3, c1}, {a2, b3, c2}}

Наконец, еще две очень популярные команды: Outer и Inner, являющиеся обобщениями тензорного произведения и свертки тензоров (в частности, скалярного произведения).

Команда Outer[f,list1,list2,...] применяет произвольную функцию f к всевозможным комбинациям элементов списков list1,list2,...:

Команда Inner[f,list1,list2,g] обобщает команду Dot[] (или ".")

In[36]:={a, b, с} . {х, у, z}
Out[36] =ax+by+cz

причем f играет роль умножения, а g - сложения:

In[37] := Inner [Min, {а, b, с}, {х, у, z}, Max]
          Inner [CircleTimes, {a, b, с}, {х, у, z}, Plus]
Out[37]=Max [Min [a, x] , Min[b, у], Min [с, z] ]

Тензорные поля

В курсе дифференциальной геометрии, тензоры используются для изучения многообразий и порождают тензорные поля, т. е. тензоры одного и того же типа, заданные на касательных пространствах к многообразию. Теперь базисы в этих линейных пространствах выбираются исходя из локальных координат, а тензорный закон, как следствие, использует матрицу Якоби перехода от одной системы координат к другой. Таким образом, чтобы получить координатное представление тензорного поля, нужно вместо базисов еь ..., е" брать локальные координаты $$x^1, \dots , x^n$$, вместо чисел $$T_{j_1 \dots j_q}^{i_1 \dots i_p}$$ - функции $$T_{j_1 \dots j_q}^{i_1 \dots i_p}(x^1, \dots , x^n) $$, а тензорный закон писать с использованием матриц Якоби.

Пример создания пакета для работы с тензорами

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

In[39] =Needs ["Tensors " " , NotebookDirectory [ ] <> "tensors.m"]

где второй аргумент - это имя файла, в котором пакет содержится. (Здесь мы предполагаем, что пакетный файл tensors.m расположен в той же папке, что и файл, с которым вы сейчас работаете, причем ваш файл или был открыт, уже находясь в указанной папке, или был сохранен в ней, иначе значение функции NotebookDirectory[] не определено.)

Правила работы в пакете Tensors

В данном разделе мы опишем возможности текущей версии пакета. Разумеется, их можно и нужно расширять.

Начало работы, создание базиса и тензоров, тензорные операции.

В пакете реализованы возможности работы как с алгебраическими тензорами (т. е. тензорами в линейном пространстве), так и с тензорными полями (в последнем случае, в текущей версии пакета, тензорные поля рассматриваются на многообразии, состоящем из одной карты, т. е. на области в $$R^n$$ ). Чтобы начать работу в пакете, следует задать или базис линейного пространства (функция makeBasis ), или локальные координаты и соответствующий канонический базис (функция makeBasisDif). Чтобы посмотреть формат интересующей вас функции пакета, следует выполнить команду ?<имя функции>. Например:

In40] := ? makeBasis

Зададим базис в трехмерном пространстве.

In[41] :=makeBasis [e, 3]

Все создаваемые базисы последовательно нумеруются. Номер текущего базиса записывается в переменную curBas, а сам базис - в список базисов basis:

In[42] := curBas
          basis[curBas] 
          val[e] 
          basis[1]
Out[42] = 1

Отметим, что на самом деле команда создает также двойственный базис. Обратите внимание, что мы пишем индексы непосредственно над и под буквами (а не справа вверху и справа внизу). Это объясняется тем, что правый верхний индекс Mathematica воспринимает как степень. Набирать такие индексы удобно с помощью сочетаний клавиш Ctrl+= для нижнего индекса и Ctrl+7 для верхнего индекса. Если же требуется поставить одновременно верхний и нижний индексы, то следует набрать одновременно Ctrl+= +5.

После того как базис фиксирован, можно создать тензор любого типа. При этом предусмотрена возможность задать абстрактный тензор с неопределенными компонентами ( makeTensor ) и возможность задать тензор с конкретными значениями компонент ( makeTensorTab ):

In[46] =makeTensor [T, 1, 1] ; makeTensorTab [W, {2, b, с}, 1, 0] ;

Чтобы посмотреть значения компонент тензора, надо выполнить команду val:

In[47] : = val[T] 
           val[W]

Out[48] = {2, b, c}

Над заданными тензорами можно выполнять тензорные операции. Реализованы линейная комбинация, тензорное произведение, свертка, перестановка индексов, симметрирование, альтернирование.

Значок тензорного произведения можно получить, набрав Esc c* Esc :

Линейная комбинация допустима только для тензоров одного типа:

Результат применения тензорных операций можно присвоить с помощью обычного оператора Set (т. е. =):

Свертка записывается просто указанием номеров индексов, по которым надо сворачивать:

Перестановка $$\tau$$ индексов указывается в виде списка, элементы которого суть $$\tau $$ -образы последовательных натуральных чисел начиная с 1. При этом список заключается в дополнительные фигурные скобки и пишется сверху или снизу, в зависимости от того, какие индексы мы собираемся переставлять:

Аналогично выглядят команды альтернирования и симметрирования:

Также часто удобно вывести конкретную компоненту тензора. Это делается с помощью явного указания списков соответствующих индексов:

Аналогично, имеется возможность породить сразу симметричный или кососимметричный тензор типа (0,q). Для этого используются команды makeSymm, makeSkew, makeSymmTab, makeSkewTab. Для симметричных и ко сосимметричных тензоров определены операции симметричного произведения $$\circ$$ ( SmallCircle, Esc sc Esc ) и, соответственно, кососимметричного произведения $$\wedge$$ ( Wedge, Esc ^ Esc ). Пакет относит к симметричным (кососимметричным) тензорам или тензоры, созданные перечисленными четырьмя операциями, или полученные в результате симметрирования или альтернирования тензоров типа (0,q) :

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

In[71]:=vall[bb]//MatrixForm

Запись тензора в тензорном базисе, вычисление значения тензора на наборе ковекторов и векторов0. Имеется возможность записать созданный тензор в базисе пространства тензоров данного типа:

In[72] : = makeTensor [T, 1, 2]; 
           toBasis[T]

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

In[76]:=val[B]

Предусмотрен также набор обратных операций, вычисляющих список компонент тензора, заданного в виде линейной комбинации базисных. Это функции toList, toListSkew, toListSymm:

In[78]:=toBasis[T]

Out[81]=basis[5]

Напомним, что тензор типа (p,q) можно представлять себе как полилинейное отображение, заданное на прямом произведении p экземпляров двойственного пространства и q экземпляров самого пространства. Чтобы вычислить значение такого отображения, в пакете имеется функция apply :

In[86]:={a,b,c}

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

Перезапуск. Чтобы отменить все созданные базисы, координаты и пр. следует выполнить команду resetTensors:

In[92] =basis [curBas] 
        resetTensors 
        basis[curBas]

Пакет Tensors перезапущен. Результаты вычислений аннулированы. 
Out[94] = basis [0]

Особенности работы с локальными координатами и тензорными полями. Если вы задали локальные координаты и соответствующий базис (команда makeBasisDif ), то появляется возможность задавать тензорные поля, компоненты которых суть абстрактные функции от координат. Соответствующий набор функций - это makeTensorField, makeSkewField и makeSymmField . Результат их работы - соответствующее тензорное поле, к которому можно применять обычные тензорные операции:

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

Приведем соответствующие примеры. Чтобы удобнее было контролировать вычисления, пакет создает ориентированный граф graphBas замен координат (базисов). Вершины этого графа - системы координат. Две вершины соединены стрелкой, если соответствующая замена координат задана явно. Этот граф удобно вывести динамически. В приведенной ниже программе предусмотрен переход к координатам, выбранным мышкой. А именно, чтобы сделать замену, нужно кликнуть в соответствующую вершину графа:

In[99] :=
      (*Создаем евклидовы координаты на плоскости*) 
     resetTensors; 
     makeBasisDif[{х, у}]; 
      (*Задаем в них постоянную билинейную форму и произвольный вектор*)

makeTensorTab[W, {а, b} , 1, 0] ;
         (*Задаем полярные координаты. При этом Записываем функции 
           перехода от них к текущим координатам, 
         т.е. функции х[r, ϕ] и у[r, ϕ]. При этом текущими координатами
            становятся полярные *) 
        makeBasisDif [ {r, ϕ , {r Cos [ϕ] , r Sin [ϕ] } ] ;
In[104] :=
       (*Программа, динамически выводящая граф замен*) 
     DynamicModule[{}, 
        Dynamic[Deploy[GraphPlot[graphBas, DirectedEdges -> True,
             VertexRenderingFunction -> 
                ({White, EdgeForm[Black], EventHandler[Disk[#, .15], 
                     { "MouseDown" : > (changeBas [curBas , #2] ) } ] , 
                  If [#2 == curBas, Red, Black],
                  EventHandler[Style[Text[Column[{coor[#2] , #2} , Center] , #l] , 
                       FontSize -> Scaled [0.03] ] , 
                     {"MouseDown" :> (changeBas[curBas, #2])}]} ), 
              ImageSize -> {300, 350}] ] ], 
        Initialization :-> Needs ["Tensors" , NotebookDirectory [] <> "tensors .m"] ]

In[105] := Dynamic [ {val [W] , apply[G, W, W] } // Simplify]
Out[105]= {val [W] , apply [G, W, W] }
In[106]:=(* Задаем еще одни координаты (u,v) . Функции перехода 
                можно задать к любым имеющимся координатам, 
             указав явно их номер. В следующем примере мы задаем
              функции перехода к декартовым координатам, 
             т.е. функции x[u,v] и y[u,v]*) 
             makeBasisDif [{u, v} , 1, {и^2 - v^2, 2 u v} ] ; 
              (*A здесь - функции перехода к полярным координатам
                  r[p,q] и ϕ [p,q]*) 
             makeBasisDif [{p, q} , 2, {р1/3 , q\3}];

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

In[108] := (*3адаем переход от 4х координат к полярным (вторым)*) 
             makeTransDif[4, 2, {r3, 3ϕ}]; 
          (*И еще переходы*)
             makeTransDif [3, 2, {Sqrt [r] Cos [ϕ / 2] , Sqrt [r] Sin [ϕ / 2] } ] ; 
             makeTransDif[2, 1, {Sqrt[х^2+у^2] , ArcTan[y/ x] } ] ;

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

Получающиеся выражения не всегда легко упрощаются. Это объективно связано с тем, что функции перехода определены не всюду. В пакете предусмотрена возможность добавлять предположения о значениях тех или иных переменных в список предположений $Assumptions, которым Mathematica пользуется при упрощениях, или убирать их оттуда. Соответствующие команды - это addAssumptions и clearAssumptions . Отметим, что предположение о вещественности заводимых координат делается автоматически.

In[111] := $Assumptions
           addAs sumptions [p>0, 0<q<π, r>0, -π / 2 < ϕ < π / 2] ;
           $As sumptions
           clearAssumptions [0 < q < π] ;
           $Assumptions
           addAs sumptions [p > 0, 0 < q < π / 6, u2-v2>0, u v > o] ;
           $As sumptions
Out[111] =
        (x | y) º Reals  (r | ϕ) º Reals  (u | v) º Reals  (p | q) º Reals
Out[113]= (p | q) º Reals  (r | ϕ) º Reals  (u | v) º Reals 
               (x | у) º Reals  p > 0  r > 0  0 < q < π  - π/2 < ϕ < π/2
Out[115]= (p | q) º Reals  (r | ϕ) º Reals  (u | v) º Reals  
               (x | y) º Reals  p > 0   r > 0  - π/2 < ϕ < π/2
Out[117] =
       (p | q) º Reals  (r | ϕ) º Reals  (u | v) º Reals  (x | у) º Reals  
      p > 0  r > 0  uv > 0  u2-v2 > 0  0 < q < π/6  -π/2 < ϕ < π/2

Отметим, что сама команда замены координат имеет вид

changeBas[<номер исходного базиса>,<номер целевого базиса>]

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

Дифференциальные формы: их разложение по базису, внешнее умножение и дифференцирование. В настоящий момент в пакете реализован ряд возможностей для работы с дифференциальными формами. Как уже упоминалось выше, можно задать кососимметрический тензор или тензорное поле типа (0,q) с помощью команд makeSkew, makeSkewTab или makeSkewField:

Дифференциальную форму можно разложить по базису. Дифференциальные формы также можно перемножить внешним образом (как описанные выше кососимметрические тензоры):

Можно задать форму с коэффициентами-функциями:

От такой формы можно взять внешний дифференциал, который, как и положено, будет дифференциальной формой:

Можно задать форму непосредственно, как вот такую форму P:

От нее тоже можно взять внешний дифференциал:

In[130]:=d[P]

однако ни P, ни d[P], в этом случае не являются тензорами:

In[131] := tensor [P] 
           skew[P]
           tensor[d[P] ] 
           tensorBas[P] 
           skewBas[P] 
           tensorBas[d[P]] 
           skewBas[d[P]]
Out[131] = False
Out[132] = False
Out[133] = False
Out[134l = True
Out[135] = True
Out[136] = True
Out[137] = True

Превратить их в тензорное поле можно так:

In[138]:=makeSkewTab[dP, toListSkew [d[P] ] ] 
         tensor[dP] 
            dP
         {1,2,3,4}
Out[139] = True

Риманова метрика, риманова и аффинная связности. В пакете предусмотрена возможность объявить рассматриваемое многообразие (в текущей версии, напомним, многообразие состоит из одной карты) римановым. При этом задается или

В обоих случаях для метрики и обратного тензора зарезервированы символы g и ig :

In[141] := resetTensors ;
           makeBasisDif[{х, у}]; 
           makeRiemann;

       Пакет Tensors перезапущен. Результаты вычислений аннулированы.

            Созданы метрические тензоры д типа (0,2) и ig типа (2,0)

In[144] :=Row [ {val [g] //MatrixForm, val[ig] //MatrixForm} , 
          Spacer[10]]

Отметим, что оба построенных тензора по определению симметричны. Аналогично:

In[145] := resetTensors ;
           makeBasisDif[{х, у}] ;

val[g] //MatrixForm

     Пакет Tensors перезапущен. Результаты вычислений аннулированы.

            Создан метрический тензор g типа
                 (0,2) и обратный к нему тензор ig типа (2,0)

При этом, конечно, метрика - это тензор, компоненты которого пересчитывают при замене координат. Так, например, можно вычислить компоненты евклидовой метрики в полярных координатах:

Если задана метрика, то можно вычислить символы Кристоффеля соответствующей симметричной римановой связности, что делается с помощью команды ? (большая греческая буква гамма, Esc G Esc ). При этом с символами Кристоффеля во многом можно обращатьсякак с тензорами типа (1,2), а именно, можно посмотреть конкретную компоненту, можно вычислять их значение в текущих координатах с помощью все той же функции val :

In[151] := resetTensors ;
           makeBasisDif[{х, у, z}];

makeCrRiemann 
            val[Г]

   Пакет Tensors перезапущен. Результаты вычислений аннулированы.

         Создан метрический тензор g типа
                 (0,2) и обратный к нему тензор ig типа (2,0)

Out[155]={{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}, 
          {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}, 
          {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}}

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

Или вот конкретная компонента:

Кроме того, предусмотрена возможность задавать абстрактную аффинную связность. При этом соответствующие функции

makeCrAffine[]

и

makeCrAffineTab[<таблица значений>>]

можно использовать как без явного указания имени связности (в этом случае ей присваивается имя Г), так и с указанием этого имени в виде последнего аргумента:

In[159] =makeCrAf f ine [ ] ; 
         val[r]

In[161] :=makeCrAf fine [AA] ; 
          val[AA]

Наконец, имеется возможность задавать метрику как тензорное поле. В этом случае компоненты метрики суть абстрактные функции от координат. Аналогичная возможность предусмотрена и для абстрактной аффинной связности. Соответствующие функции - это makeRiemannField , которая создает поля g и ig , makeCrAffineField[], которая создает связность Г , и makeCrAffineField[<имя связности>] , которая создает связность с заданным именем:

In[167] := resetTensors ;
           makeBasisDif[x, 2] ; 
           makeRiemannField; 
           makeCrRiemann;

Пакет Tensors перезапущен. Результаты вычислений аннулированы.

              Созданы метрические тензорные
                поля д типа (0,2) и ig типа (2,0)

Как работает пакет (краткое описание основных идей)

В этом разделе мы кратко опишем основные идеи и технические приемы, положенные в основу написания пакета. Сами тексты программ находятся в файле tensors.m. Этот файл тоже можно открыть в Mathematica (откройте). Возникнет специальным образом оформленное окно работы с пакетными файлами. Как и любой пакет, наш начинается оператором BeginPackage["Tensors'", {"Combinatorica'","GraphUtilities'"}], где первый аргумент - это имя пакета, а второй - список используемых пакетов, и заканчивается оператором EndPackage[].

В нашем пакете всего две клетки, первая состоит из команд вида <имя функции>::usage="текст описания.", служащих для описания пользовательских функций. Именно этот текст появляется на экране при выполнении команды ?<имя функции> . Если имеется несколько одноименных функций, то будут выведены все описания. Наконец, по команде ??<имя функции> на экран выводится весь код, задающий функцию:

In[172] := ? makeBasis

In[173] := ?? makeTensor

makeTensor [Т , р Integer, q Integer] : =

ConstantArray[nn[basis[curBas]], p + q] ;
        tensor[T] = True; type[T] = {p, q}; inBasis[T] = curBas; 
        If[TrueQ[symm[T]], (symm[T]) = . ] ; 
        If[TrueQ[skew[T] ] , (skew[T]) =.];
        If[p=0q <= l, symm[T] = True; skew[T] = True] ;)

Полный список функций пакета можно увидеть, если кликнуть в кнопку Functions в окне работы с файлом пакета.

Создание объекта с помощью описания его свойств и правил. Посмотрим, как реализован в пакете объект 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] происходит следующее:

tensor[T] = True;type[T] = {p, q}; inBasis[T] = curBas;

val[T] = <таблица значений тензора в текущей системе координат - просто символы с индексами в нужном месте и числе>;

Кроме того, написаны аналоги функций 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, в котором хранить информацию о заводимых объектах, а каждой из реализованных в пакете функций приписать необязательный аргумент - номер или имя объекта, на котором функция собирается применяться.

Следующий шаг - задание отображений между пространствами или многообразиями, реализация дифференциала отображения, переноса форм, проверки регулярности отображения и т. д.

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

Реализация этой программы, конечно, потребует труда и времени.