Crank-Nicholson method

I am trying too apply Crank-Nicholson method for the following example;

and I have applied the following scheme

I have written the following code but the error is very big and I can't see where is my mistake

clear all 
close all 
clc
% Parameters
L = 1;             % Length of the rod
T = 1;             % Total time
Nx = 10;           % Number of spatial grid points
Nt = 50;          % Number of time steps
a = 1;      % Thermal diffusivity
dx = L / (Nx);   % Spatial step size
dt = T / Nt;       % Time step size
x = linspace(0, L, Nx);  % Spatial grid
t = linspace(0, T, Nt);  % Time grid
lx = linspace(0, L, Nx);  % Spatial grid
lt = linspace(0, T, Nt);  
Exc = @(x,t) x.^2.*exp(t);
[X,T] = meshgrid(lx,lt);
figure(1)
subplot(1,2,1)
surf(X,T,Exc(X,T))
xlabel('Space')
ylabel('Time')
zlabel('Temprature')
title('The Exact solution of heat equation ')
subplot(1,2,2)
plot(lx,Exc(lx,0),'b',lx,Exc(lx,0.02),lx,Exc(lx,0.04),lx,Exc(lx,0.06))
grid on
title('The Exact at different time level')
xlabel('x')
ylabel('u')

% Initial condition
u = zeros(Nx,Nt);
% Boundary conditions (assumed zero for simplicity)
u( 1,:) = 0;
u( end,:) = exp(t(:));
%initial condition
u(2:Nx-1,1) = x(2:Nx-1).^2; 
% Source term function
f = @(x, t) x.^2-2.*exp(t);  % Example source term
% Construct the matrices A and B
Lambda = a * dt / (2 * dx^2);
A = diag((1 + 2*Lambda) * ones(Nx-2, 1)) + diag(-Lambda * ones(Nx-3, 1), 1) + diag(-Lambda * ones(Nx-3, 1), -1);
B = diag((1 - 2*Lambda) * ones(Nx-2, 1)) + diag(Lambda * ones(Nx-3, 1), 1) + diag(Lambda * ones(Nx-3, 1), -1);
% Time-stepping loop
for n = 1:Nt-1
    % Source term at time levels n and n+1
    f_n = f(x(2:end-1), t(n));
    f_np1 = f(x(2:end-1), t(n+1));
    % Right-hand side vector
    b = B * u(2:end-1,n) + 0.5 * dt * (f_n' + f_np1');
    % Solve the linear system A * u_new = b
    u(2:end-1,n+1)=b\A;
    %u(2:end-1,n+1) = inv(A)*B*u(2:end-1, n) + 0.5 *dt*inv(A)*(f_n' + f_np1');
end
u;
% Visualization
figure(3)
subplot(2,2,1)
surf(x,t,u')
colormap('pink');
title('The approximate solution by Euler forward method for \lambda =',Lambda)
xlabel('Space')
ylabel('time')
zlabel('Temperature')
grid on
subplot(2,2,2)
plot(x,u(:,1),'o',lx,Exc(lx,0),'b',x,u(:,round(Nt/4)),'o',lx,Exc(lx,lt(round(Nt/4))),'--g',x,u(:,round(Nt/2)),'o',lx,Exc(lx,lt(round(Nt/2))), x,u(:,round(3*Nt/4)),'o',lx,Exc(lx,lt(round(3*Nt/4))),':b',x,u(:,end),'o',lx,Exc(lx,lt(end)),'LineWidth',1)
%plot(x,u(1,:),'o',lx,Exc(lx,0),'b',x,u(round(Nt/4),:),'o',lx,Exc(lx,lt(round(Nt/4))),'--g',x,u(round(Nt/2),:),'o',lx,Exc(lx,lt(round(Nt/2))), x,u(round(3*Nt/4),:),'o',lx,Exc(lx,lt(round(3*Nt/4))),':b',x,u(end,:),'o',lx,Exc(lx,lt(end)),'LineWidth',1)
title('Temperature at different time level')
xlabel('x')
ylabel('T')
grid on
legend('num sol at t=0','Exact sol at t=0','num sol at t=0.02','Exact sol at t=0.02','num sol at t=0.04','Exact sol at t=0.04','num sol at t=0.06','Exact sol at t=0.06')
subplot(2,2,3)
AE=abs(Exc(lx,lt(1))-u(:,1));
AE1=abs(Exc(lx,lt(round(Nt/4)))-u(:,round(Nt/4)));
AE2=abs(Exc(lx,lt(round(Nt/3)))-u(:,round(Nt/3)));
AE3=abs(Exc(lx,lt(round(Nt/2)))-u(:,round(Nt/2)));
AE4=abs(Exc(lx,lt(round(3*Nt/4)))-u(:,round(3*Nt/4)));
AE5=abs(Exc(lx,lt(end))-u(:,end));
lt1=lt(1);
plot(x,AE,'--',x,AE1,'o',x,AE3,'b',x,AE4,'g',x,AE5,'LineWidth',1)
title('Absolute error at different time level')
grid on
legend('Error at t=0','Error at 1/4 of the time','Error at 1/2 of the time','Error at 3/4 of the time','Error at end of the time')
subplot(2,2,4)
RE=AE./Exc(lx,lt(1));
RE1=AE1./Exc(lx,lt(round(Nt/4)));
RE2=AE2./Exc(lx,lt(round(Nt/3)));
RE3=AE3./Exc(lx,lt(round(Nt/2)));
RE4=AE4./Exc(lx,lt(round(3*Nt/4)));
RE5=AE5/Exc(lx,lt(end));
plot(x,RE,x,RE1,'--',x,RE3,'*',x,RE4,'b',x,RE5,'r','LineWidth',1)
legend('Error at t=0','Error at 1/4 of the time','Error at 1/2 of the time','Error at 3/4 of the time','Error at end of the time')
grid on
title('Relative error at different time level')

Answer 1

[Torsten]

I suggest the following code:

xstart = 0;
xend = 1;
tstart = 0;
tend = 1;
nx = 10;
nt = 100;
x = linspace(xstart,xend,nx);
dx = x(2)-x(1);
t = linspace(tstart,tend,nt);
dt = t(2) - t(1);
a = 1.0;   % thermal diffusivity
lambda = a*dt/dx^2;

u = zeros(nt,nx);
u(1,:) = x.^2;

A = zeros(nx);
A(1,1) = 1.0;
for ix = 2:nx-1
    A(ix,ix-1) = -0.5*lambda;
    A(ix,ix) = 1+lambda;
    A(ix,ix+1) = -0.5*lambda; 
end
A(nx,nx) = 1.0;
dA = decomposition(A);

b = zeros(nx,1);
for it = 2:nt
  b(1) = 0.0;
  b(2:nx-1) = u(it-1,2:nx-1) + ...
              0.5*lambda*(u(it-1,1:nx-2)-2*u(it-1,2:nx-1)+u(it-1,3:nx)) + ...
              dt*(x(2:nx-1).^2 - 2)*0.5*(exp(t(it-1))+exp(t(it)));
  b(nx) = -u(it-1,nx) + (exp(t(it-1))+exp(t(it)));
  u(it,:) = (dA\b).';
end

figure(1)
plot(x,[u(10,:)-x.^2*exp(t(10));u(25,:)-x.^2*exp(t(25));u(50,:)-x.^2*exp(t(50));u(75,:)-x.^2*exp(t(75));u(100,:)-x.^2*exp(t(100))])
grid on

y0 = x.^2;
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
[T,Y] = ode15s(@(t,y)fun(t,y,x.',nx,dx,a),t,y0,options);
figure(2)
plot(x,[Y(10,:)-x.^2*exp(t(10));Y(25,:)-x.^2*exp(t(25));Y(50,:)-x.^2*exp(t(50));Y(75,:)-x.^2*exp(t(75));Y(100,:)-x.^2*exp(t(100))])
grid on

function dy = fun(t,y,x,nx,dx,a)
  dy = zeros(nx,1);
  dy(1) = 0.0;
  dy(2:nx-1) = a*(y(1:nx-2)-2*y(2:nx-1)+y(3:nx))/dx^2 + (x(2:nx-1).^2 - 2)*exp(t);
  dy(nx) = exp(t);
end

Answer 2

[Athanasios Paraskevopoulos]

The main issue with your code is in the line where you solve the linear system A⋅unew​=b. The correct approach is to solve for 𝑢𝑛𝑒𝑤unew​ using the inverse of A or directly solve the linear system using matrix division, but the syntax and the order of operations are incorrect. The correct MATLAB syntax should use the backslash operator to solve the system.

clear all 
close all 
clc

% Parameters
L = 1;             % Length of the rod
T = 1;             % Total time
Nx = 10;           % Number of spatial grid points
Nt = 50;           % Number of time steps
a = 1;             % Thermal diffusivity
dx = L / (Nx);     % Spatial step size
dt = T / Nt;       % Time step size
x = linspace(0, L, Nx);  % Spatial grid
t = linspace(0, T, Nt);  % Time grid
lx = linspace(0, L, Nx);  % Spatial grid
lt = linspace(0, T, Nt);  
Exc = @(x,t) x.^2.*exp(t);
[X,T] = meshgrid(lx,lt);

figure(1)
subplot(1,2,1)
surf(X,T,Exc(X,T))
xlabel('Space')
ylabel('Time')
zlabel('Temperature')
title('The Exact solution of heat equation ')
subplot(1,2,2)
plot(lx,Exc(lx,0),'b',lx,Exc(lx,0.02),lx,Exc(lx,0.04),lx,Exc(lx,0.06))
grid on
title('The Exact at different time levels')
xlabel('x')
ylabel('u')

% Initial condition
u = zeros(Nx,Nt);
% Boundary conditions
u(1,:) = 0;
u(end,:) = exp(t(:));
% Initial condition
u(2:Nx-1,1) = x(2:Nx-1).^2; 

% Source term function
f = @(x, t) x.^2 - 2.*exp(t);  % Example source term

% Construct the matrices A and B
Lambda = a * dt / (2 * dx^2);
A = diag((1 + 2*Lambda) * ones(Nx-2, 1)) + diag(-Lambda * ones(Nx-3, 1), 1) + diag(-Lambda * ones(Nx-3, 1), -1);
B = diag((1 - 2*Lambda) * ones(Nx-2, 1)) + diag(Lambda * ones(Nx-3, 1), 1) + diag(Lambda * ones(Nx-3, 1), -1);

% Time-stepping loop
for n = 1:Nt-1
    % Source term at time levels n and n+1
    f_n = f(x(2:end-1), t(n));
    f_np1 = f(x(2:end-1), t(n+1));
    % Right-hand side vector
    b = B * u(2:end-1,n) + 0.5 * dt * (f_n' + f_np1');
    % Solve the linear system A * u_new = b
    u(2:end-1,n+1) = A \ b;  % Corrected to use backslash operator
end

% Visualization
figure(3)
subplot(2,2,1)
surf(x,t,u')
colormap('pink');
title('The approximate solution by Crank-Nicholson method for \lambda =',Lambda)
xlabel('Space')
ylabel('time')
zlabel('Temperature')
grid on

subplot(2,2,2)
plot(x,u(:,1),'o',lx,Exc(lx,0),'b',x,u(:,round(Nt/4)),'o',lx,Exc(lx,lt(round(Nt/4))),'--g',x,u(:,round(Nt/2)),'o',lx,Exc(lx,lt(round(Nt/2))), x,u(:,round(3*Nt/4)),'o',lx,Exc(lx,lt(round(3*Nt/4))),':b',x,u(:,end),'o',lx,Exc(lx,lt(end)),'LineWidth',1)
title('Temperature at different time levels')
xlabel('x')
ylabel('T')
grid on
legend('num sol at t=0','Exact sol at t=0','num sol at t=0.02','Exact sol at t=0.02','num sol at t=0.04','Exact sol at t=0.04','num sol at t=0.06','Exact sol at t=0.06')

subplot(2,2,3)
AE = abs(Exc(lx,lt(1)) - u(:,1));
AE1 = abs(Exc(lx,lt(round(Nt/4))) - u(:,round(Nt/4)));
AE2 = abs(Exc(lx,lt(round(Nt/3))) - u(:,round(Nt/3)));
AE3 = abs(Exc(lx,lt(round(Nt/2))) - u(:,round(Nt/2)));
AE4 = abs(Exc(lx,lt(round(3*Nt/4))) - u(:,round(3*Nt/4)));
AE5 = abs(Exc(lx,lt(end)) - u(:,end));
plot(x,AE,'--',x,AE1,'o',x,AE3,'b',x,AE4,'g',x,AE5,'LineWidth',1)
title('Absolute error at different time levels')
grid on
legend('Error at t=0','Error at 1/4 of the time','Error at 1/2 of the time','Error at 3/4 of the time','Error at end of the time')

subplot(2,2,4)
RE = AE ./ Exc(lx,lt(1));
RE1 = AE1 ./ Exc(lx,lt(round(Nt/4)));
RE2 = AE2 ./ Exc(lx,lt(round(Nt/3)));
RE3 = AE3 ./ Exc(lx,lt(round(Nt/2)));
RE4 = AE4 ./ Exc(lx,lt(round(3*Nt/4)));
RE5 = AE5 ./ Exc(lx,lt(end));
plot(x,RE,x,RE1,'--',x,RE3,'*',x,RE4,'b',x,RE5,'r','LineWidth',1)
legend('Error at t=0','Error at 1/4 of the time','Error at 1/2 of the time','Error at 3/4 of the time','Error at end of the time')
grid on

title('Relative error at different time levels')

[Maryam Albuhayr]

The error is still very large, I realized this after submitting the question, this was n't the main issue.

I have fixed it as the following, but I still look for a better code to have less error as this method is well known.

clear all

close all

clc

% Parameters

L = 1; % Length of the rod

T = 1; % Total time

h=0.1; % Spatial step size

k=0.01; % Time step size % Number of time steps

a = 1; % Thermal diffusivity

Nx = round(L/h+1);

Nt = round(T/k+1);

x = linspace(0, L, Nx); % Spatial grid

t = linspace(0, T, Nt); % Time grid

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%

lx = linspace(0, L, Nx); % Spatial grid

lt = linspace(0, T, Nt);

Exc = @(x,t) x.^2.*exp(t) %x.*t;

[X,T] = meshgrid(lx,lt);

figure(1)

subplot(1,2,1)

surf(X,T,Exc(X,T))

xlabel('Space')

ylabel('Time')

zlabel('Temprature')

title('The Exact solution of heat equation ')

subplot(1,2,2)

plot(lx,Exc(lx,0),'b',lx,Exc(lx,0.02),lx,Exc(lx,0.04),lx,Exc(lx,0.06))

grid on

title('The Exact at different time level')

xlabel('x')

ylabel('u')

%%

% Initial condition

u = zeros(Nt,Nx);

% Boundary conditions (assumed zero for simplicity)

u( :,1) = 0;

u( :,end) = exp(t(:));

%initial condition

u(1,2:Nx-1) = x(2:Nx-1).^2;

% Source term function

f = @(x,t) (x.^2-2).*exp(t); % Example source term

% Construct the matrices A and B

d= a * k / ( h^2);

A = diag((2 + 2*d) * ones(Nx-2, 1)) + diag(-d * ones(Nx-3, 1), 1) + diag(-d * ones(Nx-3, 1), -1);

F=zeros(Nx-2,1);

for n=1:Nt-1

f_n = f(x(2:end-1),t(n));

f_np1 = f(x(2:end-1),t(n+1));

F =k * (f_n' + f_np1');

end

F;

B=zeros(Nx-2,1);

for n=1:Nt-1

B(1)=(2-2*d)*u(n,2)+d*u(n,3)+d*(u(n,1)+u(n+1,1));

for i=2:Nx-3

B(i)=d*u(n,i)+(2-2*d)*u(n,i+1)+d*u(n,i+2);

end

B(Nx-2)=d*u(n,Nx-2)+(2-2*d)*u(n,Nx-1)+d*(u(n,Nx)+u(n+1,Nx));

u(n+1,2:Nx-1)=A\(B+F);

end

u

% % Visualization

figure(3)

subplot(2,2,1)

surf(t,x,u')

colormap('pink');

title('The approximate solution by Euler forward method for \lambda =',d)

xlabel('Space')

ylabel('time')

zlabel('Temperature')

grid on

subplot(2,2,2)

plot(x,u(1,:),'o',lx,Exc(lx,0),'b',x,u(round(Nt/4),:),'o',lx,Exc(lx,lt(round(Nt/4))),'--g',x,u(round(Nt/2),:),'o',lx,Exc(lx,lt(round(Nt/2))),...

x,u(round(3*Nt/4),:),'o',lx,Exc(lx,lt(round(3*Nt/4))),':b',x,u(end,:),'o',lx,Exc(lx,lt(end)),'LineWidth',1)

title('Temperature at different time level')

xlabel('x')

ylabel('T')

grid on

legend('num sol at t=0','Exact sol at t=0','num sol at 1/4 of the time','Exact sol at 1/4 of the time','num sol at 1/2 of the time',...

'Exact sol at 1/2 of the time','num sol at 3/4 of the time',...

'Exact sol at 3/4 of the time','num sol at end of the time','Exact sol at end of the time')

subplot(2,2,3)

AE=abs(Exc(lx(2:end),lt(1))-u((2:end),1));

AE1=abs(Exc(lx(2:end),lt(round(Nt/4)))-u(round(Nt/4),(2:end)));

AE2=abs(Exc(lx(2:end),lt(round(Nt/3)))-u(round(Nt/3),(2:end)));

AE3=abs(Exc(lx(2:end),lt(round(Nt/2)))-u(round(Nt/2),(2:end)));

AE4=abs(Exc(lx(2:end),lt(round(3*Nt/4)))-u(round(3*Nt/4),(2:end)));

AE5=abs(Exc(lx(2:end),lt(end))-u(end,(2:end)));

plot(x(2:end),AE,'--',x(2:end),AE1,'o',x(2:end),AE3,'b',x(2:end),AE4,'g',x(2:end),AE5,'LineWidth',1)

title('Absolute error at different time level')

grid on

legend('Error at t=0','Error at 1/4 of the time','Error at 1/2 of the time','Error at 3/4 of the time','Error at end of the time')

subplot(2,2,4)

%RE=AE./Exc(lx(2:end),lt(1));

RE1=AE1./abs(Exc(lx(2:end),lt(round(Nt/4))));

RE2=AE2./abs(Exc(lx(2:end),lt(round(Nt/3))));

RE3=AE3./abs(Exc(lx(2:end),lt(round(Nt/2))));

RE4=AE4./abs(Exc(lx(2:end),lt(round(3*Nt/4))));

RE5=AE5/abs(Exc(lx(2:end),lt(end)));

plot(x(2:end),RE1,'--',x(2:end),RE3,'*',x(2:end),RE4,'b',x(2:end),RE5,'r','LineWidth',1)

legend('Error at 1/4 of the time','Error at 1/2 of the time','Error at 3/4 of the time','Error at end of the time')

grid on

title('Relative error at different time level')