====================================
Homework 3
====================================




Due extended to 23.59, Sep/30/2020 (Original due Sep/27/2020)
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Fortran and numerical integration (quadrature)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



Consider the integral

.. math::
   :nowrap:
   
   \begin{equation}
     I =  \int_{-1}^1 e^{\cos(k x)} dx,
   \end{equation}

with :math:`k = \pi` or :math:`\pi^2`. In this homework you will experiment with two different ways of computing approximate values of :math:`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 about the most basic quadrature rules in
class. For a more detailed description you can refer to `my lecture
notes`__. 

__ http://math.unm.edu/~motamed/Teaching/Fall20/HPSC/Quadrature_Motamed.pdf 

Part I: Trapezoidal rule
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

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

For a set of equi-distant grid points :math:`x_i = X_L + ih, \ \ i =
0,\ldots,n, \ \ h = \frac{X_R-X_L}{n}` on the interval :math:`[X_L, X_R]`, the composite trapezoidal rule reads:

.. math::
   :nowrap:
  
   \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 :math:`k = \pi` and :math:`k=\pi^2` and for :math:`n =
2,3,\ldots,N`, where you should pick :math:`N` so that the absolute
error is smaller than :math:`10^{-10}`. Note that here :math:`X_L =
-1` and :math:`X_R = 1`.

 * In your report, plot the error against :math:`n` using a logarithmic scale for both axes. 
 * 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 :math:`k = \pi` and
   why does it make the method converge faster than expected? For a
   careful study of the superior performance of the trapezoidal rule
   on special functions take a look at the `Euler-Maclaurin formula`__.
   


Part II: Gauss Quadrature
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

In Newton–Cotes quadrature rules, it is assumed that the value of the
integrand is known at equally spaced points. If it is possible to
change the points at which the integrand is evaluated, other methods
such as `Gauss quadrature`__ and `Clenshaw–Curtis quadrature`__  are often more suitable.


In Gauss quadrature the location of the grid-points :math:`z_i`
(usually referred to as quadrature nodes) and the quadrature weights :math:`\omega_i` are chosen so that the order of the approximation to the weighted integral   

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

is maximized, where the weight function :math:`w(z)` is assumed to be positive and
integrable (in this homework :math:`w(z) = 1`). Such "optimal" nodes turn out to be the roots of `orthogonal polynomials`__ associated with the weight function
:math:`w(z)`. For the case :math:`w(z) = 1`, we use `Legendre polynomials`__ and obtain the weights and nodes by a call to the subroutine ``lglnodes.f90`` (which is a Fortran version of Greg von Winckel's `matlab version`__). Use the subroutine ``lglnodes.f90`` from the course repository and compute the integral by the following code that implements the above formula: 
 
.. code-block:: fortran

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

Again:

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

A few notes
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 1. The errors should look something like the figure below (note that you should label the curves, I left them out on purpose).

 .. image:: quad_error.png
    :height: 600

 2. 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:

    .. code-block:: none

     $ gfortran -O3 -c gq_test.f90
     $ gfortran -O3 -c lglnodes.f90
     $ gfortran  -o gq_test.x  gq_test.o lglnodes.o 

 3. A better way of doing this is to use a `makefile`. I will talk
    about this in class. See also the `Makefile`__ section.
    
 4. It is probably a good idea to use allocatable arrays for ``x,w,f`` above

    .. code-block:: fortran
  
      allocate(x(0:n),f(0:n),w(0:n))
      ...Code here... 
      deallocate(x,f,w)

    as their size will change as you change :math:`n`. Don't forget to deallocate the arrays.

 5. Do not forget to update the top-level README file in your
    repository. You will also need to include a low-level README file
    in your HW3 subdirectory.

 6. Write up your findings neatly as a report and check it in to your
    repository. In general you should start your report with a brief description of the
    problem and continue with the procedure that you carry out and the
    methods that you use to solve the problem and obtain results. You
    then present your results (which can be in form of tables and/or
    figures, etc.) and discuss/analyze them. Finally you finish your
    report with your conclusions. Choose appropriate names for different sections of your
    report.
    

    
__ http://en.wikipedia.org/wiki/Newton-Cotes_formulas
__ http://en.wikipedia.org/wiki/Euler-Maclaurin_formula
__ http://en.wikipedia.org/wiki/Gaussian_quadrature
__ https://en.wikipedia.org/wiki/Clenshaw%E2%80%93Curtis_quadrature
__ http://people.math.sfu.ca/~cbm/aands/page_773.htm
__ https://en.wikipedia.org/wiki/Legendre_polynomials
__ http://www.mathworks.com/matlabcentral/fileexchange/4775-legende-gauss-lobatto-nodes-and-weights/content/lglnodes.m
__ http://math.unm.edu/~motamed/Teaching/Fall20/HPSC/makefiles.html




.. For this particular example we choose :math:`x_L = 0, \, x_R = 1, \, f(x) = \sin(2 \pi x)`.
.. http://amath.colorado.edu/faculty/fornberg/Docs/sirev_cl.pdf