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

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!

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.

• Hairier II is a great resource. They've even taken Van der Pol's equation as their illustrative problem. I don't know how I missed that. Thanks. BTW, I think you meant Chapter 4, Section 8. (Hairier II has only 7 chapters) Oct 16 at 5:27

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.

• Is the finite(or central) difference approximation to the Jacobian at every Newton iteration superior to Broyden's method? Oct 15 at 15:14
• Usually. With an accurate Jacobian you can expect quadratic convergence with Newton's method. Oct 15 at 16:15