0

I'm trying to implement the Backward Differentiation Formula of order 2 in Matlab and I've been stuck for more than a day. The goal is to solve the Van der Pol equation. I can't find a functioning example for this algorithm on the internet. My attempt so far:

clear; clc; close all;

u0=[1;0]; %initial value
t0=0; tf=3000; %interval
h=0.25; %step size
N = floor((tf-t0)/h)+1; %number of points
mu = 100; %Van der Pol parameter
f = @(t,y) vanderpoldemo(t,y,mu);
t = t0:h:tf;


u = zeros(2,N);
u(:,1)=u0;


%% BDF2
u(:,2)=u(:,1)+hf(t(1),u(:,1));  %estimate starting value
jacobian = @(x) [0,1; -mux(1)x(2)-1,mu(1-x(1)^2)];
G = @(x,i) [x(1);x(2)]-4/3u(:,i-1)+1/3u(:,i-2)-2/3hf(t(i),[x(1);x(2)]);
eps = 1e-07; 
%Newton iteration
for i = 3:N
    e = 1;
    x = u(:,i-1);
    while e>eps
    A = eye(2) - 2h/3jacobian(x);
    b = G(x,i);
    temp = A\b;
    x = x-temp;
    e = norm(temp);
    end
    u(:,i) = x;
end

No matter how small I make the step size, the script simply refuses to return any acceptable results. I would greatly appreciate any help.

Othman El Hammouchi
  • 734
  • 4
  • 12
  • When I tried to run the script, including making tf shorter, it was taking unreasonably long, so I suspect your while loop is simply not converging. Try examining norm(jacobian(x)), I suspect it is quite large such that your Newton iteration is likely to fail. – Ian Apr 16 '21 at 14:33
  • @Ian that was my suspicion as well, but then how can I solve this? – Othman El Hammouchi Apr 16 '21 at 14:35
  • Well, you have $x_{k+1}=x_k - A(x_k)^{-1} b(x_k)$, that's a fixed point iteration that will only converge if the norm of the Jacobian of the nonlinear map $x \mapsto x-A(x)^{-1} b(x)$ is less than 1, which means the map $x \mapsto A(x)^{-1} b(x)$ can't deviate from the identity by too much. I suspect for mu as large as you have it, this will only be true for h so small that it is impractical to run all the way to tf=3000. But then, I don't have much intuition for what that vector field $b(x)$ looks like. – Ian Apr 16 '21 at 14:42
  • 1
    A good place to start so that you have something to analyze is to cap the number of Newton iterations and give a warning when you exit the loop due to non-convergence. Then you'll at least have data to look at, but it will probably be gibberish. – Ian Apr 16 '21 at 14:47
  • Perhaps this would be more suited to Superuser Stackexchange – Filthyscrub Apr 16 '21 at 15:40
  • 1
    Cross-posted from https://stackoverflow.com/questions/67127189/implementing-bdf2-for-odes-in-matlab. The immediate problem is that something goes wrong in G, at least running it in octave. Monitoring the computed values of b, in the second or third loop the values are not small as the theory predicts. Using the formula directly and reducing the step size makes the script run correctly. – Lutz Lehmann Apr 17 '21 at 18:49
  • @LutzLehmann you are correct, using the formula directly makes the Newton root iteration behave as desired – Othman El Hammouchi Apr 17 '21 at 23:41

0 Answers0