Homework 3, due 08.00, 17/9-2015

Fortran and numerical integration (quadrature)

Consider the integral

\[\begin{equation} I = \int_{-1}^1 e^{\cos(k x)} dx, \end{equation}\]

with \(k = \pi\) or \(\pi^2\). In this homework you will experiment with two different ways of computing approximate values of \(I\). After doing this homework you will have written a Fortran program from scratch, called a subroutine that someone else has written and gained some knowledge of the accuracy of some different methods for numerical integration (quadrature). In future homework we will see how these computations can be performed in parallel.

I will talk briefly about the most basic quadrature rules in class, for a more detailed description you can take a look at Olof Runborg’s notes.

Trapezoidal rule

Recall that for a grid \(x_i = X_L + ih, \ \ i = 0,\ldots,n, \ \ h = \frac{X_R-X_L}{n}\), the composite trapezoidal rule is:

\[\begin{equation} \int_{X_L}^{X_R} f(x) dx \approx h\left(\frac{f(x_0)+f(x_n)}{2} + \sum_{i=1}^{n-1} f(x_i) \right). \end{equation}\]

Write a Fortran program that uses the composite trapezoidal rule to approximate the above integral for both \(k\) and for \(n = 2,3,\ldots,N\), where you should pick \(N\) so that the absolute error is smaller than \(10^{-10}\).

  • In your report, plot the error against \(n\) using a logarithmic scale for both axis.
  • Can you read off the order of the method from the slopes? How does this agree with theory?
  • What is special with the integrand in the case \(k = \pi\) and why does it make the method converge faster then expected? Hint: take a look at the Euler-Maclaurin formula.

The trapezoidal rule belongs to a class of quadrature called Newton-Cotes quadrature that approximates integrals using equidistant grids. The order of a composite Newton-Cotes method is typically \(s\) or math:s+1, where \(s\) is the number of points in each panel (2 for the Trapezoidal rule).

Gauss Quadrature

Another class of methods is Gauss quadrature. In Gauss quadrature the location of the grid-points (usually referred to as nodes) and weights, \(\omega_i\), are chosen so that the order of the approximation to the weighted integral

\[\begin{equation} \int_{-1}^{1} f(z) w(z) dz \approx \sum_{i=0}^n \omega_i \, f(z_i), \end{equation}\]

is maximized. Here the weight function \(w(z)\) is positive and integrable (for example \(w(z) = 1\).) For a given \(w(z)\) the nodes, \(z_i\) are the zeros of the polynomial \(\tau_n = (z-z_0)(z-z_1)\cdots(z-z_n)\) satisfying

\[\begin{equation} \int_{-1}^{1} \tau_n(z) q(z) w(z) dz = 0, \end{equation}\]

for all polynomials \(q(z)\) of degree less than \(n\). Or, equivalently, the nodes are the zeros to the degree \(n\) orthogonal polynomial associated with the weight function \(w(z)\).

Don’t worry if this sounds complicated, we will only consider the case \(w(z) = 1\) and obtain the weights and nodes by a call to the subroutine lglnodes.f90 (which is a f90 version of Greg von Winckel’s matlab version.) Use the subroutine lglnodes.f90 from the repository and compute the integral by the code below (which implements the above formula):

call lglnodes(x,w,n)
f = exp(cos(pi*pi*x))
Integral_value = sum(f*w)

Again:

  • plot the error against \(n\) using a logarithmic scale for both axis (perhaps in the same figure.)
  • For Gauss quadrature the error is expected to decrease as \(e(n) \sim C^{-\alpha n}\). Try some different \(C\) and \(\alpha\) to see if you can fit the computed error curves.

Notes

  1. The errors should look something like the figure below (note that you should label the curves, I left them out on purpose.)
_images/quad_error.png
  1. For this homework you will call a routine from another file which means that to build your executable you will have to compile two object files and then link them together:
$ gfortran -O3 -c gq_test.f90
$ gfortran -O3 -c lglnodes.f90
$ gfortran  -o gq_test.x  gq_test.o lglnodes.o
  1. A better way of doing this is to use a makefile as was discussed in class.

  2. It is probably a good idea to use allocatable arrays for x,w,f above

    allocate(x(0:n),f(0:n),w(0:n))
    ...Code here...
    deallocate(x,f,w)
    

    as their size will change as you change \(n\). Don’t forget to deallocate the arrays.