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

Brief Solution to Midterm 01

Oct 27, 2020.

1. (10 pts) The (fictional) one-and-half precision format uses 48 bits to store a binary floating point number of the form ±1.a_{1}a_{2} ⋯ a_{t}×2^{e} where a_{j} ∈ {0, 1}, −510 ≤ e ≤ 511. Find t and \underline{derive} an upper bound for relative error caused by \underline{rounding}. Express your final answer as a real number, but need not convert it to decimal expression.

Ans:

There are total 1022 different exponents (−510 ≤ e ≤ 511). It takes 10 bits to give 1022 or more different exponents (2^{10} = 1024). (2 pts) Total bits = 1 + t + 10 = 48 ⇒ t = 37 (2 pts).

Let x = ±1.a_{1}a_{2} ⋯ a_{t} … × 2^{e}. If a_{t+1} = 0, then fl_{round}(x) = ±1.a_{1}a_{2} ⋯ a_{t} × 2^{e}. A bound for the relative error is

\frac{|x − fl_{round}(x)|}{|x|} = \frac{|0.a_{t+1}a_{t+2} … |}{|1.a_{1}a_{2} … a_{t}a_{t+1} … |} × 2^{−t} ≤ 2^{−(t+1)}. (2 pts)

If a_{t+1} = 1, then fl_{round}(x) = ±(1.a_{1}a_{2} ⋯ a_{t} + 2^{−t}) × 2^{e}. The upper bound for relative error becomes

\frac{|x − fl_{round}(x)|}{|x|} = \frac{|1 − 0.a_{t+1}a_{t+2} … |}{|1.a_{1}a_{2} … a_{t}a_{t+1} … |} × 2^{−t} ≤ 2^{−(t+1)}. (2 pts)

Therefore, an upper bound for relative error caused by rounding is 2^{−38} (2 pts).

2. (10 pts) Given p_{0}, p_{1} and p_{2}, the general solution to the recursion formula p_{n} = \frac{10}{3}p_{n−1} − 3p_{n−2} + \frac{2}{3}p_{n−3} is p_{n} = c_{1}1^{n} + c_{2}2^{n} + c_{3}(\frac{1}{3})^{n} (need not show this). Find all (c_{1}, c_{2}, c_{3}) ≠ (0, 0, 0) such that the above iteration is unstable in relative error. \underline{Explain}.

Ans: Relative error = |\frac{e_{n}}{p_{n}^{exact}} |. (3 pts)

Note that

e_{n} ≈= d_{1}1^{n} + d_{2}2^{n} + d_{3}(\frac{1}{3})^{n} (3pts)

d_{1}, d_{2} and d_{3} are of O(p_{i} − fl(p_{i})), i = 0, 1, 2. Therefore d_{1}, d_{2} and d_{3} are of O(ε_{M}).(1 pts).

Therefore, relative error = \frac{O(ε_{M})1^{n}+O(ε_{M})2^{n}+O(ε_{M})( \frac{1}{3} )^{n}}{c_{1}1^{n}+c_{2}2^{n}+c_{3}( \frac{1}{3} )^{n}} . It is stable if and only if c_{2} ≠ 0. (3 pts)

3. (10 pts) The first few iteration (p_{i}, f(p_{i})), i = 0, 1, 2, 3, 4 of method of false position for some equation f(x) = 0 is given by

(0, −2), (3, 1), (∗, 2), (∗, 1), (∗, \frac{2}{9})

Find p_{5} (4 digits will do). \underline{Explain}.

Ans:

f(p_{1})f(p_{0}) < 0 ⇒ a = p_{0}, b = p_{1}, p_{2} = b − f(b)\frac{b − a}{f(b) − f(a)} = \frac{af(b) − bf(a)}{f(b) − f(a)} = 2

f(p_{2})f(p_{0}) < 0 ⇒ a = p_{0}, b = p_{2}, p_{3} = b − f(b)\frac{b − a}{f(b) − f(a)} = \frac{af(b) − bf(a)}{f(b) − f(a)} = 1

f(p_{3})f(p_{0}) < 0 ⇒ a = p_{0}, b = p_{3}, p_{4} = b − f(b)\frac{b − a}{f(b) − f(a)} = \frac{af(b) − bf(a)}{f(b) − f(a)} = \frac{2}{3}

f(p_{4})f(p_{0}) < 0 ⇒ a = p_{0}, b = p_{4}. p_{5} = b − f(b)\frac{b − a}{f(b) − f(a)} = \frac{af(b) − bf(a)}{f(b) − f(a)} = 0.6

4. (15 pts) Use any method to find a solution of \sqrt{1 + 0.9x} − \sqrt{1 − 0.8x} = 1.0 × 10^{−10} to 15 correct digits. You need to prevent loss of accuracy. Standard methods only gives you about 5 correct digits (and 1/3 partial credits).

Ans:

Apply the following identity

a^{2} − b^{2} = (a + b)(a − b)

that avoids the subtraction of two nearly identical numbers and gives

f(x) = \frac{1.7x}{\sqrt{1 + 0.9x} + \sqrt{1 − 0.8x}} − 10^{−10}. (7 pts)

Then solve f(x) = 0 by any numerical method to find the solution

x_{∗} ≈ 1.17647058823875 × 10^{−10}. (8 pts)

5. (10+5 pts) It is known that the unique solution to f(x) = x + 3 sin(x) − 0.01 = 0 is located near x = 0.

(a) Find a fixed point iteration that will converge for any x_{0} ∈ [−\frac{1}{2}, \frac{1}{2}]. Show that your method satisfies the assumptions of a relevant Theorem, but need not prove the Theorem again. You can use the numerical values of sin(\frac{1}{2}), cos(\frac{1}{2}), exp(\frac{1}{2}), etc. in your proof.

(b) Find an N (need not be optimal) such that |x_{n} − x^{∗}| < 10^{−30} for all n ≥ N with x_{0} = 0 (assuming a higher precision floating point arithmetic is used).

Ans:

(a) Direct fixed point iteration with x^{(k+1)} = g_{0}^{(k)}(x) = 0.01 − 3 sin(x^{(k)}) does not converge. Instead, a proper choice of α and g(x) = αx + (1 − α)g_{0}(x) will result in local convergence (2 pts). One could choose

α = \frac{g′_{0}(ξ)}{g′_{0}(ξ) − 1}

for some ξ near 0. Since ξ ≈ 0, g′_{0}(ξ) ≈ −3, we take α = \frac{−3}{−3−1} = \frac{3}{4}. (2 pts) Since g(x) = \frac{3}{4}(x − sin x) + 0.0025, g′(x) = \frac{3}{4}(1 − cos(x)) Therefore

0 < \frac{3}{4} (1 − cos(\frac{1}{2})) ≤ g′(x) ≤ \frac{3}{4} on [−\frac{1}{2}, \frac{1}{2}].

It follows that g is an increasing function on [−\frac{1}{2}, \frac{1}{2}],

−\frac{1}{2} < −0.012931... = g(−\frac{1}{2}) ≤ g(x) ≤ g(\frac{1}{2}) = 0.017931... < \frac{1}{2}

(thus g([−\frac{1}{2}, \frac{1}{2}]) ⊂ [−\frac{1}{2}, \frac{1}{2}] (3 pts) ) and

|g′(x)| = |\frac{3}{4}(1 − cos(x))| ≤ \frac{3}{4} = k < 1 ∀ x ∈ (−\frac{1}{2}, \frac{1}{2}) (3 pts)

(b) Correct estimate (3 pts) and correct N (2 pts). [Method 1]

|x_{n}−x_{∗}| ≤ k^{n} max\{x_{0}−a, b−x_{0}\} = (\frac{3}{4})^{n} max \{0 − (−\frac{1}{2}) , \frac{1}{2} − 0\} = (\frac{3}{4})^{n} \frac{1}{2} < 10^{−30}

⇒ n > \frac{log_{10} 2 − 30}{log_{10} \frac{3}{4}} = 237.71... ⇒ N = 238.

[Method 2]

|x_{n} − x_{∗}| ≤ \frac{k^{n}}{1 − k} |x_{1} − x_{0}| = 4 (\frac{3}{4})^{n} |g(0) − 0| = (\frac{3}{4})^{n} 0.01 < 10^{−30}

⇒ n > \frac{28}{log_{10} \frac{4}{3}} = 224.11... ⇒ N = 225.

6. (15 pts) Give a (at least) cubically convergent method to solve for e^{x} −1 = 0. Give the formula and prove that it is at least cubically convergent (locally). If you cannot do it, do the same for a locally (at least) quadratically convergent method for 1/3 partial credit.

Answer:

[Cubic] One solution is given by (there may be others)

x_{n+1} = g(x_{n}), g(x) = x − \frac{f(x)}{f′(x)} − \frac{f′′(x)}{2f′(x)} [\frac{f(x)}{f′(x)}]^{2} . (5 pts)

Check that g′(p) = g′′(p) = 0 and g′′′(p) ≠ 0. Therefore

p_{n+1} − p = g(p_{n}) − g(p) = \frac{g′′′(ξ)}{3!}(p_{n} − p)^{3}

and

lim_{n→∞} \frac{|p_{n+1} − p|}{|p_{n} − p|^{3}} = \frac{|g^{(3)}(p)|}{3!}. (10 pts)

[Quadratic]

x_{n+1} = g(x_{n}), g(x) = x − \frac{f(x)}{f′(x)}. (2 pts)

Check that g′(p) = 0 and g′′(p) ≠ 0. Therefore

p_{n+1} − p = g(p_{n}) − g(p) = \frac{g′′(ξ)}{2!}(p_{n} − p)^{2}

and

lim_{n→∞} \frac{|p_{n+1} − p|}{|p_{n} − p|^{2}} = \frac{|g′′(p)|}{2!}. (3 pts)

7. (15 pts) One way of computing π is given by the Wallis formula

\frac{π}{2} = ∏_{n=1}^{∞} (\frac{(2n)^{2}}{(2n − 1)(2n + 1)}) .

The N-term approximation is therefore given by

\frac{π_{N}}{2} = (\frac{2 · 2}{1 · 3}) (\frac{4 · 4}{3 · 5}) (\frac{6 · 6}{5 · 7}) ⋯ (\frac{2N · 2N}{(2N − 1) · (2N + 1)})

To prevent overflow too quickly, it is better to evaluate the last multiplication as ∗(2N)/(2N − 1) ∗ (2N)/(2N + 1). Find the rate of convergence of lim_{n→∞} π_{n} = π numerically. Extra points without using the limit π explicitly. Hint: The convergence is slow, try not to produce all the data points. For example, π_{100}, π_{200}, ⋯, π_{10000} should be enough to analyze. To check if an integer j is a multiple of 100 or not, you can use mod or check round(j/100) ∗ 100. Ans: Programming:

Completed program that can produce all π_{n} needed: (10 pts) Analysis: Method 1: Try semilogy and loglog plot of |π_{n} − π| vs n to determine whether |π_{n} − π| ≈ Cn^{−p} or |π_{n} − π| ≈ Cα^{−n} or something else. The results indicates that |π_{n} − π| ≈ Cn^{−p} where −p is the slope in the loglog plot and p ≈ 1. Comparing the loglog plot of |π_{n} − π| vs n, with the loglog plot of n^{−1} vs n confirms that p = 1 (5 pts). Method 2 (extra 5 pts): Proceed to find out the constants C, p or C, α. Directly try π_{n} − π ≈ Cn^{−p} and find p through

p ≈ log_{2} \frac{a_{N} − a_{2N}}{a_{2N} − a_{4N}}

Different choices of N = 300, 500, ⋯, 1000, 1500 all give consistent answer p ≈ 1. Conclusion: π_{n} − π = O(n^{−1})

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

sin(x) + \frac{2y}{1 + x} = 0.01, 5x + sin(\frac{6y}{1 + y^{2}} ) = 0.02.

Write your answer in the format of ’format long e’. Hint: the solution is near (0, 0) where sin x ≈ x, \frac{y}{1+x} ≈ y, etc. to leading order. Ans: Method 1: Let

g_{1}(x) = \frac{1}{5} (0.02 − sin (\frac{6y}{1 + y^{2}} ))

g_{2}(x) = \frac{1}{2}(0.01 − sin(x))(1 + x)

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. From the Hint, we can take x_{0} = (0, 0)^{T} and

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

Method 2: From the Hint, the linear approximation of the left hand side is (x + 2y, 5x + 6y), so we can rewrite the equation as

x + 2y = 0.01 + x + 2y − sin(x) − \frac{2y}{1 + x} ≡ h_{1}(x, y),

5x + 6y = 0.02 + 6y − sin(\frac{6y}{1 + y^{2}} ) ≡ h_{2}(x, y).

This suggests the fixed point iteration:

( \begin{matrix}x^{(k+1)} \\ y^{(k+1)} \end{matrix} ) = ( \begin{matrix}1 & 2 \\ 5 & 6 \end{matrix} )^{−1} ( \begin{matrix}h_{1}(x^{(k)}, y^{(k)}) \\ h_{2}(x^{(k)}, y^{(k)}) \end{matrix} )

Method 3: Newton’s Method. (Correct algorithm of any convergent method = 10 pts) If x_{0} = (0, 0)^{T}, then the iteration gives

x_{∗} ≈ (−4.882452175e − 03, 7.404885005e − 03)^{T}. (5 pts)