简介面向力学、土木、机械等专业本科生及考研复习者《弹性力学基础》课后习题解答以同济大学程尧舜版教材第二章为范围逐题给出推导过程适合课堂同步练习与考前查漏补缺。压缩包内仅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 与 δ 的乘积恒等式是所有化简的来源先把它们列清楚母公式化简结果出现位置δ_ij δ_jkδ_ik2.1(1)e_pqi e_ijkδ_pj δ_qk − δ_pk δ_qj2.1(2)e_ijp e_klpδ_ik δ_jl − δ_il δ_jk2.1(3)e_ijk e_ijl2δ_kl涡量与旋度互推e_ijk e_ijk6三维指标组合计数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 轴转 θ 时取值是i\j1231cosθsinθ02−sinθcosθ03001这张表是后面所有计算的唯一入口。只要它抄错一个符号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_33T_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θ把三个已知角度代进去解三元一次方程组得到系数表达式物理含义A(εa εb εc)/3面内平均正应变B(2εa − εb − εc)/3偏应变分量之一C(εb − εc)/√3偏应变分量之二代回即可还原任意角度。验证一下θ 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·rA 是与 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, realTrue) 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)这段代码有两个容易忽略的点。第一剪应变输入时必须除以 2eps矩阵里凡是(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² 开方得到剪应力。量计算式数值结果T1σ11 n1 σ12 n2 σ13 n3106.57aT2σ21 n1 σ22 n2 σ23 n3−28.03aT3σ31 n1 σ32 n2 σ33 n3−18.71a总应力T正应力 σ_nT_i n_i26.04a剪应力 τ_n√(|T|² − σ_n²)108.7aimport 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, t0eigvalsh只接受对称矩阵正好匹配应力张量不要用通用的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, realTrue) # 底边 y 0sigma_y -qtau_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), dictTrue) print(sol)solve的目标变量顺序决定返回结构dictTrue会给出键值对读起来更直观。这里两个方程式的具体系数要按题面给的 σx、σy、τxy 表达式替换代码框架不用改。解出来 A、B、C 都是 q 与 β 的比值形式与教材答案核对时要确认三件事法向取的是外法向还是内法向、x 轴方向怎么定、β 是从哪个轴量起的。这三条中任意一条反了结果都会差一个符号。4.7 题的球体一半浸在液体里边界条件按 z 的符号分两段z ≤ 0 时球面上自由σ_ij x_j 0z 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, tol1e-10): 张量恒等式左右两边作差的无穷范数 return np.abs(np.asarray(lhs) - np.asarray(rhs)).max() tol def check_invariants(T, beta, tol1e-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, tol1e-9): 协调方程数值残差eps_fn(x) 返回该点的应变矩阵 import itertools worst 0.0 for i, j, k, l in itertools.product(range(3), repeat4): 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, tol1e-8): 面力边界条件 sigma_ij n_j t_i 的残差 return np.abs(S np.asarray(n) - np.asarray(t)).max() tol def check_eigen(S, tol1e-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挂到日常脚本里每次调整公式都会被立刻拦下来。本文还有配套的精品资源点击获取