Backward Euler + Quasi Newton(Broyden) method fails to solve Van der Pol's equation(Stiff ODE)

Question

The first guess is using the forward Euler approach. The first jacobian is using finite differences. Then NR method is used to solve for the next iteration and Broyden's method is used to update the jacobian until tolerance is met. The figure shows the results that I got compared with correct results. What improvements can I make to my code so that I am moving in the right direction(at the very least). Here is the MATLAB code.

function [t,y] = bdf1_qn_broyden(fcn,tspan,Y0,h,tol)

t = (tspan(1):h:tspan(2))';     %timestep vector
n = length(t);                  %no of steps
y = zeros(n,length(Y0));        %creting output array
y(1,:) = Y0;                    %first row of o/p is initial condition
i = 2;                              
jac = zeros(length(Y0));
unity = eye(length(Y0));

while i <= n 
    H = h*unity;
    for j = 1:1:length(Y0)
        jac(:,j) = fcn(t(i),y(i-1,:)' + H(:,j)) - fcn(t(i),y(i-1,:)');
    end
    iterate = (unity - jac)\unity;
    jac = (1/h)*jac;
    ycurr = y(i-1,:)' + h*fcn(t(i-1),y(i-1,:)');
    count = 0;
    diff = 1;
    while any(diff > tol) && (count < 25)
        curreval = fcn(t(i),ycurr);
        ynext = ycurr - iterate*(ycurr - y(i-1,:)' - h*curreval);
        delta = ynext - ycurr;
        deltaeval = fcn(t(i),ynext) - curreval;
        jac = jac + (1/dot(delta,delta))*(deltaeval - jac*delta)*delta';
        iterate = (unity - h*jac)\unity;
        diff = abs(delta);
        ycurr = ynext;
        count = count +1;
    end
    if count >= 25
    disp('Iterative method failed to converge within 25 steps');
    end
    y(i,:) = ynext';
    i = i+1;    
end
end

As you might've guessed, I am new to this and my code is quite inefficient. However, I would like to atleast be accurate and slow if I'm going to be slow. Is this not possible without use of adaptive step size? Is BDF1 inept at solving such problems? Should I choose a different method? Or should I solve my implicit function better? Is there a mistake in my implementation? I don't know which way to proceed. Please point me in the right direction. Huge Thanks!

Answer 1

This is not surprising. For one, Broyden is not very stable for highly ill-conditioned nonlinear systems, and having a stiff ODE means that the nonlinear system will be ill-conditioned. This is why stiff ODE solvers don't use Broyden but use Newton iterations (well, a special form of Quasi-Newton which is Newton iteration but sometimes reusing the Jacobian a few time steps depending on the convergence).

But secondly, stiff ODE solvers know that there is a chance for nonlinear solves to diverge for many reasons (mostly a bad starting guess for the iterations), which is why all of the ones used in production-quality codes are time adaptive. If the nonlinear solve diverges then you decrease $\Delta t$ and try again, knowing that as $\Delta t \rightarrow 0$ you have $u_{n+1} \rightarrow u_n$ and so your starting guess is getting better and better for any stable extrapolant. If you are not making use of time stepping adaptivity, well then most of these codes would not work well. If you look at something like CVODE's BDF method or ode15s, you'll see that on such highly stiff problems they can have a few time steps taking tiny steps: that's essential to the method.

For more on how stiff ODE solvers are written, see https://mitmath.github.io/18337/lecture9/stiff_odes and the discussion in Hairer II Chapter 8.

Answer 2

Since performance is not a primary concern, just use the backward Euler finite difference approximation in time, full Newton method with a direct solver at each time step, and a central difference approximation to the Jacobian.

Since the Jacobian is only a 2x2 matrix, a quasi-Newton solution approach is not useful.