MAT 420, Tuesday, Oct. 16

Assignment due Tuesday, Oct. 23.

Exercise 1. Write a subroutine that solves a triangular system of linear equations. Given an n×n upper triangular matrix T and a right-hand side b, your subroutine will compute the solution x. Call your routine as follows:

call trisolve(t,n,x,b,ierr)

where on entry:

t is an n×n double-precision matrix. (Your program will ignore the elements below the diagonal.)
n gives the size of the problem (n > 0)
b is the double-precision right-hand side to be solved for.

On return:

ierr is set to 0 unless an error occurs. (One possible error is that a diagonal element of T is exactly 0.) The contents of x are undefined if the routine is unsuccessful.
x is the solution of the equation Tx=b (when ierr=0).
The contents of the other parameters are unchanged.

Example: The system of equations

x + 2y + 3z = 4
5y + 6z = 7
8z = 9

corresponds to the upper triangular matrix

1 2 3
T = 5 6
8

and right-hand side

b = 4 7 9.

Triangular systems of equations are solved by back substitution. Start with the last equation, which gives z=9/8. Then back-substitute the value of z into the second equation and solve for y. Finally, use the computed values of y and z to find x.

Your job is to generalize this approach to an n-dimensional triangular system of equations. Store your code in a file called trisolve.f90. Include only the code for trisolve. (Your main program and other code to test the subroutine will go into another file that you will not hand in.)

Your subroutine should not print anything. Any errors should be noted by setting ierr to an appropriate value and returning. Use double precision (use module precision) for declarations. Here is a code skeleton to get you started.

Exercise 2. Write a manual page, trisolve.man, that documents your routine. For an example of what a manual page should look like, see the manual page for sgesv, which is one of the LAPACK routines.

Newton's method

Newton's method is a fast and efficient way to find the roots of smooth functions, but it is not guaranteed to converge. Given f and a starting guess for the solution, x₀, Newton's method generates the following sequence of values:

x₁ = x₀ - f(x₀)/ f '(x₀)
x₂ = x₁ - f(x₁)/ f '(x₁),
and in general,
x_n+1 = x_n - f(x_n)/ f '(x_n).

For example, to find the square root of 2 using only elementary arithmetic, we can define f(x) = x² -2, f'(x)=2x and take x₀=2. Then Newton's method gives the following sequence of values:

x₀ = 2
x₁ = 1.5
x₂ = 1.4166666666667
x₃ = 1.4142156862745
x₄ = 1.4142135623747
x₅ = 1.4142135623731
x₆ = 1.4142135623731

Notice that the answer is correct to 14 significant digits after only four iterations. (Most calculators use Newton's method to compute square roots.)

Exercise 3. Write a Fortran subroutine, newton(f,df,x,tol,ierr) that finds the root of a double-precision function f(x) given its derivative df(x). (Use free source form. Do not store all the x's. On entry, x is the initial guess of the root, and on successful return, x holds the estimate of the root. Stop the iteration after 100 steps, when two successive iterates are less than tol apart, or when the absolute value of any iterate exceeds, say, the square root of the largest representable value, whichever occurs first. (Suggestion: huge(x) is an intrinsic function that returns the largest representable floating-point number whose type is the same as the type of x.)

Set the integer argument ierr to 0 if the routine converges and to a nonzero value of your choosing otherwise.

Here's some code to get you started.

Your subroutine should not print anything. The communication of your subroutine with its caller is solely through the dummy arguments. Handle errors by setting ierr to an appropriate value and returning.

To test your program, you will need to create a separate file containing code for the functions and a main program. Here is a sample program that computes the square root of 2 using Newton's method.

You will turn in only the code for your newton subroutine in a file called roots.f90.

Exercise 4. You'll need to compile three files to test your implementation of Newton's method:

precision.f90, which contains the declarations of the kind variables that you will use elsewhere;
roots.f90, which contains your implementation of Newton's method;
myfuncs.f90 (say), which contains your test functions and a main program that calls Newton's method and prints the results.

Assuming that you have three files in a single directory, write a Makefile that compiles them in the correct order, assuming that you are using gfortran as the compiler.

Submission instructions

You will use tar to package and send your homework assignments, as follows:

tar cf hw6.tar trisolve.f90 trisolve.man roots.f90 Makefile

(Although you will need other files, as mentioned in Exercise 4, you will not turn in precision.f90 or myfuncs.f90. When I test your codes, I will supply versions of these files that will call your Newton routine.)

Mail your tar archive as an attachment to

mat420hw at gmail.com

with a subject line of the form

Your_name MAT 420 HW6

Notes on grading

Both newton and trisolve will be compiled and executed by an automated testing program. The grading scheme will be as follows:

If your routine fails to compile or link, then you will be awarded at most 1/3 of the possible points. Take care not to include code for main programs in the files that you turn in.
use precision (as described in last week's exercise) to declare the kind values for double-precision variables.
If your routine compiles but gives completely wrong results, then you will be awarded at most 1/2 of the possible points. Take care to arrange that the calling sequence of your subroutines are exactly as stated in the assignment.
You will lose points to the extent that your programs do not meet the listed specifications (e.g., the return values of ierr aren't consistent with the directions).