I have implemented a backward-Euler solver in python 3 (using numpy). For my own convenience and as an exercise, I also wrote a small function that computes a finite difference approximation of the gradient so that I don't always have to determine the Jacobian analytically (if it is even possible!).
Using the descriptions provided in Ascher and Petzold 1998, I wrote this function which determines the gradient at a given point 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
I tested this function by taking a multivariate function for the pendulum and comparing the symbolic Jacobian to the numerically determined gradient for a range of points. I was pleased with the results of the test, the error was around 1e-10. When I solved the ODE for the pendulum using the approximated Jacobian, it also worked very well; I couldn't detect any difference between the two.
Then I tried testing it with the following PDE (Fisher's equation in 1D):
$$ \partial_t \mathbf{u} = \partial_x(k\partial_x\mathbf{u}) + \lambda(\mathbf{u}(C-\mathbf{u}))$$
using a finite difference discretization.
Now Newton's method blows up in the first timestep:
/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
This happens for a variety of eps values, but strangely, only when the PDE spatial step size and time step size are set so that the Courant–Friedrichs–Lewy condition is not met. Otherwise it works. (This is the behaviour you'd expect if solving with forward Euler!)
For completeness, here is the function for the Newton Method:
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
(The function la.solve is scipy.linalg.solve.)
I am confident that my backward Euler implementation is in order, because I have tested it using a function for the Jacobian and get stable results.
I can see in the debugger that newton() manages 35 iterations before the error occurs. This number remains the same for every eps I have tried.
An additional observation: when I compute the gradient with FDA and a function using the initial condition as an input and compare the two while varying the size of epsilon, the error grows as epsilon shrinks. I would expect it to be large at first, then get smaller, then larger again as epsilon shrinks. So an error in my implementation of the Jacobian is a reasonable assumption, but if so, it's so subtle that I am unable to see it. EDIT: I modified jacobian() to use forward instead of central differences, and now I observe the expected development of the error. However, newton() still fails to converge. Observing dx in the Newton iteration, I see that it only grows, there is not even a fluctuation: it nearly doubles (factor 1.9) with each step, with the factor getting progressively larger.
Ascher and Petzold do mention that difference approximations for the Jacobian do not always work well. Can an approximated Jacobian with finite differences cause instability in Newton's method? Or is the cause somewhere else? How else might I approach this problem?