Comparing implicit methods - Euler implicit, Crank Nicolson, three-time level

4 次查看（过去 30 天）
E
E 2022-3-20
回答： Jan 2022-3-20
Hello,
I'd like to compare three implicit methods- Euler implicit, Crank Nicolson, three-time level by applying to 1D advection diffusion equation.
I really need help with my codes. I can't get the right result now but I really don't know what to do.
These are parts of my codes. I'd appreciate for any help.
% this part is same for other methods
phi = zeros(n_points+1);
phi(1)=3;
phi(n_points) = 1;
Aw = zeros(n_points+1);
Ap = zeros(n_points+1);
Ae = zeros(n_points+1);
Q = zeros(n_points+1);
x(1) = 0 ;
delx = L/(n_points - 1);
for i = 2: n_points+1
x(i) = x(i-1) + delx;
end
%Euler implicit
for i = 2:n_points
Ae(i) = rho*u/(2*delx)-gamma/(delx)^2;
Aw(i) = -rho*u/(2*delx)-gamma/(delx)^2;
Ap(i) = -(Aw(i) + Ae(i))+rho/dt;
Q(n_points+1) = rho/dt*phi(i)^n_points;
end
der_phi_0 = (4*Pe/L)/(1-exp(Pe));
Aw(1) = 0;
Ap(1) = -1;
Ae(1) = 1;
Q(1) = der_phi_0*delx;
for i = 2:n_points
Ap(i) = Ap(i) - (Aw(i)*Ae(i-1)/Ap(i-1));
Q(i) = Q(i) - (Aw(i)*Q(i-1)/Ap(i-1));
end
for i = n_points-1:0:-1
Q(i)= (rho/dt)*phi(i)^n_points;
end
%Crank Nicolson
for i = 2:n_points
Ae(i) = rho*u/(4*delx)-gamma/(delx)^2;
Aw(i) = -rho*u/(4*delx)-gamma/(delx)^2;
Ap(i) = -(Aw(i) + Ae(i))+rho/dt;
Q(n_points+1) = (Aw(i)+Ae(i)+rho/dt)*phi(i)^n_points-Ae(i)*phi(i+1)^n_points-Aw(i)*phi(i-1)^n_points;
end
der_phi_0 = (4*Pe/L)/(1-exp(Pe));
Aw(1) = 0;
Ap(1) = -1;
Ae(1) = 1;
Q(1) = der_phi_0*delx;
for i = 2:n_points
Ap(i) = Ap(i) - (Aw(i)*Ae(i-1)/Ap(i-1));
Q(i) = Q(i) - (Aw(i)*Q(i-1)/Ap(i-1));
end
for i = n_points-1:0:-1
Q(i)= (Ae(i)+Aw(i)+rho/dt)*phi(i)^n_points-Ae(i)*phi(i+1)^n_points-Aw(i)*phi(i-1)^n_points;
end
% three time level
for i = 2:n_points
Ae(i) = rho*u/(2*delx)-gamma/(delx)^2;
Aw(i) = -rho*u/(2*delx)-gamma/(delx)^2;
Ap(i) = -(Aw(i) + Ae(i))+3*rho/(2*dt);
Q(n_points+1) = 2*rho/dt*phi(i)^n_points-rho/(2*dt)*phi(i)^(n_points-1);
end
der_phi_0 = (4*Pe/L)/(1-exp(Pe));
Aw(1) = 0;
Ap(1) = -1;
Ae(1) = 1;
Q(1) = der_phi_0*delx;
for i = 2:n_points
Ap(i) = Ap(i) - (Aw(i)*Ae(i-1)/Ap(i-1));
Q(i) = Q(i) - (Aw(i)*Q(i-1)/Ap(i-1));
end
for i = n_points-1:0:-1
Q(i)= (2*rho/dt)*phi(i)^n_points-(rho/(2*dt))*phi(i)^(n_points-1);
end
  1 个评论
Jan
Jan 2022-3-20
"I can't get the right result" - please explain, which problem you have with the posted code. This is more efficient than if the readers guess, what the problem is.

请先登录，再进行评论。

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

回答（1 个）

Jan
Jan 2022-3-20
A bold guess:
phi = zeros(n_points+1);
phi(1)=3;
phi(n_points) = 1;
This creates phi as square matrix. I'd assume you want a vector instead:
phi = zeros(1, n_points+1);
% ^^

类别

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

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by