【力学中的张量分析】- 力学与 Mathematica 系列 03

张量的推导与计算十分繁杂，因而使用Mathematica进行张量分析本身便是一件很自然的事，但很无奈，网上几乎没有相应的中文教程，与之相应的是，各种论坛上关于使用Mathematica进行张量分析的求助帖基本没有回应。因此，本帖系统地总结了一下Mathematica中各种张量函数的用法，并辅以部分例题，为大家展示如何使用Mathematica进行张量分析。限于篇幅，希望本帖能起到抛砖引玉的作用。

1. 张量的表示

在Mathematica中，使用多层列表表示张量。例如，二阶对称克罗内克符号为二阶张量，我们可以手动输入，或者借助内置函数产生。里奇-列维塔张量亦是如此。

{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}} // MatrixForm

Array[KroneckerDelta, {3, 3}] // MatrixForm

Table[Mod[(j - i) (k - j) (i - k), 3, -1], {i, 1, 3}, {j, 1, 3}, {k,

1, 3}] // MatrixForm

LeviCivitaTensor[3] // MatrixForm

实际上，Mathematica并不能直接区分协变张量和逆变张量，所以我们需要借助度量张量来实现。例如在极坐标下，我们已知某一矢量在自然协变基矢量下的逆变分量，求其在逆变基矢量下的协变分量。

Gx = CoordinateChartData["Polar", "Metric", {r, \[Theta]}];

Gn = CoordinateChartData["Polar", "InverseMetric", {r, \[Theta]}];

Gx // MatrixForm

Gn // MatrixForm

Pn = {1, r}

Px = Pn.Gn

Px.Gx

张量分析中，有很大一部分是对各种等式的证明，因此抽象表示张量很有必要，Mathematica也提供了这种表示方法，可以看到，这种表示方法可以指定任意维度，数据类型以及对称性等。关于张量等式的证明会在随后章节中给出。

$Assumptions =

A \[Element] Arrays[{4, 4, 4}, Reals, Symmetric[{1, 2, 3}]]

$Assumptions =

M \[Element] Matrices[{3, 3}, Reals, Antisymmetric[{1, 2}]]

$Assumptions = V \[Element] Vectors[3, Reals]

2.张量的代数运算

这里只给出张量分析中较为重要的几个代数运算，并积，缩并，内积，转置。

Clear["Global`*"];

(*并积*)

(*A\[CircleTimes]B=Subscript[A, ij]Subscript[B, mn]*)

A = {a[1], a[2], a[3]}; B = {b[1], b[2], b[3]};

TensorProduct[A, B] // MatrixForm

A\[TensorProduct]B // MatrixForm

Table[x[i, j], {i, 1, 2}, {j, 1, 2}]\[TensorProduct]Table[

y[i, j], {i, 1, 2}, {j, 1, 2}] // MatrixForm

张量的并积用TensorProduct函数来实现，或用"[TensorProduct]"符号来表示，并没有什么特殊的技巧。

Clear["Global`*"];

(*缩并*)

(*Subscript[A, ii]*)

A = Table[a[i, j], {i, 1, 3}, {j, 1, 3}];

A // MatrixForm

TensorContract[A, {{1, 2}}]

Tr[A]

张量的缩并用TensorContract函数来实现，可以具体指定缩并位置。对于矩阵来说，缩并后的结果便是矩阵的迹。

Clear["Global`*"];

(*点积*)

(*A.B=Subscript[A, ij]Subscript[B, jk]*)

A = Table[a[i, j], {i, 1, 2}, {j, 1, 2}];

B = Table[b[i, j], {i, 1, 2}, {j, 1, 2}];

Dot[A, B] // MatrixForm

A.B // MatrixForm

张量的内积用Dot函数来实现，或用"."符号表示，要注意内积是第一个张量最后一层和第二个张量第一层上的运算，因此一般是不可交换的。

Clear["Global`*"];

(*双点积*)

(*A:B=Subscript[A, ij]Subscript[B, ji]*)

A = Table[a[i, j], {i, 1, 2}, {j, 1, 2}];

B = Table[b[i, j], {i, 1, 2}, {j, 1, 2}];

TensorContract[(A\[TensorProduct]B), {{2, 3}, {1, 4}}]

张量的双点积可以使用并积函数与缩并函数一起实现，当然也可以是三点积，四点积等。例如这个例子中，双点积是定义为A矩阵和B矩阵在相邻部分的点积与远离部分的点积。

Clear["Global`*"];

(*转置*)

A = Table[i - j + k, {i, 1, 3}, {j, 1, 3}, {k, 1, 3}];

Transpose[A, {2, 1, 3}] // MatrixForm

Transpose[A] // MatrixForm

张量的转置使用Transpose函数来实现，默认参数下，在张量的前两层内转置，当然，设置参数后，Transpose函数是一种更广义的转置。

3.曲线坐标下的张量分析

本节的主要目的是推导不同坐标系下的张量方程，例如不同坐标系下弹性力学平衡方程，纳维-斯托克斯方程等。通常手动推导这类问题是相当麻烦的，比如教科书中附录会带有不同坐标系下的张量方程，但实际中，我们会遇到更多其他类似的问题却是无法查附录的。故使用Mathematica很有必要。下面先介绍梯度，散度的用法。

Clear["Global`*"];

A = Table[f[i][r, \[Theta]], {i, 1, 2}];

A // MatrixForm

Grad[A, {r, \[Theta]}, "Polar"]

Div[A, {r, \[Theta]}, "Polar"]

B = Table[f[i, j][x, y], {i, 1, 2}, {j, 1, 2}];

B // MatrixForm

Grad[B, {x, y}] // MatrixForm

Div[B, {x, y}] // MatrixForm

要注意两点，其一是这里的梯度散度均是在标准正交曲线坐标下的运算，即基矢量均为单位长度。其二是在Mathematica中，进行梯度散度计算时，增加的维度放在最外层，这意味着直接进行梯度散度运算时，是指右梯度和右散度，需要特别注意，当然，我们可以使用Transpose函数来将其化为左梯度和左散度，毕竟这才是我们最常用的。下面给出一个例子，各位读者可以好好体会一下。（C为二阶张量）

Clear["Global`*"];

c = {{Cxx[x, y], Cxy[x, y]}, {Cxy[x, y], Cyy[x, y]}};

u = {U[x, y], V[x, y]};

u.Transpose[Grad[c, {x, y}], {2, 3, 1}] // MatrixForm

(Grad[u, {x, y}]).c + c.Transpose[(Grad[u, {x, y}])] // MatrixForm

Div[Transpose@c, {x, y}] // MatrixForm

4.例题

1)证明：(a*b)*c=(a.c)b-(b.c)a

Clear["Global`*"];

$Assumptions = {a \[Element] Vectors[3, Reals],

b \[Element] Vectors[3, Reals], c \[Element] Vectors[3, Reals]};

TensorReduce[(a \[Cross]b)\[Cross]c == (a.c) b - (b.c) a]

Clear["Global`*"];

$Assumptions = True;

a = {a1, a2, a3};

b = {b1, b2, b3};

c = {c1, c2, c3};

(a \[Cross]b)\[Cross]c // Simplify // MatrixForm

(a.c) b - (b.c) a // Simplify // MatrixForm

Expand[(a \[Cross]b)\[Cross]c == (a.c) b - (b.c) a]

FullSimplify[(a \[Cross]b)\[Cross]c == (a.c) b - (b.c) a]

从上面可以看出，使用TensorReduce来证明等式本身就没多大意义，相当于验证一下等式，因此推荐写出张量分量，一步一步展示中间过程，最后作判断。这样才能对等式有更直观的了解。

2)推导球坐标系下的 不可压缩Navier-Stokes 方程（动量方程）。即：

Clear["Global`*"];

u = {ur[r, \[Theta], z], u\[Theta][r, \[Theta], z],

uz[r, \[Theta], z]};

ADV = u.Transpose[Grad[u, {r, \[Theta], z}, "Cylindrical"]];

Pre = Grad[p[r, \[Theta], z], {r, \[Theta], z}, "Cylindrical"];

Vis = Laplacian[u, {r, \[Theta], z}, "Cylindrical"];

ADV[[1]]

Pre[[1]]

Vis[[1]] // Simplify

考虑到最终完整的形式很复杂且不是重点，我们只计算了其中的微分算子，其结果跟教材一致。

3)设一点沿着半径为a的球面等速运动，运动起始位置在赤道上，速度v的方向与球的经线夹角[Alpha]保持常数，求点的轨迹（斜驶线）方程，以及到达球极点的时刻[Tau]。

在球面上建立坐标系，考虑到坐标系正交，可直接根据拉梅系数获得速度表达式。

Clear["Global`*"];

r = a;

x = r Sin[\[Theta]] Cos[\[CurlyPhi]];

y = r Sin[\[Theta]] Sin[\[CurlyPhi]];

z = r Cos[\[Theta]];

H\[CurlyPhi] =

D[{x, y, z}, \[CurlyPhi]].D[{x, y, z}, \[CurlyPhi]] // Sqrt //

Assuming[0 <= \[Theta] <= \[Pi] && a > 0, Simplify@#] &

H\[Theta] =

D[{x, y, z}, \[Theta]].D[{x, y, z}, \[Theta]] // Sqrt //

Assuming[0 <= \[Theta] <= \[Pi] && a > 0, Simplify@#] &

速度可以表示为，vφvφ=HφHφ˙φφ˙，vθvθ=HθHθ˙θθ˙，根据速度条件tan[Alpha]=-vθvθ/vφvφ,积分一次可以给出斜驶线方程，以及求出路径长度与时间。

d\[CurlyPhi] = -Tan[\[Alpha]] H\[Theta] d\[Theta] /(H\[CurlyPhi]);

ds = (H\[CurlyPhi]^2 d\[CurlyPhi]^2 + H\[Theta]^2 d\[Theta]^2) //

Sqrt // Assuming[

0 <= \[Theta] <= \[Pi] && a > 0 && d\[Theta] > 0 &&

0 < \[Alpha] < \[Pi]/2, Simplify@#] &;

s = Integrate[-ds/d\[Theta], {\[Theta], \[Pi]/2, 0}]

\[Tau] = s/v

我们可以借助Mathematica来直观感受一下斜驶线。

Clear["Global`*"];

r = 1;

\[Alpha] = \[Pi]/4;

\[Theta] = 2 ArcTan[E^(-Cot[\[Alpha]] \[CurlyPhi])];

x = r Sin[\[Theta]] Cos[\[CurlyPhi]];

y = r Sin[\[Theta]] Sin[\[CurlyPhi]];

z = r Cos[\[Theta]];

p1 = ParametricPlot3D[{ Sin[\[Theta]1] Cos[\[CurlyPhi]1],

Sin[\[Theta]1] Sin[\[CurlyPhi]1], Cos[\[Theta]1]}, {\[CurlyPhi]1,

0, 2 \[Pi]}, {\[Theta]1, 0, \[Pi]}, BoxRatios -> {1, 1, 1},

ColorFunction -> {Green, Green},

Mesh -> None,

AxesOrigin -> {0, 0, 0},

Boxed -> False, Ticks -> {{-1, 1}, {-1, 1}, {-1, 1}},

ViewPoint -> {8, -4, -6}, PlotPoints -> 80,

ImageSize -> 500];

p2 = ParametricPlot3D[{x, y, -z}, {\[CurlyPhi], 0, 10},

PlotStyle -> {Thick, Blue}];

p3 = ParametricPlot3D[{ Cos[\[Theta]1], Sin[\[Theta]1],

0}, {\[Theta]1, 0, 2 \[Pi]}, PlotStyle -> {Thick, Yellow}];

Show[p1, p2, p3]

到此，如果以上函数都比较熟悉了，那么其实就可以尝试进行各种张量分析了，比如定义协变导数，张量微积分。更进一步，微分几何与相对论的东西也可以计算了。限于篇幅，本帖不可能面面俱到。如果读者有什么疑问或者想求解的问题，可以在下方留言，下期文章会整理做出解答。

赞 (0)

相关推荐