弹性力学课后题精解：张量指标记法与Python残差校验

简介：面向力学、土木、机械等专业本科生及考研复习者，《弹性力学基础》课后习题解答以同济大学程尧舜版教材第二章为范围，逐题给出推导过程，适合课堂同步练习与考前查漏补缺。压缩包内仅1个PDF文件，1.82MB，按题号顺序排布，公式推导与文字说明完整，可打印或分屏对照教材使用。目前已有677人学习下载。

内容围绕向量与张量运算展开：从偏微分与张量乘积的计算、对称性条件下的恒等式证明，到三矢量、四矢量点乘叉乘关系的推导，均给出完整步骤；坐标变换部分结合绕z轴旋转示意图，说明矢量新旧分量及二阶张量T各分量的变换系数求法。此外还涉及张量阶数的判定、迹与单位张量点积的证明、矢量与二阶张量叉积的转置关系，以及二阶张量的对称与反对称分解、反对称部分的轴向矢量、特征值与特征矢量求解。对张量运算与坐标变换这两处难点，中间步骤较为完整。

1. 弹性力学课后题真正卡人的地方，是 δ 和 e 的指标配对

翻到《弹性力学基础》第二章的习题，题面第一行是 δ_pi δ_iq δ_qj δ_jk，第二行是 e_pqi e_ijk A_jk，参考答案一行就跳过去了。程尧舜这本同济版教材的课后题解答覆盖第二到第六章：张量代数与指标记法、坐标变换、应变张量与协调方程、应力张量与面力边界条件、线性各向同性本构，以及位移法和应力法两类边值问题。正在跟教材自学的本科生、考研复习的人都用得上，真正容易翻车的是三件事：哑标与自由标的配对顺序、绕 z 轴旋转时方向余弦矩阵的符号约定、协调方程该代哪一套形式。后面按这个顺序拆：先把 δ 和 e 的化简做成能跑的脚本，再处理坐标变换，然后分别沿应变和应力两条线往下走，最后收在一套残差校验的写法上。

2. δ 与 e 的化简：把张量恒等式写成可执行的验证脚本

第二章的题面全是哑标，答案只有一行，跳步最狠。硬啃的办法是把两条母公式背熟，再用脚本把每道题的结果跑一遍——只要残差是零，就说明指标配对没错，比对着残缺的变量名去猜要快得多。

2.1 两条母公式撑起 2.1 题

δ 的收缩和 e 与 δ 的乘积恒等式是所有化简的来源，先把它们列清楚：

2.1(1) 的思路是「相邻哑标先缩」：δ_pi δ_iq = δ_pq，δ_qj δ_jk = δ_qk，剩下 δ_pq δ_qk = δ_pk。这类题只要遵守「同一项里哑标必须成对、每对消掉后维度降 1」，眼睛扫一遍就能出结果。

2.1(2) 的关键是先算 e 与 e 的收缩。把两个 e 的公共指标按位置对齐后套母公式，得到 e_pqi e_ijk = δ_pj δ_qk − δ_pk δ_qj，再与 A_jk 收缩，两项分别给出 A_pq 和 A_qp，所以结果是 A_pq − A_qp。这个量本质上是 A 的反对称部分的两倍，后面 2.2、2.9 两题都在反复用它。

2.1(3) 同理，e_ijp e_klp 收缩成 δ_ik δ_jl − δ_il δ_jk 后与 B_ki B_lj 收缩，分别是 B_ii B_jj 和 B_ij B_ji，也就是 (tr B)² − tr(B²)。这个组合在第四章写应力不变量、第六章判 Beltrami-Michell 方程时还会再遇到。

提示：e 的指标顺序决定符号。e_pqi 和 e_piq 差一个负号，套公式前先确认题面里两个 e 的公共指标在第几位，否则结果只差符号、很难察觉。

2.2 对称性消去：2.2 题与 2.10 题的反对称部分

2.2 题给的是 a_ij = a_ji，要证 e_ijk a_jk = 0。做法是把哑标 j、k 换个名字：e_ijk a_jk = e_ikj a_kj，而 e_ikj = −e_ijk、a_kj = a_jk，于是该式等于自己的相反数，只能是零。这条「对称张量与置换符号收缩必为零」的结论值得单独记下来，它是 2.9 题和后面所有反对称运算的地基。

2.10 题给矩阵 T = [[1,2,3],[4,5,6],[7,8,9]]，先拆成对称部分和反对称部分：

• 对称部分 [T]+[T]ᵀ 的一半，等于 [[1,3,5],[3,5,7],[5,7,9]]；
• 反对称部分 [T]−[T]ᵀ 的一半，等于 [[0,−1,−2],[1,0,−1],[2,1,0]]。

反对称部分的轴向矢量需要先定约定。教材里对应的是 A_ij = e_ijk ω_k 这一套，反解就是 ω_k = ½ e_kij A_ij，代入上面三个分量算下来是 ω = −e1 + 2e2 − e3。换一套约定（比如把 Ω 写成 ½(u_j,i − u_i,j)）就会整体差一个负号，2.14 题给出的 ω = ½∇×a 与这里必须用同一套定义，核对答案前先看书里的定义式，再决定要不要补那个负号。

2.3 numpy 复算：三条式子一次跑完

手工推完之后用数值兜底，成本极低。下面这段把 2.1 题的三条式子和 e-δ 母公式全部验一遍：

import numpy as np from itertools import permutations # 构造三阶置换张量 e_ijk：偶排列 +1，奇排列 -1，其余为 0 e = np.zeros((3, 3, 3)) for perm in permutations(range(3)): inv = sum(perm[i] > perm[j] for i in range(3) for j in range(i + 1, 3)) e[perm] = -1.0 if inv % 2 else 1.0 delta = np.eye(3) # 母公式：e_pqi e_ijk = delta_pj delta_qk - delta_pk delta_qj lhs = np.einsum('pqi,ijk->pqjk', e, e) rhs = (np.einsum('pj,qk->pqjk', delta, delta) - np.einsum('pk,qj->pqjk', delta, delta)) print('e-e 恒等式残差:', np.abs(lhs - rhs).max()) rng = np.random.default_rng(0) # 2.1(2)：e_pqi e_ijk A_jk = A_qp - A_pq A = rng.normal(size=(3, 3)) r2 = np.einsum('pqi,ijk,jk->pq', e, e, A) print('2.1(2) 残差:', np.abs(r2 - (A.T - A)).max()) # 2.1(3)：e_ijp e_klp B_ki B_lj = (tr B)^2 - tr(B^2) B = rng.normal(size=(3, 3)) r3 = np.einsum('ijp,klp,ki,lj->', e, e, B, B) print('2.1(3) 残差:', abs(r3 - (np.trace(B) ** 2 - np.trace(B @ B))))

参数上要留意两处。einsum的下标串里，出现在输入但不在输出中的指标会被自动求和，'pqi,ijk,jk->pq'里的 i、j、k 都会被消掉，只留 p、q；输出串的顺序决定结果的轴排布，'pqjk'写成'pjqk'形状一样但对不上。构造 e 时用的是逆序数奇偶，这里用的是朴素双循环，3 阶只有 6 个非零元，够用。

2.4 三个高频写法错误

第一，哑标在同一项里出现三次。δ_ii δ_jj 是合法缩并，δ_ii δ_ij 里 i 就成了非法重名，必须换一个字母。

第二，等式两边自由标不一致。左边是 δ_pk，右边写成 δ_pi，即使数值上碰巧相等也是错的。

第三，忽略上下标的位置差异。直角坐标下的笛卡尔张量上下标可以统一写下标，一旦切到曲线坐标、或者要处理与 Christoffel 符号有关的量，位置就不再是装饰，这一步在教材第三章之后会越来越明显。

3. 绕 z 轴旋转的坐标变换：2.5、2.6 与 3.7 共用同一套方向余弦

2.5、2.6、3.7 三题看起来分属矢量、张量、应变三个话题，实际用的是同一个旋转矩阵，只是被作用的对象阶数不同。把方向余弦矩阵写死一次，后面三题都能直接调。

3.1 方向余弦矩阵的定义与符号

教材 2.5 题把新坐标系的基矢量对老坐标系的投影记成 β_ij，绕 z 轴转 θ 时取值是：

这张表是后面所有计算的唯一入口。只要它抄错一个符号，2.5 的矢量分量、2.6 的张量分量、3.7 的应变分量会一起错，而且错得很有规律、不容易当场发现。

3.2 矢量分量：2.5 题的二倍角捷径

一阶张量的变换是 u'_i = β_ij u_j，展开得到

• u'_1 = u_1 cosθ + u_2 sinθ
• u'_2 = −u_1 sinθ + u_2 cosθ
• u'_3 = u_3

写到这里先别急着往下算，检查一遍正交性：β 的任意两行点积为零、每行模长为 1、行列式为 1。三条都满足，说明这是一次正当的刚体转动，而不是带反射的伪旋转。很多同学在 2.6 题算出「T'_11 不守恒」之类的怪结果，回头查就是这里多写了一个负号。

3.3 二阶张量变换：2.6 题的分量展开

二阶张量按 T'_ij = β_ik β_jl T_kl 变换，把上表代进去、用二倍角化简：

• T'_11 = (T_11 + T_22)/2 + (T_11 − T_22)/2 · cos2θ + T_12 · sin2θ
• T'_12 = (T_22 − T_11)/2 · sin2θ + T_12 · cos2θ
• T'_13 = T_13 cosθ + T_23 sinθ
• T'_33 = T_33

T'_33 不变是意料之中的：绕 z 轴转动不会把 z 方向的法向应力搅进来，而 T'_11 里出现的 2θ 说明二阶张量的分量以二倍频率变化，这一点在 3.7、3.8 两题里还会再用一次。

3.4 应变张量的转动：3.7 与 3.8 题

应变是二阶对称张量，变换规律与 2.6 完全一样，只是把 T 换成 ε。把 3.7 的六个分量写出来：

ε'_x = (εx + εy)/2 + (εx − εy)/2 · cos2θ + ε_xy · sin2θ

ε'_y = (εx + εy)/2 − (εx − εy)/2 · cos2θ − ε_xy · sin2θ

ε'_xy = −(εx − εy)/2 · sin2θ + ε_xy · cos2θ

ε'_xz = ε_xz cosθ + ε_yz sinθ，ε'_yz = −ε_xz sinθ + ε_yz cosθ，ε'_z = ε_z。

这里最容易踩的坑是工程剪应变与张量剪应变差一个 2：题面里的 γ_xy 是工程量，代进公式前要除以 2 变成 ε_xy，否则结果整体偏大一倍。

3.8 题是这个变换的实际用法：在 Oxy 平面上贴三片应变片，方向分别是 0°、60°、120°，测到 εa、εb、εc，要反推任意方向的正应变。设任意方向的正应变为 ε_n = A + B cos2θ + C sin2θ，把三个已知角度代进去解三元一次方程组，得到

代回即可还原任意角度。验证一下：θ = 0 时 A + B = εa，θ = 60° 时 A − B/2 + C√3/2 = εb，θ = 120° 时 A − B/2 − C√3/2 = εc，三个点都对上，说明系数没解错。

import numpy as np def rosette_coeffs(ea, eb, ec): """由 0/60/120 度三片应变片反解 A、B、C 三个系数""" A = (ea + eb + ec) / 3.0 B = (2.0 * ea - eb - ec) / 3.0 C = (eb - ec) / np.sqrt(3.0) return A, B, C def eps_n(A, B, C, theta_deg): """任意方向的正应变，theta 为与 x 轴夹角（度）""" t = np.deg2rad(theta_deg) return A + B * np.cos(2 * t) + C * np.sin(2 * t) ea, eb, ec = 100e-6, 40e-6, -20e-6 A, B, C = rosette_coeffs(ea, eb, ec) print(eps_n(A, B, C, 0), eps_n(A, B, C, 60), eps_n(A, B, C, 120))

函数的输入是三个实测应变，单位保持一致即可（这里用微应变）；输出是三个系数，量纲与输入相同。回代 0°、60°、120° 应当原样复现输入值，这是检查公式有没有写错的最快办法。实际布片时如果三片不是严格的 0/60/120，需要按真实角度重新列方程，别硬套上面的系数。

3.5 用不变量做交叉验证

对任意 θ 跑一遍变换，再检查两个量：tr(T') 与 tr(T) 相等，det(T') 与 det(T) 相等。转动是正交变换，这两条必须成立。三行代码就能挂进流程里：

def rot_z(theta): """绕 z 轴转 theta 弧度，返回 x'_i = beta_ij x_j 形式的方向余弦矩阵""" c, s = np.cos(theta), np.sin(theta) return np.array([[c, s, 0.0], [-s, c, 0.0], [0.0, 0.0, 1.0]]) T = np.array([[12.0, 3.0, -2.0], [3.0, 5.0, 1.0], [-2.0, 1.0, -3.0]]) for deg in (0, 15, 30, 60, 90): b = rot_z(np.deg2rad(deg)) Tp = b @ T @ b.T # T'_ij = beta_ik beta_jl T_kl assert np.allclose(np.trace(Tp), np.trace(T)) assert np.allclose(np.linalg.det(Tp), np.linalg.det(T))

b @ T @ b.T就是 T'_ij = β_ik β_jl T_kl 的矩阵写法，两次乘法分别吃掉两个变换系数。这类断言一旦挂上，任何一处符号写反都会被立刻拦住，比事后盯着一长串三角函数找错快得多。

4. 应变张量、协调方程与位移场重建

第三章的题分两类：一类是从位移场出发求应变，一类是从应变反推位移是否存在。前者是求导，后者是可积性问题，考的是协调方程。

4.1 位移梯度拆成对称与反对称两部分

3.2 题给的位移场是 u = A·r，A 是与 r 无关的二阶常张量。对 u 求梯度直接得到 ∇u = A，然后拆成对称与反对称两部分：

def strain_and_spin(A): """输入位移梯度 A_ij = u_i,j，返回应变张量和对数转动张量""" eps = 0.5 * (A + A.T) # 对称部分，应变 om = 0.5 * (A - A.T) # 反对称部分，刚体转动 return eps, om def axial_vector(om, e): """由反对称张量取轴向矢量：omega_k = 0.5 * e_kij * om_ij""" return 0.5 * np.einsum('kij,ij->k', e, om)

输入 A 是位移梯度而不是位移本身，这一步很容易搞混：位移场 u = A·r 里的 A 求完导就是它自己，所以 ∇u 直接等于 A；如果位移场写成别的形式，得先老老实实求偏导。

输出里的 ε 是对称张量，六个独立分量；Ω 是反对称张量，三个独立分量，正好对应刚体转动的三个自由度。3.6 题讨论面积变化率、3.5 题讨论两条微线段夹角的变化，用到的都是 ε；而刚体转动部分不影响任何长度和角度，只能通过 Ω 表现出来。

4.2 协调方程：3.9 题的常数约束

应变是由位移求导得来的，六个分量之间必须满足可积条件，也就是圣维南协调方程。直角坐标下标量形式可以统一写成

ε_ij,kl + ε_kl,ij − ε_ik,jl − ε_jl,ik = 0

它对 i、j、k、l 共 81 个组合成立，但真正独立的只有 6 个。3.9 题给了一组含参数 a、b 的应变分量，要求判断是否可能发生，做法就是把这些分量代进去，看能不能让残差恒为零。手工展开容易漏项，交给符号计算更稳：

import sympy as sp x, y, z, a, b = sp.symbols('x y z a b', real=True) coords = (x, y, z) # 应变张量（注意：剪应变要用张量分量，即工程剪应变的一半） eps = sp.Matrix([ [a * y**2, 0, (a * x**2 + b * y**2) / 2], [0, a * x**2 * y, (a * y**2 + b * z**2) / 2], [(a * x**2 + b * y**2) / 2, (a * y**2 + b * z**2) / 2, a * x * y], ]) def compat_residual(eps, coords): """返回所有非零的协调方程残差""" out = {} for i in range(3): for j in range(3): for k in range(3): for l in range(3): r = (sp.diff(eps[i, j], coords[k], coords[l]) + sp.diff(eps[k, l], coords[i], coords[j]) - sp.diff(eps[i, k], coords[j], coords[l]) - sp.diff(eps[j, l], coords[i], coords[k])) r = sp.simplify(r) if r != 0: out[(i, j, k, l)] = r return out res = compat_residual(eps, coords) sols = set() for expr in res.values(): sols |= set(sp.solve(expr, (a, b))) print(sols)

这段代码有两个容易忽略的点。第一，剪应变输入时必须除以 2，eps矩阵里凡是(a*x**2 + b*y**2)/2这类写法的都是工程剪应变转换过来的；输入的 γ_yz = a y² + b z²、γ_xz = a x² + b y²、γ_xy = 0。第二，sp.solve返回的是嵌套解集，用集合展开去重后才好读。跑完能看到残差只在 a、b 不同时为零时才被消掉，所以 a = b = 0 是唯一可能的情形。

4.3 3.10 题：让符号计算给出常数关系

3.10 题是反过来的问法：给一组含 A0、A1、B0、B1、C0、C1、C2 的应变分量，要求确定各常数之间的关系，使协调方程成立。方法完全一样，只是把上面代码里的 eps 换成题给的表达式，最后对七个常数求解。得到的约束只会落在 A1、B1、C1、C2 上，A0、B0、C0 保持自由，这与解答里「其余三个常数可以是任意的」的说法一致。用符号求解的好处是它不会漏项——81 个组合里只要有一个残差没被消掉，solve就会把它带出来。

4.4 从应变反推位移：均匀应变下的通式

3.11、3.12 两题讨论的是均匀应变（应变张量与坐标无关）时的位移一般表达式。结论是位移可以拆成三块：任意的刚体平移 u0、任意的刚体转动 ω0 × (r − r0)、以及由应变积分出来的变形部分：

u = u0 + ω0 × (r − r0) + ε · (r − r0)

3.12 题给的是 ε = a e1e1 + b e2e2 + c e3e3 这类只有正应变、且各自是坐标函数的情形，变形部分积分出来是 ∇[½(a x² + b y² + c z²)] 的形式。数值复核时可以验证一件事：把求出的位移场代回几何方程，应当原样还原给定的应变分量，同时 Ω 只贡献刚体转动、不影响任何长度。

5. 应力张量：斜截面、主应力与面力边界条件

第四章的题从「一点的应力状态」出发，先算斜截面上的三个量，再求主应力和不变量，最后落到面力边界条件。这三步在有限元前后处理里都能直接对应到代码。

5.1 斜截面上的总应力、正应力与剪应力

4.1 题给的是 σx = 50a、σy = 0、σz = −30a、τyz = −75a、τzx = 80a、τxy = 50a，法线方向余弦是 (1/2, 1/2, √2/2)。计算分三步：先求应力矢量 T_i = σ_ij n_j，再求正应力 σ_n = T_i n_i，最后用 |T|² − σ_n² 开方得到剪应力。

import numpy as np def traction(S, n): """给定应力矩阵与单位法向，返回总应力、正应力、剪应力""" T = S @ n # T_i = sigma_ij n_j Tn = float(T @ n) # 正应力 tau = np.sqrt(max(float(T @ T) - Tn ** 2, 0.0)) return T, Tn, tau a = 1.0 S = np.array([[50.0, 50.0, 80.0], [50.0, 0.0, -75.0], [80.0, -75.0, -30.0]]) * a n = np.array([0.5, 0.5, np.sqrt(2) / 2]) T, Tn, tau = traction(S, n) print(T, np.linalg.norm(T), Tn, tau)

S必须是完整的三阶对称矩阵，题面给的六个应力分量按 σx、σy、σz 放对角线，剪应力填到对应位置，注意 τyz 与 τzy 数值相同。n要归一化，否则 σ_n 会被法向长度放大。最后那个max(..., 0.0)是防浮点误差把开方项弄成微小负数的，工程计算里很常见。

5.2 主应力、不变量与八面体应力

4.9、4.10、4.11 三题都围绕特征值展开。主应力是应力张量的特征值，特征矢量是主方向，三个不变量分别是 I1 = σx + σy + σz、I2 = σxσy + σyσz + σzσx − τxy² − τyz² − τzx²、I3 = det(σ)。数值上用对称矩阵的特征值求解器最省事：

def stress_invariants(S): """返回三个主应力和三个不变量""" principal = np.linalg.eigvalsh(S) # 对称矩阵，实特征值升序 I1 = np.trace(S) I2 = 0.5 * (np.trace(S) ** 2 - np.trace(S @ S)) I3 = np.linalg.det(S) return principal, (I1, I2, I3) def octahedral(principal): """八面体正应力与剪应力""" s1, s2, s3 = principal s0 = (s1 + s2 + s3) / 3.0 t0 = np.sqrt((s1 - s2) ** 2 + (s2 - s3) ** 2 + (s3 - s1) ** 2) / 3.0 return s0, t0

eigvalsh只接受对称矩阵，正好匹配应力张量，不要用通用的eig，否则复数特征值会让结果没法读。八面体应力公式里三个主应力之差的平方和除以 3，等价于 √(2/3) 乘上偏应力第二不变量，两个写法都可以，核对答案时注意系数。

4.11 题是个值得动手的例子：σ11 = σ22 = σ33 = 0，三个剪应力都等于 σ。手算特征多项式得到特征值 2σ 和 −σ（二重根），对应主方向之一是 n = (e1 + e2 + e3)/√3，另外两个主方向在与 n 垂直的平面内任取。用上面的函数跑一遍，eigvalsh会给出 [−σ, −σ, 2σ]，顺序是从小到大，别被顺序误导。

5.3 面力边界条件：从 4.4 到 4.8

面力边界条件的统一形式是 σ_ij n_j = t_i，自由面上 t_i = 0，受法向压力 p 的面上 t_i = −p n_i。4.6 题讨论的是曲面 f(x, y, z) = 0，外法向由梯度给出：

n = ∇f / |∇f|，代入边界条件后写成 σ_ij f_j = p f_i 的指标形式。

4.4 题的三角柱体是个典型例子。它有两段边界：底边 y = 0 上受均匀压力 q，斜面上自由。底边的条件是 σ_y = −q、τ_xy = 0；斜面上两个分量都得为零，把应力表达式代进去就得到一个关于 A、B、C 的线性方程组。手算容易在斜面法向的方向余弦上翻车，用符号求解稳一点：

import sympy as sp q, beta, A, B, C = sp.symbols('q beta A B C', real=True) # 底边 y = 0：sigma_y = -q，tau_xy = 0（后者已自动满足） eq1 = sp.Eq(-(A + B), -q) # 具体表达式按题面代入 # 斜面 y = x*tan(beta)：外法向 (sin(beta), cos(beta), 0)，两个分量分别为零 n1, n2 = sp.sin(beta), sp.cos(beta) eq2 = sp.Eq((A * sp.sin(beta) + A * sp.cos(beta) + C) * n1, 0) eq3 = sp.Eq((A * sp.sin(beta) - B * sp.cos(beta)) * n2, 0) sol = sp.solve([eq1, eq2, eq3], (A, B, C), dict=True) print(sol)

solve的目标变量顺序决定返回结构，dict=True会给出键值对，读起来更直观。这里两个方程式的具体系数要按题面给的 σx、σy、τxy 表达式替换，代码框架不用改。解出来 A、B、C 都是 q 与 β 的比值形式，与教材答案核对时要确认三件事：法向取的是外法向还是内法向、x 轴方向怎么定、β 是从哪个轴量起的。这三条中任意一条反了，结果都会差一个符号。

4.7 题的球体一半浸在液体里，边界条件按 z 的符号分两段：z ≤ 0 时球面上自由，σ_ij x_j = 0；z > 0 时受液体压力 ρgz，σ_ij x_j = −ρgz x_i。写法上直接套 4.6 的梯度形式即可，因为球的方程 x² + y² + z² = a² 的梯度正比于矢径。

4.8 题讨论的是静水应力状态 σ_ij = σ δ_ij。代进平衡方程后可以发现体力是有势的，且势函数就是 σ 本身；表面上的面力则退化成 t_i = σ n_i，方向永远与法向一致。这类题的价值在于让你意识到：体积力与面力在静水应力下没有区别，都是同一个球张量的不同表现。

6. 用残差断言把三十多道题压成五个检查

习题做多了会发现，真正需要反复核对的只有五类量：指标缩并的结果、旋转后的不变量、协调方程的残差、特征值、以及边界残差。把这五类各写一个断言函数，后面无论改公式还是换参数，跑一遍就知道有没有写错。

import numpy as np def check_identity(lhs, rhs, tol=1e-10): """张量恒等式：左右两边作差的无穷范数""" return np.abs(np.asarray(lhs) - np.asarray(rhs)).max() < tol def check_invariants(T, beta, tol=1e-10): """坐标变换后迹与行列式不变""" Tp = beta @ T @ beta.T return (abs(np.trace(Tp) - np.trace(T)) < tol and abs(np.linalg.det(Tp) - np.linalg.det(T)) < tol) def check_compat(eps_fn, coords, tol=1e-9): """协调方程数值残差：eps_fn(x) 返回该点的应变矩阵""" import itertools worst = 0.0 for i, j, k, l in itertools.product(range(3), repeat=4): h = 1e-4 # 用二阶中心差分近似 eps_ij,kl val = (eps_fn(np.array(coords) + h * np.eye(3)[k] + h * np.eye(3)[l])[i, j] - eps_fn(np.array(coords) + h * np.eye(3)[k] - h * np.eye(3)[l])[i, j] - eps_fn(np.array(coords) - h * np.eye(3)[k] + h * np.eye(3)[l])[i, j] + eps_fn(np.array(coords) - h * np.eye(3)[k] - h * np.eye(3)[l])[i, j]) / (4 * h * h) worst = max(worst, abs(val)) return worst < tol def check_boundary(S, n, t, tol=1e-8): """面力边界条件 sigma_ij n_j = t_i 的残差""" return np.abs(S @ np.asarray(n) - np.asarray(t)).max() < tol def check_eigen(S, tol=1e-8): """特征值分解自洽：S v = lambda v""" w, v = np.linalg.eigh(S) return np.abs(S @ v - v * w).max() < tol

这几个函数都不长，但覆盖了教材第二到第六章的主要验算点。check_compat用的是二阶中心差分，步长取 1e-4 是精度与舍入误差之间的折中，步长再小会被浮点相减吃掉有效位，再大截断误差就上来了；如果是符号表达式，还是用第 4 章的 sympy 版本更靠谱，数值差分只适合快速排查。

使用时的顺序也有讲究。先跑check_identity确认指标和符号都没错，再跑check_invariants确认变换矩阵是正交的、行列式为 1，然后才去看特征值和边界残差——前两步没过，后面的结果没有意义。

有一个坑必须单独提：这份解答是扫描件转出来的，上下标、希腊字母和运算符在 OCR 阶段丢得很厉害，2.6 题里那几个分数和 2.10 题的矩阵符号都需要对照教材原文重新辨认。稳妥的读法是把解答当成「最终答案的提示」，推导过程自己补一遍，最后用上面这套残差断言做裁判。凡是断言过不去的，先怀疑自己抄错了指标顺序或者漏了一个负号，而不是先怀疑答案。把assert挂到日常脚本里，每次调整公式都会被立刻拦下来。

本文还有配套的精品资源，点击获取