Can You Help Me Solve a 1D Diffusion Equation with a FTCS Scheme?

[Juliousceasor]

I have a 1_D diffusion equation

dc/dt = D*d^2c/dx^2-Lc

where L,D = constants

I am trying to solve the equation above by following b.c. by FTCS scheme

-D*dc/dx = J0*delta(t); where delta(t)= dirac delta function ----(upper boundary)

I have written the code for it

but i just don't get the satisfactory answers. Could anyone help?

% --- Define constants and initial condition
clc;
clear;

L = 0.02; % length of domain in x direction
I0 = 0.0000829*24*60*60; % Bq m^-2 day^-1
L1 = (0.693/53.2); % Decay constant day^-1
tmax = 120; % end time
nx = 90;% number of nodes in x direction
nt = 121; % number of time steps
dx = L/(nx-1);
dt = tmax/(nt-1);
alpha= 2.5*10^-13*24*3600;
r = alpha*dt/dx^2; rr = 1 - 2*r-L1*dt;
v = 2.5*10^-6; % Convection velocity m day^-1
J0 = I0/sqrt(L1*alpha); % total inventory of Be-7 in soil
% --- Create arrays to save data for export
x = linspace(0,L,nx);
t = linspace(0,tmax,nt);
U = zeros(nx,nt);

% --- Set IC and BC
U(:,1)= 0;

% --- Loop over time steps

for m= 2:nt
U(1,m) = J0; %--- Upper boundary
U(end,m) = 0; %--- Lower boundary

for i= 2:nx-1


U(i,m) = r*U(i-1,m-1)+ rr*U(i,m-1)+ r*U(i+1,m-1);


end
end

 

[hunt_mat]

Can you post the results? What do they look like? What don't you like about them?

 

[Juliousceasor]

The numerical results are giving me very high values. I also have the analytical solution to it. I have posted a figure(please see the attachment). Blue dots in the figure is the numerical solution and the very small red line (which is hardly visible) is the analytical solution. I just do not understand what wrong with my Program. Could you help?

 

[hunt_mat]

Trying to understand the method from someones MATLAB code is hard, I use MATLAB a great deal, it's a wonderful piece of kit. I believe it has an inbuilt parabolic equation solver, if that would be of any help for you.

Can you write you algorithm out (preferably in LaTeX) so I can see what you're doing please.

I will say one thing now though, you have the lines:
U(1,m) = J0; %--- Upper boundary
U(end,m) = 0; %--- Lower boundary
inside your double loop. It is good programming practice to keep these out of the looping as for one it cuts down on computational time and secondly it may inadvertently affect your algorithm.

 

[Juliousceasor]

First boundary condition is discretized as

D*dU/dx = D* (U(i+1,m)-U(i-1,m))/(2*dx)

The main differential equation is dicrretized as

d^2U/dx^2= (U(i+1,m)-2*U(i,m)+U(i-1,m))/(dx^2)--------(1)

dU/dt = (U(i,m+1)-U(i,m))/dt------------------------(2)

Using discretizations (1) and (2) in the main differential equation and reaaranging it we get,

U(i,m+1) = r*U(i-1,m)+ rr*U(i,m)+ r*U(i+1,m); (this is the main loop of the method)

here r = D*dt/dx^2; rr = 1 - 2*r-L*dt;

I hope this helps

 

[hunt_mat]

I have a better method, I don't have time now but I will post one later this evening.