我在 python 3 中实现了一个反向欧拉求解器（使用 numpy）。为了我自己的方便和作为练习，我还编写了一个小函数来计算梯度的有限差分逼近，这样我就不必总是分析地确定雅可比行列式（如果可能的话！）。

使用Ascher 和 Petzold 1998中提供的描述，我编写了这个函数来确定给定点 x 的梯度：

def jacobian(f,x,d=4):
    '''computes the gradient (Jacobian) at a point for a multivariate function.

    f: function for which the gradient is to be computed
    x: position vector of the point for which the gradient is to be computed
    d: parameter to determine perturbation value eps, where eps = 10^(-d).
        See Ascher und Petzold 1998 p.54'''

    x = x.astype(np.float64,copy=False)
    n = np.size(x)
    t = 1 # Placeholder for the time step
    jac = np.zeros([n,n])
    eps = 10**(-d)
    for j in np.arange(0,n):
        yhat = x.copy()
        ytilde = x.copy()
        yhat[j] = yhat[j]+eps
        ytilde[j] = ytilde[j]-eps
        jac[:,j] = 1/(2*eps)*(f(t,yhat)-f(t,ytilde))
    return jac

我通过对钟摆采用多元函数并将符号雅可比行列式与一系列点的数值确定梯度进行比较来测试此函数。我对测试结果很满意，误差在 1e-10 左右。当我使用近似雅可比求解摆的 ODE 时，它的效果也很好；我无法发现两者之间的任何区别。

然后我尝试使用以下 PDE（一维费舍尔方程）对其进行测试：

$$\partial_t \mathbf{u} = \partial_x(k\partial_x\mathbf{u}) + \lambda(\mathbf{u}(C-\mathbf{u}))$$

使用有限差分离散化。

现在牛顿的方法在第一个时间步中爆炸了：

/home/sfbosch/Fisher-Equation.py:40: RuntimeWarning: overflow encountered in multiply
  du = (k/(h**2))*np.dot(K,u) + lmbda*(u*(C-u))
./newton.py:31: RuntimeWarning: invalid value encountered in subtract
  jac[:,j] = 1/(2*eps)*(f(t,yhut)-f(t,yschlange))
Traceback (most recent call last):
  File "/home/sfbosch/Fisher-Equation.py", line 104, in <module>
    fisher1d(ts,dt,h,L,k,C,lmbda)
  File "/home/sfbosch/Fisher-Equation.py", line 64, in fisher1d
    t,xl = euler.implizit(fisherode,ts,u0,dt)
  File "./euler.py", line 47, in implizit
    yi = nt.newton(g,y,maxiter,tol,Jg)
  File "./newton.py", line 54, in newton
    dx = la.solve(A,b)
  File "/usr/lib64/python3.3/site-packages/scipy/linalg/basic.py", line 73, in solve
    a1, b1 = map(np.asarray_chkfinite,(a,b))
  File "/usr/lib64/python3.3/site-packages/numpy/lib/function_base.py", line 613, in asarray_chkfinite
    "array must not contain infs or NaNs")
ValueError: array must not contain infs or NaNs

对于各种 eps 值都会发生这种情况，但奇怪的是，只有当 PDE 空间步长和时间步长设置为不满足 Courant-Friedrichs-Lewy 条件时。否则它工作。（如果使用正向 Euler 求解，这是您所期望的行为！）

为了完整起见，这里是牛顿法的函数：

def newton(f,x0,maxiter=160,tol=1e-4,jac=jacobian):
    '''Newton's Method.

    f: function to be evaluated
    x0: initial value for the iteration
    maxiter: maximum number of iterations (default 160)
    tol: error tolerance (default 1e-4)
    jac: the gradient function (Jacobian) where jac(fun,x)'''

    x = x0
    err = tol + 1
    k = 0
    t = 1 # Placeholder for the time step
    while err > tol and k < maxiter:
        A = jac(f,x)
        b = -f(t,x)
        dx = la.solve(A,b)
        x = x + dx
        k = k + 1
        err = np.linalg.norm(dx)
    if k >= maxiter:
        print("Maxiter reached. Result may be inaccurate.")
        print("k = %d" % k)
    return x

（函数 la.solve 是 scipy.linalg.solve。）

我相信我的反向欧拉实现是有序的，因为我已经使用雅可比函数对其进行了测试并获得了稳定的结果。

我可以在调试器中看到 newton() 在错误发生之前管理 35 次迭代。对于我尝试过的每个 eps，这个数字都保持不变。

另一个观察结果：当我使用 FDA 和一个使用初始条件作为输入的函数计算梯度并在改变 epsilon 大小的同时比较两者时，误差会随着 epsilon 的缩小而增长。我希望它一开始会很大，然后变小，然后随着 epsilon 的缩小再次变大。因此，在我的雅可比实现中的错误是一个合理的假设，但如果是这样，它是如此微妙以至于我无法看到它。编辑：我修改了 jacobian() 以使用正向而不是中心差异，现在我观察到错误的预期发展。但是， newton() 仍然无法收敛。观察牛顿迭代中的 dx，我看到它只会增长，甚至没有波动：每一步它几乎翻倍（因子 1.9），并且因子逐渐变大。

Ascher 和 Petzold 确实提到雅可比行列式的差分近似并不总是很好。具有有限差分的近似雅可比行列式会导致牛顿方法的不稳定性吗？还是其他原因？我还能如何解决这个问题？