13
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ "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." Please clarify. Are you saying that you run the same problem with an exact Jacobian and the solution converges to the exact solution of the PDE? That is important information. $\endgroup$ – David Ketcheson May 16 '15 at 6:15
  • $\begingroup$ @DavidKetcheson Yes, that is what I'm saying. Apologies if my terminology is incorrect or incomplete. (I suppose I should also have said "I get stable and expected results.") $\endgroup$ – Stephen Bosch May 16 '15 at 8:03
3
$\begingroup$

More of a long comment than anything else:

Using the descriptions provided in Ascher and Petzold 1998, I wrote this function which determines the gradient at a given point x:

Look at the code for SUNDIALS' difference quotient approximation to get a better idea of what you should do in an implementation. Ascher and Petzold is a good book to get started, but SUNDIALS is actually used in production work, and thus has been better tested. (Also, SUNDIALS is related to DASPK, which Petzold worked on.)

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?

Empirically, approximate Jacobians can cause convergence failures in Newton's method. I don't know that I'd characterize them as "instability"; in some cases, it's just not possible to achieve the desired error tolerances in the termination criteria. In other cases, it could manifest as an instability. I'm almost certain there's a more quantitative result on this phenomenon in Higham's numerical methods book, or Hairer and Wanner's discussion of W-methods.

Or is the cause somewhere else? How else might I approach this problem?

It depends where you think the error might be. If you're extremely confident in your implementation of backward Euler, I wouldn't start there. Experience has made me paranoid in my implementations of numerical methods, so if it were me, I'd start by coding up a few really basic test problems (a couple nonstiff and stiff linear problems, the heat equation with a centered finite difference approximation, things like that) and I'd use the method of manufactured solutions to assure myself that I know what the solution will be, and what I should be comparing against.

However, you've already done some of that:

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.

That would be the next thing I would test: use an analytical Jacobian. After that, you might also look at the extremal eigenvalues of your finite difference Jacobian on the off chance you're in the unstable region of backward Euler. Looking at extremal eigenvalues of your analytical Jacobian as a basis for comparison might give you some insight. Assuming those all check out, the problem is probably in the Newton solve.

$\endgroup$
  • $\begingroup$ Thanks for the thoughtful analysis (plus the SUNDIALS hint and the alternate sources). My professor suggested setting lambda = 0, arguing that the FDA of the PDE then becomes linear, so we would expect the FDA Jacobian to equal the analytical Jacobian. When I do this, it manages three timesteps, newton() hitting maxiter each time, before finally blowing up as before. $\endgroup$ – Stephen Bosch May 23 '15 at 15:41
  • $\begingroup$ He also said it was not common practice to use approximated Jacobians for solving PDEs and suggested that it might be having trouble due to the many degrees of freedom (without providing an explanation, although having just looked at Hairer and Wanner's discussion of W-methods, I can see that it is probably not trivial). $\endgroup$ – Stephen Bosch May 23 '15 at 15:45
  • 1
    $\begingroup$ Your professor's statement is somewhat surprising, given the amount of literature on the subject, for instance this famous review by Knoll and Keyes. I probably should have cited this paper in my answer, as the sources in the bibliography may also be of some help in diagnosing your problems. $\endgroup$ – Geoff Oxberry May 24 '15 at 1:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.