# Numerical Treatment of Differential Equations

M 383G, Unique #56470 and CAM 386K, Unique #64325, Spring 2011

### Instructor

Todd Arbogast
Office: RLM 11.162, Phone: 471-0166
Office: ACE 5.334, Phone: 475-8628
E-Mail: arbogast@ices.utexas.edu
Office hours: M and W from 10:00-11:30 in RLM 11.162.
Also, unless the instructor has pressing business, he will be available to help students who find him in one of his offices.

### Teaching Assistant

Yoonsang Lee
E-Mail: ylee@math.utexas.edu

Syllabus

### Homework, Projects, and Exams

Finals week office hours: M 5/9 and Tu 5/10, 9-10:30 a.m.

Final Exam. Wednesday, May 11, 9-12:00 noon, in RLM 7.114. The exam is comprehensive, but emphasizes parabolic and hyperbolic equations. You may use a page of notes on this exam.

Notes for Hyperbolic PDEs. hyperbolicNotes.pdf

Problem Set 12. Not to be handed in.
1. For the linear hyperbolic equation, find the CFL condition so that the Lax-Wendroff scheme is stable.
2. Iserles Exercise 17.8 (μ is aΔt/h).
3. Consider the Riemann problem for Burgers equation ut + (f(u))x = 0, where f(u) = ½u2.
For ε>0, verify that the solution w of
wt + wwx = εwxx
is
w(x,t) = W(x-st), s=½(uR+uL)
W(y) = uR + ½(uL-uR)[1-tanh((uL-uR)y/4ε)],
and that as ε tends to 0+, w(x,t) tends to u(x,t).
4. Consider ut + (f(u))x = 0, where f(u) = u3.
(a) Give an example of a smooth initial condition u(x,0)=u0(x) so that a shock forms.
(b) Give an example of a smooth initial condition u(x,0)=u0(x) so that a rarefaction forms.
(c) Find the Rankin-Hugoniot shock speed for the Riemann problem with uL=0.8 and uR=0.2.

Problem Set 11. Due W 5/4.
1. For 0 < a* ≤ a ≤ a* < ∞, consider the parabolic equation
ut - ∇· a∇ u = f    x ∈ Ω, t > 0,
u = g    x ∈ ∂Ω, t > 0,
u = u0    x ∈ Ω, t = 0.
Suppose that we approximate it by using backward Euler in time and a Galerkin method in space. Derive a stability bound for the Galerkin approximation using the solution itself (i.e., u-g) as test function.
2. Consider the parabolic equation above with a=1 and g=0. Approximate using a Crank-Nicolson time approximation
( (un+1 - un)/Δt , v) + ½(∇ (un+1 + un) , ∇ v ) = ½((fn+1 + fn) , v)
Derive an error estimate for the Galerkin approximation.

Notes for Parabolic PDEs. parabolicNotes.pdf

Exam 2. Monday, April 25, 5-7:00 p.m., in RLM 10.176, on Part II: Approximation of Elliptic Partial Differential Equations.

The exam should take about 1.5 hours, but I will give you 2 hours. You may use a page of notes on this exam.

90 A
80 A-
70 B+
60 B
50 B-
40 C+

Problem Set 10. Due M 4/18.
1. Show that Qk,k on a rectangle E is affine equivalent to the reference rectangle R=[-1,1]d, and that there is a universal constant C such that the following estimate holds:
infp in Qk,k ||u-p||m,2,EC hk+1-m |u|k+1,2,E
where 0 ≤ m ≤ k and h is the largest side length of E.
2. Write a code to demonstrate the above result for k=1 and both m=0,1, when d=2. Use a fixed function u(x,y) on [0,1]2 and a rectangular grid of spacing h=1/N, as h goes to zero. Let p be the nodal interpolant at the four corners of each grid element. You can use a tensor product composite two-point Gauss approximate integration over each grid element to compute the integrals.

Problem Set 9. Due F 4/8.
1. Larsson & Thomee Problem 5.10.
2. Iserles Exercise 9.9 (Iserles describes the basis, he means Q2,2).
3. Iserles Exercise 9.10.
4. Let the element domain E be a tetrahedron in 3-D, the shape functions P be the polynomials of degree less than or equal to 1, and the nodal values be either the values at the corners or the values at the centers of the faces. Prove that either choice gives a finite element. Describe how the first can be joined to form a finite element space in H1, but the second cannot.

Problem Set 8. Due F 4/1.
1. Larsson & Thomee Problem 5.6 (p. 74).
2. Prove that the abstract Galerkin method (using a finite dimentional approximating space) has a unique solution by using only the coercivity assumption to show uniqueness.
3. Iserles Exercise 9.11 (p. 203).
4. Prove that the element described in Iserles Equation (9.33) on page 199 is unisolvent. That is, if the degrees of freedom vanish for a cubic polynomial, then that cubic is identically zero

Problem Set 7. Due F 3/25.
1. Larsson & Thomee Problem 5.1 (p. 73).

2. For c≥0 a constant, write a finite element code to solve the problem
- u'' + cu = f    in (0,1)
u(0) = uL   and   u(1) = uR
Use piecewise continuous linear functions on a uniform mesh. Let the user select c, uL, uR, and the number of subintervals. Assume that f(x) is a given function subroutine. Approximate integrals of f(x) times the basis functions using the trapezoidal rule
∫ab F(x)dx ≈ ½(b-a)[F(b)+F(a)]
(Note that the basis function will be zero at one side of the interval! Nevertheless, it works fine.)

3. Test the convergence order of your code with error measured as
(∫ (u'-uh')2)½
You can choose c and u(x), and then define uL, uR, and f(x) from them to get a problem with a known true solution. Does the code work for c=0 and for c large (maybe 10 or 100)?

4. Test the convergence order of your code with error measured in the L2 norm.

5. What happens if you try c<0?

Problem Set 6. Due F 3/11.
Using your code from Problem Set 5 (or a friend's, if your code did not work), in this problem set, you will investigate the eigenvalues and eigenvectors of the discrete approximation to the Laplace operator.

In all runs, solve the Poisson problem
- uxx - u_yy = f    in (0,1)2
u = 0     for x=0,1 or y=0,1

(i.e., g=0) using a uniform grid. Recall that you solve
AU = F
and that the eigenvalues are
λa,b = (4/h2){sin2(aπh/2) + sin2(bπh/2)}
for 1≤a≤m and 1≤b≤m, with corresponding eigenvectors
va,bk(i,j) = sin(a i πh)sin(b j πh)

1. Use m=9 and h=1/10. Set f so that F=λa,bva,b, and show that the solution gives back the eigenvectors for the 6 cases where a=2,5,8 and b=3,7.

2. Use the exact eigenvalues and eigenfunctions
λa,b = (a2 + b2)&pi2
va,b(x,y) = sin(a π x)sin(b π y)

to set f so that F=λa,bva,b, solve the problem, and report the errors ||u-va,b||. In this case, use h=1/10,1/20, 1/40, 1/80, so you can determine the convergence rate for the reduction of the error. Just use the 4 cases a=1,9 and b=1,9.

Problem Set 5. Due F 3/4.
1. On a uniform grid, find an approximation to f''(x0) using only points to the right (bigger than) x0. The approximation should be second order in the grid spacing scale h, and it should use as few points as possible.
2. Show that on a uniform grid x-1, x0, x1, we have

(au')' = (1/h2) ( a1/2 (u1 - u0) - a-1/2 (u0 - u-1) ) + O(h2)
where a1/2=a((x1 + x0)/2).
3. Write a code to implement the standard 5-point finite difference method for the problem
au - div(d grad u) = f     for 0 < x < Lx, 0 < y < Ly
u = g
on the boundary
using uniform rectangular grids. Assume that a(x,y) and d(x,y) are given functions, such that a ≥ 0 and d(x,y) ≥ d* > 0.

(a) Write the complete code. You do not need to write your own linear solver. For example, you may solve the linear system by using the lapack routines (see http://www.netlib.org). These routines solve various linear algebra problems. They are written in Fortran, but can be called by C or C++ routines with only a bit of work. The routine you call is called dgbsv.f (it solves a general banded system of linear equations AX=B). You link in the routines (if they are installed) by using the appropriate compiler flags -llapack and possibly -lblas).
(b) By testing your code on some known problems, verify that your code is O(h2) for uniform grids. This is easily accomplished by choosing a solution u(x,y) and coefficients a(x,y) and d(x,y). Then define the functions
f(x,y) = a(x,y) u(x,y) - div(d grad u) = a(x,y) u(x,y) - ( dx ux + d uxx + dy uy + d uyy )
g(x,y) = u(x,y)
Thus you need to provide functions for
u, ux, uxx, uy, uyy, a, d, dx, dy
The error should be
sqrt( sumi,j (u(xi,yj) - ui,j)2 / (nx * ny) )
where nx and ny are the number of points in each coordinate direction.

Exam 1. Thursday, February 24, 7-9:00 or Friday, February 25, 2-4:00 on Part I: Approximation of Ordinary Differential Equations.

The exam should take about 1.5 hours, but I will give you 2 hours. You may use a page of notes on this exam.

To accomodate everyone, you should take the exam during one of the two time slots:

1. Thursday, February 24, 7-9:00 in RLM 12.166
2. Friday, February 25, 2-4:00 in RLM 11.176
90 A
80 A-
70 B+
60 B
50 B-
40 C+

Problem Set 4. Due F 2/18.
1. Compute the domain of linear stability for the following RK2 scheme:

y_{n+1} = y_n + (h/2)[f(t_n,y_n) + f(t_{n+1}, y_n + hf(t_n,y_n))] .
2. Stability tests.
stability.C    newton.C    newton.h    backEuler.C    backEuler.h
Compile the code using the linux command
c++ -o stability stability.C newton.C backEuler.C
Run the program by typing "stability".
The code solves 2-D problems of the form
y' = f(t,y)
with the Euler and backward Euler methods.
(a) Test problem 2 solves the linear problem y' = Ay where the matrix A is

and so has the solution

If, initially,

(with exactly 6 3's in 0.333333), up to a small amount of rounding error, only the first component of the solution exists.
(i) Solve the problem for 0 <= t <= 1 and h=0.1. Describe your results and explain them in terms of stability.
(ii) Solve the problem for 0 <= t <= 10 and h=0.01. Describe your results and explain them in terms of stability.
(b) Test problem 3 solves the nonlinear problem

over the time interval 0 <= t <= 1.
(i) Solve the problem with h=0.05. Describe your results and explain them in terms of linearized stability.
(ii) Solve the problem with h=0.2. Describe your results and explain them in terms of linearized stability.

Problem Set 3. Due M 2/14.
1. Iserles Exercises 3.1, 3.3, 3.6, 3.8.
2. Write a code to implement the standard 4th order Runge-Kutta (RK4) algorithm, and verify that it is 4th order.
3. Write a code to implement the Runge-Kutta-Fehlberg (RK4/5) error controlling method:

yn+1 = yn + (25.0/216.0)K1 + (1408.0/2565.0)K3 + (2197.0/4104.0)K4 - (1.0/5.0)K5
Yn+1 = yn + (16.0/135.0)K1 + (6656.0/12,825.0)K3 + (28,561.0/56,430.0)K4 - (9.0/50.0)K5 + (2.0/55.0)K6
where
Z1 = y_n
K1 = h f(t_n , Z1)

Z2 = y_n + (1.0/4.0)K1
K2 = h f(t_n + h/4 , Z2)

Z3 = y_n + (3.0/32.0)K1 + (9.0/32.0)K2
K3 = h f(t_n + 3h/8 , Z3)

Z4 = y_n + (1932.0/2197.0)K1 - (7200.0/2197.0)K2 + (7296.0/2197.0)K3
K4 = h f(t_n + 12h/13 , Z4)

Z5 = y_n + (439.0/216.0)K1 - 8K2 + (3680.0/513.0)K3 - (845.0/4104.0)K4
K5 = h f(t_n + h , Z5)

Z6 = y_n - (8.0/27.0)K1 + 2K2 - (3544.0/2565.0)K3 + (1859.0/4104.0)K4 - (11.0/40.0)K5
K6 = h f(t_n + h/2 , Z6)
and the local truncation error is estimated as
LTE = | Yn+1 - yn+1 | / h
For time step control, if E <= tolerance, then accept the step. Otherwise, let
delta = min( 4, max( 0.1, 0.84*sqrt(sqrt(tolerance/E)) ) )
h = h*delta
and recompute the step.
The coefficients in a C++ header file: rkfCoefs.h

My multistep code: multistep.C and multistep.h
An example of nonconvergence: nonconvergence.C

Problem Set 2. Due F 2/4 (actually M 2/7)
1. Iserles Exercises 2.3, 2.6.
2. Construct the interpolating polynomial for the function f(x)=x+x1/2 at the points x=1, 4, and 9. Find a bound for the error on the interval [1,9].
3. Construct the Adams-Bashforth method of order 4.
4. Multistep methods.
(a) Construct an explicit multistep method with local truncation error O(h3) with

ρ(w) = w(w-1/3)(w-1) = w3 - 4/3 w2 + 1/3 w
(b) Write a code to implement your method (use Euler's method to get y1 and y2 using many small steps).
(c) Verify that your code gives an O(h3) solution. Since Euler is O(h), to get accurate y0 and y1, you will need to use about 1/h2 small steps for each call to Euler, which solves between time 0 and h, and then between time h and 2h (thus the Euler error is about h2 times h = h3).

Problem Set 1. Due F 1/28.
1. Iserles Exercises 1.1, 1.3, 1.5, 1.7
2. Euler's method for the system of equations

y'(t) = f(t,y(t)),  a < t ≤ b
y(0) = ya
where the system has size d.
(a) Write a C++ or Fortran or Matlab (or other computer) subroutine to implement Euler's method. Be sure to allow the function f to be given as a separate subroutine. Also let the user specify the initial condition ya, the times a and b, and either the time step size or the number of steps to take. You should also return the solution at the final time.
(b) Write a main program and a function for f to test your code. You will need to call the Euler subroutine several times to get a graph of the solution over time.
(c) Verify that Euler is O(h) order accurate by testing it with a couple of examples for which the solution is known. [Hint: Plot log(|error|) versus log(1/h), and verify that a line of slope -1 results.]