Solve a partial differential equations using the explicit method

Hello,

I have two codes. Can anyone help me to check them and tell me what's wrong with these codes.

the first code, I'm not getting the error results neither its plot

the second code, the exacct solution and the error are not appearing, and the results of the explict equation are doubtful.

I attached the problems and the codes

Answer 1

[Abraham Boayue]

%% Cp2ex2nd

% 1. Space steps

beta = sqrt(0.03); % sqrt(0.003); sqrt(1/878);

xa = 0;

xb = 1;

dx = 1/40;

N = (xb-xa)/dx ;

x = xa:dx:xb;

%2.Time steps

ta = 0;

tb = 0.5;

dt = 1/3300;

M = (tb-ta)/dt ;

t = ta:dt:tb;

%3. Controling Parameters

R1 = (beta/dx)^2*dt;

R2 = dt/(2*dx);

%4. Define equations A, B , C and phi(x,t)

A = @(x,t) (50/3)*(x-0.5+4.95*t);

B = @(x,t) (250/3)*(x-0.5+0.75*t);

C = @(x,t) (500/3)*(x-0.375);

phi = @(x,t)(0.1*exp(-A(x,t)) +0.5*exp(-B(x,t))+exp(-C(x,t)))./...

(exp(-A(x,t)) +exp(-B(x,t))+exp(-C(x,t)));

% 5 Initial and boundary conditions

f = @(x) phi(x,t(1)); % Initial condition

g1 = @(t)phi(x(1),t); % Left boundary condition

g2 = @(t)phi(x(N+1),t); % Right boundary condition

% 6 Implementation of the explicit method

u = zeros(N+1,M+1);

u(2:N,1) = f(x(2:N)); % Put in the initial condition

u(1,:) = g1(t); % The boundary conditions, g1 and g2 at

u(N+1,:) = g2(t); % x = 0 and x = 1

for i = 2:N

u(i,2) = R1*(u(i+1,1)-2*u(i,1)+u(i-1,1))-R2*(u(i+1,1)-u(i-1,1))*u(i,1)+u(i,1);

end

for j=2:M % Time Loop

for i= 2:N % Space Loop

u(i,j+1) = R1*(u(i+1,j)-2*u(i,j)+u(i-1,j))-R2*(u(i+1,j)-u(i-1,j))*u(i,j)+u(i,j);

end

end

% Make some plots

T= 0:.1:M;

V = [];

for i= 1: length(T)

P = find(t==T(i));

V = [V P];

end

figure

subplot(131)

for i = 1:length(V)

hold on

plot(x,u(:,V(i)),'*-','linewidth',2.5,'DisplayName',sprintf('t = %1.4f',t(V(i))))

end

legend('-DynamicLegend','location','bestoutside');

a = ylabel('Heat');

set(a,'Fontsize',14);

a = xlabel('x');

set(a,'Fontsize',14);

a=title(['Using The Explicit Method - Lamda =' num2str(R1)]);

set(a,'Fontsize',16);

grid;

% disp(u(:,V)'); % Each row corresponds to a particular value of t and

% Each column corresponds to a particular value of x

% Implement the exact solution and compare it to the exact solution

Exact = @(x,t) 0.5*erfc((0.5*(x-t))./(beta*sqrt(t)))+0.5*exp(x/beta^2).*...

erfc((0.5*(x+t))./(beta*sqrt(t)));

subplot(132)

for j = 1:length(V)

hold on

plot(x,Exact(x,t(V(j))),'*-','linewidth',2.5,'DisplayName',sprintf('t = %1.4f',t(V(j))))

end

legend('-DynamicLegend','location','bestoutside');

a = ylabel('Heat');

set(a,'Fontsize',14);

a = xlabel('x');

set(a,'Fontsize',14);

a=title(' Exact solution');

set(a,'Fontsize',16);

grid;

% 6 Create the error function

[X,T] = meshgrid(x,t);

Exact2 = 0.5*erfc((0.5*(X-T))./(beta*sqrt(T)))+0.5*exp(X/beta^2).*...

erfc((0.5*(X+T))./(beta*sqrt(T)));

Error =abs(Exact2'-u);

subplot(133)

for j = 1:length(V)

hold on

plot(x,Error(:,(V(j))),'*-','linewidth',2.5,'DisplayName',sprintf('t = %1.4f',t(V(j))))

end

legend('-DynamicLegend','location','bestoutside');

a = ylabel('Heat');

set(a,'Fontsize',14);

a = xlabel('x');

set(a,'Fontsize',14);

a=title(' Error ');

set(a,'Fontsize',16);

grid;

%% Cp2Exo12nd

beta = sqrt(0.03); % sqrt(1/878);

xa = 0;

xb = 1;

dx = 1/40;

N = (xb-xa)/dx ;

x = xa:dx:xb;

%2.Time steps

ta = 0;

tb = 0.5;

dt = 1/3300;

M = (tb-ta)/dt ;

t = ta:dt:tb;

%3. Controling Parameters

R1 = (beta/dx)^2*dt;

R2 = dt/(2*dx);

% 5 Initial and boundary conditions

f = @(x) 0; % Initial condition

g1 = @(t)1; % Left boundary condition

g2 = @(t)0;

% 6 Implementation of the explicit method

u = zeros(N+1,M+1);

u(2:N,1) = f(x(2:N)); % Put in the initial condition

u(1,:) = g1(t); % Insert the left boundary condition at beginning

% Perform the inner computations

for j=2:M % Time Loop

for i= 2:N % Space Loop

u(i,j+1) = R1*(u(i+1,j)-2*u(i,j)+u(i-1,j))-R2*(u(i+1,j)-u(i-1,j))+u(i,j);

end

end

% Insert the right boundary condition at end

for j = 1:M

U = u(1,j)-2*dx*g2(t(j));

u(1,j+1) = R1*(u(2,j)-2*u(1,j)+U)-R2*(u(2,j)-U)+u(1,j);

end

% Make some plots

T= 0:.1:M;

V = [];

for j= 1: length(T)

P = find(t==T(j));

V = [V P];

end

figure

subplot(131)

for j = 1:length(V)

hold on

plot(x,u(:,V(j)),'*-','linewidth',2.5,'DisplayName',sprintf('t = %1.4f',t(V(j))))

end

legend('-DynamicLegend','location','bestoutside');

a = ylabel('Heat');

set(a,'Fontsize',14);

a = xlabel('x');

set(a,'Fontsize',14);

a=title(['Using The Explicit Method - Lamda =' num2str(R1)]);

set(a,'Fontsize',16);

grid;

% disp(u(:,V)'); % Each row corresponds to a particular value of t and

% Each column corresponds to a particular value of x

% Implement the exact solution and compare it to the exact solution

Exact = @(x,t) 0.5*erfc((0.5*(x-t))./(beta*sqrt(t)))+0.5*exp(x/beta^2).*...

erfc((0.5*(x+t))./(beta*sqrt(t)));

subplot(132)

for j = 1:length(V)

hold on

plot(x,Exact(x,t(V(j))),'*-','linewidth',2.5,'DisplayName',sprintf('t = %1.4f',t(V(j))))

end

legend('-DynamicLegend','location','bestoutside');

a = ylabel('Heat');

set(a,'Fontsize',14);

a = xlabel('x');

set(a,'Fontsize',14);

a=title(' Exact solution');

set(a,'Fontsize',16);

grid;

% 6 Create the error function

[X,T] = meshgrid(x,t);

Exact2 = 0.5*erfc((0.5*(X-T))./(beta*sqrt(T)))+0.5*exp(X/beta^2).*...

erfc((0.5*(X+T))./(beta*sqrt(T)));

Error =abs(Exact2'-u);

subplot(133)

for j = 1:length(V)

hold on

plot(x,Error(:,(V(j))),'*-','linewidth',2.5,'DisplayName',sprintf('t = %1.4f',t(V(j))))

end

legend('-DynamicLegend','location','bestoutside');

a = ylabel('Heat');

set(a,'Fontsize',14);

a = xlabel('x');

set(a,'Fontsize',14);

a=title(' Error ');

set(a,'Fontsize',16);

grid;

[Hana Bachi]

@Abraham Boayue

thank you for your answer.

I have few questions, for %% Cp2Exo12nd:

1- why did you put sqrt(0.03) instead of sqrt(1/878). While sqrt(0.33)= 0.1732 and sqrt(1/878)=0.3374

Also I have a boundary condition that I need it to implemented, which is:

dU/dX(1,t)=0

dU/dx(1,t)= Ui+1 -U1-1)/2Dx=0

so Ui+1 - Ui-1=0 .......... Ui+1=Ui-1

2- In %% Cp2ex2nd, for the analytical solution

the equaton is phi = @(x,t)(0.1*exp(-A(x,t)) +0.5*exp(-B(x,t))+exp(-C(x,t)))./ (exp(-A(x,t)) +exp(-B(x,t))+exp(-C(x,t)));

and not Exact = @(x,t) 0.5*erfc((0.5*(x-t))./(beta*sqrt(t)))+0.5*exp(x/beta^2).* erfc((0.5*(x+t))./(beta*sqrt(t))); this is why the error is large . I edited it but the error still out of the tolerance, even the graphs are differents as it shows in the attached graph.

[Abraham Boayue]

Recall that it is mentioned in the paper that an accurate approximation of solution is affordable for small beta.... This means that beta is the controling factor in order for the solution to be correct to some extend. If you use 1/878 for beta squre, you will get an unstable solution. Try it and see. Hence, beta must be a small number for the solution to converge.

[Abraham Boayue]

I'm not sure if I understand your second question. Can you send me an email like you did previously? What wrong with boundary conditions? The expressions you wrote not so clear.

[Hana Bachi]

Yes, I got unstable resoluts for 1/878 and I was wondring why.

What about dU/dX(1,t)=0. I think this condition is messing in the code

and also the equation of the exact solution for the second code should be Phi(x,t) not the one from the first problem. I got a results after modifying it but they don't look accurate escpecailly with that high error

[Hana Bachi]

I Just sent you a file

Answer 2

[Abraham Boayue]

What about dU/dX(1,t)=0. I think this condition is messing in the code.

No, it is not. Remember the last attachment I emailed you where I derived the discretized FDM equations of the PDE? The boundary condition was incorporated in the equations that I placed in the rectangular box. It is in this section of the code:

for j = 1:M

U = u(1,j)-2*dx*g2(t(j));

u(1,j+1) = R1*(u(2,j)-2*u(1,j)+U)-R2*(u(2,j)-U)+u(1,j);

end

[Hana Bachi]

Thank you for making it clearer.

did you check the figure where I changed the exact equation of the second problem? the error is large