Numerical Analysis I, Fall 2017 (http://www.math.nthu.edu.tw/˜ wangwc/)

Final Exam

Jan 12, 2018.

1. (15 pts) Use any desingularization method to evaluate the improper integral

∫_{0}^{∞} \frac{1}{1 + x^{3}} dx

with a composite quadrature rule of 4th order accuracy. Demonstrate numerically that the result is indeed 4th order accurate.

Answer: \underline{Method 1.} Divide the integral into ∫_{0}^{1} and ∫_{1}^{∞}. Replace x = s^{−q} in the second integral to obtain

∫_{1}^{∞} \frac{1}{1 + x^{3}} dx = q ∫_{0}^{1} \frac{s^{2q−1}}{s^{3q} + 1} ds.

To make 4th order accuracy, we take q = 1. (5 pts)

Apply Simpson’s rule or Gaussian quadrature, then we can demonstrate numerically the result

log_{2} \frac{I(100) − I(200)}{I(200) − I(400)} ≈ 4. (10 pts)

\underline{Method 2.} Use change of variable t = \frac{1}{1+x}.

C lagrange. (Extra 3 pts)

2. (15 pts) Consider discretizing

(∂_{x}^{2} + ∂_{y}^{2})u(x, y) = f(x, y), (x, y) ∈ (0, 5) × (0, 1) u = 0, on the boundary of (0, 5) × (0, 1)

with uniformly spaced grids 0 = x_{0} < x_{1} < ⋯ < x_{5N} = 5, 0 = y_{0} < y_{1} < ⋯ < y_{N} = 1, x_{i} − x_{i−1} = y_{j} − y_{j−1} = h = \frac{1}{N}, using second order centered finite difference method. The following two possible ordering of the unknowns u_{i,j},

(i, j) = (1, 1) → (1, 2) → ⋯ → (1, N − 1) → (2, 1) → (2, 2) → ⋯,  (1)

and

(i, j) = (1, 1) → (2, 1) → ⋯ → (5N − 1, 1) → (1, 2) → (2, 2) → ⋯,  (2)

result in two linear systems A^{(1)}u^{(1)} = f^{(1)} and A^{(2)}u^{(2)} = f^{(2)}. For large N, which linear system is faster to solve using Gaussian elimination without pivoting, combined with backward substitution? or are they the same? Explain.

Answer: For the 1st ordering, the leading order term of number of operations is 5N^{4}. (7 pts) For the 2nd ordering, the leading order term of number of operations is 125N^{4}. (7 pts) Thus, the 1st one is faster to solve the system by Gaussian elimination. (1 pt)

3. (7+7 pts)

(a) Perform the required row interchanges for the following linear system using scaled partial pivoting:

\begin{aligned}& x_{2}+ & x_{3} = & 4 \\ x_{1}− & 2x_{2}− & x_{3} = & 5 \\ x_{1}− & x_{2}+ & x_{3} = & 6 \end{aligned} ,  (3)

Give all details and the final linear system. Need not solve it. Answer:

[ \begin{matrix}0 & 1 & 1 & | & 4 \\ 1 & −2 & −1 & | & 5 \\ 1 & −1 & 1 & | & 6 \end{matrix} ] ⟶^{E_{1}↔E_{3}} [ \begin{matrix}1 & −1 & 1 & | & 6 \\ 1 & −2 & −1 & | & 5 \\ 0 & 1 & 1 & | & 4 \end{matrix} ] ⟶ [ \begin{matrix}1 & −1 & 1 & | & 6 \\ 0 & −1 & −2 & | & −1 \\ 0 & 1 & 1 & | & 4 \end{matrix} ] (4 pts)

⟶^{E_{2}↔E_{3}} [ \begin{matrix}1 & −1 & 1 & | & 6 \\ 0 & 1 & 1 & | & 4 \\ 0 & −1 & −2 & | & −1 \end{matrix} ] ⟶ [ \begin{matrix}1 & −1 & 1 & | & 6 \\ 0 & 1 & 1 & | & 4 \\ 0 & 0 & −1 & | & 3 \end{matrix} ] (3 pts)

(b) Find P, L, U in the factorization P A = LU for the matrix

A = ( \begin{matrix}0 & 0 & −1 & 1 \\ −1 & −1 & 2 & 0 \\ 1 & 1 & −1 & 2 \\ 1 & 2 & 0 & 2 \end{matrix} )  (4)

where P is the permutation matrix corresponding to partial pivoting. Answer: Finding U:

U = ( \begin{matrix}−1 & −1 & 2 & 0 \\ 0 & 1 & 2 & 2 \\ 0 & 0 & 1 & 2 \\ 0 & 0 & 0 & 3 \end{matrix} ) (3 pts)

Finding P:

P = ( \begin{matrix}0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 \end{matrix} ) (2 pts)

Finding L:

L = ( \begin{matrix}1 & 0 & 0 & 0 \\ −1 & 1 & 0 & 0 \\ −1 & 0 & 1 & 0 \\ 0 & 0 & −1 & 1 \end{matrix} ) (2 pts)

4. (15 pts) True or False? Explain.

Suppose that A = M − N is nonsingular. Then the iteration Mx^{(k+1)} = b + Nx^{(k)} converges for any x^{(0)} if M is nonsingular and T = M^{−1}N satisfies ∥T ∥ < 1 for some matrix norm ∥ · ∥.

Answer: True. Let x be the exact solution and E_{k} := x^{(k)} − x. Then

∥E_{k}∥ = ∥x^{(k)} − x∥ = ∥(T x^{(k−1)} + M^{−1}b) − (T x + M^{−1}b)∥ = ∥T x^{(k−1)} − T x∥ = ∥T E_{k−1}∥

≤ ∥T ∥∥E_{k−1}∥ ≤ … ≤ ∥T ∥^{k}∥E_{0}∥ →^{∥T ∥<1} 0

as k → ∞. (15 pts)

5. (15 pts) Write down the definition of ’condition number’. Derive a related inequality that gives an estimate of relative error of the solution of a linear system.

Answer: See Definition 7.28. for the definition. (3 pts) See Theorem 7.27 for the inequality. Correct inequality. (3 pts) Details of derivations:∥b∥ ≤ ∥A∥∥x∥, ∥Δx∥ ≤ ∥A^{−1}∥∥Δb∥. (9 pts)

6. (15 pts) Consider the integral equation

u(x_{i}) = f(x_{i}) + ∫_{0}^{1} K(x_{i}, t) u(t) dt, i = 0, ⋯, M,  (5)

where 0 = x_{0} < x_{1} < ⋯ < x_{M} = 1, with x_{j} − x_{j−1} = h = \frac{1}{M}, f(x) = x^{2}, K(x, t) = −\frac{1}{3}e^{−|x−t|} and ∫_{0}^{1} dt is evaluated using trapezoidal rule.

Form the linear system, find an convergent iterative method, explain why your method converges. Need not solve it.

Answer: The Jacobi iteration will converge. Consider

∥T_{j}∥_{∞} = max_{i=1,... M+1} ∑_{j=1}^{M+1} |(T_{j})_{ij}| ≈ max_{i=1,... M+1} \frac{\frac{1}{3}}{\frac{h}{3} + 1} ∫_{0}^{1} e^{−|x_{i}−t|} dt

= max_{i=1,... M+1} \frac{1}{h + 3}(2 − e^{−x_{i}} − e^{x_{i}−1})
≈ max_{i=1,... M+1} \frac{1}{3}(2 − e^{−x_{i}} − e^{x_{i}−1})

≈ max_{x∈[0,1]} \frac{1}{3}(2 − e^{−x} − e^{x−1})
= \frac{2}{3} (1 − \frac{1}{\sqrt{e}}) < 1.

Find the linear system. (5 pts) Find an convergent iterative method. (5 pts) Explain why your method converges. (5 pts)

7. (15 pts) Find the polynomial p(x) that minimizes ∫_{0}^{1}(sin x − p(x))^{2}dx among all odd polynomials (p(−x) = −p(x)) with degree less than or equal to 5. Derive the linear system and need not solve it numerically.

Answer: Let p(x) = a_{1}x + a_{3}x^{3} + a_{5}x^{5}. Then the minimization problem results in the linear system

\frac{1}{3}a_{1} + \frac{1}{5}a_{3} + \frac{1}{7}a_{5} = ∫_{0}^{1} x sin x dx (5 pts)

\frac{1}{5}a_{1} + \frac{1}{7}a_{3} + \frac{1}{9}a_{5} = ∫_{0}^{1} x^{3} sin x dx (5 pts)

\frac{1}{7}a_{1} + \frac{1}{9}a_{3} + \frac{1}{11}a_{5} = ∫_{0}^{1} x^{5} sin x dx. (5 pts)

8. (8 pts) Find the constants a, b, c so that, the quadrature

∫_{−h}^{h} f(x)dx ≈ h · (af(−h) + bf(0) + cf(h)) has highest degree of precision, where f is

sufficiently smooth. Under the assumption (need not prove this assumption) that the error of this quadrature is of the form

∫_{−h}^{h} f(x)dx − h · (af(−h) + bf(0) + cf(h)) = Kf^{(n)}(ξ)h^{p},

find K, n and p.

Ans: Details for finding a, b, c. (1 pt) a = \frac{1}{3} (1 pt), b = \frac{4}{3} (1 pt), c = \frac{1}{3} (1 pt), DOP= 3. Details for finding K, n, p. (1 pt) K = −\frac{1}{90} (1 pt), n = 4 (1 pt), p = 5 (1 pt).

9. (8 pts) Use any method to solve the nonlinear system of equations

\begin{aligned}1x_{1} + 2x_{2} + 0.03 ∗ sin(x_{1} + x_{2}) & = 4 \\ 5x_{1} + 6x_{2} + 0.07 ∗ cos(x_{1} − x_{2}) & = 8 \end{aligned}

Write your answer in the format of ’format long e’ and hand in your code.

Answer: Let

g_{1}(x) = −2x_{2} − 0.03 sin(x_{1} + x_{2}) + 4

g_{2}(x) = −\frac{5}{6}x_{1} − \frac{0.07}{6} cos(x_{1} − x_{2}) + \frac{8}{6}

and

G(x) = (g_{1}(x), g_{2}(x))^{t}
Ḡ(x) = αx + (I − α)G(x)

where α is a 2 × 2 matrix and I is the identity matrix.

Consider

x_{n+1} = Ḡ(x_{n}) = αx_{n} + (I − α)G(x_{n})
x_{∗} = Ḡ(x_{∗}) = αx_{∗} + (I − α)G(x_{∗})

⇒^{MV T} x_{n+1} − x_{∗} = α(x_{n} − x_{∗}) + (I − α) [ \begin{matrix}\frac{∂g_{1}}{∂x_{1}} (z_{n,1}) & \frac{∂g_{1}}{∂x_{2}} (z_{n,1}) \\ \frac{∂g_{2}}{∂x_{1}} (z_{n,2}) & \frac{∂g_{2}}{∂x_{2}} (z_{n,2}) \end{matrix} ] (x_{n} − x_{∗}). − −(1)

for some z_{n,1}, z_{n,2} on the straight line path connecting x_{n} and x_{∗}. Since the nonlinear term is small, so the solution is near that of the linear part.

1x_{1} + 2x_{2} = 4
5x_{1} + 6x_{2} = 8
⇒ x_{1} = −2, x_{2} = 3.

Therefore, we choose initial point x_{0} = (−2, 3)^{t} and approximate

z_{n,1} ≈ z_{n,2} ≈ x_{0}. − −(2)

Combine (1) and (2), then

x_{n+1} − x_{∗} ≈ α(x_{n} − x_{∗}) + (I − α)DG(x_{0})(x_{n} − x_{∗}) = (α + (I − α)DG(x_{0}))(x_{n} − x_{∗}).

To accelerate the convergence, we choose α such that

α + (I − α)DG(x_{0}) = 0
⇒ α = (DG(x_{0}) − I)^{−1}DG(x_{0}).

The approximation solution is (−1.97005844315465e + 00, 2.97238830860496e + 00)^{t}.

Another possible iteration: x_{k+1} = [ \begin{matrix}1 & 2 \\ 5 & 6 \end{matrix} ]^{−1} [ \begin{matrix}4 − 0.03... \\ 8 − 0.07... \end{matrix} ].

Code+Answer. (8 pts, either all pts or no pts, partial credits only for typo in inputting data, outputting answer, or something else.) C language. (Extra 2 pts)