Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

1D diffusion equation

  1. Jun 30, 2011 #1
    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 dont get the satisfactory answers. Could anyone help?

    % --- Define constants and initial condition

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

  2. jcsd
  3. Jul 4, 2011 #2


    User Avatar
    Homework Helper

    Can you post the results? What do they look like? What don't you like about them?
  4. Jul 5, 2011 #3
    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?

    Attached Files:

  5. Jul 5, 2011 #4


    User Avatar
    Homework Helper

    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.
  6. Jul 5, 2011 #5
    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
  7. Jul 5, 2011 #6


    User Avatar
    Homework Helper

    I have a better method, I don't have time now but I will post one later this evening.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook