{
"cells": [
{
"cell_type": "markdown",
"id": "53fdd565",
"metadata": {},
"source": [
"# NYU Paris - Numerical Analysis"
]
},
{
"cell_type": "markdown",
"id": "a6d022bb",
"metadata": {},
"source": [
"## Final Exam\n",
"\n",
"- Submit this notebook on Brightspace before 6:00 PM\n",
"\n",
"- The exam comprises six independent exercises. In each exercise, the\n",
" cells may depend on the previous cells.\n",
"\n",
"- To facilitate correction of your code,\n",
" do not change the signatures of the functions to implement.\n",
"\n",
"- Run the cell below to import necessary packages and\n",
" load a macro used for unit tests thoughout the notebook."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e738e4bb",
"metadata": {},
"outputs": [],
"source": [
"using ForwardDiff\n",
"using LinearAlgebra\n",
"using Polynomials\n",
"using SpecialFunctions\n",
"using TestImages\n",
"using LaTeXStrings\n",
"using Plots\n",
"\n",
"macro mark(bool_expr)\n",
" return :( print($bool_expr ? \"✔️\" : \"❌\"))\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "0e9ebb0d",
"metadata": {},
"source": [
"### 1. Floating point arithmetic\n",
"\n",
"Throughout this exercise,\n",
"$(\\bullet)_\\beta$ denotes base $\\beta$ representation.\n",
"True or false? (+1/0/-1)\n",
"1. It holds that\n",
" $\n",
" (1.0111)_2 + (0.1001)_2 = (10)_2.\n",
" $\n",
"1. It holds that\n",
" $\n",
" (1000)_5 \\times (0.003)_5 = 3.\n",
" $\n",
"1. It holds that\n",
" $\n",
" (0.\\overline{4})_5 = 4\n",
" $\n",
"1. In base 16, all the natural numbers from 1 to 1000 can be represented using 3 digits.\n",
"\n",
"1. Machine multiplication according to the IEEE754 standard is a commutative operation.\n",
"\n",
"1. Machine addition according to the IEEE754 standard is an associative operation.\n",
"\n",
"1. Only finitely many real numbers can be represented exactly in the `Float64` format.\n",
"\n",
"1. The spacing (in absolute value) between successive `Float64` numbers is strictly increasing.\n",
"\n",
"1. In Julia, `eps(Float32)` is the smallest positive number representable in the `Float32` format.\n",
"\n",
"1. In Julia, `nextfloat(1f0) - 1f0` equals the machine epsilon for the `Float32` format.\n",
"\n",
"1. Assume that $x \\in \\real$ belongs to $\\mathbf F_{64}$. Then $2x \\in \\mathbf F_{64}$.\n",
"\n",
"1. The real number $\\sqrt{2}$ can be represented exactly in the single-precision format."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "aaa7cbcc",
"metadata": {
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"answers_q1 = zeros(Int, 12)\n",
"\n",
"# Use -1 for false, 1 for true\n",
"answers_q1[1] = 0\n",
"answers_q1[2] = 0\n",
"answers_q1[3] = 0\n",
"answers_q1[4] = 0\n",
"answers_q1[5] = 0\n",
"answers_q1[6] = 0\n",
"answers_q1[7] = 0\n",
"answers_q1[8] = 0\n",
"answers_q1[9] = 0\n",
"answers_q1[10] = 0\n",
"answers_q1[11] = 0\n",
"answers_q1[12] = 0"
]
},
{
"cell_type": "markdown",
"id": "dde04a88",
"metadata": {},
"source": [
"### 2. Exercise on interpolation and approximation"
]
},
{
"cell_type": "markdown",
"id": "d7bd2802",
"metadata": {},
"source": [
"1. Find the quadratic polynomial $p(x) = ax^2 + b x + c$ that goes through the points $(0, 1)$, $(1, 3)$ and $(2, 7),$\n",
" and then plot the polynomial $p$ together with the data points on the same graph."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "25bc7ff6",
"metadata": {},
"outputs": [],
"source": [
"# Calculate a, b and c here\n",
"\n",
"@mark begin p(x) = a*x^2 + b*x + c; p(0) == 1 end\n",
"@mark begin p(x) = a*x^2 + b*x + c; p(1) == 3 end\n",
"@mark begin p(x) = a*x^2 + b*x + c; p(2) == 7 end\n",
"\n",
"# Create the plot here"
]
},
{
"cell_type": "markdown",
"id": "7b48172d",
"metadata": {},
"source": [
"2. Find the function $f$ of the form $f(x) = a \\sin(x) + b \\cos(x) + c \\sin(2x)$ that goes through the points $(x_1, y_1)$, $(x_2, y_2)$ and $(x_3, y_3)$,\n",
" where the data is given below.\n",
" Plot the function $f$ together with the data points on the same graph."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "46a58f5a",
"metadata": {},
"outputs": [],
"source": [
"x = [1, 2, 3]\n",
"y = [2.8313730233698577, -0.6797987415765313, -2.11828048333995]\n",
"\n",
"# Calculate a, b and c here\n",
"\n",
"@mark begin f(z) = a*sin(z) + b*cos(z) + c*sin(2z); f(x[1]) == y[1] end\n",
"@mark begin f(z) = a*sin(z) + b*cos(z) + c*sin(2z); f(x[2]) == y[2] end\n",
"@mark begin f(z) = a*sin(z) + b*cos(z) + c*sin(2z); f(x[3]) == y[3] end\n",
"\n",
"# Create the plot here"
]
},
{
"cell_type": "markdown",
"id": "905931ce",
"metadata": {},
"source": [
"3. Find the function $g$ of the form $g(x) = a \\sin(x) + b \\cos(x)$\n",
" such that the sum of squared residuals\n",
" $$\n",
" \\sum_{i=1}^{n} |g(x_i) - y_i|^2, \\qquad n = 10,\n",
" $$\n",
" is minimized for the data given below.\n",
" Plot on the same graph the function $g$ and the data points."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bba82caa",
"metadata": {},
"outputs": [],
"source": [
"x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]\n",
"y = [0.46618950237090034, 0.9365762460801551, 0.5907432672662498,\n",
" -0.21370874419278404, -0.8869010982313607, -0.6906040605442342,\n",
" 0.1784931250350807, 1.024036713535387, 0.7837248688922458,\n",
" -0.1544083454499539]\n",
"\n",
"# Calculate a and b and create the plot here"
]
},
{
"cell_type": "markdown",
"id": "56d0ae3b",
"metadata": {},
"source": [
"### 3. Exercise on numerical integration\n",
"\n",
"Our goal in this exercise is to write a program to calculate integrals of the form\n",
"$$\n",
"I[f] := \\int_{-\\infty}^{\\infty} f(x) \\, e^{-x^2} \\, d x.\n",
"$$\n",
"To this end, we will use Hermite polynomials,\n",
"which are orthogonal polynomials for the following inner product:\n",
" $$\n",
" \\langle f, g \\rangle := \\int_{-\\infty}^{\\infty} f(x) \\, g(x) \\, e^{-x^2} \\, d x\n",
" $$\n",
"1. Using that Hermite polynomials satisfy the recurrence relation\n",
" $$\n",
" H_{k + 1}(x) = 2x H_k(x) - H_k'(x), \\qquad H_0(x) = 1,\n",
" $$\n",
" write a function `hermite(n)` which returns the Hermite polynomial of degree $n$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4da66024",
"metadata": {},
"outputs": [],
"source": [
"function hermite(n)\n",
" # Write your code here\n",
"end\n",
"\n",
"@mark hermite(2) == Polynomial([-2, 0, 4])\n",
"@mark hermite(3) == Polynomial([0, -12, 0, 8])\n",
"@mark hermite(4) == Polynomial([12, 0, -48, 0, 16])\n",
"@mark degree(hermite(10)) == 10"
]
},
{
"cell_type": "markdown",
"id": "fde9bb30",
"metadata": {},
"source": [
"2. Write a function `integrate_monomial(n)`, which return the value of the integral\n",
" $$\n",
" \\int_{-\\infty}^{\\infty} x^n e^{-x^2} \\, dx =\n",
" \\frac{1}{2} \\Bigl((-1)^n + 1\\Bigr) \\Gamma\\left( \\frac{n+1}{2} \\right),\n",
" $$\n",
" where $\\Gamma$ is the Gamma function,\n",
" available as `gamma` in the `SpecialFunctions` library."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dfcc14cd",
"metadata": {},
"outputs": [],
"source": [
"function integrate_monomial(n)\n",
" # Write your code here\n",
"end\n",
"\n",
"@mark integrate_monomial(0) ≈ √π\n",
"@mark integrate_monomial(2) ≈ √π/2\n",
"@mark integrate_monomial(3) ≈ 0\n",
"@mark integrate_monomial(4) ≈ 3√π/4\n",
"@mark integrate_monomial(5) == 0"
]
},
{
"cell_type": "markdown",
"id": "af088001",
"metadata": {},
"source": [
"3. Write a function `get_nodes_and_weights(n)` which returns the nodes and weights of the Gauss-Hermite quadrature with $n$ nodes.\n",
"\n",
" **Hints:**\n",
" - Remember that the nodes of the quadrature are given by the roots of $H_{n}$.\n",
" - In order to construct Lagrange polynomials, you may find it useful to use the `fromroots` function."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0fc23392",
"metadata": {},
"outputs": [],
"source": [
"function get_nodes_and_weights(n)\n",
" # Write your code here\n",
" return nodes, weights\n",
"end\n",
"\n",
"@mark length(get_nodes_and_weights(1)) == 2\n",
"@mark length(get_nodes_and_weights(3)[1]) == 3\n",
"@mark length(get_nodes_and_weights(5)[2]) == 5"
]
},
{
"cell_type": "markdown",
"id": "48eb41cf",
"metadata": {},
"source": [
"4. Write a function `integrate_hermite(f, n)`, which returns an approximation of $I[f]$ obtained by Gauss-Hermite integration with $n$ nodes."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dfdc0f2d",
"metadata": {},
"outputs": [],
"source": [
"function integrate_hermite(f, n)\n",
" # Write your code here\n",
"end\n",
"\n",
"@mark integrate_hermite(x -> 1, 10) ≈ integrate_monomial(0)\n",
"@mark abs(integrate_hermite(x -> x, 10)) <= 1e-10\n",
"@mark integrate_hermite(x -> x^2, 10) ≈ integrate_monomial(2)\n",
"@mark integrate_hermite(x -> x^4, 10) ≈ integrate_monomial(4)\n",
"@mark integrate_hermite(x -> x^6, 10) ≈ integrate_monomial(6)\n",
"@mark abs(integrate_hermite(sin, 10)) <= 1e-10"
]
},
{
"cell_type": "markdown",
"id": "308f1595",
"metadata": {},
"source": [
"5. Set $n = 4$, and calculate numerically the degree of precision of the corresponding quadrature rule."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c482007d",
"metadata": {},
"outputs": [],
"source": [
"# Write your code here"
]
},
{
"cell_type": "markdown",
"id": "efd8b98a",
"metadata": {},
"source": [
"6. Set $f(x) = \\cos(x)$, and plot the integration error as a function of $n$.\n",
" You may use that\n",
" $$\n",
" \\int_{-\\infty}^{\\infty} \\cos(x) \\, e^{-x^2} \\, d x = \\frac{\\sqrt{\\pi}}{\\sqrt[4]{e}}\n",
" $$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a964494b",
"metadata": {},
"outputs": [],
"source": [
"# Write your code here"
]
},
{
"cell_type": "markdown",
"id": "50f6b33d",
"metadata": {},
"source": [
"### 4. Exercise on linear systems\n",
"\n",
"This exercise focuses on solving the Euler-Bernoulli beam equation in one dimension,\n",
"with homogeneous Dirichlet boundary conditions:\n",
"$$\n",
"u''''(x) = 1, \\qquad u(0) = u'(0) = u'(1) = u(1) = 0.\n",
"$$\n",
"This equation models the deflection of a beam clamped at both ends,\n",
"under a uniform load.\n",
"A discretization of this equation on a uniform grid using the finite difference method leads to the following linear system:\n",
"$$\n",
"\\begin{pmatrix}\n",
" 6 & -4 & 1 \\\\\n",
" -4 & 6 & -4 & 1 \\\\\n",
" 1 & -4 & 6 & -4 & 1 \\\\\n",
" & 1 & -4 & 6 & -4 & 1 \\\\\n",
" & & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots \\\\\n",
" & & & 1 & -4 & 6 & -4 & 1 \\\\\n",
" & & & & 1 & -4 & 6 & -4 \\\\\n",
" & & & & & 1 & -4 & 6 \\\\\n",
"\\end{pmatrix}\n",
"\\begin{pmatrix}\n",
" u_2 \\\\\n",
" u_3 \\\\\n",
" u_4 \\\\\n",
" u_5 \\\\\n",
" \\vdots \\\\\n",
" u_{n-4} \\\\\n",
" u_{n-3} \\\\\n",
" u_{n-2}\n",
"\\end{pmatrix}\n",
"=\n",
"h^4\n",
"\\begin{pmatrix}\n",
" 1 \\\\\n",
" 1 \\\\\n",
" 1 \\\\\n",
" 1 \\\\\n",
" \\vdots \\\\\n",
" 1 \\\\\n",
" 1 \\\\\n",
" 1\n",
"\\end{pmatrix},\n",
"\\quad\n",
"h := \\frac{1}{n},\n",
"\\quad\n",
"x_i := ih.\n",
"$$\n",
"where $h$ is the spacing between the discretization points, and $(u_2, u_3, \\cdots, u_{n-3}, u_{n-2})$ are the values of the unknown function $u$ at the points $(x_2, x_3, \\cdots, x_{n-3}, x_{n-2})$.\n",
"\n",
"1. Write a function `build_matrix(n)`, which returns the matrix of the linear system.\n",
"\n",
" **Hint:** the function `diagm` is useful here."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0e66494d",
"metadata": {},
"outputs": [],
"source": [
"function build_matrix(n)\n",
" # Write your code here\n",
"end\n",
"\n",
"@mark size(build_matrix(20)) == (17, 17)\n",
"@mark build_matrix(20)[5, 3] ≈ 1\n",
"@mark build_matrix(20)[5, 4] ≈ -4\n",
"@mark build_matrix(20)[5, 5] ≈ 6\n",
"@mark build_matrix(20)[5, 6] ≈ -4\n",
"@mark build_matrix(20)[5, 7] ≈ 1\n",
"@mark build_matrix(20)[5, 8] ≈ 0"
]
},
{
"cell_type": "markdown",
"id": "c67b4993",
"metadata": {},
"source": [
"2. Write a function `build_rhs(n)`, which returns the right-hand side of the linear system."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2bf19861",
"metadata": {},
"outputs": [],
"source": [
"function build_rhs(n)\n",
" # Write your code here\n",
"end\n",
"\n",
"@mark length(build_rhs(20)) == 17\n",
"@mark build_rhs(20)[end] == 1 / 20^4\n",
"@mark build_rhs(20)[1] == 1 / 20^4"
]
},
{
"cell_type": "markdown",
"id": "e5de6676",
"metadata": {},
"source": [
"3. Write a function `forward_substitution(L, y)`\n",
" that solves the linear system $\\mathsf L \\mathbf x = \\mathbf y$,\n",
" in the case where `L` is a lower-triangular matrix,\n",
" by using forward substitution."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a789548e",
"metadata": {},
"outputs": [],
"source": [
"function forward_substitution(L, y)\n",
" # Write your code here\n",
"end\n",
"\n",
"@mark begin n = 10; L = [1 0; 2 1]; forward_substitution(L, [1; 0]) ≈ [1; -2] end\n",
"@mark begin n = 10; L = LowerTriangular(ones(n, n)); y = fill(2, n); forward_substitution(L, L*y) ≈ y end\n",
"@mark begin n = 10; L = LowerTriangular(randn(n, n)); y = randn(n); forward_substitution(L, L*y) ≈ y end\n",
"@mark begin n = 20; L = LowerTriangular(2ones(n, n)); y = rand(n); forward_substitution(L, L*y) ≈ y end"
]
},
{
"cell_type": "markdown",
"id": "10eeaf06",
"metadata": {},
"source": [
"4. The successive over-relaxation method is a splitting method for solving linear systems of the form $\\mathsf A \\mathbf x = \\mathbf b$.\n",
" It is based on the iteration\n",
" $$\n",
" \\mathsf M \\mathbf x_{k+1} = \\mathsf N \\mathbf x_{k} + \\mathbf{b}, \\qquad \\text{where} \\quad \\mathsf M = \\frac{\\mathsf D}{\\omega} + \\mathsf L \\quad \\text{and} \\quad \\mathsf N = \\mathsf M - \\mathsf A.\n",
" \\tag{SOR}\n",
" $$\n",
" \n",
" Here $\\omega \\in (0, 2)$ is a parameter,\n",
" $\\mathsf D$ is diagonal matrix containing the diagonal of $\\mathsf A$,\n",
" and $\\mathsf L$ is a strictly lower triangular matrix containing the strict lower triangular part of $\\mathsf A$,\n",
" not including the diagonal.\n",
" Write a function `sor(A, b, ω)` that\n",
" implements this iteration, without using Julia's `\\` and `inv` functions.\n",
" Initialize the iteration with $\\mathbf x_0 = \\mathbf 0$,\n",
" and use as stopping criterion that\n",
" $$\n",
" \\lVert \\mathsf A \\mathbf x - \\mathbf b \\rVert\n",
" \\leq \\varepsilon \\lVert \\mathbf b \\rVert,\n",
" \\qquad \\varepsilon := 10^{-10}.\n",
" $$\n",
" If after `maxiter` iterations,\n",
" a solution that satisfies this stopping criterion has not been found,\n",
" return `nothing`.\n",
"\n",
" **Hints**:\n",
" - The functions `Diagonal` and `LowerTriangular` are useful to construct the matrices $\\mathsf D$ and $\\mathsf L$.\n",
" - In order to solve (SOR) at each iteration,\n",
" use the function `forward_substitution` you wrote above."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "735e42a1",
"metadata": {},
"outputs": [],
"source": [
"function sor(A, b; ω = 1, ε = 1e-10, maxiter = 10^5)\n",
" # Write your code here\n",
"end\n",
"\n",
"# Test code\n",
"n = 20\n",
"X = qr(randn(n, n)).Q\n",
"A = X*Diagonal(rand(n) .+ 5)*X'\n",
"b = randn(n)\n",
"@mark begin ε = 1e-10; sol = sor(A, b; ω = 1.5, ε = ε); norm(A*sol - b) ≤ ε*norm(b) end\n",
"@mark begin ε = 1e-10; sol = sor(A, b; ω = 1.5, ε = ε); norm(A*sol - b) ≥ 1e-15*norm(b) end\n",
"@mark begin ε = 1e-12; sol = sor(A, b; ω = 1.5, ε = ε); norm(A*sol - b) ≤ ε*norm(b) end\n",
"@mark begin ε = 1e-12; sor(A, b; ω = 2, ε = ε) == nothing end\n",
"@mark begin ε = 1e-12; sor(A, b; ω = 1, ε = ε) ≈ A\\b end"
]
},
{
"cell_type": "markdown",
"id": "43e63cfc",
"metadata": {},
"source": [
"5. Use the relaxation method implemented in the previous item\n",
" with parameters $\\omega = 1.5$ and $\\varepsilon = 10^{-8}$ to solve the linear system of this exercise with $n = 40$.\n",
" Compare on a graph the solution you obatined with the exact solution,\n",
" which is given by\n",
" $$ u(x) = \\frac{1}{24} x^2(x - 1)^2. $$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1f5ba663",
"metadata": {
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"u(x) = 1/24 * x^2 * (x - 1)^2\n",
"# Write your code here"
]
},
{
"cell_type": "markdown",
"id": "badd118f",
"metadata": {},
"source": [
"### 5. Exercise on nonlinear equations\n",
"\n",
"We wish to find all the solutions to the following transcendental equation for $x \\in [-5, -5]$.\n",
"$$\n",
"x = 5 \\sin(\\pi x)\n",
"\\tag{1}\n",
"$$\n",
"\n",
"1. Plot the functions $x \\mapsto x$ and $x \\mapsto 5 \\sin(\\pi x)$ on the same graph,\n",
" for $x$ in the range $[-5, 5]$,\n",
" and count the number of solutions"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a8060fe7",
"metadata": {},
"outputs": [],
"source": [
"# Write your code here"
]
},
{
"cell_type": "markdown",
"id": "7803e926",
"metadata": {},
"source": [
"2. Using the bisection method, calculate precisely the only solution in the interval $[\\frac{1}{2}, 1]$."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "969f0622",
"metadata": {},
"outputs": [],
"source": [
"function bisection(f, a, b, ε = 1e-10)\n",
" # Write your code here\n",
"end\n",
"\n",
"# Calculate solution here"
]
},
{
"cell_type": "markdown",
"id": "6e92a82a",
"metadata": {},
"source": [
"3. Write a function `newton_raphson(f, x₀; maxiter = 100, ε = 1e-12)` that returns the result of the Newton-Raphson iteration applied to the equation $f(x) = 0$ and initialized at `x₀`,\n",
" or `nothing` if a solution is not found after `maxiter` iterations.\n",
" Use the `ForwardDiff` library to compute derivatives,\n",
" and use the following stopping criterion\n",
" $$\n",
" |f(x_k)| ≤ \\varepsilon.\n",
" $$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a2f62f9d",
"metadata": {},
"outputs": [],
"source": [
"function newton_raphson(f, x₀; maxiter = 100, ε = 1e-12)\n",
" # Write your code here\n",
"end\n",
"\n",
"@mark newton_raphson(x -> x^2 - 2, 1) ≈ √2\n",
"@mark newton_raphson(x -> x^2 - 2, -1) ≈ -√2\n",
"@mark newton_raphson(x -> x^3 - 2, 1) ≈ cbrt(2)\n",
"@mark newton_raphson(x -> x^2 + 1, 1.5) == nothing\n",
"@mark newton_raphson(x -> x^2 + 1, 1) == nothing"
]
},
{
"cell_type": "markdown",
"id": "2219fc5f",
"metadata": {},
"source": [
"4. Using the Newton-Raphson method you implemented,\n",
" calculate all the solutions to the transcendental equation (1)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ce44d968",
"metadata": {},
"outputs": [],
"source": [
"# Write your code here."
]
},
{
"cell_type": "markdown",
"id": "d4bcb668",
"metadata": {},
"source": [
"5. We now consider another approach,\n",
" called the *secant method*.\n",
" An iteration of this method is given by\n",
" $$\n",
" x_{k+2} = \\frac{x_{k} f(x_{k+1}) - x_{k+1}f(x_k)}{f(x_{k+1}) - f(x_{k})}.\n",
" $$\n",
" Write a function `secant(f, x₀, x₁; maxiter = 100, ε = 1e-12)`\n",
" that returns the result of the secant iteration applied to the equation $f(x) = 0$, and initialized with `x₀` and `x₁`,\n",
" or `nothing` if a solution is not found after `maxiter` iterations.\n",
" Use the same stopping criterion as above."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8819b0c1",
"metadata": {},
"outputs": [],
"source": [
"function secant(f, x₀, x₁; maxiter = 100, ε = 1e-12)\n",
" # Write your code here\n",
"end\n",
"\n",
"@mark secant(x -> x^2 - 2, 1., 2.) ≈ √2\n",
"@mark secant(x -> x^2 - 2, -1., -2.) ≈ -√2\n",
"@mark secant(x -> x^2 + 1, 1., 2.) == nothing"
]
},
{
"cell_type": "markdown",
"id": "48435545",
"metadata": {},
"source": [
"6. Use the secant method you implemented to solve the transcendental equation\n",
"$$x = e^{-x}$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bb25e735",
"metadata": {},
"outputs": [],
"source": [
"# Write your code here."
]
},
{
"cell_type": "markdown",
"id": "f6721f3b",
"metadata": {},
"source": [
"### 6. Exercise on eigenvalue problems"
]
},
{
"cell_type": "markdown",
"id": "cb8d2606",
"metadata": {},
"source": [
"Our goal in this exercise is to implement,\n",
"without using the functions `eigen` and `svd` from the `LinearAlgebra` library,\n",
"a simple algorithm for calculating the dominant *singular values* and associated *singular vectors* of a matrix,\n",
"and to use this algorithm for image compression.\n",
"Recall that any matrix $\\mathsf A \\in \\mathbf R^{n \\times n}$ admits a *singular value decomposition* of the form\n",
"$$\n",
"\\mathsf A = \\mathsf U \\mathsf \\Sigma \\mathsf V^\\top,\n",
"$$\n",
"where $\\mathsf U \\in \\mathbf R^{n \\times n}$ and $\\mathsf V \\in \\mathbf R^{n \\times n}$ are orthogonol matrices,\n",
"and $\\mathsf \\Sigma \\in \\mathbf R^{n \\times n}$ is a diagonal matrix.\n",
"The diagonal entries of $\\mathsf \\Sigma$ are called *singular value*,\n",
"while the columns of $\\mathsf U$ and $\\mathsf V$ are the associated left and right *singular vectors*, respectively.\n",
"We restrict our attention to square matrices, for simplicity.\n",
"\n",
"1. Write a function `myEigen(B, p, niter)` which returns the `p` dominant eigenvalues,\n",
" together with the associated eigenvectors,\n",
" of a symmetric matrix `B`.\n",
" To this end, implement the simultaneous iteration seen in class.\n",
" You can use the `myQR` function given below in order to calculate reduced QR decompositions."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c1dce274",
"metadata": {},
"outputs": [],
"source": [
"function myQR(B)\n",
" Q, R = qr(B)\n",
" return Matrix(Q), R\n",
"end\n",
"\n",
"function myEigen(B, p, niter)\n",
" # Write your code here\n",
"end\n",
"\n",
"@mark begin n = 2; A = randn(n, n); myEigen(A'A, n, 100)[1] ≈ reverse(eigen(A'A).values) end\n",
"@mark begin n = 4; A = randn(n, n); myEigen(A'A, n, 100)[1] ≈ reverse(eigen(A'A).values) end\n",
"@mark begin A = randn(5, 5); q, r = qr(A); B = q*Diagonal(1:5)*q'; myEigen(B, 3, 100)[1] ≈ [5; 4; 3] end"
]
},
{
"cell_type": "markdown",
"id": "6ec0399d",
"metadata": {},
"source": [
"2. Write a function `mySVD(B, p, niter)`\n",
" that returns the `p` dominant singular values of square matrix `B` (in a vector `σs`),\n",
" together with the associated left and right singular vectors (in matrices `Up` and `Vp`).\n",
" To this end,\n",
" notice that\n",
" $$\n",
" \\mathsf A \\mathsf A^\\top = \\mathsf U \\mathsf \\Sigma^2 \\mathsf U^\\top, \\qquad\n",
" \\mathsf A^\\top \\mathsf A = \\mathsf V \\mathsf \\Sigma^2 \\mathsf V^\\top.\n",
" $$\n",
" Therefore, the left singular vectors of $\\mathsf A$\n",
" are the eigenvectors of $\\mathsf A \\mathsf A^\\top$,\n",
" while the right singular vectors of $\\mathsf A$\n",
" are the eigenvectors of $\\mathsf A^\\top \\mathsf A$,\n",
"\n",
" **Hints**:\n",
"\n",
" - In order to calculate singular vectors,\n",
" apply the function `myEigen` you wrote above to the matrices $\\mathsf A \\mathsf A^\\top$\n",
" (in order to obtain left singular vectors),\n",
" and $\\mathsf A^\\top \\mathsf A$ (in order to obtain right singular vectors).\n",
"\n",
" - Once you have calculated the left and right singular vectors\n",
" associated with the `p` dominant singular values,\n",
" the singular values themselves can be obtained by extracting the diagonal from the matrix\n",
" $$\n",
" \\Sigma_p = \\mathsf U_p^\\top \\mathsf B \\mathsf V_p.\n",
" $$\n",
" Here $\\mathsf U_p$ and $\\mathsf V_p$ are the matrices containing as columns the left and right singular vectors associated with the `p` dominant singular values,\n",
" respectively."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "991a28c2",
"metadata": {},
"outputs": [],
"source": [
"function mySVD(B, p, niter)\n",
" # Write your code here\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "1dfe8108",
"metadata": {},
"source": [
"3. Write a function to test your implementation of `mySVD`.\n",
" Proceed as follows:\n",
" - Create a 10 x 10 matrix `B` with random entries using the `randn` function.\n",
" - Calculate its full singular value decomposition using `mySVD`.\n",
" - Check that `U` and `V` are orthogonal,\n",
" and that `U*Diagonal(σs)*V'` is close to `B`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "81b070a6",
"metadata": {},
"outputs": [],
"source": [
"# Write you code here"
]
},
{
"cell_type": "markdown",
"id": "e0113e32",
"metadata": {},
"source": [
"4. The singular value decomposition is very useful for compressing matrices.\n",
" The idea of matrix compression based on SVD is the following:\n",
" given $p \\leqslant n$, a matrix $B \\in \\mathbf R^{n\\times n}$\n",
" can be approximated by\n",
" $$\n",
" \\widetilde {\\mathsf B} := \\mathsf U_p \\mathsf \\Sigma_p \\mathsf V_p^\\top,\n",
" $$\n",
" where $\\Sigma_p \\in \\mathbf R^{p \\times p}$ is a diagonal matrix containing the $p$ dominant singular values of $\\mathsf B$ on its diagonal,\n",
" and where $\\mathsf U_p \\in \\mathbf R^{n \\times p}$ and $\\mathsf V_p \\in \\mathbf R^{n \\times p}$ are rectangular matrices containing the associated left and right singular vectors,\n",
" respectively.\n",
"\n",
" Since a grayscale image can be represented by a matrix containing the intensity values of all the pixels,\n",
" this approach for compressing matrices can be used for compressing grayscale images.\n",
" Use this method, i.e. calculate $\\widetilde {\\mathsf B}$, for $p \\in \\{5, 10, 20, 30\\}$,\n",
" in order to compress the image `fabio_gray_256` given below,\n",
" and plot the compressed image for these values of `p`."
]
},
{
"cell_type": "markdown",
"id": "01a449fc",
"metadata": {},
"source": [
" **Remarks**:\n",
" - (For information only) In practice, instead of storing the full matrix $\\widetilde {\\mathsf B}$, which contains $n^2$ entries,\n",
" we can store only the matrices $\\mathsf U_p$, $\\mathsf \\Sigma_p$ and $\\mathsf V_p$,\n",
" which together contain only $2np + p^2$ entries.\n",
" If $p \\ll n$,\n",
" then the memory required to store these matrices is much smaller than the memory required to store the initial matrix $\\mathsf B$.\n",
"\n",
" - A function for drawing images based on the matrix of pixel intensity values is provided below."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "612d3d33",
"metadata": {},
"outputs": [],
"source": [
"# Do not change the code in this cell\n",
"A = testimage(\"fabio_gray_256\")\n",
"\n",
"# Convert image to matrix of Float64\n",
"M = Float64.(A)\n",
"\n",
"# Function to plot a grayscale image from the matrix\n",
"# containing the intensity values of all the pixels\n",
"function plot_matrix(B)\n",
" plot(Gray.(B), ticks=false, showaxis=false)\n",
"end\n",
"\n",
"plot_matrix(M)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "56bda9b6",
"metadata": {},
"outputs": [],
"source": [
"p = 5\n",
"# Write you code here"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5801459e",
"metadata": {},
"outputs": [],
"source": [
"p = 10\n",
"# Write you code here"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b028396c",
"metadata": {},
"outputs": [],
"source": [
"p = 20\n",
"# Write you code here"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a3d591f0",
"metadata": {},
"outputs": [],
"source": [
"p = 30\n",
"# Write you code here"
]
}
],
"metadata": {
"jupytext": {
"encoding": "# -*- coding: utf-8 -*-",
"main_language": "julia"
}
},
"nbformat": 4,
"nbformat_minor": 5
}