Assignment 3

Due April 1, 1997

Class questions and responses

This assignment is very imposing to look at. It is actually much easier than it appears. Almost all of the really nasty parts are provided for you in the form of routines referenced later on. The major part is to implement the algorithm given in class. This can be done in a number of ways including ones completely unacceptable in an advanced numerical analysis course (which this is not), e.g., form the matrix Ah as a dense matrix.

Part 1

You will be implementing the Rayleigh-Ritz-Galerkin (RRG) method for a boundary value problem in one dimension:

-(p(x)ux)x + r(x)u = f(x) in (0,1);
u(0) = u(1) = 0.
The functions p(x) and r(x) satisfy the following criteria in (0,1):

p(x) > pmin > 0 and r(x) > 0
Additionally, you will have a symmetric, positive definite operator.

The variational problem is to find uh in S(k,x,z) such that

a(uh,vh) = f(vh) for all vh in S(k,x,z),
where the bilinear form is given by

a(u,v) = Int(0,1) p(x)Du(x)Dv(x) + r(x)u(x)v(x) dx
and the finear functional is given by

f(v) = Int(0,1) f(x)u(x) dx.

Using the algorithm given in class, construct

Ahc = k
for C0-piecewise linear basis functions (i.e., S(2,x,1)) and classical cubic splines (i.e., S(4,x,1)). Use a uniform mesh for x and numerical quadrature to evaluate the integrals. The solution u(x) is given by

u(x) = sin(pi*x).

In completing the assignment, you may use the following routines:

Routine Language Description
GQRULE (Fortran) Returns the G-L points in [0,1] and the weights.
SPDCMP (Fortran) Factors a symmetric, positive definite matrix stored in lower profile format.
SPRSLV (Fortran) Solves Ax=b after A has been factored into LDLt by SPDCMP.
BSPLVD (Fortran) Evaluates B-splines and derivatives.
BSPLVB (Fortran) Evaluates B-splines (used by BSPLVD).
BVALUE (Fortran) Evaluates first j derivatives of a B-spline.
INTERV (Fortran) Computes an interval for use by BVALUE.
Makefile (make) Constructs a small executable.
Main (Fortran or C) A skeleton of a main program.
a3.tar (tar) A tar file with many files. ONLY UNPACK THIS FILE ONCE SINCE IT DESTROYS YOUR PREVIOUS WORK.
You can unpack the tar file with the command

tar xvf a3.tar
Use gtar if you get weird error messages (use gtar anyway since it is a superior program). WARNING: tar is a very destructive program. It will overwrite your work of the same file name(s) without so much as a warning. Use it very carefully or copy your files to a backup directory before using it.

Should there not be a Fortran routine and you are using C, you can always call the Fortran routine from C using the name followed by an underscore. (On some machines, you will need to call the Fortran routine using its name in capital letters and no underscore.) Remember that Fortran always passes a pointer to a variable, not a value.

You can use the program f2c to translate any of these from Fortran to C (f2c is on the mscf cluster of machines and can also be found on netlib). You must follow the directions at the top of a translated file to get the right libraries and include file.

You should produce a modular piece of software. In particular, you should have a file with four functions in it: p, r, f, and u. These should take arguments of x and mpys, where x is a point and mpys is the cumulative multiply count.

There should be an assembly phase and a solution phase. These should be completely independent of each other. Do the integration using a (k-1)-point Gauss-Legendre quadrature rule in each interval (xi,xi+1). Solve the linear algebra problems using SPRCMP and SPRSLV.

Debugging hints

First, write the right hand side assembly phase, debugging on the simplest cases when f is identically one and then a linear function. Second, write the matrix assembly phase, debugging on the simplest cases when p, r, and f are identically zero or one. Third, combine with the solution phase, debugging the complete program on simple problems with u(x) already in the subspace: except for error due to the numerical quadrature, yor program should reproduce the classical solution.

There is an example in SPDCMP. Be sure to understand it since the mu values are correct, but tricky. Also, all of the routines keep adding to the mpys variable. You must set this to zero in your main program before calling any other routine.

Part 2

Each of you must choose a unique p(x) and r(x). As you pick p(x) and r(x), they will be added to the table below. In each case, print out the multiply count and the approximate pointwise and L2 errors measured using the 8-point Gauss-Legendre quadrature rule in each interval. Estimate the rate of convergence with respect to N for each (k,z) combination that you try (you only have to look at (2,1) and (4,1), but you can earn bonus points by looking at others). Interesting cases to run are N=4,8,16,32,64 (where N is the number of elements). Which methods seem to produce the smallest error for the lowest cost?

Here are the coefficients chosen by the class (in order of acceptance so far):

Whop(x)r(x) Code
Sue Donim 1 0 Fortran or C
Ben Childers 1 + x2 x Fortran
Sean Townsend 2 + x2 ex Fortran
Karen Walters 1-x2 (x+1)-1 Fortran
Jonathan Hu sin(x) xex Fortran
Chris Yung 1+ln(2+x) ln(2+x)
HaoHua Tu (1+sin(pi*x))*pi-2 1+2sin(pi*x)

Part 3 (Bonus points)

You can get more than the usual high score on this assignment.

Cheers,
Craig C. Douglas

This is http://www.ms.uky.edu/~ccd/ma625/a3.html.
Back to http://www.ms.uky.edu/~ccd/ma625.