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

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

[Jan]

"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.

Answer 1

[Jan]

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);
%           ^^