{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b18627f2",
   "metadata": {},
   "source": [
    "### 1. Exercise on Interpolation\n",
    "\n",
    "We aim to implement a numerical method to approximately solve the Euler-Bernoulli beam equation with homogeneous Dirichlet boundary conditions:\n",
    "\n",
    "$$\n",
    "  u\\in C^4([0,1]),\\quad\\left\\{\\begin{aligned}  u''''(x) &= \\varphi(x) \\qquad \\forall\\, x\\in(0,1),\\\\\n",
    "  u(0) &= u'(0) = u'(1) = u(1) = 0, \\end{aligned}\\right.\n",
    "$$\n",
    "where $\\varphi(x) = (2\\pi)^4\\cos(2\\pi x)$.\n",
    "In order to solve the equation approximately, we approximate the right-hand side $\\varphi$ by its interpolating polynomial $\\widehat \\varphi$, and then we solve the equation exactly with the right-hand side $\\widehat \\varphi$ instead of $\\varphi$.\n",
    "\n",
    "1. Write a function `fit_values_and_slopes(u₀, up₀, u₁, up₁)` which returns the unique polynomial $p$ of degree 3 such that\n",
    "   $$\n",
    "   p(0) = u_0, \\qquad p'(0) = up_0, \\qquad p(1) = u_1, \\qquad p'(1) = up_1.\n",
    "   $$\n",
    "1. Write a function `approx(n)` implementing the approach described above for solving the PDE. The function should return a polynomial approximation of the solution based on an interpolation of **degree** $n$ of the right-hand side at equidistant points between 0 and 1, inclusive.\n",
    "\n",
    "1. Write a function `estimate_error(n)` that approximates the error,\n",
    "   in $L^\\infty$ norm,\n",
    "   between the exact and approximate solutions.\n",
    "   Note that the exact solution is given by\n",
    "   $$\n",
    "      \\varphi(x) = \\cos(2\\pi x) - 1.\n",
    "   $$\n",
    "\n",
    "1. Plot the error for $n$ in the range $\\llbracket 5, 50 \\rrbracket$. Use a logarithmic scale for the $y$ axis.\n",
    "\n",
    "*Hints:*\n",
    "- The Wikipedia page on *cubic Hermite splines* may be useful for the first exercise.\n",
    "- You can use the `fit` function from the `Polynomials.jl` library to obtain the interpolating polynomial:\n",
    "\n",
    "    ```julia\n",
    "    p = fit(x, y)\n",
    "    ```\n",
    "\n",
    "    where `x` are the interpolation nodes, and `y` are the values of the function to interpolate.\n",
    "\n",
    "- To calculate the analytical solution with a polynomial right-hand side, notice that all solutions are polynomials, and without boundary conditions, the solution is unique modulo $\\mathbf{P}_3$.\n",
    "\n",
    "- You can use the `integrate` function from the `Polynomials.jl` library, which calculates an antiderivative of a polynomial:\n",
    "\n",
    "    ```julia\n",
    "    P = integrate(p)\n",
    "    ```\n",
    "\n",
    "- Use the `BigFloat` format to limit rounding errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfd34e75",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LaTeXStrings\n",
    "using Polynomials\n",
    "using Plots\n",
    "\n",
    "function fit_values_and_slopes(u₀, up₀, u₁, up₁)\n",
    "    # We look for polynomials p(x) = a₀ + a₁ x + a₂ x² + a₃ x³\n",
    "    A = [1 0 0 0; 0 1 0 0; 1 1 1 1; 0 1 2 3]\n",
    "    α = A\\[u₀; up₀; u₁; up₁]\n",
    "    return Polynomial(α)\n",
    "end\n",
    "\n",
    "# Test our code\n",
    "p = fit_values_and_slopes(-1, -1, 1, 1)\n",
    "plot(p, xlims=(0, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0455862",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "φ(x) = (2π)^4 * cospi(2*x)\n",
    "u(x) = cospi(2*x) - 1\n",
    "\n",
    "function approx(n)\n",
    "    X = LinRange{BigFloat}(0, 1, n + 1)\n",
    "    Y = φ.(X)\n",
    "    p = fit(X, Y)\n",
    "    uh = integrate(integrate(integrate(integrate(p))))\n",
    "    ∂uh = derivative(uh)\n",
    "    uh -= fit_values_and_slopes(uh(0), ∂uh(0), uh(1), ∂uh(1))\n",
    "    return uh\n",
    "end\n",
    "\n",
    "plot(approx(3), xlims=(0, 1), label=\"Exact solution\")\n",
    "plot!(u, xlims=(0, 1), label=\"Approximate solution\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fa2f4fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "function estimate_error(n)\n",
    "    un = approx(n)\n",
    "    x_fine = LinRange(0, 1, 1000)\n",
    "    un_fine, u_fine = un.(x_fine), u.(x_fine)\n",
    "    return maximum(abs.(u_fine - un_fine))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d7bde23",
   "metadata": {},
   "outputs": [],
   "source": [
    "ns = 5:50\n",
    "errors = estimate_error.(ns)\n",
    "plot(ns, errors, marker = :circle, label=L\"$L^{\\infty}$ Error\")\n",
    "plot!(yaxis=:log, lw=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a4a0325",
   "metadata": {},
   "source": [
    "### 2. Exercise on numerical integration\n",
    "\n",
    "Our goal in this exercise is to write a program in order to calculate integrals of the form\n",
    "$$\n",
    "I[f] := \\int_0^{\\infty} f(x) e^{-x} \\, d x\n",
    "$$\n",
    "To this end, we will use Laguerre polynomials,\n",
    "which are orthogonal polynomials for the following inner product:\n",
    "   $$\n",
    "    \\langle f, g \\rangle := \\int_0^{\\infty} f(x) g(x) e^{-x} \\, d x\n",
    "   $$\n",
    "1. Using that Laguerre polynomials satisfy the recurrence relation\n",
    "   $$\n",
    "       L_{k + 1}(x) = \\frac{(2k + 1 - x)L_k(x) - k L_{k - 1}(x)}{k + 1}, \\qquad L_0(x) = 1, \\qquad L_1(x) = 1-x,\n",
    "   $$\n",
    "   write a function `laguerre(n)` which returns the Laguerre polynomial of degree $n$.\n",
    "\n",
    "1. Write a function `get_nodes_and_weights(n)` which returns the nodes and weights of the Gauss-Laguerre quadrature    with $n$ nodes. In order to construct Lagrange polynomials, you may find it useful to use the `fromroots` function.\n",
    "You will need to use that\n",
    "   $$\n",
    "   \\int_0^{\\infty} x^n e^{-x} \\, dx = n!\n",
    "   $$\n",
    "\n",
    "1. Write a function `integrate_laguerre(f, n)`, which returns an approximation of $I[f]$ obtained by Gauss-Laguerre integration with $n$ nodes.\n",
    "\n",
    "1. Set $n = 3$, and calculate numerically the degree of precision of the corresponding quadrature rule.\n",
    "\n",
    "1. Set $f(x) = \\sin(x)$, and plot the integration error as a function of $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd50802e",
   "metadata": {},
   "outputs": [],
   "source": [
    "function laguerre(n)\n",
    "    if n == 0\n",
    "        return Polynomial([1])\n",
    "    elseif n == 1\n",
    "        return Polynomial([1, -1])\n",
    "    else\n",
    "        k = n-1\n",
    "        x = Polynomial([0, 1])\n",
    "        return ((2k + 1 - x) * laguerre(k) - k*laguerre(k-1))/(k+1)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3af72ebf",
   "metadata": {},
   "outputs": [],
   "source": [
    "function get_nodes_and_weights(n)\n",
    "    nodes = roots(laguerre(n))\n",
    "    weights = zero(nodes)\n",
    "    for i in 1:n\n",
    "        ℓ = fromroots(nodes[1:end .!= i])\n",
    "        ℓ /= ℓ(nodes[i])\n",
    "        weights[i] = factorial.(0:n-1)'ℓ.coeffs\n",
    "    end\n",
    "    return nodes, weights\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a111240",
   "metadata": {},
   "outputs": [],
   "source": [
    "function integrate_laguerre(f, n)\n",
    "    nodes, weights = get_nodes_and_weights(n)\n",
    "    return f.(nodes)'weights\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0070e085",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 5\n",
    "for i in 1:15\n",
    "    correct = integrate_laguerre(x -> x^i, n) ≈ factorial(i)\n",
    "    println(\"f = x^$i, Rule exact? \", correct)\n",
    "    if !correct\n",
    "        println(\"Degree of precision = \", i - 1)\n",
    "        break\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdf9e1f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "ns = 1:20\n",
    "f(x) = sin(x)\n",
    "In = 1/2\n",
    "Ih = integrate_laguerre.(f, ns)\n",
    "plot(ns, abs.(Ih .- In), yscale=:log10, xlabel=L\"n\", ylabel=\"Error\")\n",
    "scatter!(ns, abs.(Ih .- In))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8230b37",
   "metadata": {},
   "source": [
    "### 3. Exercise on linear systems\n",
    "\n",
    "We consider in this exercise Poisson's equation in the domain $\\Omega = (0, 2) \\times (0, 1)$,\n",
    "equipped with homogeneous Dirichlet boundary conditions:\n",
    "$$\n",
    "    \\begin{aligned}\n",
    "        - \\Delta f(x, y) &= b(x, y), \\qquad x \\in \\Omega, \\\\\n",
    "        f(x) &= 0, \\quad \\qquad x \\in \\partial \\Omega.\n",
    "    \\end{aligned}\n",
    "$$\n",
    "The right-hand side is\n",
    "$$\n",
    "    b(x, y) = \\sin(4\\pi x) + \\sin(2\\pi y).\n",
    "$$\n",
    "A number of methods can be employed in order to discretize this partial differential equation.\n",
    "After discretization, a finite-dimensional linear system of the form $\\mathsf A \\mathbf x = \\mathbf b$ is obtained.\n",
    "A Julia function for calculating the matrix $\\mathsf A$ and the vector $\\mathbf b$ using the finite difference method is given to you below,\n",
    "as well as a function to plot the solution.\n",
    "The goal of this exercise is to solve the linear system using the conjugate gradient method.\n",
    "\n",
    "1. Based on the conjugate gradient algorithm given in the lecture notes,\n",
    "   write a function `conjugate_gradient(A, b)` to solve the linear system $\\mathsf A \\mathbf x = \\mathbf b$.\n",
    "   Use a vector of zeros for $\\mathbf x_0$, and the stopping criterion\n",
    "   $$\n",
    "   \\lVert \\mathsf A \\mathbf x_k - \\mathbf b \\rVert \\leq \\varepsilon \\lVert \\mathbf b \\rVert.\n",
    "   $$\n",
    "\n",
    "1. Make a graph illustrating the evolution of the residual $\\lVert \\mathbf r_k \\rVert$, the error $\\lVert \\mathbf x_k - \\mathbf x_* \\rVert$,\n",
    "     and the error $\\lVert \\mathbf x_k - \\mathbf x_* \\rVert_{\\mathsf A}$ as a function of $k$.\n",
    "     Choose appropriate scales for the graph (linear-linear, log-log, linear-log, or log-linear)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35fc530f",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "using Plots\n",
    "import SparseArrays\n",
    "\n",
    "# Domain size\n",
    "Lx, Ly = 2, 1\n",
    "\n",
    "# Number of discretization points along x and y, including the boundary points\n",
    "nx, ny = 101, 101\n",
    "\n",
    "function discretize(nx, ny)\n",
    "    hx, hy = Lx/(nx - 1), Ly/(ny - 1)\n",
    "    Dxx = (1/hx^2) * SymTridiagonal(2ones(nx-2), -ones(nx-3))\n",
    "    Dyy = (1/hy^2) * SymTridiagonal(2ones(ny-2), -ones(ny-3))\n",
    "    A = kron(Dxx, I(ny-2)) + kron(I(nx-2), Dyy)\n",
    "    xgrid = Lx/(nx-1) * (1:nx-2)\n",
    "    ygrid = Ly/(ny-1) * (1:ny-2)\n",
    "    x_2d = reshape([x for y in ygrid, x in xgrid], (nx-2)*(ny-2))\n",
    "    y_2d = reshape([y for y in ygrid, x in xgrid], (nx-2)*(ny-2))\n",
    "    b = sin.(4π*x_2d) + sin.(2π*y_2d)\n",
    "    return SparseArrays.SparseMatrixCSC(A), b\n",
    "end\n",
    "\n",
    "function plot_solution(f)\n",
    "    f = reshape(f, ny-2, nx-2)\n",
    "    z = [zeros(nx)'; zeros(ny-2) f zeros(ny-2); zeros(nx)']  # Add boundary\n",
    "    xgrid = Lx/(nx-1) * (0:nx-1)\n",
    "    ygrid = Ly/(ny-1) * (0:ny-1)\n",
    "    heatmap(xgrid, ygrid, z, c=:viridis, levels=50)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "716509d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "function conjugate_gradients(A, b)\n",
    "    n = length(b)\n",
    "    x = ones(n)\n",
    "    d = r = A*x - b\n",
    "    ε = 1e-10\n",
    "    xs = [x]\n",
    "    while √(r'r) ≥ ε * √(b'b)\n",
    "        ω = d'r / (d'A*d)\n",
    "        x -= ω*d\n",
    "        r = A*x - b\n",
    "        β = r'A*d/(d'A*d)\n",
    "        d = r - β*d\n",
    "        push!(xs, x)\n",
    "    end\n",
    "    xs\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2604a85",
   "metadata": {},
   "outputs": [],
   "source": [
    "A, b = discretize(nx, ny)\n",
    "x_ = A\\b\n",
    "xs = conjugate_gradients(A, b)\n",
    "plot_solution(xs[end])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c56535b",
   "metadata": {},
   "outputs": [],
   "source": [
    "rs = [A*x - b for x in xs]\n",
    "es = [x - x_ for x in xs]\n",
    "plot(xlabel=L\"k\", yaxis=:log)\n",
    "plot!([√(r'r) for r in rs], label=L\"\\Vert Ax - b \\Vert\")\n",
    "plot!([√(e'r) for (e, r) in zip(es, rs)], label=L\"\\Vert x - x_{*} \\Vert_A\")\n",
    "plot!([√(e'e) for e in es], label=L\"\\Vert x - x_{*} \\Vert\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03659acf",
   "metadata": {},
   "source": [
    "### 4. Exercise on nonlinear equations\n",
    "\n",
    "Data from a biological experiment are collated in the cell below.\n",
    "We wish to fit to this data a model of the form\n",
    "$$\n",
    "    R = \\frac{\\beta_1 C}{\\beta_2 + C},\n",
    "$$\n",
    "where $R$ is the rate, $C$ is the concentration,\n",
    "and $\\beta_1$ and $\\beta_2$ are parameters to be determined.\n",
    "To this end, we shall minimize the sum of square residuals:\n",
    "$$\n",
    "    J(\\mathbf \\beta) := \\sum_{i=1}^{N} \\lvert r_i(\\mathbf \\beta) \\rvert^2, \\qquad r_i(\\mathbf \\beta) := R_i - \\frac{\\beta_1 C_i}{\\beta_2 + C_i}\n",
    "$$\n",
    "\n",
    "1. Find the optimal $\\mathbf \\beta$ using the chord method, applied to the nonlinear equation $\\nabla J(\\beta) = 0$.\n",
    "\n",
    "1. Find the optimal $\\mathbf \\beta$ using the Newton-Raphson method.\n",
    "\n",
    "1. In applications, it is often desirable solve nonlinear least squares problems of the type considered here without calculating second derivatives, i.e. without resorting to the Newton-Raphson method. The most famous algorithm to this end is the so-called Gauss-Newton method.\n",
    "\n",
    "   The idea of the Gauss-Newton method is the following:\n",
    "   given an approximation $\\beta_k$ of the optimal vector of parameters $\\beta$,\n",
    "   we approximate each residual $r_i(\\beta)$ by its linearization around $\\beta_k$,\n",
    "   which is given by\n",
    "   $$\n",
    "      \\widetilde r_i(\\beta) = r_i(\\beta_k) + \\nabla r_i(\\beta_k)^\\top (\\beta - \\beta_k).\n",
    "   $$\n",
    "   We then calculate $\\beta_{k+1}$ as the minimizer of\n",
    "   $$\n",
    "   \\widetilde J(\\mathbf \\beta) = \\sum_{i=1}^{N} \\lvert \\widetilde r_i(\\mathbf \\beta) \\rvert^2.\n",
    "   $$\n",
    "   This is now a **linear** least squares problem,\n",
    "   the solution of which is given by the normal equations\n",
    "   $$\n",
    "   \\beta_{k+1} - \\beta_{k} = (\\mathsf A^\\top \\mathsf A)^{-1} \\mathsf A^\\top \\mathbf b,\n",
    "   $$\n",
    "   where\n",
    "   $$\n",
    "   \\mathsf A =\n",
    "   \\begin{pmatrix}\n",
    "   \\nabla r_1(\\beta_k)^\\top \\\\\n",
    "   \\vdots \\\\\n",
    "   \\nabla r_N(\\beta_k)^\\top\n",
    "   \\end{pmatrix},\n",
    "   \\qquad\n",
    "   \\mathbf b =\n",
    "   -\\begin{pmatrix}\n",
    "   r_1(\\beta_k) \\\\\n",
    "   \\vdots \\\\\n",
    "   r_N(\\beta_k)\n",
    "   \\end{pmatrix}\n",
    "   $$\n",
    "\n",
    "*Hints*:\n",
    "\n",
    "- You may use `ForwardDiff` in order to compute derivatives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83305a26",
   "metadata": {},
   "outputs": [],
   "source": [
    "concentrations = [0.038, 0.194, 0.425, 0.626, 1.253, 2.500, 3.740]\n",
    "rates = [0.050, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8197eb1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "using ForwardDiff\n",
    "r(i, β) = rates[i] - β[1] * concentrations[i] / (β[2] + concentrations[i])\n",
    "r(β) = [r(i, β) for i in 1:length(rates)]\n",
    "J(β) = r(β)'r(β)\n",
    "∇J(β) = ForwardDiff.gradient(J, β)\n",
    "HJ(β) = ForwardDiff.hessian(J, β)   \n",
    "\n",
    "# Chord method\n",
    "β = [0.3; 0.5]\n",
    "α = 1000\n",
    "maxiter = 100_000\n",
    "for i in 1:maxiter\n",
    "    β -= ∇J(β) / α\n",
    "    if norm(∇J(β)) < 1e-12\n",
    "        println(\"Converged in $i iterations!\")\n",
    "        break\n",
    "    end\n",
    "end\n",
    "println(β, \" \", norm(∇J(β)))\n",
    "\n",
    "# Newton-Raphson\n",
    "β = [0.3; 0.5]\n",
    "α = 1000\n",
    "maxiter = 100_000\n",
    "for i in 1:maxiter\n",
    "    β -= HJ(β)\\∇J(β)\n",
    "    if norm(∇J(β)) < 1e-12\n",
    "        println(\"Converged in $i iterations!\")\n",
    "        break\n",
    "    end\n",
    "end\n",
    "println(β, \" \", norm(∇J(β)))\n",
    "\n",
    "# Gauss-Newton\n",
    "β = [0.3; 0.5]\n",
    "α = 1000\n",
    "maxiter = 100_000\n",
    "for i in 1:maxiter\n",
    "    A = ForwardDiff.jacobian(r, β)\n",
    "    β -= (A'A)\\(A'r(β))\n",
    "    if norm(∇J(β)) < 1e-6\n",
    "        println(\"Converged in $i iterations!\")\n",
    "        break\n",
    "    end\n",
    "end\n",
    "println(β)"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "encoding": "# -*- coding: utf-8 -*-",
   "main_language": "julia"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}