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)=bA;
%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’)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)=bA;
%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’) 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)=bA;
%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’) mathematics, nonlinear, differential equations MATLAB Answers — New Questions
