# Math 4083 (Numerical Analysis II) - Spring 2000

Instructor:

John Albert
Office: PHSC 1004, Ext. 5-3782.
Office hours: Mon and Thurs, 11:30-12:30; Wed 2:30-3:30 (or by appointment)
E-mail: jalbert@ou.edu

## Assignment 1 (Due Monday, Jan. 24)

• Problems 6.1 (p. 250): #3(a,d), 6, 7(c).
• Computer Problems 6.1 (p. 252): #1. Specifically, solve equation (8) on page 248 using bge.m, for n = 1, 2, 3, ... , and answer the following questions: (1) what is the first value of n for which the computed solution is not close to the correct solution? (2) For this value of n, find the condition number of A (use the "cond" function in Matlab). (3) For this value of n, use gepivot.m to compute the solution. Is it better than the solution found using bge.m?

## Assignment 2 (Due Monday, Jan. 31)

• Problems 6.2 (p. 266): #12, 14, 19.
• Computer Problems 6.2 (p. 270): #4, 7.

## Assignment 3 (Due Wednesday, Feb. 9)

• Problems 6.3 (p. 278): #2, 4.
• Computer Problem: write an M-file called penta.m, which implements the pseudocode on p. 278 for solving pentadiagonal systems. It would be a good idea to test your code on an example: try the one on p. 240.
• Problems 6.4 (p. 297-8): #1(b), 4.

## Assignment 4 (Due Wednesday, Feb. 16)

• Problems 6.5 (p. 315): #1. (Do by hand for a half an hour, or for how long you feel like doing it, whichever is longer.)
• Computer Problems 6.5 (p. 316): #2.

## Assignment 5 (due Fri. Mar. 3)

• Problems 13.1 (p. 529): #3
• (5 pts.) Suppose you use the explicit scheme for solving the initial-value problem for the heat equation, on the interval 0 <= x <= 1, for 0 <= t <= 10, with h = 10^(-4). How many points will your computational grid have to contain, in order to solve the problem with stability?
• Run heatexpl.m with initial data f(x) = sin(pi*x), taking (a) h = 0.1, k = 0.005125, M = 200; and (b) h = 0.1, k = 0.006, and M = 171. Print out graphs of the solutions, and explain what you see.
• Run the same two cases as above with heatimpl.m. Also try (c) h=0.1, k=0.1, and M=10. Compare the errors observed in these three runs.
• Carry out the same kind of stability analysis for the fully implicit scheme as we did in class for the explicit scheme. The correct conclusion should be that the scheme is stable for all values of h and k.
• Modify the code for heatimpl.m so that it performs the "Alternative Version of the Crank-Nicholson Method" described on pp. 527-528 of the text. Try this code on the three runs (a), (b), and (c) described above, and compare its performance to that of the fully implicit code.

## Assignment 6 (due Fri. Mar. 24)

• On pp. 540-541 of the text is the matrix equation Au=b which one must solve to find the solution of Laplace's equation on the square. Write down what the matrix A would look like if (a) you reordered the unknowns according to lexicographic ordering, so that u is the transpose of [u22, u23, u24, u32, u33, u34, u42, u43, u44] (b) you reordered the unknowns according to checkerboard ordering, so that u is the transpose of [u22, u42, u33, u24, u44, u32, u23, u43, u44].
• Modify the program laplacesq.m so that it solves the matrix equation Au=b by the SOR method rather than by the Gauss-Seidel method. Run the modified program on a test case and compare it with laplacesq.m.

## Assignment 7 (due Wed. Apr. 5)

• Try laxwend.m and wendrff.m, with initial data f(x)=hat(x), using (i) n=50, k=0.024, M=400; (ii) n=50, k=0.02, M=400; (iii) n=50, k=0.01, M=400. Explain what you see.
• Two finite-difference schemes for the equation u_t = cu_x are given by (i) v(x,t+k)=v(x,t-k)+(ck/h) (v(x+h,t)-v(x-h,t)) and (ii) v(x,t+k)=v(x,t)+(ck/2h) (v(x+h,t)-v(x-h,t)). Do von Neumann stability analyses for each of these schemes.

## Assignment 8 (due Fri. Apr. 21)

• Use FixedRK.m to solve the equation dx/dt=-5x with initial value x(0)=1, on the interval 0 <= x <= 5, starting at h=0.1 and increasing the value of h gradually. Plot the stable solutions on a single graph to show how their accuracy worsens as you approach the region of instability. Report h_0, the value of h at which instability first appears. Do this for each of the possible orders of the method; i.e., orders one through five. Thus you should turn in 5 plots and 5 values for h_0. (Remember, for Euler's order-one method, instability sets in at h=h_0=2/5; so for this method a good choice of values of h for your plot would be h=0.1, h=0.2, h=0.25, h=0.4.)
• Prove that when the classical 4th order Runge-Kutta method is applied to the problem dx/dt = a x, the formula for advancing the solution will be:

x(t+h) - x(t) = [ah + (1/2) a^2 h^2 + (1/6) a^3 h^3 + (1/24) a^4 h^4] * x(t)

• Prove that the local truncation error in the preceding problem is O(h^5). (Hint: the exact solution of the problem is x(t)=C*exp(at), where C is a constant. Show that the exact value of x(t+h)-x(t) differs from the right-hand side of the preceding formula by a term of size O(h^5).)
• Use the method of undetermined coefficients to derive the fourth-order Adams-Bashforth formula:

x(t+h)=x(t) + (h/24)[55 f(x(t),t) -59 f(x(t-h),t-h) + 37 f(x(t-2h,t-2h)) -9 f(x(t-3h),t-3h)]

• Use the method of undetermined coefficients to derive the fourth-order Adams-Moulton formula:

x(t+h)=x(t) + (h/24)[9 f(x(t+h,t+h) + 19 f(x(t),t) - 5 f(x(t-h),t-h) + f(x(t-2h),t-2h)]

## Assignment 9 (due at final exam)

• Determine whether each of the following methods is stable, consistent, and/or convergent.

(i) x(t+h)-x(t-h) = 2 h f(x(t),t)

(ii) x(t+h)-x(t-h) = (h/3) [7 f(x(t),t) - 2 f(x(t-h),t-h) + f(x(t-2h),t-2h)]

(iii) x(t+h)-x(t) = (h/24)[9 f(x(t+h,t+h) + 19 f(x(t),t) - 5 f(x(t-h),t-h) + f(x(t-2h),t-2h)]

• Use FixedRK with k=4, h=0.01 to solve the initial-value problem

d(x_1)/dt = (x_1)^2 + x_2 + t

d(x_2)/dt = - x_1 + (x_2)^2 + t^2

with initial data x_1 = 0.43 and x_2 = -0.69 at t=(-1). Find the values of x_1 and x_2 at t=1.

• Find the region of absolute stability for the "Implicit Euler's method", given by

x(t+h)-x(t) = h f(x(t+h),t+h).

## M-files

Below are links to M-files which we have used in our class so far. Here is how to use them.

• `bge.m          Basic Gaussian elimination. `
• `gepivot.m      Gaussian elimination with scaled partial pivoting. `
• `tridiag.m      Solves a tridiagonal system. `
• `jacobi.m       Solves a linear system by Jacobi iteration. `
• `heatexpl.m     Solves the heat equation by explicit forward differencing. `
• `start.m        Subroutine for initializing heatexpl.m. `
• `hat.m          Defines an initial function which can be used in heatexpl.m. `
• `heatimpl.m     Solves the heat equation using a fully implicit scheme. `
• `tridiag2.m     A function which solves tridiagonal systems; for use in heatimpl.m `
• ```laplacesq.m    Solves Laplace's equation on the square with Dirichlet boundary data,
using a 5-point difference scheme. ```
• ```laxwend.m      Solves the wave equation for left-moving waves using a Lax-Wendroff
difference scheme. ```
• ```wendrff.m      Solves the wave equation for left-moving waves using Wendroff's
implicit method. ```
• ```FixedRK.m      Solves an ordinary differential equation using a Runge-Kutta
method with fixed step size. Uses the files RKStep.m and f1.m (see below).
Written by Charles Van Loan, available at
http://www.cs.cornell.edu/cv/Books/SCMV/index.htm ```
• `RKStep.m       Used in FixedRK.m (see above). `
• ```FixedAB.m      Solves an ordinary differential equation using an Adams-Bashforth
method with fixed step size. Uses the files ABStart.m, ABStep.m and f1.m
(see below). Written by Charles Van Loan.   ```
• `ABStart.m      Used in FixedAB.m (see above). `
• `ABStep.m       Used in FixedAB.m (see above). `
• ```f1.m           The above codes require you to input (in single quotes) the name of a function
which defines the right-hand side of the ODE you are solving.  You can use this
m-file to define that function. ```