Skip to main content

NNESMO: Session 5 – Stability and numerical diffusion

NNESMO: Session 5 – Stability and numerical diffusion

Dr Simon A Mathias

Department of Earth Sciences

Durham University

Contents

This session’s learning objectives

At the end of this session you should be able to:

  • Derive multiple finite difference solutions for the advection diffusion equation.
  • Speculate about stability criteria by assessing the sign of finite difference coefficients.
  • Calculate numerical diffusion due to different finite difference schemes using Taylor expansions.

The advection diffusion equation (ADE)

Most of the problems considered up till now can be thought of as representing special forms of the so-called advection diffusion equation (ADE)

___(1)

where is the dependent variable, is distance, is time and and are constant parameters.

We specifcally solved the ADE when considering the previous pollution transport scenario. Recall it was found that the ADE exhibited some numerical instabilities when high values of were applied. The reason that the ADE is challenging to solve in this context is due to its composite hyperbolic and parabolic nature.

When , the ADE reduces to a parabolic partial differential equation (PDE) identical to the groundwater flow equation, chemical diffusion equation and heat conduction equation considered previously. A parabolic equation is characterized by having all time derivatives of first-order and all spatial derivatives of second-order. Such equations are relatively straightforward to solve using finite difference methods.

When , the ADE reduces to a hyperbolic PDE, similar to the wave equation. Hyperbolic equations are characterised by having both time and space derivatives of the same order. Hyperbolic equations can be solved analytically using the so-called “Method of Characteristics” (MOC).

When both and are non-zero, the ADE cannot be solved using MOC and it is also challenging to solve using finite differences. Situations when both and are non-zero frequently come up in Earth Science. This is a particular problem when seeking to solve diffusion type equations, which have spatially varying diffusion coefficients. For example, consider the situation:

In this session we will examine why such difficulties occur and explore how some of the consequences can be dealt with.

For the next few examples, the following initial and boundary conditions are considered:

where and represent specified values of and is the length of the domain.

For the special case where , the above system of equations can be solved exactly to obtain

To evaluate the analytical solution, type in the following code in a single script-file and save it in an appropriate place as “ADEStudy.m”.

 function ADEStudy
 %A sequence of different subfunctions to explore stability and
 %numerical diffusion associated with finite difference solutions of
 %the advection diffusion equation.
 %Define model parameters
 L=2; %m
 a=1; %m/day
 b=0.01; %m^2/day
 uI=0; %mg/l
 u0=1; %mg/l
 %Define location of points for plotting
 x=linspace(0,L,50);
 %Define values of time to be used for plotting
 t=[1];
 %Evaluate the exact solution
 uA=AnalSol(x,t,a,b,uI,u0);
 %Plot results
 figure(1)
 hold off
 plot(x,uA,'ko')
 xlabel('Distance (m)')
 ylabel('Concentration (mg/l)')
 title('Concentration profile after one day')
 %**********************************************************************
 function u=AnalSol(x,t,a,b,uI,u0)

 %Subfunction containing the analytical solution
 [x,t]=ndgrid(x,t);
 Pe=a*x/b;
 F1=erfc((x-a*t)./2./sqrt(b*t))/2;
 F2=exp(Pe).*erfc((x+a*t)./2./sqrt(b*t))/2;
 F=F1+F2;
 ind=Pe>700;
 F(ind)=F1(ind);
 u=(u0-uI)*F+uI;

The outline of the remainder of this session is as follows. First we will derive a set of finite difference approximations using Taylor expansions. Following on from this we will develop four different finite difference solutions to the ADE along with associated stability criteria. We will then use the Taylor expansions to derive expressions for the numerical diffusion associated with each of the finite difference schemes. Finally, we will develop a higher order time integration scheme using ode45 and compare this with the earlier first-order explicit and implicit time-stepping schemes.

Taylor expansions and finite difference

Consider a set of discrete points in time and space, , where and represent the i-th point in space and n-th point in time, respectively. It is assumed that these discrete points are separated by uniform intervals in both space and time of and , respectively. The values of a function, , at the locations of these discrete points can be written as .

It is possible to derive finite difference approximations for the derivatives of with respect to and by considering the Taylor expansion:

from which it can be seen that:

It is also possible to state that

from which it can be seen that:

Interestingly, it can also be seen that:

Rearranging the above set of equations leads to the following set of finite difference approximations:

___(2a)

___(2b)

___(2c)

___(2d)

___(2e)

ETCS – Explicit time-stepping with central differences in space

First we will develop an Euler explicit time-stepping scheme for the ADE problem described above and use a central difference approximation for the term. Ignoring all the truncation errors associated with the Taylor expansions, the ADE can then be written as

___(3a)

which can be rearranged to get

where

Hereafter this scheme will be referred to as the ETCS scheme.

For numerical stability of explicit time-stepping schemes to be ensured it is generally found that the coefficients, , and , need always to be positive. Given that , , and should also always be positive, it follows that there are two stability criteria for this scheme:

where and are often referred to as the Courant number and Peclet number, respectively.

To evaluate the ETCS scheme add the following subfunction to ADEStudy.m:

 %**********************************************************************

 function u=FDSol(xI,tI,a,b,uI,u0,Pe,Cr)

 %Subfunction containing the finite difference solutions
 %Determine the space-step and time-step using a specified
 %Courant and Peclet number
 dx=Pe*b/a;
 dt=Cr*dx^2/b;
 %Determine number of nodes to solve for
 Nx=ceil((xI(end)-xI(1))/dx)+1;
 Nt=ceil(tI(end)/dt)+1;
 %Determine the location of solution points in x and t
 x=[0:dx:dx*(Nx-1)]+xI(1);
 t=[0:dt:dt*(Nt-1)];
 %Initialise the solution matrix
 u=zeros(Nx,Nt);
 %Apply initial condition
 u(:,1)=uI;
 %Apply boundary conditions
 u(1,:)=u0;
 u(Nx,:)=uI;
 %Determine A, B and C coefficients
 A=a*dt/2/dx+b*dt/dx^2;
 B=1-2*b*dt/dx^2;
 C=-a*dt/2/dx+b*dt/dx^2;
 %Determine Jacobian matrix
 J=spdiags(ones(Nx,1)*[A B C],[-1 0 1],Nx,Nx);
 %Apply boundary conditions
 J(1,:)=0;
 J(1,1)=1;
 J(Nx,:)=0;
 J(Nx,Nx)=1;
 %Solve problem
 for n=1:Nt-1
 u(:,n+1)=J*u(:,n);
 end
 %Interpolate solution to plotting points
 u=interp2(t,x',u,tI,xI');

Note that the values of and are determined by considering a specified Peclet and Courant number. Also note that the solution is then interpolated to the values of and specified by the input arguments of the sub-function. This removes the need to plot all the solution points solved for, which can be excessive depending on the Peclet number and Courant number specified.

It is also interesting to see that the , and coefficients represent the values of the first, second and third diagonal of the Jacobian matrix, , defined in this case by the formula

To evaluate the ETCS scheme, replace:

 %Plot results
 figure(1)
 hold off
 plot(x,uA,'k')
 xlabel('Distance (m)')
 ylabel('Concentration (mg/l)')
 title('Concentration profile after one day')
in the main function of ADEStudy.m with:
 %Specify a Peclet number and Courant number
 Pe=2;
 Cr=0.5;
 %Evaluate the finite difference solution
 uFD=FDSol(x,t,a,b,uI,u0,Pe,Cr);
 %Plot results
 figure(1)
 hold off
 plot(x,uA,'k')
 hold on
 plot(x,uFD)
 xlabel('Distance (m)')
 ylabel('Concentration (mg/l)')
 title('Concentration profile after one day')
 legend('Analytical solution','Finite difference solution')

Run the script-file and compare the results from the analytical solution and the finite difference solution. Experiment by varying the specified value of and . Do the pre-determined stability criteria work? Does stability also ensure accuracy?

ITCS – Implicit time-stepping with central differences in space

Here we will develop an Euler implicit time-stepping scheme and use a central difference approximation for the term. Ignoring all the truncation errors associated with the Taylor expansions, the ADE can then be written as

___(3b)

which can be rearranged to get

where

Hereafter this scheme will be referred to as the ITCS scheme.

For numerical stability of implicit time-stepping schemes to be ensured it is generally found that the coefficients, and need always to be negative whereas needs to be positive. Given that , , and should always be positive, it follows that there remains still one stability criterion:

To evaluate and compare the ITCS scheme to the ETCS scheme and the analytical solution we need to make the following modifications to ADEStudy.m.

First we need to add a term called Scheme as an input argument to the FDSol subfunction. This can be done by modifying the start of the subfunction to read as follows.

 function u=FDSol(xI,tI,a,b,uI,u0,Pe,Cr,Scheme)

The term, Scheme, can be set to ‘ITCS’ or ‘ETCS’. This in turn can be used to dictate which expressions should be used to calculate the , and coefficients. To enable this flexibility replace the code:

 %Determine A, B and C coefficients
 A=a*dt/2/dx+b*dt/dx^2;
 B=1-2*b*dt/dx^2;
 C=-a*dt/2/dx+b*dt/dx^2;

in the FDSol subfunction with:

 %Determine A, B and C coefficients
 switch Scheme
 case 'ETCS'
 %Explicit time-stepping central difference in space
 A=a*dt/2/dx+b*dt/dx^2;
 B=1-2*b*dt/dx^2;
 C=-a*dt/2/dx+b*dt/dx^2;
 case 'ITCS'
 %Implicit time-stepping central difference in space
 A=-a*dt/2/dx-b*dt/dx^2;
 B=1+2*b*dt/dx^2;
 C=a*dt/2/dx-b*dt/dx^2;
 end

The above modification allows the use of Scheme to determine whether the user requires the ETCS or the ITCS scheme. Based on that information, the switch function can select the appropiate set of expressions for the , and coefficients. Read about switch in the help documentation to learn more about switch and how it can be used.

Again it should be noted that the , and coefficients represent the values of the first, second and third diagonals of a tridiagonal Jacobian matrix, . But for the ITCS scheme, is defined by the formula

MATLAB’s left divide operator can be used to solve such a system of linear equations. For example, x=A\B solves the system of linear equations A*x=B.

To implement the above equation replace

 %Solve problem
 for n=1:Nt-1
 u(:,n+1)=J*u(:,n);
 end
with
 %Solve problem
 switch Scheme(1:2)
 case 'ET'
 %Explicit time-stepping
 for n=1:Nt-1
 u(:,n+1)=J*u(:,n);
 end
 case 'IT'
 %Implicit time-stepping
 for n=1:Nt-1
 u(:,n+1)=J\u(:,n);
 end
 end

Note that we only need to use the first two characters of the string contained in Scheme, because the above code is only concerned with whether the scheme involves implicit time-stepping (i.e., ‘IT’) or explicit time-stepping (i.e., ‘ET’).

To compare the analytical solution with ETCS and ITCS schemes collectively, replace the following code within the main function

 %Evaluate the finite difference solution
 uFD=FDSol(x,t,a,b,uI,u0,Pe,Cr);

with

 %Evaluate the explicit time-stepping with central differences in space scheme
 uFD(:,1)=FDSol(x,t,a,b,uI,u0,Pe,Cr,'ETCS');
 %Evaluate the implicit time-stepping with central differences in space scheme
 uFD(:,2)=FDSol(x,t,a,b,uI,u0,Pe,Cr,'ITCS');

and update the legend appropriately.

Run the modified script-file and compare the results from the analytical solution and the finite difference solutions. Experiment by varying the specified value of and .

ITBS – Implicit time-stepping with backward differences

If we consider an Euler implicit time-stepping scheme and use a backward difference approximation for the term, ignoring all the truncation errors associated with the Taylor expansions, the ADE can then be written as

___(3c)

which can be rearranged to get

where

from which it can be seen that this scheme is unconditionally stable.

Hereafter this scheme will be referred to as the ITBS scheme.

ETBS – Explicit time-stepping with backward differences

If we consider an Euler explicit time-stepping scheme and use a backward difference approximation for the term, ignoring all the truncation errors associated with the Taylor expansions, the ADE can then be written as

___(3d)

which can be rearranged to get

where

from which it can be seen that this scheme is conditionally stable providing

Hereafter this scheme will be referred to as the ETBS scheme.

Classroom Assignment – Part 1

Modify the ADEStudy.m file further to also include and compare the ITBS and ETBS schemes. Generate a single figure overlaying the plots of concentration against distance using each of the different schemes. How do they compare with the results from the analytical solution?

Numerical diffusion

It is apparent that the stability criteria for the different schemes listed above do not necessarily to lead to accurate solutions. If you set and , you will see that all the finite difference solutions overestimate the amount of diffusion taking place, with the exception of the ETCS scheme.

It is possible to determine exactly how much the different schemes are overestimating the diffusion by consideration of the truncation errors in the Taylor expansions provided earlier.

The ETCS scheme

Consider again the ETCS scheme (i.e., Eq. (3a)). Substituting Eqs. (2a), (2b) and (2e) into Eq. (3a) whilst retaining the truncation errors leads to

___(4a)

Now consider Eq. (1), it follows that

After some further rearranging (and recalling that and are constants) it can be further shown that

___(4b)

which on substitution into Eq. (4a) leads to

from which it can be seen that the ETCS scheme gives rise to a numerical diffusion, , of magnitude

which is negative, hence explaining why ETCS appears to underestimate the amount of diffusion predicted by the analytical solution.

The ITCS scheme

Consider again the ITCS scheme (i.e., Eq. (3b)). Substituting Eqs. (2a), (2b) and (2d) into Eq. (3b) and then substituting Eq. (4b), whilst retaining the truncation errors, leads to

from which it can be seen that the ITCS scheme gives rise to a numerical diffusion, , of magnitude

The ITBS scheme

Consider again the ITBS scheme (i.e., Eq. (3c)). Substituting Eqs. (2a), (2c) and (2d) into Eq. (3c) and then substituting Eq. (4b), whilst retaining the truncation errors, leads to

from which it can be seen that the ITBS scheme gives rise to a numerical diffusion, , of magnitude

The ETBS scheme

Consider again the ETBS scheme (i.e., Eq. (3d)). Substituting Eqs. (2a), (2c) and (2e) into Eq. (3d) and then substituting Eq. (4b), whilst retaining the truncation errors, leads to

from which it can be seen that the ETBS scheme gives rise to a numerical diffusion, , of magnitude

Classroom Assignment – Part 2

1) Modify the ADEStudy.m file further to also evaluate the analytical solution with using the four different values of associated with the four finite difference schemes discussed above. Plot the results as dot markers over the finite difference results. How effective are the expressions for at predicting the diffusion associated with each of the finite difference schemes? Which finite difference scheme do you think is best and why?

2) Solve the above problem using ode45 for the time integration and try using central differences and then backward differences for the term. Compare your results with those from the other finite difference schemes. Also, derive expressions for the numerical diffusion associated with the two new schemes and check these work by evaluating the analytical solution with the additional diffusion terms (as in part 2 of the classroom assignment above).


Published with MATLAB® 7.13