如何获取scipy.optimize.fsolve求解方程组的各解精度

scipy.optimize.fsolve本身不会直接返回每个解的误差或不确定度，但可以通过它返回的infodict结合雅可比矩阵的信息来估算，具体步骤如下：

关键思路 #

fsolve在full_output=1时返回的infodict包含了求解过程的核心数据，我们可以利用这些数据重构雅可比矩阵，再通过协方差矩阵计算每个解的标准误差（即不确定度）。协方差矩阵的近似公式为：
(J^T J)^{-1} * (fvec^T fvec)/(n - p)
其中：

• J是雅可比矩阵
• fvec是求解结束时的残差向量
• n是方程个数，p是未知数个数

结合示例代码实现 #

首先修正你示例里的拼写错误（full_outupt改为full_output），然后添加误差估算的代码：

import numpy as np
from scipy.optimize import fsolve
from scipy.linalg import inv

def func(x):
    return [x[0] * np.cos(x[1]) - 4,
            x[1] * x[0] - x[1] - 5]

# 开启完整输出，获取求解过程数据
root, infodict, ier, msg = fsolve(func, [1, 1], full_output=1)

# 从infodict提取所需数据
fjac = infodict['fjac']  # 雅可比矩阵QR分解的正交矩阵转置
r = infodict['r']        # QR分解的上三角矩阵
fvec = infodict['fvec']  # 最终残差向量
n_eq = len(fvec)         # 方程数量
n_vars = len(root)       # 未知数数量

# 重构雅可比矩阵J：J = fjac.T @ r（fsolve内部用QR分解存储J，J = Q@R，fjac是Q^T）
J = fjac.T @ r[:n_vars, :]

# 计算协方差矩阵
residual_sq_sum = np.sum(fvec ** 2)
cov_matrix = inv(J.T @ J) * (residual_sq_sum / (n_eq - n_vars))

# 每个解的标准误差是协方差矩阵对角线元素的平方根
std_errors = np.sqrt(np.diag(cov_matrix))

print("求解结果:", root)
print("每个解的标准误差:", std_errors)

补充说明 #

• 当方程组是超定（方程数>未知数）时，这个误差估算方法更准确；如果是正定（方程数=未知数），残差理论上应该很小，此时可以简化协方差矩阵为(J^T J)^{-1}乘以残差的平均平方。
• 你也可以通过调整fsolve的xtol（参数相对误差阈值）和ftol（残差相对误差阈值）参数，直接控制求解的精度，参数值越小，求解精度越高，但计算时间也会更长。

内容的提问来源于stack exchange，提问作者oscarcapote