Crank-Nicholson method
11 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
I am trying too apply Crank-Nicholson method for the following example;
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1708066/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1708071/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1708076/image.png)
and I have applied the following scheme
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1706916/image.png)
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')
0 Commenti
Risposta accettata
Torsten
il 1 Giu 2024
Modificato: Torsten
il 2 Giu 2024
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
0 Commenti
Più risposte (1)
Athanasios Paraskevopoulos
il 1 Giu 2024
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')
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!