本文揭示了在实现多变量牛顿-拉夫逊法时一个典型但隐蔽的矩阵运算错误——将残差向量与雅可比逆矩阵的乘法顺序颠倒，导致迭代公式数学失真，进而引发振荡不收敛；同时提供修正代码、数值稳定性建议与工程实践优化方案。
本文揭示了在实现多变量牛顿-拉夫逊法时一个典型但隐蔽的矩阵运算错误——将残差向量与雅可比逆矩阵的乘法顺序颠倒，导致迭代公式数学失真，进而引发振荡不收敛；同时提供修正代码、数值稳定性建议与工程实践优化方案。
牛顿-拉夫逊法求解非线性方程组的核心迭代公式为：
$$ \mathbf{x}_{n+1} = \mathbf{x}_n - \mathbf{J}^{-1}(\mathbf{x}_n)\, \mathbf{F}(\mathbf{x}_n) $$
其中 $\mathbf{F} = [f_1, f_2]^\top$ 是列向量形式的残差，$\mathbf{J}$ 是 $2\times2$ 雅可比矩阵，$\mathbf{J}^{-1}\mathbf{F}$ 表示对线性系统 $\mathbf{J}\,\Delta\mathbf{x} = -\mathbf{F}$ 求解校正量 $\Delta\mathbf{x}$。
而原代码中关键一行存在根本性数学错误：
result = -np.matmul(np.array([f1(x0, y0), f2(x0, y0)]), np.linalg.inv(jacobian_mat(x0, y0))) + np.array([x0, y0])
此处将 1×2 行向量（[f1, f2]）左乘 2×2 矩阵逆，等价于计算：
$$
\mathbf{F}^\top \mathbf{J}^{-1}
$$
而非正确的：
$$
\mathbf{J}^{-1} \mathbf{F}
$$
由于矩阵乘法不可交换，且 $\mathbf{F}^\top \mathbf{J}^{-1} \neq (\mathbf{J}^{-1} \mathbf{F})^\top$（除非 $\mathbf{J}$ 对称），该操作在数学上无明确意义，实际计算出的是错误的更新方向——这正是高精度（如 TOLERANCE=0.001）下出现振荡、无法收敛的根本原因。低容差（如 0.01）下偶然“成功”，实为误差掩盖了算法缺陷，并非真正进入收敛域。
✅ 正确做法是：将残差构造为列向量，用 np.linalg.solve 直接求解线性系统，避免显式求逆（更稳定、更高效）：
def compute_newton_raphson(x0, y0):
x, y = x0, y0
for i in range(ITER_LIMIT):
F = np.array([[f1(x, y)], [f2(x, y)]]) # 2×1 列向量
J = jacobian_mat(x, y) # 2×2 矩阵
try:
delta = np.linalg.solve(J, -F).flatten() # 解 J·δ = -F → δ 为 2×1 → flatten 为 (2,)
except np.linalg.LinAlgError:
print(f"Iteration {i}: Singular Jacobian at ({x:.4f}, {y:.4f})")
return None
x_new, y_new = x + delta[0], y + delta[1]
print(f"Iter {i}: x={x_new:.6f}, y={np.rad2deg(y_new):.6f}°")
if max(abs(x_new - x), abs(y_new - y)) <= TOLERANCE:
return (x_new, y_new)
# 防发散：若步长过大或变量溢出，提前终止
if abs(x_new) > 1e4 or abs(y_new) > 100*np.pi:
print(f"Iteration {i}: Divergence detected.")
return None
x, y = x_new, y_new
print("Max iterations reached.")
return None? 关键改进说明：
- ✅ 使用
np.linalg.solve(J, -F)替代@ np.linalg.inv(J)：数值更稳定，避免病态矩阵求逆失败； - ✅ 显式构造
F为2×1列向量（[[f1],[f2]]），确保维度匹配； - ✅ 增加
LinAlgError异常捕获，防止雅可比奇异导致崩溃； - ✅ 加入发散检测（如
|x| > 1e4或|y| > 100π），避免无限循环； - ✅ 初始值策略优化：
y是角度（弧度制），不应在[-100, 100]内随机整数采样（跨度超30个周期）。推荐改用：y0 = np.random.uniform(-np.pi, np.pi) # 合理角度范围 x0 = np.random.uniform(-50, 50) # 根据物理量纲设定
? 其他工程建议：
-
避免手动实现：对于生产级应用，优先使用
scipy.optimize.root(method='hybr')（本质为改良的 Powell Hybrid 方法）或method='krylov'，它们内置了阻尼、回溯线搜索与自动缩放，鲁棒性远超手写 NR； -
雅可比验证：可用有限差分近似
scipy.optimize.approx_fprime对比解析雅可比，确保导数无误； -
收敛判据升级：除步长范数外，建议同步检查残差范数
np.linalg.norm(F) ，双保险； -
参数归一化：若
x（单位 mm）与y（单位 rad）量纲差异巨大，应对变量做预缩放（如x̃ = x/10,ỹ = y），提升雅可比条件数。
综上，牛顿法失效往往不源于“函数难解”，而始于一个微小但致命的线性代数误用。修正矩阵乘法语义、采用数值稳健的求解器、辅以合理的初值与监控机制，即可让算法在 1e-6 级精度下稳定收敛。
