从张量迹到行列式：用Python验证连续介质力学中的关键性质

连续介质力学中的张量运算常常让初学者望而生畏。那些复杂的数学符号和抽象概念，如何在实践中验证它们的性质？本文将带你用Python的NumPy和SymPy库，通过代码实现张量迹、行列式、逆矩阵和正交性等核心概念的验证，让理论变得触手可及。

1. 环境准备与基础概念

在开始之前，我们需要准备好Python环境。推荐使用Jupyter Notebook进行交互式编程，这样可以即时看到每一步的计算结果。

import numpy as np
import sympy as sp
from IPython.display import display, Math

张量在连续介质力学中扮演着重要角色。二阶张量可以表示为3×3矩阵，而张量的迹(Tr)定义为矩阵对角线元素之和。让我们先创建一个随机张量：

A = np.random.rand(3,3)
print("随机张量A:\n", A)
print("张量A的迹:", np.trace(A))

注意：NumPy的trace函数直接计算矩阵的迹，这与数学定义Tr(A)=A₁₁+A₂₂+A₃₃完全一致。

2. 张量迹的性质验证

张量迹具有几个重要性质，我们可以通过编程来验证这些数学定理。

2.1 迹的不变性

张量的迹是一个不变量，与坐标系选择无关。我们可以通过坐标变换来验证这一点：

# 创建一个随机正交矩阵Q
Q,_ = np.linalg.qr(np.random.rand(3,3))
A_transformed = Q @ A @ Q.T  # 相似变换

print("原始迹:", np.trace(A))
print("变换后迹:", np.trace(A_transformed))

2.2 迹的线性性质

迹运算具有线性性质：Tr(A+B)=Tr(A)+Tr(B)。我们可以轻松验证：

B = np.random.rand(3,3)
print("Tr(A+B):", np.trace(A+B))
print("Tr(A)+Tr(B):", np.trace(A)+np.trace(B))

2.3 迹与矩阵乘法的关系

对于矩阵乘法，Tr(AB)=Tr(BA)。让我们用代码验证：

print("Tr(AB):", np.trace(A @ B))
print("Tr(BA):", np.trace(B @ A))

更一般地，对于三个矩阵的乘积，循环置换不改变迹：

C = np.random.rand(3,3)
print("Tr(ABC):", np.trace(A @ B @ C))
print("Tr(BCA):", np.trace(B @ C @ A))
print("Tr(CAB):", np.trace(C @ A @ B))

3. 张量的行列式及其性质

行列式是张量的另一个重要不变量，在连续介质力学中常用于体积变化的分析。

3.1 行列式计算

NumPy提供了直接计算行列式的函数：

det_A = np.linalg.det(A)
print("det(A):", det_A)

对于符号计算，我们可以使用SymPy：

a11, a12, a13, a21, a22, a23, a31, a32, a33 = sp.symbols('a11 a12 a13 a21 a22 a23 a31 a32 a33')
A_sym = sp.Matrix([[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]])
display(Math(r'\det(A) = ' + sp.latex(A_sym.det())))

3.2 行列式的性质验证

行列式具有以下重要性质，我们可以通过编程验证：

1. 乘积的行列式 ：det(AB)=det(A)det(B)

print("det(AB):", np.linalg.det(A @ B))
print("det(A)*det(B):", np.linalg.det(A)*np.linalg.det(B))

2. 转置的行列式 ：det(Aᵀ)=det(A)

print("det(A):", np.linalg.det(A))
print("det(A.T):", np.linalg.det(A.T))

3. 逆矩阵的行列式 ：det(A⁻¹)=1/det(A)

A_inv = np.linalg.inv(A)
print("det(A⁻¹):", np.linalg.det(A_inv))
print("1/det(A):", 1/np.linalg.det(A))

3.3 特殊情况的验证

当行列式为零时，矩阵是奇异的。我们可以构造一个秩亏矩阵来验证：

# 构造一个秩为2的矩阵
C = np.array([[1,2,3], [4,5,6], [7,8,9]])
print("det(C):", np.linalg.det(C))  # 应该接近0

4. 张量的逆与正交性

4.1 张量逆的计算

在连续介质力学中，张量的逆经常出现在本构关系中。我们可以用两种方法计算逆矩阵：

1. 直接计算法 ：

A_inv_np = np.linalg.inv(A)
print("NumPy逆矩阵:\n", A_inv_np)

2. 伴随矩阵法 （验证数学公式A⁻¹=adj(A)/det(A)）：

def adjugate(matrix):
    cofactor = np.zeros((3,3))
    for i in range(3):
        for j in range(3):
            minor = np.delete(np.delete(matrix, i, axis=0), j, axis=1)
            cofactor[i,j] = (-1)**(i+j) * np.linalg.det(minor)
    return cofactor.T

A_adj = adjugate(A)
A_inv_adj = A_adj / np.linalg.det(A)
print("伴随矩阵法逆矩阵:\n", A_inv_adj)

4.2 正交张量的验证

正交张量在连续介质力学中描述刚体旋转，满足QᵀQ=I。我们可以验证一个矩阵是否正交：

def is_orthogonal(Q, tol=1e-6):
    return np.allclose(Q @ Q.T, np.eye(3), atol=tol)

# 创建一个旋转矩阵（绕z轴旋转θ角度）
theta = np.pi/4
Q = np.array([[np.cos(theta), -np.sin(theta), 0],
              [np.sin(theta), np.cos(theta), 0],
              [0, 0, 1]])

print("Q是正交矩阵吗?", is_orthogonal(Q))
print("det(Q):", np.linalg.det(Q))  # 旋转矩阵行列式为+1

4.3 正交变换的性质

正交变换保持向量长度和内积不变，这是刚体运动的重要特性：

v1 = np.array([1, 0, 0])
v2 = np.array([0, 1, 0])

# 变换后的向量
v1_t = Q @ v1
v2_t = Q @ v2

print("原始内积:", v1 @ v2)
print("变换后内积:", v1_t @ v2_t)
print("原始长度:", np.linalg.norm(v1))
print("变换后长度:", np.linalg.norm(v1_t))

5. 张量的特征值与正定性

在连续介质力学中，正定张量常出现在应力应变关系中。

5.1 特征值计算

# 创建一个对称正定矩阵
C = A @ A.T  # 对于任意实矩阵A，AᵀA和AAᵀ都是对称半正定的

eigvals = np.linalg.eigvals(C)
print("C的特征值:", eigvals)

5.2 正定性验证

正定矩阵的所有特征值都大于零：

def is_positive_definite(matrix):
    if not np.allclose(matrix, matrix.T):
        return False
    eigvals = np.linalg.eigvals(matrix)
    return np.all(eigvals > 0)

print("C是正定矩阵吗?", is_positive_definite(C))

5.3 张量的加性分解

任意二阶张量可以分解为球量和偏量部分：

def spherical_deviatoric_decomposition(S):
    sph = np.trace(S)/3 * np.eye(3)
    dev = S - sph
    return sph, dev

S_sph, S_dev = spherical_deviatoric_decomposition(A)
print("球量部分:\n", S_sph)
print("偏量部分:\n", S_dev)
print("偏量迹:", np.trace(S_dev))  # 应该为0

6. 可视化张量变换

为了更直观理解张量变换，我们可以可视化向量在变换前后的变化。

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 创建单位球面上的点
theta = np.linspace(0, 2*np.pi, 30)
phi = np.linspace(0, np.pi, 30)
x = np.outer(np.cos(theta), np.sin(phi))
y = np.outer(np.sin(theta), np.sin(phi))
z = np.outer(np.ones(30), np.cos(phi))

# 应用张量变换
def transform_surface(X, Y, Z, T):
    shape = X.shape
    coords = np.vstack([X.ravel(), Y.ravel(), Z.ravel()])
    transformed = T @ coords
    return transformed[0].reshape(shape), transformed[1].reshape(shape), transformed[2].reshape(shape)

# 绘制原始和变换后的球面
fig = plt.figure(figsize=(10, 5))
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_surface(x, y, z, alpha=0.5)
ax1.set_title('单位球面')

ax2 = fig.add_subplot(122, projection='3d')
x_t, y_t, z_t = transform_surface(x, y, z, A)
ax2.plot_surface(x_t, y_t, z_t, alpha=0.5)
ax2.set_title('变换后的椭球面')

plt.tight_layout()
plt.show()

7. 实际应用案例

让我们考虑一个连续介质力学中的实际例子：变形梯度张量F的性质验证。

# 假设一个简单的变形梯度
F = np.array([[1.1, 0.2, 0], 
              [0.05, 0.9, 0], 
              [0, 0, 1]])

# 右Cauchy-Green张量
C = F.T @ F
# 左Cauchy-Green张量
b = F @ F.T

print("det(F):", np.linalg.det(F))  # 体积变化率
print("C的特征值:", np.linalg.eigvals(C))
print("b的特征值:", np.linalg.eigvals(b))

提示：在有限变形理论中，C和b的特征值反映了材料在不同方向上的拉伸比。

通过这种实践性的编程验证，抽象的数学概念变得具体而直观。在Jupyter Notebook中尝试修改这些代码，观察不同张量变换的效果，是理解连续介质力学张量运算的最佳方式之一。