exams/computerexam2122s.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MATH50003 Numerical Analysis (2021–2022) Computer-based Exam\n",
    "\n",
    "Instructions for uploading and downloading:\n",
    "\n",
    "1. Rename the file to include your CID.\n",
    "2. You have 30 mins to download the exam beginning at 3pm on 18 March.\n",
    "2. You have 1 hour to complete the exam beginning at 3:30pm on 18 March.\n",
    "3. Deadline is 5pm on 18 March to upload the completed Jupyter notebook (`.ipynb`) to Blackboard. \n",
    "Late submissions received before 7pm will be capped at 40%.\n",
    "5. Once uploaded, re-download the file before the final submission time to confirm it is correct.\n",
    "You are allowed to upload additional submissions but only last valid upload before 5pm will be used.\n",
    "6. If uploading via Blackboard fails you may e-mail the UG Office: maths.exams@imperial.ac.uk\n",
    "\n",
    "Instructions for the exam:\n",
    "\n",
    "1. For each problem, replace the `# TODO` to complete the question.\n",
    "The unit tests are provided to help you test your answers.\n",
    "3. Problems are marked A/B/C to indicate difficulty (\"A\" being most difficult).\n",
    "Partial credit will be awarded for reasonable attempts even if the tests\n",
    "are not passed. A and B questions are worth 12 marks while C questions are worth 10 marks.\n",
    "3. If you have technical queries please email s.olver@imperial.ac.uk. Any other queries\n",
    "should be sent to the UG Office: maths.exams@imperial.ac.uk\n",
    "4. You may use existing code from anywhere\n",
    "but you are **REQUIRED** to cite the source if it is not part of the module material,\n",
    "ideally by including a weblink in a comment. \n",
    "5. You **MUST NOT** ask for help online or\n",
    "communicate with others within or outside the module.\n",
    "Failure to follow these rules will be considered misconduct.\n",
    "\n",
    "\n",
    "\n",
    "You should use the following packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:30.955108Z",
     "iopub.status.busy": "2022-04-06T10:51:30.518176Z",
     "iopub.status.idle": "2022-04-06T10:51:31.189263Z",
     "shell.execute_reply": "2022-04-06T10:51:31.188740Z"
    }
   },
   "outputs": [],
   "source": [
    "using LinearAlgebra, Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Note `SetRounding` is not necessary.)\n",
    "\n",
    "**WARNING** It may be necessary to restart the kernel if issues arise. Remember to reload the packages\n",
    "when you do so.\n",
    "\n",
    "## 1. Numbers\n",
    "\n",
    "**Problem 1.1 (C)** \n",
    "Implement the function `issub` that determines whether a `Float16` is a sub-normal number.\n",
    "DO NOT use the inbuilt routine `issubnormal`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:31.218175Z",
     "iopub.status.busy": "2022-04-06T10:51:31.217489Z",
     "iopub.status.idle": "2022-04-06T10:51:31.715015Z",
     "shell.execute_reply": "2022-04-06T10:51:31.714611Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "issub (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function issub(x::Float16)\n",
    "    # TODO: return `true` if `x` is a sub-normal float. Otherwise return `false`\n",
    "    ## SOLUTION\n",
    "    σ = 15\n",
    "    abs(x) < 2.0^(1-σ)\n",
    "    ## END\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:31.717033Z",
     "iopub.status.busy": "2022-04-06T10:51:31.716528Z",
     "iopub.status.idle": "2022-04-06T10:51:32.301958Z",
     "shell.execute_reply": "2022-04-06T10:51:32.301562Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m\n",
       "  Expression: issub(reinterpret(Float16, 0x03ff))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@test issub(Float16(0))\n",
    "@test issub(nextfloat(Float16(0)))\n",
    "@test issub(prevfloat(Float16(0)))\n",
    "@test !issub(Float16(1))\n",
    "@test !issub(reinterpret(Float16,0b0000010000000000))\n",
    "@test issub(reinterpret(Float16,0b0000001111111111))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Differentiation\n",
    "\n",
    "**Problem 2.1 (C)** Use second-order finite differences\n",
    "with an appropriately chosen $h$ to approximate the second derivative of\n",
    "$$\n",
    "f(x) = \\cos(x^2)\n",
    "$$\n",
    "at $x = 0.1$ to 5 digits accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:32.304073Z",
     "iopub.status.busy": "2022-04-06T10:51:32.303526Z",
     "iopub.status.idle": "2022-04-06T10:51:32.406306Z",
     "shell.execute_reply": "2022-04-06T10:51:32.405774Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fd2 (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function fd2(x)\n",
    "    # TODO: implement a second-order finite-difference rule \n",
    "    # to approximate f''(x)\n",
    "    # for f(x) = cos(x^2)\n",
    "    # with step-size h chosen to get sufficient accuracy\n",
    "    ## SOLUTION\n",
    "    h = cbrt(eps())\n",
    "    f = x -> cos(x^2)\n",
    "    (f(x + h) - 2f(x) + f(x - h)) / h^2\n",
    "    ## END\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:32.408273Z",
     "iopub.status.busy": "2022-04-06T10:51:32.407758Z",
     "iopub.status.idle": "2022-04-06T10:51:32.596051Z",
     "shell.execute_reply": "2022-04-06T10:51:32.594639Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m\n",
       "  Expression: abs(fd2(0.1) + 2 * sin(0.1 ^ 2) + 4 * 0.1 ^ 2 * cos(0.1 ^ 2)) ≤ 1.0e-5\n",
       "   Evaluated: 2.2397068676072163e-7 ≤ 1.0e-5"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@test abs(fd2(0.1) + 2*sin(0.1^2) + 4*0.1^2*cos(0.1^2)) ≤ 1E-5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem 2.2 (A)** Consider a 2nd order version of a dual number:\n",
    "$$\n",
    "a + b ϵ_1 + c ϵ_2\n",
    "$$\n",
    "such that\n",
    "$$\n",
    "\\begin{align*}\n",
    "ϵ_1^2 &= ϵ_2 \\\\\n",
    "ϵ_2^2 &= ϵ_1 ϵ_2 =  0.\n",
    "\\end{align*}\n",
    "$$\n",
    "Complete the following implementation supporting `+` and `*` (and\n",
    "assuming `a,b,c` are `Float64`). Hint: you may need to work out on paper\n",
    "how to multiply `(s.a + s.b ϵ_1 + s.c ϵ_2)*(t.a + t.b ϵ_1 + t.c ϵ_2)` using the\n",
    "relationship above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:32.598655Z",
     "iopub.status.busy": "2022-04-06T10:51:32.597940Z",
     "iopub.status.idle": "2022-04-06T10:51:32.881464Z",
     "shell.execute_reply": "2022-04-06T10:51:32.881014Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "* (generic function with 406 methods)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import Base: *, +, ^\n",
    "struct Dual2\n",
    "    a::Float64\n",
    "    b::Float64\n",
    "    c::Float64\n",
    "end\n",
    "\n",
    "function +(s::Dual2, t::Dual2)\n",
    "    ## TODO: Implement Dual2(...) + Dual2(...), returning a Dual2\n",
    "    ## SOLUTION\n",
    "    Dual2(s.a + t.a, s.b + t.b, s.c + t.c)\n",
    "    ## END\n",
    "end\n",
    "\n",
    "function +(s::Dual2, c::Real)\n",
    "    ## TODO: Implement Dual2(...) + c, returning a Dual2\n",
    "    ## SOLUTION\n",
    "    Dual2(s.a + c, s.b, s.c)\n",
    "    ## END\n",
    "end\n",
    "\n",
    "function *(c::Number, s::Dual2)\n",
    "    ## TODO: Implement c * Dual2(...), returning a Dual2\n",
    "    ## SOLUTION\n",
    "    Dual2(c * s.a, c * s.b, c * s.c)\n",
    "    ## END\n",
    "end\n",
    "\n",
    "function *(s::Dual2, t::Dual2)\n",
    "    ## TODO: Implement Dual2(...) * Dual2(...), returning a Dual2\n",
    "\n",
    "    ## SOLUTION\n",
    "    # we deduce (s.a + s.b*ϵ_1 + s.c*ϵ_2)*(t.a + t.b*ϵ_1 + t.c*ϵ_2) == \n",
    "    # s.a * t.a + (s.a*t.b + s.b*t.a)*ϵ_1 + (s.a*t.c + s.b*t.b + s.c*t.a)*ϵ_2\n",
    "    Dual2(s.a * t.a, s.a * t.b + s.b * t.a, s.a * t.c + s.b * t.b + s.c * t.a)\n",
    "    ## END\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:32.883524Z",
     "iopub.status.busy": "2022-04-06T10:51:32.883025Z",
     "iopub.status.idle": "2022-04-06T10:51:32.936487Z",
     "shell.execute_reply": "2022-04-06T10:51:32.936109Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m\n",
       "  Expression: f(Dual2(x, 1.0, 0.0)) == Dual2(f(x), 3 * x ^ 2 + 2, (6x) / 2)\n",
       "   Evaluated: Dual2(1.201, 2.03, 0.30000000000000004) == Dual2(1.201, 2.03, 0.30000000000000004)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = x -> x*x*x + 2x + 1\n",
    "x = 0.1\n",
    "@test f(Dual2(x,1.,0.)) == Dual2(f(x), 3x^2+2, 6x / 2)\n",
    "\n",
    "# This has computed the first and second derivatives as\n",
    "# as f(x) + f'(x)*ϵ_1 + f''(x)/2*ϵ_2\n",
    "# == (x^3 + x) + (3x^2+1)*ϵ_1 + 6x/2*ϵ_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Structured Matrices\n",
    "\n",
    "**Problem 3.1 (C)** Complete the implementation of `LowerTridiagonal` which represents a banded matrix with\n",
    "bandwidths $(l,u) = (2,0)$ by storing only its diagonal, sub-diagonal, and second-sub-diagonal as vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:32.938372Z",
     "iopub.status.busy": "2022-04-06T10:51:32.937868Z",
     "iopub.status.idle": "2022-04-06T10:51:33.444751Z",
     "shell.execute_reply": "2022-04-06T10:51:33.444332Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m\n",
       "  Expression: LowerTridiagonal(d, dl, dl2) == diagm(0 => d, -1 => dl, -2 => dl2)\n",
       "   Evaluated: [-1.3060650083779026 0.0 … 0.0 0.0; -1.7190655458441395 -0.32549142470133213 … 0.0 0.0; … ; 0.0 0.0 … 0.7486471167904769 0.0; 0.0 0.0 … 0.11458311309250767 -1.6333698695008252] == [-1.3060650083779026 0.0 … 0.0 0.0; -1.7190655458441395 -0.32549142470133213 … 0.0 0.0; … ; 0.0 0.0 … 0.7486471167904769 0.0; 0.0 0.0 … 0.11458311309250767 -1.6333698695008252]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import Base: getindex,  size, *\n",
    "\n",
    "struct LowerTridiagonal <: AbstractMatrix{Float64}\n",
    "    d::Vector{Float64}   # diagonal entries of length n\n",
    "    dl::Vector{Float64}  # sub-diagonal entries of length n-1\n",
    "    dl2::Vector{Float64} # second-sub-diagonal entries of length n-2\n",
    "end\n",
    "\n",
    "size(L::LowerTridiagonal) = (length(L.d),length(L.d))\n",
    "\n",
    "function getindex(L::LowerTridiagonal, k::Int, j::Int)\n",
    "    d, dl, dl2 = L.d, L.dl, L.dl2\n",
    "    # TODO: return L[k,j].\n",
    "    # If `k == j` then it should be equal to `d[k]`.\n",
    "    # If `k == j+1` then it should be equal to `dl[j]`.\n",
    "    # If `k == j+2` then it should be equal to `dl2[j]`.\n",
    "    # Otherwise, it should return 0.0\n",
    "    ## SOLUTION\n",
    "    if k == j\n",
    "        d[k]\n",
    "    elseif k == j+1\n",
    "        dl[j]\n",
    "    elseif k == j+2\n",
    "        dl2[j]\n",
    "    else\n",
    "        0.0\n",
    "    end\n",
    "    ## END\n",
    "end\n",
    "\n",
    "n = 10\n",
    "d, dl, dl2 = randn(n), randn(n-1), randn(n-2)\n",
    "@test LowerTridiagonal(d, dl, dl2) == diagm(0 => d, -1 => dl, -2 => dl2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem 3.2 (B)** Complete the implementation of `*` for a `LowerTridiagonal` matrix\n",
    "so that it takes $O(n)$ operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:33.447295Z",
     "iopub.status.busy": "2022-04-06T10:51:33.446792Z",
     "iopub.status.idle": "2022-04-06T10:51:33.755476Z",
     "shell.execute_reply": "2022-04-06T10:51:33.755036Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m\n",
       "  Expression: LowerTridiagonal(d, dl, dl2) * x ≈ diagm(0 => d, -1 => dl, -2 => dl2) * x\n",
       "   Evaluated: [1.0426642047679044, -1.805807111400234, -0.47593613900239584, 0.4063813641392797, 0.32988846510003056, 0.6728956798190676, -0.31723859650332453, 1.067084915281285, 1.178928708164844, 0.9993355798370727] ≈ [1.0426642047679044, -1.805807111400234, -0.4759361390023958, 0.4063813641392797, 0.32988846510003056, 0.6728956798190676, -0.31723859650332453, 1.0670849152812851, 1.178928708164844, 0.9993355798370727]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function *(L::LowerTridiagonal, x::AbstractVector)\n",
    "    ## TODO: Return L*x but computed in O(n) operations\n",
    "    ## SOLUTION\n",
    "    n,m = size(L)\n",
    "    b = zeros(n) # returned vector\n",
    "\n",
    "    for j = 1:n, k = j:min(j+2,n)\n",
    "        b[k] += L[k,j]*x[j]\n",
    "    end\n",
    "    b\n",
    "    ## END\n",
    "end\n",
    "\n",
    "n = 10\n",
    "d, dl, dl2 = randn(n), randn(n-1), randn(n-2)\n",
    "x = randn(n)\n",
    "@test LowerTridiagonal(d, dl, dl2)*x ≈ diagm(0 => d, -1 => dl, -2 => dl2)*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Decompositions\n",
    "\n",
    "**Problem 4.1 (C)** Approximate $\\exp x$ by a cubic polynomial by minimising\n",
    "the least squares error when sampled at $n$ evenly spaced points in $[0,1]$,\n",
    "that is, $x_k = (k-1)/(n-1)$,\n",
    "returning the coefficients in the monomial basis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:33.757421Z",
     "iopub.status.busy": "2022-04-06T10:51:33.756924Z",
     "iopub.status.idle": "2022-04-06T10:51:33.837496Z",
     "shell.execute_reply": "2022-04-06T10:51:33.837061Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "expfit (generic function with 1 method)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function expfit(n)\n",
    "    ## TODO: return the coefficients [c_0,c_1,c_2,c_3] of the polynomial\n",
    "    # c_0 + c_1*x + c_2*x^2 + c_3*x^3 that minimises the L^2 error at `n`\n",
    "    # evenly spaced samples\n",
    "    ## SOLUTION\n",
    "    x = range(0,1; length=n)\n",
    "    V = x .^ (0:3)'\n",
    "    V \\ exp.(x)\n",
    "    ## END\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:33.839432Z",
     "iopub.status.busy": "2022-04-06T10:51:33.838851Z",
     "iopub.status.idle": "2022-04-06T10:51:34.814138Z",
     "shell.execute_reply": "2022-04-06T10:51:34.813615Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m\n",
       "  Expression: abs((c₀ + c₁ * x + c₂ * x ^ 2 + c₃ * x ^ 3) - exp(x)) ≤ 0.001\n",
       "   Evaluated: 0.00021335436497960103 ≤ 0.001"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c₀,c₁,c₂,c₃ = expfit(1000)\n",
    "    x = 0.1\n",
    "    @test abs(c₀ + c₁*x + c₂*x^2 + c₃*x^3 - exp(x)) ≤ 1E-3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Problem 4.2 (A)** Complete the function `lq(A)` that\n",
    "returns a LQ decomposition, that is, `A = LQ` where  `L` is lower triangular and `Q` is an orthogonal\n",
    "matrix. You may assume that `A`\n",
    "is a square `Matrix{Float64}`. Hint: think of how a Householder reflection\n",
    "can be constructed such that, for $𝐱 ∈ ℝ^n$,\n",
    "$$\n",
    "𝐱^⊤ Q = \\|𝐱\\|𝐞_1^⊤.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:34.816027Z",
     "iopub.status.busy": "2022-04-06T10:51:34.815532Z",
     "iopub.status.idle": "2022-04-06T10:51:34.906421Z",
     "shell.execute_reply": "2022-04-06T10:51:34.905955Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "lq (generic function with 1 method)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function lq(A)\n",
    "    m,n = size(A)\n",
    "    m == n || error(\"not square\")\n",
    "    ## TODO Create Q and L such that A = L*Q, Q'Q == I and L is lower triangular\n",
    "    ## SOLUTION\n",
    "    L = copy(A)\n",
    "    Q = Matrix(1.0I, n, n)\n",
    "    for k = 1:n-1\n",
    "        y = L[k, k:end]\n",
    "        y[1] -= norm(y)\n",
    "        w = y / norm(y)\n",
    "        Qₖ = I - 2 * w * w'\n",
    "        L[k:end, k:end] = L[k:end, k:end] * Qₖ\n",
    "        Q[k:end, :] = Qₖ * Q[k:end, :]\n",
    "    end\n",
    "    L,Q\n",
    "    ## END\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:34.908926Z",
     "iopub.status.busy": "2022-04-06T10:51:34.908092Z",
     "iopub.status.idle": "2022-04-06T10:51:36.070385Z",
     "shell.execute_reply": "2022-04-06T10:51:36.069952Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m\n",
       "  Expression: L ≈ tril(L)\n",
       "   Evaluated: [3.7416573867739413 0.0 -1.1102230246251565e-16; 9.621404708847278 2.32992949004287 5.551115123125783e-17; 1.6035674514745462 -0.6131393394849659 0.22941573387056163] ≈ [3.7416573867739413 0.0 0.0; 9.621404708847278 2.32992949004287 0.0; 1.6035674514745462 -0.6131393394849659 0.22941573387056163]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [1.0 2 3; 1 4 9; 1 1 1]\n",
    "L,Q = lq(A)\n",
    "@test Q'Q ≈ I\n",
    "@test L*Q ≈ A\n",
    "@test L ≈ tril(L) # it is acceptable to have small non-zero entries in L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Singular Value Decomposition\n",
    "\n",
    "**Problem 5.1 (B)** Implement `pseudoinv` that returns the pseudo-inverse $A^+$\n",
    "for an arbitrary square matrix, assuming that any singular value less than\n",
    "$10^{-15}$ is in fact exactly zero. DO NOT use the inbuilt routine `pinv`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:36.072321Z",
     "iopub.status.busy": "2022-04-06T10:51:36.071814Z",
     "iopub.status.idle": "2022-04-06T10:51:36.158140Z",
     "shell.execute_reply": "2022-04-06T10:51:36.157764Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "pseudoinv (generic function with 1 method)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function pseudoinv(A)\n",
    "    m,n = size(A)\n",
    "    m == n || error(\"A must be square\")\n",
    "    tol = 1E-15 # threshold below which we assume a singular value is zero\n",
    "    ## TODO: construct and return the pseudo inverse of A\n",
    "    ## SOLUTION\n",
    "    U,σ,V = svd(A)\n",
    "    r = 0\n",
    "    for k = 1:length(σ)\n",
    "        if σ[k] > tol\n",
    "            r += 1\n",
    "        end\n",
    "    end\n",
    "    V[:,1:r] * Diagonal(inv.(σ[1:r])) * U[:,1:r]'\n",
    "    ## END\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:36.160022Z",
     "iopub.status.busy": "2022-04-06T10:51:36.159527Z",
     "iopub.status.idle": "2022-04-06T10:51:37.513823Z",
     "shell.execute_reply": "2022-04-06T10:51:37.513210Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m\n",
       "  Expression: A * A⁺ * A ≈ A\n",
       "   Evaluated: [1.0000000000000018 2.0000000000000004 2.999999999999999; 4.000000000000003 5.000000000000001 6.0; 7.000000000000005 8.000000000000004 9.000000000000004] ≈ [1 2 3; 4 5 6; 7 8 9]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [1 2 3; 4 5 6; 7 8 9]\n",
    "A⁺ = pseudoinv(A)\n",
    "@test A⁺*A*A⁺ ≈ A⁺\n",
    "@test A*A⁺*A ≈ A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Differential Equations\n",
    "\n",
    "**Problem 6.1 (B)** Complete the function `mathieu(n)` that returns a length-$n$ `Vector{Float64}`\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "u_1 \\\\\n",
    "⋮ \\\\\n",
    "u_n\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "such that $u_k$ approximates the solution to the time-evolution equation\n",
    "$$\n",
    "\\begin{align*}\n",
    "u(0) &= 0 \\\\\n",
    "u'(0) &= 1 \\\\\n",
    "u''(t) &= cos(t) u(t)\n",
    "\\end{align*}\n",
    "$$\n",
    "at the point $t_k = (k-1)/(n-1)$ using the Forward Euler method, by first recasting the problem\n",
    "as a vector ODE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:37.516968Z",
     "iopub.status.busy": "2022-04-06T10:51:37.516082Z",
     "iopub.status.idle": "2022-04-06T10:51:37.624504Z",
     "shell.execute_reply": "2022-04-06T10:51:37.624072Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mathieu (generic function with 1 method)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function mathieu(n)\n",
    "    # TODO: return a Vector{Float64} approximating the solution to the ODE\n",
    "    # at n evenly spaced points between 0 and 1.\n",
    "    ## SOLUTION\n",
    "    t = range(0, 1; length = n)\n",
    "    A = t -> [0 1; cos(t) 0]\n",
    "    h = step(t)\n",
    "\n",
    "    U = zeros(2, n) # each column is a time-slice\n",
    "    U[:, 1] = [0.0, 1.0] # initial condition\n",
    "    for k = 1:n-1\n",
    "        U[:, k+1] = (I + h * A(t[k])) * U[:, k]\n",
    "    end\n",
    "\n",
    "    U[1,:]\n",
    "    ## END\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-04-06T10:51:37.626521Z",
     "iopub.status.busy": "2022-04-06T10:51:37.626024Z",
     "iopub.status.idle": "2022-04-06T10:51:37.871877Z",
     "shell.execute_reply": "2022-04-06T10:51:37.871180Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m\n",
       "  Expression: abs(u[end] - 1.148783704310448) ≤ 2.0e-5\n",
       "   Evaluated: 4.409274623640158e-6 ≤ 2.0e-5"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u = mathieu(100_000)\n",
    "@test u isa Vector\n",
    "@test abs(u[1]) ≤ 1E-15\n",
    "# this compares against the exact formula\n",
    "@test abs(u[end] - 1.148783704310448) ≤ 2E-5"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.7.0",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}