Crank-Nicholson method

33 views (last 30 days)
Maryam Albuhayr
Maryam Albuhayr on 30 May 2024
Edited: Torsten on 2 Jun 2024
Ran in:
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')

Sign in to answer this question.

Accepted Answer

Torsten
Torsten on 1 Jun 2024
Edited: Torsten on 2 Jun 2024
Ran in:
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

More Answers (1)

Athanasios Paraskevopoulos
Ran in:
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')
  1 Comment
Maryam Albuhayr
Maryam Albuhayr on 1 Jun 2024
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')

Sign in to comment.

Categories

Find more on Pie Charts in Help Center and File Exchange

Tags

Products


Release

R2024a

Community Treasure Hunt

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

Start Hunting!