Crank-Nicolson scheme for system of differential equations

I am currently working on the numerical time integration of a system of differential equations

function dydt = odeSystem(t,y)

dydt(1,1) = y(2);

dydt(2,1) = -y(1);

end

I want to implicitly find the solution using Crank-Nicolson and a Newton-Raphson scheme. I have the following code:

if Method == 'CNM' | Method == 'All'

for i5 = 1:length(x)-1 % Over Time

% Forward Euler for initial guess

k1 = odeSystem(x(i5),y(:,i5));

y_guess = y(:,i5) + k1*h ;

% Newton Raphson to estimate n+1

error = 1;

tolerance = 1e-3;

y_new = y_guess;

y_old = y(:,i5);

count = 0;

while error >= tolerance

y_new;

f = y_new - (y_old + h*odeSystem(x(i5+1),y_new));

J = [1 -h;

h 1];

y_iter = y_new - J\f;

error = max(abs(y_iter - y_new)); % (local maximum) absolute error

y_new = y_iter;

count = count+1;

end

y(:,i5+1) = y_new;

% Overwrite k1 and add k2

k1 = odeSystem(x(i5),y(:,i5));

k2 = odeSystem(x(i5+1),y(:,i5+1));

y(:,i5+1) = y(:,i5) + (k1+k2)*(h/2);

count

end

subplot(212)

hold on

plot(x,y(1,:),'r')

plot(x,y(2,:),'b')

plot(T,Y(:,1),'r--') % Solution according to ODE45

plot(T,Y(:,2),'b--') % Solution according to ODE45

axis([0 100 -1 1])

%legend('Crank Nicolson 1','Crank Nicolson 2','Exact 1','Exact 2')

title(['Crank Nicolson with h = ',num2str(h)])

grid on;

hold off

end

I use initial conditions

y(:,1) = [0.2;0]; % Initial conditions

, a time step that is h = 0.125 for a total time of x = 0:h:100 seconds. My Crank-Nicolson scheme is blowing up, as the solution (2 pure sin/cos waves) grows in amplitude. In this same script I also performed Heun's method and they are remarkably alike.

I know that for this specific system, this scheme should be unconditionally stable, thus the exploding of the solution indicatest a mistake in my implementation.

Do you know where my mistake is and how to fix this problem?

Answer 1

[Torsten]

  y0(1,1) = 0.2;
  y0(2,1) = 0;
  x = 0:0.125:100;
  h = x(2)-x(1);
  y = zeros(2,numel(x));
  y(:,1) = y0;
  for i4 = 1:numel(x)-1
     % Forward Euler for initial guess
     k1 = odeSystem(x(i4),y(:,i4));
     y_guess = y(:,i4) + k1*h ;
     % Newton Raphson to estimate n+1
     error = 1;
     tolerance = 1e-6;
     y_new = y_guess;
     y_old = y(:,i4);
     while error >= tolerance
        f  = y_new - y_old - h/2*(odeSystem(x(i4),y_old)+odeSystem(x(i4+1),y_new));
        J = [1 -h/2;
             h/2 1];
        y_iter = y_new - J\f;
        error = max(abs(y_iter - y_new));                   % (local maximum) absolute error
        y_new = y_iter;
     end
     y(:,i4+1) = y_new;
  end
  plot(x,y(1,:),x,y(2,:))

function dydt = odeSystem(t,y)
  dydt(1,1) = y(2);
  dydt(2,1) = -y(1);
end

[Felix Tieleman]

Thank you for your help!

[Torsten]

f is always the definition of the method you try to implement.

Euler-implicit:

f = y_new - y_old - h*odeSystem(x(i4+1),y_new)

corresponds to

(y_new - y_old)/h = f(x_new,y_new)

Crank-Nicolson:

f  = y_new - y_old - h/2*(odeSystem(x(i4),y_old)+odeSystem(x(i4+1),y_new));

corresponds to

(y_new - y_old)/h = (f(x_old,y_old)+f(x_new,y_new))/2