Numerical Treatment of Differential Equations

Instructor

Teaching Assistant

General Information

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 u_t + (f(u))_x = 0, where f(u) = ½u².
For ε>0, verify that the solution w of
   w_t + ww_x = εw_xx
is
   w(x,t) = W(x-st), s=½(u_R+u_L)
   W(y) = u_R + ½(u_L-u_R)[1-tanh((u_L-u_R)y/4ε)],
and that as ε tends to 0⁺, w(x,t) tends to u(x,t).
4. Consider u_t + (f(u))_x = 0, where f(u) = u³.
   (a) Give an example of a smooth initial condition u(x,0)=u₀(x) so that a shock forms.
   (b) Give an example of a smooth initial condition u(x,0)=u₀(x) so that a rarefaction forms.
   (c) Find the Rankin-Hugoniot shock speed for the Riemann problem with u_L=0.8 and u_R=0.2.

Problem Set 11. Due W 5/4.
1. For 0 < a_* ≤ a ≤ a^* < ∞, consider the parabolic equation
   u_t - ∇· a∇ u = f    x ∈ Ω, t > 0,
   u = g    x ∈ ∂Ω, t > 0,
   u = u₀    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
   ( (uⁿ⁺¹ - uⁿ)/Δt , v) + ½(∇ (uⁿ⁺¹ + uⁿ) , ∇ v ) = ½((fⁿ⁺¹ + fⁿ) , 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.

Grading scale:
   90 A
   80 A-
   70 B+
   60 B
   50 B-
   40 C+

Problem Set 10. Due M 4/18.
1. Show that Q_k,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:
    inf_{p in Qk,k} ||u-p||_m,2,E ≤ C h^k+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]² 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 Q_2,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 H¹, 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
    ∫_a^b 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'-u_h')²)^½
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 L² 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
   - u_xx - u__yy = f    in (0,1)²
   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/h²){sin²(aπh/2) + sin²(bπh/2)}
for 1≤a≤m and 1≤b≤m, with corresponding eigenvectors
   v^a,b_k(i,j) = sin(a i πh)sin(b j πh)

1. Use m=9 and h=1/10. Set f so that F=λ^a,bv^a,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 = (a² + b²)&pi²
   v_a,b(x,y) = sin(a π x)sin(b π y)
to set f so that F=λ_a,bv_a,b, solve the problem, and report the errors ||u-v_a,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''(x₀) using only points to the right (bigger than) x₀. 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, x₀, x₁, we have

(au')' = (1/h²) ( a_1/2 (u₁ - u₀) - a_-1/2 (u₀ - u_-1) ) + O(h²)

where a_1/2=a((x₁ + x₀)/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 < L_x, 0 < y < L_y
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(h²) 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) - ( d_x u_x + d u_xx + d_y u_y + d u_yy )
g(x,y) = u(x,y)

Thus you need to provide functions for

u, u_x, u_xx, u_y, u_yy, a, d, d_x, d_y

The error should be

sqrt( sum_i,j (u(x_i,y_j) - u_i,j)² / (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

Grading scale:
   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.
Set-up: Download the code
    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:

y_n+1 = y_n + (25.0/216.0)K1 + (1408.0/2565.0)K3 + (2197.0/4104.0)K4 - (1.0/5.0)K5
Y_n+1 = y_n + (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 = | Y_n+1 - y_n+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+x^1/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(h³) with

ρ(w) = w(w-1/3)(w-1) = w³ - 4/3 w² + 1/3 w

(b) Write a code to implement your method (use Euler's method to get y₁ and y₂ using many small steps).
(c) Verify that your code gives an O(h³) solution. Since Euler is O(h), to get accurate y₀ and y₁, you will need to use about 1/h² 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 h² times h = h³).

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) = y_a

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 y_a, 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.]
(d) Test your code on

y'(t) = -1/y,  0 < t
y(0) = 1

for which the solution is

y(t) = (1 - 2t)^1/2

What goes wrong?

My code: euler.C and euler.h

---