Crank-Nicolson scheme for system of differential equations

16 views (last 30 days)
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?

Accepted Answer

Torsten
Torsten on 12 May 2022
Edited: Torsten on 12 May 2022
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
  2 Comments
Torsten
Torsten on 12 May 2022
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

Sign in to comment.

More Answers (0)

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!