Finite differencing II
Return to the simplest equation, the diffusion equation
which we previously descretized as
and then solved to provide an explicit scheme that defines the future in terms of the present
Consider the following proposition
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 equations for unknowns.
For instance, rearranging the finite difference formulation,
and cleaning up the notation so that super-scripts refer to time and subscripts refer to space
which can be written for each point as
with the normal "issues" at the boundary.
This much more easily written if you know a bit about matrix vector multiplication. Define , and the matrix ,
Producing the very interesting situation where
Implementation in Fortran 90
First, you'll need to create the matrix . 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 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 contains 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 n=size(a) 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 .
- 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 (the diffusion constant) of 1.0, and is 1. Set a fixed boundary on both sides to zero. Use the descretization discussed above.
- 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?
- 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.
Other parts of finite differencing
The Lab Classes on Finite Differencing are split up into the following parts: