MATH567-Homework 3 Solved

1.   Linear two-point boundary value problems: Finite differences

(a)    Implement a general function in whatever language you like for numerically solving the general two-point boundary value problem (TPBVP)

u00 = p(x)u0 + q(x)u + r(x), x ∈ [0,1], u(0) = α, u(1) = β

using centered, second-order finite differences. Your routine should be called fd2tpbvp, and should be entirely general (i.e. it should take as input functions representing p(x), q(x), & r(x), values for α & β, and the number of interior discretization points m). Send me an electronic copy of the code and include a listing of the code in your write-up.

(b)   Use your function from part (a) to numerically approximate the TPBVP

u00 = 2tan(x)u0 + 2, u(0) = u(1) = 0,

which has the exact solution u(x) = (x − 1)tan(x). Compare the approximate solution to the exact solution using m = 10,20,40,80,160 equally spaced interior points over [0,1] (i.e. xj = jh, j = 1,...,m, with h = 1/(m + 1)) and verify numerically that the approximate solution is second-order accurate.

2.   Fictitious point method for Robin boundary conditions. Consider the TPBVP

u00 = p(x)u0 + q(x)u + r(x), x ∈ [a,b],

with mixed boundary conditions

u(a) = α and β1u(b) + β2u0(b) = β3.

Suppose we discretize the equation with m+1 equally-spaced subintervals of spacing h. Using the fictitious point method described in in class, derive the following second-order accurate finite difference approximation for um+1:

 ,

where u(b) ≈ um+1, u(b − h) ≈ um, pm+1 = p(b), qm+1 = q(b), and rm+1 = r(b).

3.   Neumann-Neumman boundary conditions. Consider the 1D Poisson equation with NeumannNeumann boundary conditions

                                                       u00(x) = f(x), a ≤ x ≤ b, u0(a) = σ0, u0(b) = σ1.                                     (1)

As discussed in class this equation has an infinite number of solutions if 

(the so-called compatibility condition), otherwise there is no solution. If σ0 = σ1 = 0 (i.e. the zero-flux condition) then it can be shown that solutions to the Poisson equation are unique up to an additive constant. If we discretize this BVP at m + 2 equally-spaced points across [a,b] and use second-order accurate FD formulas in the interior and the second-order accurate fictitious point method at the boundary, we arrive at the linear system

 

b

where h = (b − a)/(m + 1), fj = f(xj), and xj = a + jh, j = 0,...,m + 1.

As discussed in class (and by simple inspection), the matrix A in this equation is singular since every row sums to zero. This means that the vector of length m + 2

                                                                           e = [1     1     ...     1]T

is an eigenvector corresponding to a zero eigenvalue: Ae = 0. Furthermore, it is the only eigenvector (up to a non-zero multiplicative constant) with this property, i.e. it is the only vector in the null-space of A. Another property of A is that the vector of length m + 2

w  

is an eigenvector of AT corresponding to a zero eigenvalue: ATw = 0 or wTA = 0T. It is also the only eigenvector (up to a non-zero multiplicative constant) with this property, i.e. it is the only vector in the null-space of AT.

A result from linear algebra (known as the Fredholm Alternative) is that the system Au = b has an infinite number of solutions if wTb = 0 (which guarantees b is in the column space of A). This condition can be viewed as a discrete version of the continuous compatibility condition .

This is all beautiful theory, but in practice what we want to do is actually compute one of the solutions to the BVP. This homework problem and the next one discuss two approaches (with this problem being the most general of the two) for computing a solution with a “direct” method. It should be noted that effective “iterative” methods can also be used. The overall goal of any of these methods is to make the problem unique by adding on an additional constraint. In the continuous problem, especially for the zero-flux case, this is typically accomplished by specifying the “total amount” of u in the domain, which mathematically we write as

                                                                                                                                  (3)

For the zero-flux case, this fixes the arbitrary constant in the solution to be equal to U. This same approach is followed in the discrete problem as well, with the discrete version of this integral taking the form

 

where we have used the fact that h = (b − a)/(m + 1). Note that this is just the Trapezoidal rule approximation to (3).

The discrete problem (2) can now be converted into the following equality-constrained quadratic programming problem:

min J(u) = 1 uTAu − uTb

                                                                                           2                                                                          (4)

subject to wTu = (m + 1)U

The method of Lagrange multipliers can be used to solve this problem. The idea is to form the Lagrangian

 ,

and find where the gradient is zero (in general one looks for the saddle points of L). Here the gradient is defined as

  .

Applying this to the Lagrangian gives

∇L(u,λ) = Au − b + λw = 0,

which is a linear system of m+2 equations and m+3 unknowns, u and λ. Using the constraint wTu = (m + 1)U gives one extra equation and leads to the (m + 3)-by-(m + 3) linear system of equations

                                                                     .                                                 (5)

This is an example of a more general type of linear system called a “saddle point system”.

(a)    Show that the linear system (5) has a unique solution regardless of b.

Hint: One way to proceed is to show that the only solution to (i) Au + λw = 0 and (ii) wTu = 0 (i.e. the system (5) with the right hand side set to zero) is u = 0 and λ = 0. One way to show this is to consider what happens when multiplying the first equation by wT. This should give you the result that Au = 0, which means u = αe for some α. You can then use the equation wTu = 0 to show that α = 0.

(b)   Show that if wTb = 0 in (5) then λ = 0. This means that the solution u from (5) solves the original system (2) and thus gives an approximate solution to the BVP (6).

Note: If wTb 6= 0 then the solution u from (5) does not solve (2), however, u may still approximately solve the BVP (6). This will be reflected in the size of λ, since kAu − bk = |λ|kwk. If |λ| = O(h2) then u will solve the BVP up to O(h2) since this is the truncation error of the FD method.

(c)    Solve the BVP (6) numerically with a = 0, b = 2π, m = 99, f(x) = −4cos(2x), and σ0 = σ1 = 0 using the augmented linear system (5). Set U = 0 and check your answer against the true solution u(x) = cos(2x). Plot the error in the approximate solution and report the relative two-norm of the error. Also report the value of λ. What does this value tell you about your solution in regards to solving the system (2).

(d)   Repeat part (c) but now for the function f(x) = x, σ0 = −π2, and σ1 = π2. Plot the numerical solution you obtain and report the value of λ. Does it seem reasonable that this approximately solves (2)? Explain.

4. Neumann-Neumman boundary conditions and the Discrete Cosine Transform. Consider again the Neumann-Neumann Poisson equation, but now with strictly zero-flux boundary conditions

                                                         u00(x) = f(x), a ≤ x ≤ b, u0(a) = 0, u0(b) = 0.                                       (6)

In this problem, we will look at different way to solve the corresponding linear system (2) that results from discretizing this problem with a second-order accurate FD formula. Here we again let h = (b − a)/(m + 1), xj = a + jh, uj ≈ u(xj), and fj = f(xj) for j = 0,...,m + 1. In this approximation, the fictitious point idea has been applied, which used the assumption

                                   and      ,            (7)

where u−1 ≈ u(x0 − h) and um+2 ≈ u(xm+1 + h). The method that we will use is based on the Discrete Cosine Transform (DCT).

(a)   Consider the vectors u and f = b in (2) of length m+2. Let uˆ and ˆf denote the DCT of these vectors, respectively, i.e.

 ,

for k = 0,...,m + 1. Now, we can express the entries of u and b as

 ,

Substitute these expressions into the linear system (2) to show that row j of this systems simplifies to

  ,

in the case of j = 0 and j = m + 1 use (7).

(b)   Use the results from part (a) to show that the DCT of the solution to (2) for k = 1,...,m + 1 is given by

 .

However, for k = 0, the value of ˆuk is undefined. For this case, note that if fˆ0 = 0 (or

“numerically zero”) then ˆu0 can be chosen arbitrarily. Show how the condition fˆ0 = 0 corresponds to the discrete compatibility condition wTb = wTf = 0 discussed in the previous problem.

(c)    Put parts (a) and (b) together to explain how to obtain the solution to (2) when the right hand satisfies wTf = 0. Also explain how one makes the solution unique by fixing the arbitrary constant to U.

(d)   The Matlab codes dct.m and idct.m on the course webpage compute the DCT and inverse DCT (these codes use the FFT to do the computation and are thus ‘fast’, although not as fast as they could be). Use these codes (or the equivalent codes in another language) and the procedure outlined in (a)–(c) to solve the problem from 4(c). Check your answer against the true solution u(x) = cos(2x). Plot the error in the approximate solution and report the relative two-norm of the error. Compare your solution from this problem to your solution from 3(c) above. You should find that they are the same (up to rounding errors).

(e)   Extra credit (5 points): Why should the solution computed using the techniques from this problem be mathematically equivalent to the solution using the techniques from problem

3.