我正在尝试使用残差的雅可比矩阵相对于守恒变量的向量来研究数值离散化方案的稳定性。

对于中心差分格式离散的简单扩散方程，RHS 可以写为：

$RHS=\nu\frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{\Delta x^2}$

这个方程很有趣，因为雅可比矩阵有一个解析形式，可以用带状矩阵表示法写成：

$J=\nu/dx^2 (+1,-2,+1)$

基本上，我要做的是将数值近似雅可比行列式的特征值与解析解的特征值进行比较。

但是，我遇到了一些麻烦，因为结果没有意义。与精确的雅可比相比，数值近似的雅可比呈现完全不同的特征值。我怀疑问题出在我试图计算 Frechet 导数以计算雅可比行列式的方式上。下面的代码片段显示了我是如何这样做的：

# Residue function of a simple centered three point scheme.

def frhs(um1,uo1,up1,dx,nu):
    return nu*( (up1 - 2.0 * uo1 + um1)/dx**2.0 )

# Computes the derivative of the residues with respect to the solution vector.

eps = 0.00001

drhs_du = np.zeros((nx-1,nx-1))  # In order to take the eigenvalues, this shall be a matrix.

e = np.zeros(nx-1)

# This loop computes the jacobian matrix according to http://www.netlib.org/math/docpdf/ch08-04.pdf

for i in range(1,nx-1):
    for j in range(1,nx-1):

        e[j] = 1.0

        drhs_du[i,j] = ( frhs(un[i-1]+eps*e[j], un[i]+eps*e[j], un[i+1]+eps*e[j],dx,nu) - frhs(un[i-1], un[i],un[i+1],dx,nu) )/eps

        e[j] = 0.0

完整代码见： https ://github.com/lmarmotta/disc_stability/blob/master/1dDiffusionEq.py

我做错什么了吗？

你们如何在数字上近似雅可比？对于这些类型的情况，还有其他算法可以近似雅可比矩阵吗？

谢谢 ！