Crank-Nicholson method

47 次查看（过去 30 天）
Maryam Albuhayr
Maryam Albuhayr 2024-5-30
编辑： Torsten 2024-6-2
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')

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2024-6-1
编辑：Torsten 2024-6-2
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

更多回答（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 个评论
Maryam Albuhayr
Maryam Albuhayr 2024-6-1
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')

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Pie Charts 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by