{ "cells": [ { "cell_type": "markdown", "id": "8eee5db0", "metadata": {}, "source": [ "---\n", "format:\n", " html: default\n", "---\n", "\n", "# A2: Stability & simple iterations {.unnumbered}\n", "\n", "Assignment 2: Floating point numbers, conditioning, and solving nonlinear equations in 1d
\n", "Published: Tue 22 Sept\n", "\n", "
\n", "Note:\n", "\n", "
\n", "\n", "
\n", " \n", " Download (.ipynb)\n", " \n", "
" ] }, { "cell_type": "markdown", "id": "8bfaaf1b", "metadata": {}, "source": [ "The assignment is worth **100 pts**. The symbols ✍️, 💻 indicate theory ($\\LaTeX$), practice (code), respectively." ] }, { "cell_type": "code", "execution_count": null, "id": "773dbca1", "metadata": {}, "outputs": [], "source": [ "Name = \"YOUR NAME HERE\"" ] }, { "cell_type": "markdown", "id": "58b5f49d", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 24, "id": "8732dea5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "order_table (generic function with 1 method)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# | code-fold: true\n", "\n", "using Plots, LaTeXStrings, Printf\n", "\n", "# simple iteration x_{n+1} = g(x_n) (from Lecture 4)\n", "function simple_iteration(g, x1; N=100, tol=1e-10, flag=false)\n", " x = Vector{Float64}(undef, N)\n", " x[1] = x1\n", " \n", " for n ∈ 2:N\n", " x[n] = g(x[n-1])\n", " \n", " if !isfinite(x[n])\n", " return x[1:n], false \n", " end\n", "\n", " # compute relative error\n", " rel_err = abs(x[n] - x[n-1]) / max(abs(x[n]), 1e-12)\n", " \n", " if rel_err < tol\n", " return x[1:n], true \n", " end\n", " end\n", " \n", " if flag\n", " @warn \"Fixed-point iteration did not converge within $N iterations.\"\n", " end\n", " return x, false\n", "end\n", "\n", "# prints the errors e_n = |x_n - ξ| together with \n", "# e_n / e_{n-1}^α (\\to μ if the order of convergence is α) \n", "# log(e_n) / log(e_{n-1}) (\\to order of convergence)\n", "function order_table(x, ξ; α=1)\n", " e = @. abs(x - ξ)\n", " @printf(\"%4s %22s %12s %16s %20s\\n\", \"n\", \"x_n\", \"e_n\", \"e_n/e_{n-1}^α\", \"log e_n/log e_{n-1}\")\n", " @printf(\"%4d %22.16f %12.3e\\n\", 1, x[1], e[1])\n", " for n ∈ 2:length(x)\n", " @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]))\n", " end\n", "end" ] }, { "cell_type": "markdown", "id": "9598f32d", "metadata": {}, "source": [ "
\n", "\n", "## A: Conditioning & Stability {.unnumbered}" ] }, { "cell_type": "markdown", "id": "21b42c60", "metadata": {}, "source": [ "***Exercise 1. (Harmonic numbers)***\n", "\n", "The harmonic numbers are defined as $H_N := \\sum_{k=1}^N \\frac1k$. \n", "\n", "(i) 💻 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? \n", "\n", "(ii) ✍️ Which approximation is more accurate? Explain why.\n", "\n", "*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```." ] }, { "cell_type": "code", "execution_count": null, "id": "45be8256", "metadata": {}, "outputs": [], "source": [ "y = 16.695311365859851815399118939540451884249869752373080462785\n", "\n", "# your code here" ] }, { "cell_type": "markdown", "id": "8b8c666a", "metadata": {}, "source": [ "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "e9fd78c1", "metadata": {}, "source": [ "***Exercise 2. Conditioning vs stability*** \n", "\n", "Let $f(x) := \\sqrt{x+1} - \\sqrt{x}$ for $x > 0$.\n", "\n", "(i) ✍️ Compute the relative condition number $\\kappa_f(x)$ and show that $\\kappa_f(x) \\leq \\frac12$ for all $x > 0$. \n", "\n", "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "(ii) 💻 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)$." ] }, { "cell_type": "code", "execution_count": null, "id": "2d262b29", "metadata": {}, "outputs": [], "source": [ "# your code here" ] }, { "cell_type": "markdown", "id": "0dca633c", "metadata": {}, "source": [ "(iii) ✍️ Explain why your results in (ii) do not contradict (i).\n", "\n", "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "(iv) ✍️💻 Find a mathematically equivalent formula for $f$ which avoids subtractive cancellation and repeat (ii) using your new formula.\n", "\n", "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "547c6f8c", "metadata": {}, "outputs": [], "source": [ "# your code here" ] }, { "cell_type": "markdown", "id": "9807762e", "metadata": {}, "source": [ "
\n", "\n", "\n", "## B: Simple iteration {.unnumbered}" ] }, { "cell_type": "markdown", "id": "3b8b1c08", "metadata": {}, "source": [ "***Exercise 3.*** \n", "\n", "Let $f(x) := x^3 + 4x^2 - 10$ and consider\n", "\n", "\\begin{align}\n", " g_1(x) &:= x - x^3 - 4x^2 + 10, \\nonumber\\\\ \n", " g_2(x) &:= \\frac12 \\sqrt{10 - x^3}, \\nonumber\\\\ \n", " g_3(x) &:= \\sqrt{\\frac{10}{4+x}}. \\nonumber\n", "\\end{align}\n", "\n", "(i) ✍️ Show that $f$ has a root $\\xi \\in [1,2]$. \n", "\n", "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "(ii) ✍️ Show that $\\xi$ is a fixed point of each of $g_1, g_2$, and $g_3$. \n", "\n", "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "(iii) 💻 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? " ] }, { "cell_type": "code", "execution_count": 23, "id": "bd5d1a39", "metadata": {}, "outputs": [], "source": [ "# your code here" ] }, { "cell_type": "markdown", "id": "119cbbf5", "metadata": {}, "source": [ "(iv) ✍️ 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]$. \n", "\n", "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "\n", "(v) ✍️ 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]$. \n", "\n", "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "(vi) ✍️💻 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?\n", "\n", "
Answer. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "481d3e13", "metadata": {}, "outputs": [], "source": [ "# your code here" ] }, { "cell_type": "markdown", "id": "26b7b965", "metadata": {}, "source": [ "
\n", "\n", "\n", "## C. Reflection {.unnumbered}\n", "\n", "**Goal:** Reflect on your learning strategies, track your workload, and help us improve future course material. \n", "**Time:** take 5-10 mins (this shouldn't take too long)\n", "\n", "Estimate the approximate hours you spent: \n", "* **Assignment:** ____ \n", "* **Reading notes/materials:** ____ \n", "* **Exercises:** ____ \n", "\n", "(Reminder: you should be spending around 8 hours per week on this course outside of the 4 hours of lectures)\n", "\n", "
Any other comments: \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "Take a moment to reflect on your learning over the past two weeks:\n", "\n", "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?\n", "\n", "
Comments: \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
\n", "\n", "(Optional) Any changes that would improve the lectures, exercises, or assignments for the future? You can also suggest things anonymously \n", " \n", " here\n", " (opens in a new tab)\n", " \n", "\n", "
Comments: \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.12.6", "language": "julia", "name": "julia-1.12" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.12.6" } }, "nbformat": 4, "nbformat_minor": 5 }