Finite differencing II

From Interactive System for Ice sheet Simulation
Jump to: navigation, search



Return to the simplest equation, the diffusion equation

 \frac{\partial u(x,t) }{\partial t} = \frac{\partial ^2 u(x,t)}{\partial x^2},

which we previously descretized as

 \frac{u(x,t+\Delta t) - u(x,t)}{\Delta t} = \frac{u(x+\Delta x,t) - 2u(x,t) + u(x-\Delta x,t)}{\Delta x^2}.

and then solved to provide an explicit scheme that defines the future in terms of the present

u(x,t + \Delta t) = u(x,t) + \Delta t \frac{u(x+\Delta x,t) - 2u(x,t) +
u(x-\Delta x,t)}{\Delta x^2}

Consider the following proposition

 \frac{u(x,t+\Delta t) - u(x,t)}{\Delta t} = \frac{u(x+\Delta x,t+\Delta t) - 2u(x,t +\Delta t) + u(x-\Delta x,t +\Delta t)}{\Delta x^2},

demanding that you find the future in terms of the future! Is this possible?

A System of Linear Equations

Sure it can be solved! One only needs to observe that the each degree of freedom corresponds to a linear equation. None of the equations can be solved by itself, but the system of equations, one for each degree of freedom (node) in the problem, can be solved because there are N equations for N unknowns.

For instance, rearranging the finite difference formulation,

u(x,t+\Delta t)  - 
\Delta t \frac{u(x-\Delta x,t+\Delta t) - 2u(x,t + \Delta t) + u(x+\Delta x,t + \Delta t)}{\Delta x^2}=u(x,t)

and cleaning up the notation so that super-scripts refer to time and subscripts refer to space

u_j^{n+1} - \frac{\Delta t}{\Delta x ^2} \left( u_{j-1}^{n+1} - 2 u_j^{n+1} + u_{j+1}^{n+1} \right) = u_j^n.

which can be written for each point j as

u_1^{n+1} - \frac{\Delta t}{\Delta x ^2} \left( u_{-1}^{n+1} - 2 u_1^{n+1} + u_{2}^{n+1} \right) = u_1^n,   , 
u_2^{n+1} - \frac{\Delta t}{\Delta x ^2} \left( u_{1}^{n+1} - 2 u_2^{n+1} + u_{3}^{n+1} \right) = u_2^n,   , 

with the normal "issues" at the boundary.

This much more easily written if you know a bit about matrix vector multiplication. Define  \alpha = \frac{\Delta t}{\Delta x^2}, and the matrix \mathbf{A},

\mathbf{A} = \begin{bmatrix}
1 & 0 & 0 & 0 && \cdots & 0\\
-\alpha & 1+2\alpha & -\alpha && \cdots & &0 \\
0&-\alpha & 1+2\alpha & -\alpha & \cdots && 0 \\
 & \ddots & \ddots &    \ddots  &&        & \\
\vdots & \cdots & -\alpha & 1+2\alpha & -\alpha  & \cdots & \vdots \\
 &&& \ddots & \ddots &    \ddots & \\

0 & \cdots && 0 & -\alpha & 1+2\alpha & -\alpha \\
0 & \cdots && 0 & 0 &0 & 1 \\
\end{bmatrix}   , 

Producing the very interesting situation where

\mathbf{A} \mathbf{u^{n+1}} = \mathbf{u^{n}}
The Stencil for the most common explicit method for the parabolic equation.
The implicit method stencil.

Implementation in Fortran 90

Sparse Matrices

First, you'll need to create the matrix \mathbf{A}. We note that this matrix is `sparse' or mostly zeros. This is an important properties as it greatly reduces the storage requirements and can lead to efficient implementations of the algorithms which solve linear systems.

At this point we will simply observe that the matrix is tridiagonal, meaning non-zeros entries occur along the diagonal, the diagonal above the diagonal, and the diagonal below the diagonal. This is easy to see by inspection of the matrix for \mathbf{A} written above.

We will return to some of the subtleties of sparse matrices next week. For now, let's just look at how to get the job done in Fortran. You'll be making use of a subroutine which solves linear systems with a tridiagonal matrix. The code for that follows

module solvers
      subroutine tridiag(a,b,c,x,y)
        !*FD Tridiagonal solver. All input/output arrays should have the 
        !*FD same number of elements.
        real,dimension(:),intent(in) ::  a !*FD Lower diagonal; a(1) is ignored.
        real,dimension(:),intent(in) ::  b !*FD Centre diagonal
        real,dimension(:),intent(in) ::  c !*FD Upper diagonal; c(n) is ignored.
        real,dimension(:),intent(out) :: x !*FD Unknown vector
        real,dimension(:),intent(in) ::  y !*FD Right-hand side 
        real,dimension(size(a)) :: aa
        real,dimension(size(a)) :: bb
        integer :: n,i
        aa(1) = c(1)/b(1)
        bb(1) = y(1)/b(1)
        do i=2,n
           aa(i) = c(i)/(b(i)-a(i)*aa(i-1))
           bb(i) = (y(i)-a(i)*bb(i-1))/(b(i)-a(i)*aa(i-1))
        end do 
        x(n) = bb(n)
        do i=n-1,1,-1
           x(i) = bb(i)-aa(i)*x(i+1)
        end do
  end subroutine tridiag
end module solvers

I'll leave it to you to use this idea to create the full \mathbf{A}.


  1. Solve the one dimensional diffusion problem on a lattice of 500 sites, with an initial condition of zeros at all sites, except the central site which starts at the value 1.0. Use k(x) (the diffusion constant) of 1.0, and \Delta x is 1. Set a fixed boundary on both sides to zero. Use the descretization discussed above.
  2. Repeat the previous exercise using the Crank Nicolson method, which averages in the time domain. The stencil for this method appears below. Does this allow for a larger time step?  
 \frac{u_j^{n+1} - u_j^{n}}{\Delta t} = \frac{1}{2} \left(\frac{u_{j+1}^{n+1} - 2u_j^{n+1} + u_{j-1}^{n+1}}{\Delta x^2}+\frac{u_{j+1}^{n} - 2u_j^{n} + u_{j-1}^{n}}{\Delta x^2}\right).\,  
  • Vary the grid spacing and time step in both the models above. Examine the height of the central peak. Does it vary? Does this mean the simulation is not working? Consider the integral under the curve.
The Crank-Nicolson implicit method stencil.

Other parts of finite differencing

The Lab Classes on Finite Differencing are split up into the following parts: