[Solved] MATH567-Homework 3

1. Linear two-point boundary value problems: Finite differences
   • Implement a general function in whatever language you like for numerically solving the general two-point boundary value problem (TPBVP)

u⁰⁰ = p(x)u⁰ + 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.

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

u⁰⁰ = 2tan(x)u⁰ + 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. x_j = 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

u⁰⁰ = p(x)u⁰ + q(x)u + r(x), x [a,b],

with mixed boundary conditions

u(a) = and ₁u(b) + ₂^u0(b) = ₃.

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 u_m+1:

,

where u(b) u_m+1, u(b h) u_m, p_m+1 = p(b), q_m+1 = q(b), and r_m+1 = r(b).

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

u⁰⁰(x) = f(x), a x b, u⁰(a) = ₀, u⁰(b) = ₁. (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 (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), f_j = f(x_j), and x_j = 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 A^T corresponding to a zero eigenvalue: A^Tw = 0 or w^TA = 0^T. 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 A^T.

A result from linear algebra (known as the Fredholm Alternative) is that the system Au = b has an infinite number of solutions if w^Tb = 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) = 1uTAu uTb

2 (4)

subject to w^Tu = (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 w^Tu = (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.

• 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) w^Tu = 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 w^T. This should give you the result that Au = 0, which means u = e for some . You can then use the equation w^Tu = 0 to show that = 0.

• Show that if w^Tb = 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 w^Tb 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(h²) then u will solve the BVP up to O(h²) since this is the truncation error of the FD method.

• Solve the BVP (6) numerically with a = 0, b = 2, m = 99, f(x) = 4cos(2x), and ₀ = ₁ = 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).
• Repeat part (c) but now for the function f(x) = x, ₀ = ², and ₁ = ². 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

u⁰⁰(x) = f(x), a x b, u⁰(a) = 0, u⁰(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), x_j = a + jh, u_j u(x_j), and f_j = f(x_j) for j = 0,,m + 1. In this approximation, the fictitious point idea has been applied, which used the assumption

and , (7)

where u₁ u(x₀ h) and u_m+2 u(x_m+1 + h). The method that we will use is based on the Discrete Cosine Transform (DCT).

• 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).

• 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 u_k is undefined. For this case, note that if f₀ = 0 (or

numerically zero) then u₀ can be chosen arbitrarily. Show how the condition f₀ = 0 corresponds to the discrete compatibility condition w^Tb = w^Tf = 0 discussed in the previous problem.

• Put parts (a) and (b) together to explain how to obtain the solution to (2) when the right hand satisfies w^Tf = 0. Also explain how one makes the solution unique by fixing the arbitrary constant to U.
• The Matlab codes 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).
• 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.