Soving PDE using explicit method

I need to solve a PDE using the explict method. I created a code by the plot didn't look god.So I could'nt move forward for the following question.

Could anyone check it and tell me what's wrong?

Answer 1

[Davide Masiello]

Hi, I have edited your code as shown below

%% Case 1. Numerical solution of the PDE Ut=-Ux+Uxx-U using the explicit FDM
%  Written by Viridiana Salazar
%  Edited by Davide Masiello
clear,clc

%% Parameters
B2 = 1/878;

%% Space steps
Lx = 1;                                                                     % Max value of x
dx = 1/40;                                                                 % Space step
N  = Lx/dx+1;                                                               % Number of space steps
x  = 0:dx:Lx;

%% Time steps
tf  = 0.5;                                                                   % Final time
dt  = 0.9*0.5*dx^2;                                                          % Time step such that dt/dx^2<=0.5
t   = 0:dt:tf;
M   = length(t);                                                               % Number of time of steps

%% Initial and Boundary Conditions
U       = zeros(N,M);                                                       % Create a Matrix of NxM
U(:,1)  = 0;                                                                % Initial conditions
U(1,:)  = 1;                                                                % Boundary conditions at the left side

%% Explicit Method
for j = 1:M-1                                                               %Time Loop
   U(2:end-1,j+1) = -(U(2:end-1,j)-U(1:end-2,j))*dt/(dx)+B2*(U(3:end,j)-2*U(2:end-1,j)+U(1:end-2,j))*dt/(dx^2)+U(2:end-1,j);
   U(end,j)       = U(end-1,j);
end

%% Analytical Solution
Ur      = zeros(N,M);
Ur(:,1) = 0;                                                                % Initial conditions
Ur(1,:) = 1;                                                                % Boundary conditions at the left side
for j = 1:M-1                                                               %Time Loop
   Ur(:,j+1) = 0.5*erfc((x-t(j+1))/(2*B2*t(j+1)^0.5))+0.5*exp(x/B2).*erfc((x+t(j+1))/(2*B2*t(j+1)^0.5));
end 

%% Error
Error = abs(Ur-U);

%% Plots
figure(1)
for idx = 1:50:length(t)
    
    subplot(3,1,1)
    plot(x,U(:,idx));
    hold on
    xlabel('x-values')
    ylabel('U-values')
    
    subplot(3,1,2)
    plot(x,Ur(:,idx));
    hold on
    title('Analytical Solution')
    xlabel('x-values')
    ylabel('U-values')
    hold on
    
    subplot(3,1,3)
    plot(x,Error(:,idx));
    title('Error')
    xlabel('x-values')
    ylabel('Error')
    hold on
    
end

EXPLANATION

First, you were mis-calculating B, because you wrote B = (1/878)^1/2, which first raises the content of the brackets to the power of 1 and then divides the result by 2. I avoided any mistake by just defining B^2 directly.

This was your only actual mistake, the rest of the code was fine.

However, I have changed a couple of things.

1. I have set the value dt to be computed automatically from the value of dx (using stability criterion). In this way, you can change the value of dx at your discretion without having to manually compute the new dt.
2. I have eliminated the inner loop in the PDE solving, because you can code spatial discretization making intelligent use of indexing (see code).
3. I have approximated the first order spatial derivative with a first order accurate finite difference using an upwind scheme, which is notoriously better for convective terms.

With these changes, you get the following plot

So numerical and analytical solutions do not match and there's a large error.

However, I think there might be some problem with how the analytical solution is formulated.

In fact, given the presence of the second-order spatial derivative (i.e., diffusion), I would expect a behaviour like the one yielded by the numerical integration, while the analytical solution seems to only advect the initial profile.

Could you double check the correctness of the analytical solution?

[Hana Bachi]

Hello Davide,

The equation of the analytical solution is at it shows in ( i) in the file named CP2.

About checking the analytical solution, I'm waiting to see the results of others in my group,

Answer 2

[Hana Bachi]

Hello davide,

Thank you for edditing my code, and sepcially where you put B^2, that helps me to avoid mis-calculating it.

I have a question about dt, why did you write it this way (dt = 0.9*0.5*dx^2)? what are the 0.9 and 0.5?

Thank you

[Hana Bachi]

Acctually I added the analytical solution, the caluculation of the error, and the plots of them. But while runing the code I got a message says:

Error using ^ (line 52)

Incorrect dimensions for raising a matrix to a power. Check that the matrix is square and the power is a scalar. To perform

elementwise matrix powers, use '.^'.

Error in Cp2Exo1nd (line 31)

Ur(i,j) = 0.5*erfc((x(i)-t(j))/(2*B2*t^0.5))+0.5*exp(x(i)/t(j)^2)*erfc((x(i)+t(j))/(2*B2*t(j)^0.5))

[Davide Masiello]

Hi @Hana Bachi, I have edited the answer so now it includes also the analytical solution.

Regarding your question:

The stability criterion is

but of course you can't specify an inequaity on matlab, so I have just set to be slightly less than , e.g.