Physics 6730 Assignment 1

For each of the following exercises create the specified files with your answer(s). Submit your homework using the course submit utility.

Exercise 1 Illustrating Gauss-Laguerre Quadrature

  P(x) = e-x(x5 - 4x4 + 5x3 - 2x2 + x - 1)

Here are the Gauss-Laguerre abscissas and weights for n=3

        x_i         c_i
       0.415774558  .7110930104
       2.294280360  .2785177333
       6.289945083  .0103892566

Please give your calculation, explanation, and results in a text file glag3.txt.

Hint: Remember that with Laguerre polynomials, the factor e^-x is part of the integration weight here, so it isn't included in the rhs expression: c₁ f(x₁) + c₂ f(x₂) + c₃ f(x₃).

Exercise 2 Understanding Gauss Quadrature

1. Divide the quintic polynomial coefficient in Exercise 1 above by the Laguerre polynomial (Maple) expand(LaguerreL(3,x)) to get quotient and remainder polynomials. (Hint: the Maple quo and rem functions help here.) Express the result in the form f(x) = Q(x)L(3,x) + R(x), where Q(x) and R(x) are polynomials. What is the degree of Q(x)? What is the degree of R(x)?
2. Expand the quotient Q(x) in terms of Laguerre polynomials.
3. Using the orthogonality of Laguerre polynomials, show that (Maple notation)

      int(exp(-x)*f(x),x=0..infinity) = int(exp(-x)*R(x),x=0..infinity) 
   

4. Finally, show that the integral of e^-xR(x) is given exactly by three-point Gauss-Laguerre quadrature.

Exercise 3 Program of the week: numerical integration

The error function erf(x) is frequently used in statistics. It is the integral of a normalized Gaussian distribution. In Maple notation

  erf(x) = 2/sqrt(Pi) * int((exp(-t^2), t=0..x)

You should carry out the integral using two numerical methods:

1. Eight-interval composite Simpson's rule
2. Eight-point Gauss-Legendre quadrature

For this exercise you will submit the following files:

File	Purpose of procedure
erf.cc	Supervises the integration, reads parameters, writes results
simp.cc	Does composite Simpson's rule integration
Makefile	Rule for compiling the above.

The procedure for doing the 8-point Gauss-Legendre integration is provided in ~p6730/examples/a01/gleg8.cc. You may simply copy it. (There is also an 8-pt Gauss-Laguerre program in the same directory called glag8.f.)

All computations are to be done in double precision where possible.

1. For convenience in this assignment you should put two program units in your file erf.cc. a main program that supervises the calculation and a function subprogram normdist that calculates the integrand. The latter name is passed to the integration routines.

   So we can check your code, the main program should

   1. FIRST prompt for a value of x (upper limit of integration range)
   2. SECOND prompt for the number of Simpson's intervals n
   3. call for the n-interval Simpson's rule evaluation,
   4. call for the 8-point Gauss-Legendre integration,
   5. write the result of both evaluations. Write 8 digits after the decimal point. See cout_format to review controlling cout precision.
   Input is on standard input and consists of the value x and the value n.
2. Your function normdist (only a few lines long) This procedure supplies the integrand for the composite Simpson's rule and Gauss-Legendre subroutines. It should return the value

         2/sqrt(Pi) * exp(-t^2)
   

   Please don't forget to #include <cmath> to get the correct exp.

3. Your procedure simp

   You will use this procedure to evaluate the composite Simpson's rule integral from zero to x. However, it is supposed to be general enough to handle any integrand, any (even) number of intervals, and any end points, so you can use it for other applications.

   This procedure should be coded as a function-type subprogram, so that it can be called from the main program as follows:

      asimp = simp(n,a,b,normdist);
   

   The result of this call is the evaluation of the integral of the function normdist over the range a to b by means of an n-interval (n+1-point) Simpson's rule. (The first point is at x = a, the last point at x = b, with n equal intervals (b-a)/n between.)

   You should design your subprogram so that it will integrate any function that is provided to it. See the handout to see how this is done. That is the reason for setting up the call to simp so that the name of the function that evaluates the integrand is passed as a dummy argument.

4. The supplied program gleg8. This subprogram behaves like your composite Simpson's rule procedure, but the number of points is fixed (8). It can handle any integrand. The calling sequence is

      a = 0; b = x;
      agleg = gleg8(a,b,normdist);
   

5. Your Makefile

   Copy the program gleg8.cc into the directory where you keep your files for this assignment and include gleg8.o in the list of object files to be linked to make the executable. Remember that make is smart enough to know that if gleg8.o doesn't exist, but gleg8.cc exists, it can be generated with g++, and so it will compile it automatically, without your having to say so explicitly.

Exercise 4 Running your program and interpreting your results.

(a) Run your program, using n = 8 and x = 2. Compare the results of the two numerical methods with the "exact" result from Maple. Be sure to use at least 10 digit accuracy in the Maple evalf.

(b) Estimate the error for your Simpson's rule result. Use the Richardson method to do that.

Do this by repeating the calculation with n = 4. We then have two equations in two unknowns:

     I = S(4) + k h^4       = S(4) + 16*e
     I = S(8) + k (h/2)^4   = S(8) + e

(c) Does the error estimate in (b) give an reasonable accounting (order of magnitude) of the difference between I and the Maple answer?

Please put your calculations and results in a file erf.txt.