A2: Stability & simple iterations

Assignment 2: Floating point numbers, conditioning, and solving nonlinear equations in 1d
Published: Tue 22 Sept

Note:

Complete the following and submit to Canvas before Wed 30 Sept, 23:59.
Late work will receive 0%.
Each assignment is worth the same.
Please get in contact with the TA/instructor in plenty of time if you need help.
Before submitting your work, make sure to check everything runs as expected. Click Kernel > Restart Kernel and Run All Cells. Submit your notebook in .pdf format by clicking File > Export Notebook As > PDF or by printing to PDF in the browser.
Feel free to add more cells to experiment or test your answers.
I encourage you to discuss the course material and assignment questions with your classmates. However, unless otherwise explicitly stated on the assignment, you must complete and write up your solutions on your own.
The use of GenAI is prohibited as outlined in the course syllabus. If suspected of cheating, you may be asked to complete a written or oral exam on the content of this assignment.

The assignment is worth 100 pts. The symbols ✍️, 💻 indicate theory (\(\LaTeX\)), practice (code), respectively.

Code

Name = "YOUR NAME HERE"

The following cell contains some code that you may find useful: simple_iteration is the function from Lecture 4 and order_table(x, ξ; α) prints the errors \(e_n := |x_n - \xi|\) together with \(\frac{e_{n}}{e_{n-1}^\alpha}\) and \(\frac{\log e_{n}}{\log e_{n-1}}\) (which are useful when estimating the asymptotic error constant and the order of convergence). You may need to ignore the last few rows of the table once \(e_n\) is of the order of machine precision.

Code

using Plots, LaTeXStrings, Printf

# simple iteration x_{n+1} = g(x_n) (from Lecture 4)
function simple_iteration(g, x1; N=100, tol=1e-10, flag=false)
    x = Vector{Float64}(undef, N)
    x[1] = x1
    
    for n ∈ 2:N
        x[n] = g(x[n-1])
        
        if !isfinite(x[n])
            return x[1:n], false 
        end

        # compute relative error
        rel_err = abs(x[n] - x[n-1]) / max(abs(x[n]), 1e-12)
        
        if rel_err < tol
            return x[1:n], true 
        end
    end
    
    if flag
        @warn "Fixed-point iteration did not converge within $N iterations."
    end
    return x, false
end

# prints the errors e_n = |x_n - ξ| together with 
#   e_n / e_{n-1}^α          (\to μ if the order of convergence is α)  
#   log(e_n) / log(e_{n-1})  (\to order of convergence)
function order_table(x, ξ; α=1)
    e = @. abs(x - ξ)
    @printf("%4s %22s %12s %16s %20s\n", "n", "x_n", "e_n", "e_n/e_{n-1}^α", "log e_n/log e_{n-1}")
    @printf("%4d %22.16f %12.3e\n", 1, x[1], e[1])
    for n ∈ 2:length(x)
        @printf("%4d %22.16f %12.3e %16.6f %20.6f\n", n, x[n], e[n], e[n]/e[n-1]^α, log(e[n])/log(e[n-1]))
    end
end

order_table (generic function with 1 method)

A: Conditioning & Stability

Exercise 1. (Harmonic numbers)

The harmonic numbers are defined as \(H_N := \sum_{k=1}^N \frac1k\).

💻 Compute \(H_{N}\) for \(N = 10^7\) in Float32 (single precision) by summing forwards (\(k = 1, 2, \dots, N\)) and backwards (\(k = N, N-1, \dots, 1\)). What are the relative errors of these two approximations?
✍️ Which approximation is more accurate? Explain why.

Hint: You may use y defined below as the “exact” value of \(H_{10^7}\). You can use Float32(1) or 1f0 for \(1\) in Float32. Use typeof( x ) to check the type of variable x.

Code

y = 16.695311365859851815399118939540451884249869752373080462785

# your code here

Answer.

Exercise 2. Conditioning vs stability

Let \(f(x) := \sqrt{x+1} - \sqrt{x}\) for \(x > 0\).

✍️ Compute the relative condition number \(\kappa_f(x)\) and show that \(\kappa_f(x) \leq \frac12\) for all \(x > 0\).

Answer.

💻 Evaluate \(f(x)\) (using the formula above) for \(x = 10^2, 10^4, \dots, 10^{16}\) in Float64 and compute the relative errors. You may use sqrt(big(x)+1) - sqrt(big(x)) (which uses higher precision BigFloat arithmetic) as the “exact” value of \(f(x)\).

Code

# your code here

✍️ Explain why your results in (ii) do not contradict (i).

Answer.

✍️💻 Find a mathematically equivalent formula for \(f\) which avoids subtractive cancellation and repeat (ii) using your new formula.

Answer.

Code

# your code here

B: Simple iteration

Exercise 3.

Let \(f(x) := x^3 + 4x^2 - 10\) and consider

\[\begin{align} g_1(x) &:= x - x^3 - 4x^2 + 10, \nonumber\\ g_2(x) &:= \frac12 \sqrt{10 - x^3}, \nonumber\\ g_3(x) &:= \sqrt{\frac{10}{4+x}}. \nonumber \end{align}\]

✍️ Show that \(f\) has a root \(\xi \in [1,2]\).

Answer.

✍️ Show that \(\xi\) is a fixed point of each of \(g_1, g_2\), and \(g_3\).

Answer.

💻 Run the simple iterations \(x_{n+1} = g_j(x_n)\) with \(x_1 = 1.5\) for \(j = 1,2,3\). Which of these iterations appear to converge?

Code

# your code here

✍️ Show that \(g_3 : [1,2] \to [1,2]\) and that \(|g_3'(x)| \leq L < 1\) for all \(x \in [1,2]\). Explain why this means that \(x_{n+1} = g_3(x_n)\) converges to \(\xi\) for all \(x_1 \in [1,2]\).

Answer.

✍️ Show that the contraction mapping theorem can not be applied to \(g_2\) on \([1,2]\) but that it can be applied to \(g_2\) on \([1, 3/2]\).

Answer.

✍️💻 What are the order of convergence and the asymptotic error constants of the iterations \(x_{n+1} = g_2(x_n)\) and \(x_{n+1} = g_3(x_n)\)? Compare your theoretical values with numerical estimates (you may use order_table). Which iteration gains more digits of accuracy per iteration and approximately how many?

Answer.

Code

# your code here

C. Reflection

Goal: Reflect on your learning strategies, track your workload, and help us improve future course material.
Time: take 5-10 mins (this shouldn’t take too long)

Estimate the approximate hours you spent:
* Assignment: ____
* Reading notes/materials: ____
* Exercises: ____

(Reminder: you should be spending around 8 hours per week on this course outside of the 4 hours of lectures)

Any other comments:

Take a moment to reflect on your learning over the past two weeks:

What was the easiest topic or task and why? Where did you struggle, and how did you overcome those challenges? Which resources (e.g. lecture notes, practice exercises, peer discussions, …) helped you the most? How did you manage your time and workload (e.g., all at once, spread out over several days, or working with peers)? Looking back, what would you do differently if you were starting these two weeks over today?

Comments:

(Optional) Any changes that would improve the lectures, exercises, or assignments for the future? You can also suggest things anonymously here (opens in a new tab)

Comments: