实现Burgers方程时出现Singular Matrix错误求助

问题描述 #

我在实现Burgers方程以模拟压力波传播时，遭遇了Singular matrix错误，无法定位问题成因。作为Python和Jupyter Notebook新手，希望得到针对性的分析与排查建议。

报错信息 #

23 J = eval(jac)(x,n,fnon,F,*fnonargs)
     25 
---> 26 delta = np.linalg.solve(J,-F)
     27 x = x + delta
     29 F = eval(fnon)(x,*fnonargs)
    111 def _raise_linalgerror_singular(err, flag):
--> 112     raise LinAlgError("Singular matrix")

LinAlgError: Singular matrix

相关代码 #

我曾怀疑Jacobian函数存在问题，但它在之前的模型中运行正常，函数代码如下：

import numpy as np
import copy  # 原代码未导入copy，补充后可避免运行报错

def Jacobian(x,n,fnon,F0,*fnonargs):
  J = np.zeros((n,n), dtype=np.float64)
  h = 10e-8

  for k in range(0,n):
    xb = copy.copy(x)
    xb[k] = xb[k] + h

    F = eval(fnon)(xb,*fnonargs)

    for i in range(0,n):
        J[i,k] = (F[i] - F0[i]) / h

  return J

排查建议 #

• 调整数值微分步长：当前h=10e-8过小，易引发浮点精度损失，尝试改为1e-6或1e-5，观察Jacobian矩阵是否仍奇异。
• 校验残差函数实现：奇异矩阵通常源于残差函数梯度线性相关，检查Burgers方程的离散逻辑（如有限差分格式、边界条件）是否正确，确保残差函数能准确反映方程的物理意义。
• 更换迭代初始值：初始猜测值可能处于残差函数的奇异区域，尝试用均匀场或简单解析解作为初始条件，看是否能避免矩阵奇异。
• 替换eval()调用：eval(fnon)既不安全又易引发隐式错误，直接传入函数对象替代函数名字符串，比如将fnon参数改为函数本身。
• 检查矩阵秩：在np.linalg.solve前添加print(np.linalg.matrix_rank(J))，若秩小于n，说明存在线性相关的行/列，可进一步定位问题残差项。
• 验证边界条件：过度约束或错误的边界条件会导致残差系统出现冗余方程，进而让Jacobian矩阵奇异，确认边界条件的离散是否符合物理场景。

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