Central Difference Scheme

[robby991]

Hello, I have a question on the central difference scheme for differential equations. If we are evaluating at a point xi, where xi = 0, We first evaluate from xi to x(i+1), then what happens xi to x(i-1) that does not exist since it is beyond x=0? I am confused as to how this works. Any suggestions would be great, thank you.

 

[AlephZero]

You have to use the boundary conditions to "invent" a value for x(i-1).

For example if you know the slope at x(i) at the boundary, you can use the slope to form a difference equation between x(i+1) and x(i-1).

For a central difference in time, you can assume the values of x at time = -h is are same as the initial conditions at at time 0, if the errors caused by that assumption will decay as you integrate forwards in time.

All finite difference schemes have good and bad features. This issue can be one of the bad features of central difference, if it messes up the numerical solution.

 

[robby991]

Thank you, I understand now. I am trying to implement this in Matlab for solving the diffusion equation. However, I am still a little confused as to if I completely understand how to evaluate this. My initial conditions and boundary conditions are the following:

[B]IC:[/B]
S(x,0)=0   0 <= x < d
D(d,0)=S0

[B]BC[/B]

dS/dt(x=0) = 0
S(d,t)=S0

This is what I have so far for the code:

Length = 1E-6;                     %1 micron grid length
GridNum = 60;
TotalTime = 3;                     %run for 3 seconds
Dt= .2;                            
M=30;
S0 = 1E-3;
D = .1;                             %diffusion coefficient for Glucose

X = linspace(0,Length,GridNum);    %Defines the grid from 0 to L (1 micron)

Dx = X(2)-X(1);                    %Defines the grid spacing


%******Define initial condition******%

S(:,1)=0;                          %defines the initial substrate concentration at x=0
S(:,GridNum)=S0;                   %defines the initial substrate concentration at x=d
Time(1)=0;                         %t=0 at i=1

for n = 2:M
       
    for j=2:GridNum-1
        
        %%Define second derivative of concentration on the points between
        %%x=0 and x=d.
                
        DS2(j,n)=(S(GridNum-1,n-1)-2.*S(GridNum, n-1)+S(GridNum+1,n-1))./Dx.^2;
        
        %%define second derivative using Neumann boundary condition
        
        DS2(1)=2*(S(2,n-1)-S(1,n-1))./Dx.^2;
               
        %%Time step the concentration
        
        S(:,n)=S(:,n)+Dt.*DC2;
        Time(n)=Time(n)+Dt;
        
    end
    
    clf       %clear contents of figure
    plot(X,S(:,n));
    
end

Please let me know if there is a fundamental error in this. The code does not run so I assume there is a problem in that respect also. Thank you very much.

 

[Leperous]

For (4th-order) accurate time-stepping that only requires one set of initial data, try a Runge-Kutta scheme. It's easier than it looks, trust me :)

Also, I don't know how Matlab array index notation works, but it looks like you're not initializing S everywhere, and what is DC2?

 

[robby991]

I am sorry. That DC2 should be DS2. Also, it does look like I am defining or initializing S incorrectly. Matlab is giving me these errors:

? Error using ==> plus
Matrix dimensions must agree.

Error in ==> Substrate_Diffusion at 34
        S(:,n)=S(:,n)+Dt.*DS2;

I will look into the method you suggested. However, how does it differ from what I am currently doing? Any other suggestions about my code would be great. Thanks.

 

[AlephZero]

For (4th-order) accurate time-stepping that only requires one set of initial data, try a Runge-Kutta scheme. It's easier than it looks, trust me :)

Runge-Kutta with automatic timestep control is a good "ODE solver for dummies", if you know nothing about how the equation behaves and you don't care how long it takes to get a solution. That's why it is a good "default" solver in Matlab. It has one big advantange, namely that assuming your code calculates the derivatives correctly, it is very unlikely you will ever get a solution that look believable, but is actually wrong. (It's fairly easy to tell if the solution is wrong by several orders of magnitude, but harder to tell if it is wrong by say 20% or 30%).

However when you do know something about the ODE you are solving, R-K is usually a terrible choice compared with a finite difference scheme that has the right properties to match the equation you want to solve. For the diffusion equation, a good implementation of a central difference based scheme with fixed space points and variable time steps will run "real-life" problems literally hundreds or thousands of times faster than R-K, because the conditional stability of R-K limits the time step size.

 

[robby991]

Yes I feel as I have a handle on the system as I am aware of the boundary conditions I need and the Eulers step method to solve. However, I think there is an issue in the time stepping along the y axis. The x-axis seems to be coded ok.

 

[robby991]

I worked on my code a little more and I think I understand what is going on a little better interms of how the solution runs through space and time. I think I correctly defined my mesh. However I arbitrarily choose the inputs, so I don't know if there is a stability problem or other issue when the algorithm runs. Can someone please check and make sure I am evaluating the central difference scheme correctly. The only error I am getting in MatLab is:

? Error using ==> plus
Matrix dimensions must agree.

Error in ==> Substrate_Diffusion at 39
         S(:,n)=S(:,n)+(D*Dt*DS2);

My new code is:

Length = 1E-6;                     %1 micron grid length
nx = 60;                           %number of grid poitns
TotalTime = 3;                     %run for 3 seconds
nt=30;                             %number of Time Steps
S0 = 1E-3;                         %initial substrate concentration
D = .1;                            %diffusion coefficient 

                                   
X = linspace(0,Length,nx);         %Defines the grid from 0 to L (1 micron)
Y = linspace(0,TotalTime,30);      %Defines the grid from 0 to total simulation time
Dx = X(2)-X(1);                    %Defines the grid spacing in x
Dt = Y(2)-Y(1);                    %Defines grid spacing in y (time)
S = zeros(nx,nt);

%******Define initial condition******%

S(:,1)=0;                          %defines the initial substrate concentration at x=0
S(:,nx)=S0;                   %defines the initial substrate concentration at x=d
Time(1)=0;                         %t=0 at i=1

for n = 2:nt
    
    S(:,n)=S(:,n-1);
    Time(n)=Time(n-1);
                
    for j=2:nx-1
        
        %%Define second derivative of concentration on the points between
        %%x=0 and x=d.
                
        DS2(j,n)=(S(j-1,n-1)-2.*S(j,n-1)+S(j+1,n-1))./Dx.^2;
        
        %%define second derivative using Neumann boundary condition
        
        DS2(1)=2*(S(2,n-1)-S(1,n-1))./Dx.^2;
               
        %%Time step the concentration
        
         S(:,n)=S(:,n)+(D.*Dt.*DS2);
         Time(n)=Time(n)+Dt;
        
    end
    
    clf       %clear contents of figure
    plot(X,S(:,n));
    
end

Thanks so much.

 

[Edgardo]

Please describe exactly what you are trying to program.

As far as I understood you are trying to solve the diffusion equation:

S_t = D * S_xx


where S(x,t) is the concentration at point x at time t.
S_t is the time derivative, S_xx is the second derivative in space and
D is the diffusion coefficient.

You then use the explicit Euler scheme for the time step:
S(x,t+dt) = S(x,t) + dt*S_t(x,t)
= S(x,t) + dt*D*S_xx(x,t)

(Note: Choose a time step small enough otherwise you get
oscillations)


The second derivative S_xx is approximated by the [1 -2 1] stencil:
S_xx(x,t+dt) = [S(x-dx,t) - 2*S(x,t) + S(x+dx,t)]/dx^2

Plugging this into the equation from above yields:

S(x,t+dt) = S(x,t) + dt*S_t(x,t)
= S(x,t) + dt*D*S_xx(x,t)
= S(x,t) + dt*D*[S(x-dx,t) - 2*S(x,t) + S(x+dx,t)]/dx^2
= S(x,t) + D*dt/dx^2*[S(x-dx,t) - 2*S(x,t) + S(x+dx,t)]

--------------------------
I have a couple of questions regarding your code:


1) At the beginning you write:
%******Define initial condition******%
%defines the initial substrate concentration at x=0
%defines the initial substrate concentration at x=d
S(:,1)=0;
S(:,nx)=S0;

=> So, S has as first argument time and as second argument space:
S(t,x)

However, later on when you define DS(j,n) you write:
DS2(j,n)=(S(j-1,n-1)-2.*S(j,n-1)+S(j+1,n-1))./Dx.^2;

Now, the arguments are switched, i.e. the first argument is space
and the second argument is time:
S(x,t)

I gather this from your loop description at the beginning:
for j=2:nx-1 (space?)
for n = 2:nt (time?)




2) You write:
%%define second derivative using Neumann boundary condition
DS2(1)=2*(S(2,n-1)-S(1,n-1))./Dx.^2;

Shouldn't it be DS2(1,:) instead of just DS2(1)?




3) You write:
%%Time step the concentration

S(:,n)=S(:,n)+(D.*Dt.*DS2);
Time(n)=Time(n)+Dt;

First question:
Here, you use the explicit Euler scheme, but shouldn't
it be:
S(:,n) = S(:,n-1)+(D.*Dt.*DS2)

Or is this the reason why you write S(:n)=S(:n-1)
directly beneath the first loop definition?

Second question:
Why doesn't DS2 have an argument?


4) You write:
%%define second derivative using Neumann boundary condition
DS2(1)=2*(S(2,n-1)-S(1,n-1))./Dx.^2;
Shouldn't DS2 have two arguments, i.e. DS2(1,:)?

 

[robby991]

Thank you Edgardo. Yes, you defined the problem nicely. I went through your concerns and made some changes per your suggstions.

Now, the arguments are switched, i.e. the first argument is space
and the second argument is time:

I believe I correctly fixed this.

First question: 
Here, you use the explicit Euler scheme, but shouldn't 
it be: 
S(:,n) = S(:,n-1)+(D.*Dt.*DS2)

Or is this the reason why you write S(:n)=S(:n-1)
directly beneath the first loop definition?

I changed this as it was confusing, and implemented the way you suggested.

It seems my Matlab code is running fine, however there is just a fundamental mistake in the implementation of the algorthim. I think I am having problems implementing the Neumann boundary condition. Does it need two lines of code like I have it? Currently, I run the code and it generates a graph with a flat line. x-axis from 0 to 1E-6, and y-axis 0 to 1.

[CODE]
Length = 1E-6;                     %1 micron grid length
nx = 60;                           %number of grid poitns
TotalTime = 3;                     %run for 3 seconds
nt=30;                             %number of Time Steps
S0 = 1E-3;                         %initial substrate concentration
%D = 3E-10;                            %diffusion coefficient 

                                   
X = linspace(0,Length,nx);         %Defines the grid from 0 to L (1 micron)
Y = linspace(0,TotalTime,30);      %Defines the grid from 0 to total simulation time
Dx = X(2)-X(1);                    %Defines the grid spacing in x
Dt = Y(2)-Y(1);                    %Defines grid spacing in y (time)
S = zeros(nx,nt);

%******Define initial condition******%

S(:,1)=0;                          %defines the initial substrate concentration at x=0
S(:,nx)=S0;                   %defines the initial substrate concentration at x=d
Time(1)=0;                         %t=0 at i=1

for n = 2:nt
                    
    for j=2:nx-1
        
        %%Define second derivative of concentration on the points between
        %%x=0 and x=d.
                
        DS2(n,j)=(S(n-1,j-1)-2.*S(n,j-1)+S(n+1,j-1))./Dx.^2;
        
        %%define second derivative using Neumann boundary condition
        
        DS2(:,1)=2*(S(2,j-1)-S(1,j-1))./Dx.^2;
        DS2(:,nx)=2*(S(nx-1,n)-S(nx,n))./Dx^2;   
        
        %%Time step the concentration
        
         S(:,n)=S(:,n-1)+(Dt.*DS2(n,j));
         Time(n)=Time(n-1)+Dt;
        
    end
    
    clf       %clear contents of figure
    plot(X,S(:,n));
    
end

[/CODE]

 

[Edgardo]

1) You wrote:

S = zeros(nx,nt);

Here, you have defined the first argument as space and the second argument as time: S(space,time)
However, later on you define it vice versa as S(time,space).
(I find the names of your variables confusing, for example you use n for time
and j for space. I suggest using t for time and x for space.)

2) You wrote:

Time(1)=0;                         %t=0 at i=1
...
Time(n)=Time(n-1)+Dt;

Time is a vector. Don't you have to initialize it as you have done e.g. with S,X and Y?

3) You wrote:

for n = 2:nt                    
    for j=2:nx-1        
        %%Define second derivative of concentration on the points between
        %%x=0 and x=d.                
        DS2(n,j)=(S(n-1,j-1)-2.*S(n,j-1)+S(n+1,j-1))./Dx.^2;

Again, there is a problem with the variables.
Here, you are calculating the second derivative with respect to time (n is your time variable?) instead of the the second derivative in space.


4) You wrote:

%%define second derivative using Neumann boundary condition        
        DS2(:,1)=2*(S(2,j-1)-S(1,j-1))./Dx.^2;
        DS2(:,nx)=2*(S(nx-1,n)-S(nx,n))./Dx^2;

Do you also have the Neumann boundary condition for the right boundary?

 

[robby991]

Here, you have defined the first argument as space and the second argument as time: S(space,time)
However, later on you define it vice versa as S(time,space). 
(I find the names of your variables confusing, for example you use n for time 
and j for space. I suggest using t for time and x for space.)

In reply to (1) and (3), I matched the space and time arguments for the functions in the entire code. I am new to MatLab coding and this didn't occur to me right away, thanks for pointing it out.

In regards to the Time() argument in my code, I honestly am confused as to if I need to do this and how it should be coded. I already initialized and defined the "time" domain (t) long the y-axis, and is included in the S() function in the "for" loop, so do I actually need to step the "Time". I would like to hear your suggestions on this.

%%define second derivative using Neumann boundary condition        
        DS2(:,1)=2*(S(2,j-1)-S(1,j-1))./Dx.^2;
        DS2(:,nx)=2*(S(nx-1,n)-S(nx,n))./Dx^2;

I know I need a Neumann boundary for the x=0 point because this point is unknown and must be calculated as I have given, taking into account the x-1 point. For the boundary at x=nx, the boundary is given as SO, which I have defined in the beginning of the code. Is it necessary to define it with the DS2(:,nx) function I used? I am very confused at this.

Here is my code. Again, I am very confused as to the Time() function I am using and defining the boundaries. I have ignored the diffusion coefficient for now

The plot for this code gives all zeros, with x-axis ranging from 0 to 1E-6, and the y-axis from 0 to 1E-3. The x-axis is correct as it has been defined as the length, however the y-axis should not stop at 1E-3. This is S0 which is the righthand boundary at x=nx.

    Length = 1E-6;                     %1 micron grid length
nx = 60;                           %number of grid poitns
TotalTime = 3;                     %run for 3 seconds
nt=30;                             %number of Time Steps
S0 = 1E-3;                         %initial substrate concentration 
 
                                   
X = linspace(0,Length,nx);         %Defines the grid from 0 to L (1 micron)
Y = linspace(0,TotalTime,nt);      %Defines the grid from 0 to total simulation time
S = zeros(nt,nx);
Dx = X(2)-X(1);                    %Defines the grid spacing in x
Dt = Y(2)-Y(1);                    %Defines grid spacing in y (time)


%******Define initial condition******%

S(:,1)=0;                          %defines the initial substrate concentration at x=0
S(:,nx)=S0;                   %defines the initial substrate concentration at x=d
Time(1,:)=0;                         %t=0 at i=1

for t = 2:nt
                    
    for x = 2:nx-1
        
        %%Define second derivative of concentration on the points between
        %%x=0 and x=d.
                
        DS2(t,x)=(S(t-1,x-1)-2.*S(t-1,x)+S(t-1,x+1))./Dx.^2;
        
        %%define second derivative using Neumann boundary condition
        
        DS2(:,1)=2*(S(2,x-1)-S(1,x-1))./Dx.^2;
        DS2(:,nx)=2*(S(t,nx-1)-S(t,nx))./Dx^2;   
        
        %%Time step the concentration
        
         S(:,x)=S(:,x-1)+(Dt.*DS2(t,x));
         Time(t)=Time(t-1)+Dt;
        
    end
    
    clf       %clear contents of figure
    plot(X,S(nt,:));
    
end

 

[Edgardo]

I've managed to create an animation for the diffusion/heat equation with Dirichlet boundary conditions (click on the attachment). For this I've written a program in http://www.gnu.org/software/octave/" (a Matlab clone) and saved it as "heatEquation_dirichlet_pics.m" (see attachment):
(I've done this in Linux with the console version of Octave. Trying to run the file with the graphical version QtOctave resulted in errors.)

% This solves the one dimensional diffusion equation S_t = S_xx
% using the explcit Euler method

% This prevents displaying the curve for each frame
set(0, 'defaultfigurevisible', 'off');

Length = 1;                     % Define length of domain
n = 101;                        % n-1 is the number of intervals
				% n is the number of grid points
				% Note: Changing n will also change the effect
				% of the time step, i.e. how big the time step will be
				% For a big n, e.g. n=501 the change from one frame to the next one
				% will be less than for n=101
				
			 
 
X = linspace(0,Length,n);         %Defines the grid from 0 to Length 

for i=1:n			% creates a vector S and initiates it with zeros
	S(i)=0;
end		
	  
Dx = X(2)-X(1);                 %Defines the grid spacing in x


% You have to be careful with this setting. 
% If you choose Dt too big you get numerical errors!
Dt = 0.8*Dx^2;


%******Define initial condition******%
for i=2:n-1
	%S(i)=100*X(i)*(1-X(i))*sin(4*pi*X(i))^2;
	S(i)=23*X(i)*sin(2*pi*X(i))^2;	
end

% Dirichlet boundary conditions
S(1) = 0;
S(n) = 0;




% plot initial condition *******************************************************************
   clf;   
   plot(X,S);
   
   % http://www.network-theory.co.uk/docs/octave3/octave_160.html
   title ("Solution of diffusion equation S_t = S_{xx}");
   xlabel ("x");
   ylabel ("S(x)");
   
   
   %note: You have to place "axis([0, 1, 0, 25])" here but don't ask me why
   % If you place it elsewhere the axis range in y-direction will be set automatically
   % which you don't want for your animation
   axis([0, 1, 0, 20]);
   
   %http://www.krizka.net/2009/11/06/creating-animations-with-octave/ 
   %Creates png files with 5 digits as label
   
   filename=sprintf('output/%05d.png',0);
   print(filename);






%*******Time Evolution*************%
max_t = 20;
for t=1:max_t         % time step


   for j=2:n-1    % Index of vector S(j)
	S(j) = S(j) + Dt/Dx^2*(S(j-1)-2*S(j)+S(j+1));  
   end
   
   
   clf;   
   plot(X,S);
   
   % http://www.network-theory.co.uk/docs/octave3/octave_160.html
   title ("Solution of diffusion equation S_t = S_{xx}");
   xlabel ("x");
   ylabel ("S(x)");
   
   
   %note: You have to place "axis([0, 1, 0, 25])" here but don't ask me why
   % If you place it elsewhere the axis range in y-direction will be set automatically
   % which you don't want for your animation
   axis([0, 1, 0, 20]);
   
   %Creates png files with 5 digits as label
   filename=sprintf('output/%05d.png',t);
   print(filename);
   
end


   % To create a gif file showing the animation do the following:
   % Go into the folder 'output' and type in the console
   % convert -delay 100 -loop 0 *png output_dirichlet.gif

I have also plotted the total heat as a function of time. You will notice that the total heat decreases with time since I've chosen the boundaries to be zero (cold ends).

 

[Edgardo]

I've also written an Octave program for the diffusion/heat equation with Neumann boundary conditions.
(see attachment: heatEquation_neumann_pics.m).

For the Neumann boundary conditions I have chosen the derivative to be zero at the ends. This corresponds to an isolating boundary. You will notice that the temperature at the boundary does not stay zero (as opposed to the solution with the Dirichlet boundary conditions).

Also, the total heat should be constant since no heat is flowing out of our domain (isolating boundary).
To show this I've created a plot "total_heat vs time". It is not constant but fluctuates due to numerical errors.
The total heat fluctuates with -7 and +7 around the value 575. This corresponds to a relative error of 1.2 percent.

 

[robby991]

Thanks for the post Edgardo. I have one question. Why did you use a sin function as your initial condition?

for i=1:n
	%S(i)=100*X(i)*(1-X(i))*sin(4*pi*X(i))^2;
	S(i)=23*X(i)*sin(2*pi*X(i))^2;	
end

I don't understand the use/need for the sin function, could you explain it to me? Should I be using this in my code for my application?

 

[Edgardo]

I just used the sin function to have an interesting looking initial heat distribution.

Besides, for the Dirichlet boundary condition I only changed the inner points (the for loop goes from j=2 to n-1). So I assumed that the equation
S(x,t+dt) = S(x,t) + D*dt/dx^2*[S(x-dx,t) - 2*S(x,t) + S(x+dx,t)]
is valid for $$x \in (0,1)$$.

%*******Time Evolution*************%
...
   for j=2:n-1    % Index of vector S(j)
	S(j) = S(j) + Dt/Dx^2*(S(j-1)-2*S(j)+S(j+1));  
   end
...

Whereas for the Neumann boundary condition I changed every point in the domain, so here I assumed that the equation
S(x,t+dt) = S(x,t) + D*dt/dx^2*[S(x-dx,t) - 2*S(x,t) + S(x+dx,t)]
is valid for $$x \in [0,1]$$.

You can also have the case where you have mixed boundary conditions, for example the left boundary has a constant value (Dirichlet boundary condition) and right boundary has a given derivative (Neumann boundary condition).

I should also mention that to run my code you should create a folder "output" in which the png files are placed. Do the attached ".m files" also work in Matlab?

 

[robby991]

I got my code working in Matlab for a solution to the diffusion equation using the central difference scheme with mixed boundary conditions, one dirichlet boundary condition and one Neumann. It is as follows:

numx = 100;         %number of grid points in space
numt = 2000;        %number of time steps to be iterated over 
dx = 1/(numx - 1);  %Define grid spacing in time
dt = 0.00005;       %Define grid spacing in time
Length = 1E-6;      %length of grid
Ds = .5;

x = linspace(0,Length,numx);   %vector of x values, to be used for plotting
S = zeros(numx,numt);          %initialize everything to zero

%specify initial conditions

t(1) = 0;      %t=0

S(1,numx) = 0; %S=0 at x=1 -- dirichlet boundary  condition

S(:,1) = sin(pi*x/1);  

%iterate central difference equation 

for j=1:numt  
    
    %Time stepping
       
    t(j+1) = t(j) + dt;
    
    %2nd Derivative Substrate
   
    for i=2:numx-1
      S(i,j+1) = S(i,j) + (dt/dx^2)*Ds*(S(i+1,j) - 2*S(i,j) + S(i-1,j)); 
      S(1,1) = 2*(S(i,2)-S(i,1))./dx.^2;             %Neumann Boundary Condition
           
   end
end


plot(x,S(:,numt));

I am having trouble understanding the need to the sin function, sin(pi*x/L), and I was wondering if someone can clarify this for me. I noticed that if removed, the plot becomes zero. I am just very confused on how it integrates into the solution (I found it on the internet in somone elses code and in Edgardos recent post). I have also seen exponential functions. Shouldn't the differential equations be solved by the Euler stepping alone? How do we know which sinusodial function is the best to use and why?

Also, currently, the length of the diffusion medium is 1 micron (1E-6) and the diffusion coefficient is .5. Could someone explain the relations of these parameters to the time and space stepping (dx and dt)? If I reduce the diffusion coefficient how will that effect the solution?

Ultimately, the diffusion equation I have will be integrated into other equations, so I am wondering what this sin() business is all about. Thank you.

Edgardo, I tried running your program in Matlab and it didn't work. There are some syntax commands that are different. I got some ideas from your code however.

 

[robby991]

Also, I noticed that everyones code plots the "nt" column or row. WHat is the significance of plotting only these values out of the whole matrix? Thanks.

 

[Edgardo]

I've noticed that the following part of my code is wrong:

%*******Time Evolution*************%
...
   for j=2:n-1    % Index of vector S(j)
	S(j) = S(j) + Dt/Dx^2*(S(j-1)-2*S(j)+S(j+1));  
   end
...

It is wrong because I am changing the "old" vector S already when I am calculating the "new" vector S.
This is not allowed since we should calculate a totally new vector from the old vector where the old
vector must not be changed.

The correct code is:

%*******Time Evolution*************%
...
   for j=2:n-1    % Index of vector S(j)
	S_new(j) = S_old(j) + Dt/Dx^2*(S_old(j-1)-2*S_old(j)+S_old(j+1));  
   end
...
   S_old = S_new;

I've created a new animation (see attachment) and also got something very interesting. When you choose the time step too large you get huge oscillations. The keyword here is "CFL number". Here is a reference on the "CFL number": http://www.cfd-online.com/Forums/main/1592-cfl-number-insight.html" .

Oscillations occur for: Dt=0.8*Dx^2
No oscillations for: Dt=0.4*Dx^2

 

[Edgardo]

I am having trouble understanding the need to the sin function, sin(pi*x/L), and I was wondering if someone can clarify this for me. I noticed that if removed, the plot becomes zero.

You can choose any arbitrary function, not just the sin function. It really depends on your problem and what you want to have as initial "heat" distribution (or concentration distribution). What do you mean with the plot becoming zero?

Also, currently, the length of the diffusion medium is 1 micron (1E-6) and the diffusion coefficient is .5. Could someone explain the relations of these parameters to the time and space stepping (dx and dt)? If I reduce the diffusion coefficient how will that effect the solution?

Think about the meaning of the diffusion coefficient. What does a higher/lower value mean?
Note: If you change the diffusion coefficient you also have to be careful with the time step $$\Delta t$$ otherwise you could get these wild oscillations. See also my last post on the CFL number.

 

[robby991]

Oscillations occur for: Dt=0.8*Dx^2
No oscillations for: Dt=0.4*Dx^2

This phenomenon can be more clearly understood when you incorporate a diffusion coefficient, say D, into the code. The requirement for Dt with a diffusion coefficient becomes:

 D*Dt/(Dx^2)<1

The idea is you want to choose your time and space stepping so that Dt isn't too large compared to Dx. If the time stepping is too large in comparison with the space stepping, information is being pushed through the system at a faster pace than what the diffusion coefficient and partial derivatives are moving, hence the oscillations and instability.

You can choose any arbitrary function, not just the sin function. It really depends on your problem and what you want to have as initial "heat" distribution (or concentration distribution).

I understand this now.

I am having problems I think implementing the Neumann Boundary conditions still. Also, I am unsure where they should be going, into the inner for loop or out of. If I place them inside the inner for loop the peak of the sin function after the code runs is 14 x 10^247 which must be incorrect, and out of the inner for loop gives about .25. In addition to the placement, I don't think I am indexing it correctly and I am having a real hard time figuring out why.

clear all;

numx = 20;                  %number of grid points in space
numt = 1000;                %number of time steps to be iterated over 
tmax = .00000005;
Length = 1E-6;              %length of grid
Ds = 3E-6;                  %requirement Ds(dt)/dx^2 < .5

x = linspace(0,Length,numx)';   %vector of x values, to be used for plotting
t = linspace(0,tmax,numt);    %vector of t values, to be used for plotting
S = zeros(numx,numt);          %initialize everything to zero

dx = x(2)-x(1);                %Define grid spacing in time
dt = t(2)-t(1);                %Define grid spacing in time

%specify initial conditions

t(1) = 0;      %1st t position = 0

S(:,1) = sin(pi*x/Length);  

%iterate central difference equation 

for j=1:numt  
    
    %Time stepping, t(1) = 0 from IC
       
    t(j+1) = t(j) + dt;
     
     %2nd Derivative Substrate
   
    for i=2:numx-1
      S(i,j+1) = S(i,j) + (dt/dx^2)*Ds*(S(i+1,j) - 2*S(i,j) + S(i-1,j)); 
    end
    
    S(1,1) = 2*(S(i,j)-S(i-1,j))./dx.^2;           %Neumann Boundary Condition
    S(numx,1) = 2*(S(numx-1,j)-S(numx,j))./dx.^2;  %Neumann Boundary Condition
  
    plot(x,S(:,numt));    
end

 

[Edgardo]

This part is wrong:

  S(1,1) = 2*(S(i,j)-S(i-1,j))./dx.^2;           %Neumann Boundary Condition
  S(numx,1) = 2*(S(numx-1,j)-S(numx,j))./dx.^2;  %Neumann Boundary Condition

There is no index j on the left hand side and on the right hand side it ought to be j-1.
Furthermore, it should be "inner point" minus "outer point" if I'm not mistaken.
Also, why is there no dt and Ds anymore?Try the following:

  S(1,j)=S(1,j-1)+dt*Ds*2*(S(2,j-1)-S(1,j-1))./dx.^2;     %Neumann Boundary Condition
  S(numx,j)=S(numx,j-1)+dt*Ds*2*(S(numx-1,j-1)-S(numx,j-1))./dx.^2;

1) How did you derive your expression for the Neumann boundary condition?
2) What are the dots for, e.g. dx.^2?

 

[robby991]

I derived the Neumann BCs the following way. As far as I am concerned, the math is sound. But my feeling is I am not implementing it correctly. This in combination with the placement of the BCs in the loop I am unsure of.

The diffusion equation dS/dt = Ds (d^2S/dx^2)

Which has the following Central Difference Solution

S(i,j+1) = S(i,j) + (dt/dx^2)*Ds*(S(i+1,j) - 2*S(i,j) + S(i-1,j)); 

I would like Neumann Boundary conditions at each end. That is, dS/dx = 0.

To obtain an equation for this, we look at the central difference approximation to dS/dx, which is:

(S i+1 - S i-1) / dx^2

For the case where i =1 (first node in the x direction), we get the following:

[S(2)-S(0)] / dx^2 = 0. Therefore S(0) = S(2).

Thus, the second derivative central difference now becomes, 

d^2S / dx^2 = 2(S(2) - S(1)) / dx^2

This is implemented in MatLab as a boundary condition of 

S(1,1) = 2*(S(i,j)-S(i-1,j))./dx.^2;

Simularly, at a i = numx, S(numx+1) = S(numx-1), so we have

S(numx,1) = 2*(S(numx-1,j)-S(numx,j))./dx.^2;

Therefore, my boundary conditions are the following:

S(1,1) = 2*(S(i,j)-S(i-1,j))./dx.^2;
S(numx,1) = 2*(S(numx-1,j)-S(numx,j))./dx.^2;

The period after the vaiables tells Matlab to do this math on every element of that variable. Without the period Matlab tries to treat the variable array as a matrix and do matrix algebra. I think this is necessary?