Finite difference examples in Fortran¶

This section discusses implementations of the finite difference approximations discussed in the section Finite difference approximations. The examples can be found in the files differentiate_vX.f90 in the repository. Here we consider an example with a first order derivative.

First we approximate the derivative using a three point second order accurate stencil on an equidistant grid. Then we use Bengt Fornberg’s subroutine weights1 to find the approximation which uses all the function values at all the grid-points to achieve an approximation of maximal order. We will find that the this second approach does not work well on an equidistant grid. In the third example we find that the maximal order of accuracy finite difference stencil (usually called spectral finite differences) gives excellent! results on graded Legendre-Gauss-Lobatto grids. These grids are also discussed in Homework 3, due 23.59, 23/9-2016 in the context of numerical integration.

Example 1, a low order finite difference method¶

To this end consider the function $$f(x) = \exp(\cos \pi x)$$ and its derivative $$f'(x) = -\pi \sin (\pi x) \exp(\cos \pi x)$$ on the domain $$x\in[-1,1]$$.

First we consider the approximation (derived in the same way as in the finite difference section)

\begin{eqnarray} (f')_i = \left\{ \begin{array}{ll} \frac{-1.5 f_{i}+2f_{i+1}-0.5 f_{i+2}}{h} & i = 0\\ \frac{0.5 f_{i+1}-f_{i-1}}{h} & i = 1,\ldots,n-1\\ \frac{1.5 f_{i}-2f_{i-1} + f_{i-2}}{h} & i = n \end{array} \right. \end{eqnarray}

on the equidistant grid

$$x_i = -1 + i h, \, i = 0,\ldots,n, \ \ h = \frac{2}{n}.$$

As we could find the exact derivative we can now try out our finite difference approximation for a sequence of smaller and smaller $$h$$ to see if the error decreases as we expect. As this means that we will have different number of grid-points for different $$h$$ we declare the array x to be a single dimension allocatable array (see line 6). The difference stencils on the other hand will be the same independent of the number of grid points and can thus be stored in a fixed size array diff_weights(1:3,1:3). Note that in Fortran the default indexation starts with one so an equivalent declaration would be diff_weights(3,3). We also declare the arrays holding the function, the approximation to the derivative and the exact derivative as allocatable.

The code then loops over all $$n=4,5,\ldots,400$$ and allocates arrays of size $$n+1$$ (x(0:n)) (line 10), computes the grid spacing, the grid, the function and the exact derivative (lines 12-18). When the function is known on the grid we set up the array that holds the finite difference weights (note that we bake in the $$h$$ dependence and thus have to change it for every loop) and compute the finite difference approximation to $$f'(x)$$ (lines 23-45). Finally, on line 46 we compute the maximum error on the grid and print it to screen.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 program differentiate implicit none real(kind = 8), parameter :: pi = acos(-1.d0) integer :: n,i,j real(kind = 8) :: h,diff_weights(1:3,1:3) real(kind = 8), dimension(:), allocatable :: x,f,df,df_exact do n = 4,100 ! Allocate memory for the various arrays allocate(x(0:n),f(0:n),df(0:n),df_exact(0:n)) ! Set up the grid. h = 2.d0/dble(n) do i = 0,n x(i) = -1.d0+dble(i)*h end do ! The function and the exact derivative f = exp(cos(pi*x)) df_exact = -pi*sin(pi*x)*exp(cos(pi*x)) ! Set up weights for a finite difference stencil ! using three gridpoints ! In the interior we use a centered stencil diff_weights(1:3,1) = (/-0.5d0,0.d0,0.5d0/) ! To the left we use a biased stencil diff_weights(1:3,2) = (/-1.5d0,2.d0,-0.5d0/) ! To the right we use a biased stencil too diff_weights(1:3,3) = (/1.5d0,-2.d0,0.5d0/) ! scale by 1/h diff_weights = diff_weights/h ! Now differentiate ! To the left is a special case df(0) = 0.d0 do j = 1,3 df(0) = df(0) + diff_weights(j,2)*f(j-1) end do ! Now interior points do i = 1,n-1 df(i) = sum(diff_weights(1:3,1)*f((i-1):(i+1))) end do ! Finally, special case to the right df(n) = 0.d0 do j = 1,3 df(n) = df(n) + diff_weights(j,3)*f(n-(j-1)) end do write(*,'(I3,2(ES12.4))') n,h,maxval(abs(df-df_exact)) ! Deallocate the arrays deallocate(x,f,df,df_exact) end do end program differentiate 

The results from the output of the above program is displayed the figure below. As can be seen the slope of the error is the same as a $$h^{2} \sim n^{-2}$$ indicating that the implementation may be correct.

Example 2, maximum order on an equidistant grid¶

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 ! This program differentiates functions ! Using a standard second order finite difference stencil ! and a difference of maximal width computed by ! Bengt Fornberg's weights subroutine. program differentiate implicit none real(kind = 8), parameter :: pi = acos(-1.d0) integer :: n,i,j real(kind = 8) :: h,diff_weights(1:3,1:3) real(kind = 8), dimension(:), allocatable :: x,f,df,df_exact ! Variables for weights integer :: nd,m real(kind = 8) :: z real(kind = 8), dimension(:,:), allocatable :: c m = 1 do n = 4,100 ! Allocate memory for the various arrays allocate(x(0:n),f(0:n),df(0:n),df_exact(0:n)) allocate(c(0:n,0:m)) ! Set up the grid. h = 2.d0/dble(n) do i = 0,n x(i) = -1.d0+dble(i)*h end do ! The function and the exact derivative f = exp(cos(pi*x)) df_exact = -pi*sin(pi*x)*exp(cos(pi*x)) !!!! SKIPPING SOME REDUNDANT CODE HERE (SEE ABOVE)!!!! ! Finally, special case to the right df(n) = 0.d0 do j = 1,3 df(n) = df(n) + diff_weights(j,3)*f(n-(j-1)) end do write(*,'(I3,2(ES12.4))',advance='no') n,h,maxval(abs(df-df_exact)) ! Now find the difference stencil that uses all the points nd = n do i = 0,n df(i) = 0.d0 z = x(i) call weights1(z,x,n,nd,m,c) do j = 0,nd df(i) = df(i) + c(j,1)*f(j) end do end do write(*,'(ES12.4)') maxval(abs(df-df_exact)) ! Deallocate the arrays deallocate(x,f,df,df_exact,c) end do end program differentiate 

Example 3, maximum order on a Legendre-Gauss-Lobatto grid¶

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 ! ! This program differentiates functions ! Using a standard second order finite difference stencil ! and a difference of maximal width computed by ! Bengt Fornberg's weights subroutine. ! This third version adds differentiation on a ! Legendre-Gauss-Lobatto grid. program differentiate implicit none !!!! SKIPPING SOME REDUNDANT CODE HERE (SEE ABOVE)!!!! do n = 4,100 ! Allocate memory for the various arrays allocate(x(0:n),f(0:n),df(0:n),df_exact(0:n),w(0:n)) !!!! SKIPPING SOME REDUNDANT CODE HERE (SEE ABOVE)!!!! ! The above approach is plagued by numerical ! roundoff due to ill-conditioning. ! A much better grid distribution is the LGL nodes. call lglnodes(x,w,n) f = exp(cos(pi*x)) df_exact = -pi*sin(pi*x)*exp(cos(pi*x)) ! Now find the difference stencil that uses all the points nd = n do i = 0,n df(i) = 0.d0 z = x(i) call weights1(z,x,n,nd,m,c) do j = 0,nd df(i) = df(i) + c(j,1)*f(j) end do end do write(*,'(ES12.4)') maxval(abs(df-df_exact)) ! Deallocate the arrays deallocate(x,f,df,df_exact,c,w) end do end program differentiate