{
 "cells": [
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "# Computer session 2\n",
    "\n",
    "The goal of this computer class is to get a good feel of the Newton method and its\n",
    "variants. In a (maybe) surprising way, we actually start with the dichotomy method\n",
    "in the one-dimensional case.\n",
    "\n",
    "## The dichotomy method in the one-dimensional case\n",
    "\n",
    "When trying to solve the equation $\\phi(x) = 0$ in the one-dimensional case, the\n",
    "most naive method, which actually turns out to be quite efficient, is the dichotomy\n",
    "method. Namely, starting from an initial pair $(a_L , a_R ) \\in \\mathbb{R}^2$ with $a_L < a_R$ such\n",
    "that $\\phi(a_L)\\phi(a_R)<0$, we set $b ∶=\\frac{a_L+a_R}{2}$. If $\\phi(b) = 0$, the algorithm stops. If\n",
    "$\\phi(a_L)\\phi(b) < 0$ we set $a_L\\to a_L$ and $a_R \\to b$. In this way, we obtain a linearly\n",
    "converging algorithm. In particular, it is globally converging.\n",
    "\n",
    "\n",
    "Write a function `Dichotomy(phi,aL,aR,eps)` that take sas argument\n",
    "a function `phi`, an initial guess `aL,aR` and a tolerance `eps` and that runs the dichotomy\n",
    "algorithm. Your argument should check that the condition $\\phi(a_L)\\phi(a_R) < 0$ is satisfied,\n",
    "stop when the function `phi` reaches a value lower than `eps` and return the number\n",
    "of iteration. Run your algorithm on the function $f = tanh$ with initial guesses\n",
    "$a_L = −20$ , $a_R = 3$.\n",
    "\n"
   ]
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-18T16:19:14.314484Z",
     "start_time": "2025-03-18T16:19:13.728014Z"
    }
   },
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "def dichotomy(phi, aL, aR, eps=1e-8):\n",
    "    iter = 0\n",
    "    b = (aL + aR) / 2\n",
    "    while (aR - aL) / 2 > eps and phi(b) != 0:\n",
    "        b = (aL + aR) / 2\n",
    "        if phi(aL) * phi(b) < 0:\n",
    "            aR = b\n",
    "        else:\n",
    "            aL = b\n",
    "        iter += 1\n",
    "    return b, iter\n",
    "\n",
    "\n",
    "f = lambda x: np.tanh(x)\n",
    "aL, aR = -20, 3\n",
    "print(dichotomy(f, aL, aR))"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-2.3283064365386963e-09, 31)\n"
     ]
    }
   ],
   "execution_count": 1
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Solving one-dimensional equation with the Newton and the secant method\n",
    "\n",
    "We work again in the one-dimensional case with a function φ we want to find the\n",
    "zeros of.\n",
    "\n",
    "### Newton method\n",
    "\n",
    "Write a function `Newton(phi,dphi,x0,eps)` that takes, as arguments, a function\n",
    "`phi`, its derivative `dphi`, an initial guess `x0` and a tolerance `eps`\n",
    "and that returns an approximation of the solutions of the equation $\\phi(x) = 0$. The\n",
    "tolerance criterion should again be that $|\\phi| ≤\\text{\\texttt{eps}}$. Your\n",
    "algorithm should return an error message in the following cases:\n",
    "1. If the derivative is zero (look up the `try` and `except` commands in Python).\n",
    "2. If the method diverges.\n",
    "\n",
    "Apply this code to the minimisation of $x\\mapsto \\ln(e^x + e^{−x})$, with initial condition `x0=1.8`.\n",
    "Compare this with the results of Exercise 3.10."
   ]
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-18T16:19:17.447647Z",
     "start_time": "2025-03-18T16:19:17.442560Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def Newton(phi, dphi, x0, eps=1e-10):\n",
    "    iter = 0\n",
    "    while abs(phi(x0)) > eps:\n",
    "        try:\n",
    "            x0 -= phi(x0) / dphi(x0)\n",
    "        except ZeroDivisionError:\n",
    "            return np.inf, \"Derivative is zero\"\n",
    "        iter += 1\n",
    "        if iter > 10000 or phi(x0) == np.inf:\n",
    "            return np.inf, \"Method diverges\"\n",
    "    return x0, iter\n",
    "\n",
    "\n",
    "f = lambda x: np.log(np.exp(x) + np.exp(-x))\n",
    "x0 = 1.8\n",
    "df = lambda x: (np.exp(x) - np.exp(-x)) / (np.exp(x) + np.exp(-x))\n",
    "print(Newton(f, df, x0))"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(inf, 'Method diverges')\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/tp/_ld5_pzs6nx6mv1pbjhq1l740000gn/T/ipykernel_25957/3868809151.py:14: RuntimeWarning: overflow encountered in exp\n",
      "  f = lambda x: np.log(np.exp(x) + np.exp(-x))\n"
     ]
    }
   ],
   "execution_count": 2
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "### Secant method\n",
    "\n",
    "Write a function `Secant(phi,x0,x1,eps)` that takes, as arguments, a function `phi`, two initial positions `x0`, `x1` and a tolerance\n",
    "`eps` and that returns an approximation of the solutions of the equation\n",
    "$\\phi(x) = 0$. The tolerance criterion should again be that $|\\phi| ≤\\text{\\texttt{eps}}$. Apply this code to the minimisation\n",
    "of $x\\mapsto \\ln(e^x + e^{−x})$, with initial conditions `x0=1`, `x1=1.9`, then\n",
    "`x0=1`, `x1=2.3` and\n",
    "`x0=1`, `x1=2.4`. Compare with the results of Exercise 3.10."
   ]
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-18T16:19:19.649523Z",
     "start_time": "2025-03-18T16:19:19.456149Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def Secant(phi, x0, x1, eps=1e-8):\n",
    "    iter = 0\n",
    "    while abs(phi(x1)) > eps:\n",
    "        x1, x0 = x1 - phi(x1) * (x1 - x0) / (phi(x1) - phi(x0)), x0\n",
    "        iter += 1\n",
    "        if iter > 10000:\n",
    "            return np.inf, \"Method diverges\"\n",
    "    return x0, iter\n",
    "\n",
    "\n",
    "f = lambda x: np.log(np.exp(x) + np.exp(-x))\n",
    "xx = [(1, 1.9), (1, 2.3), (1, 2.4)]\n",
    "\n",
    "for x0, x1 in xx:\n",
    "    print(Secant(f, x0, x1))"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(inf, 'Method diverges')\n",
      "(inf, 'Method diverges')\n",
      "(inf, 'Method diverges')\n"
     ]
    }
   ],
   "execution_count": 3
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Combining dichotomy and the Newton method\n",
    "\n",
    "A possibility to leverage the advantages of dichotomy (the global convergence of\n",
    "the method) and of the Newton method (the quadratic convergence rate) is to\n",
    "combine both: start from an initial interval `[aL,aR]` of length `InitialLength`\n",
    "with $\\phi(a_L)\\phi(a_R)<0$ and fix a real\n",
    "number $s \\in [0; 1]$. Run the dichotomy algorithm until the new interval is of length\n",
    "`s*InitialLength`. From this point on, apply the Newton method.\n",
    "\n",
    "Implement this algorithm with `s = 0.1`. Include a possibility to switch\n",
    "back to the dichotomy method if, when switching to the Newton method, the new\n",
    "iterate falls outside of the computed interval `[aL,aR]`. Apply this to the minimisation\n",
    "of the function $f ∶ x\\mapsto \\ln(e^x + e^{−x})$ with an initial condition that made the Newton\n",
    "method diverge. What can you  say about the number of iterations?"
   ]
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-18T16:44:41.592150Z",
     "start_time": "2025-03-18T16:44:41.584318Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def DichotomyNewton(phi, dphi, aL, aR, s=0.1, eps=1e-10):\n",
    "    iter = 0\n",
    "    inital_length = aR - aL\n",
    "    while (aR - aL) >= s * inital_length:\n",
    "        b = (aL + aR) / 2\n",
    "        if phi(aL) * phi(b) < 0:\n",
    "            aR = b\n",
    "        else:\n",
    "            aL = b\n",
    "        iter += 1\n",
    "    x0 = (aL + aR) / 2\n",
    "    while abs(phi(x0)) > eps:\n",
    "        try:\n",
    "            x0 -= phi(x0) / dphi(x0)\n",
    "        except ZeroDivisionError:\n",
    "            return np.inf, \"Derivative is zero\"\n",
    "        iter += 1\n",
    "        if iter > 10000 or phi(x0) == np.inf:\n",
    "            return np.inf, \"Method diverges\"\n",
    "    return x0, iter\n",
    "\n",
    "\n",
    "f = lambda x: np.log(np.exp(x) + np.exp(-x))\n",
    "df = lambda x: (np.exp(x) - np.exp(-x)) / (np.exp(x) + np.exp(-x))\n",
    "print(DichotomyNewton(f, df, -20, 3))"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(inf, 'Method diverges')\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/tp/_ld5_pzs6nx6mv1pbjhq1l740000gn/T/ipykernel_25957/1578277506.py:23: RuntimeWarning: overflow encountered in exp\n",
      "  f = lambda x: np.log(np.exp(x) + np.exp(-x))\n"
     ]
    }
   ],
   "execution_count": 22
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Solving an optimisation problem using the Newton method\n",
    "\n",
    "An island (denoted by a point $I$ below) is situated 2 kilometers from the shore (its projection on the shore\n",
    "is a point $P$). A guest staying at a nearby hotel $H$ wants to go from the hotel to the\n",
    "island and decides that he will run at 8km/hr for a distance $x$, before swimming at\n",
    "speed 3km/hr to reach the island.\n",
    "\n",
    "![illustration of the problem](./images/optiNewton.png)\n",
    "\n",
    "Taking into account the fact that there are 6 kilometers between the hotel $H$ and $P$,\n",
    "how far should the visitor run before swimming?\n",
    "\n",
    "Model the situation as a minimisation problem, and solve it numerically.\n",
    "Compare the efficiency of the dichotomy method and of the Newton algorithm."
   ]
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-18T17:43:43.061916Z",
     "start_time": "2025-03-18T17:43:43.042625Z"
    }
   },
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "\n",
    "u = lambda x: np.sqrt((6 - x) ** 2 + 4)\n",
    "\n",
    "\n",
    "def objective_function(x):\n",
    "    return x / 8 + u(x) / 3\n",
    "\n",
    "\n",
    "def gradient(x):\n",
    "    return 1 / 8 + (6 - x) / (3 * u(x))\n",
    "\n",
    "\n",
    "def hessian(x):\n",
    "    return (12 * u(x) - (2 * x - 12) ** 2 / 12 * u(x)) / 36 * u(x) ** 2\n",
    "\n",
    "\n",
    "# Newton's method for optimization\n",
    "def newton_method(initial_guess, tolerance=1e-6, max_iterations=100):\n",
    "    x = initial_guess\n",
    "    iterations = 0\n",
    "    while iterations < max_iterations:\n",
    "        grad = gradient(x)\n",
    "        hess = hessian(x)\n",
    "        if np.abs(grad) < tolerance:\n",
    "            break\n",
    "        x -= grad / hess\n",
    "        iterations += 1\n",
    "    return x, iterations\n",
    "\n",
    "\n",
    "# Dichotomy method for optimization\n",
    "def dichotomy_method(aL, aR, eps=1e-6, max_iterations=1000):\n",
    "    iterations = 0\n",
    "    x0 = (aL + aR) / 2\n",
    "    grad = gradient(x0)\n",
    "    hess = hessian(x0)\n",
    "    while abs(grad) > eps:\n",
    "        try:\n",
    "            x0 -= grad / hess\n",
    "        except ZeroDivisionError:\n",
    "            return np.inf, \"Derivative is zero\"\n",
    "        iterations += 1\n",
    "        if iterations > max_iterations or grad == np.inf:\n",
    "            return np.inf, \"Method diverges\"\n",
    "        grad = gradient(x0)\n",
    "        hess = hessian(x0)\n",
    "    return x0, iterations\n",
    "\n",
    "\n",
    "# Initial guess for Newton's method\n",
    "initial_guess_newton = 4\n",
    "\n",
    "# Run Newton's method\n",
    "optimal_point_newton, iterations_newton = newton_method(initial_guess_newton)\n",
    "print(f\"Optimal point (Newton): {optimal_point_newton}\")\n",
    "print(f\"Objective function value at optimal point (Newton): {objective_function(optimal_point_newton)}\")\n",
    "print(f\"Number of iterations (Newton): {iterations_newton}\")\n",
    "\n",
    "# Initial interval for dichotomy method\n",
    "aL, aR = 0, 6\n",
    "\n",
    "# Run dichotomy method\n",
    "optimal_point_dichotomy, iterations_dichotomy = dichotomy_method(aL, aR)\n",
    "print(f\"Optimal point (Dichotomy): {optimal_point_dichotomy}\")\n",
    "print(f\"Objective function value at optimal point (Dichotomy): {objective_function(optimal_point_dichotomy)}\")\n",
    "print(f\"Number of iterations (Dichotomy): {iterations_dichotomy}\")"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal point (Newton): 0.9299901531755377\n",
      "Objective function value at optimal point (Newton): 1.9329918821224974\n",
      "Number of iterations (Newton): 100\n",
      "Optimal point (Dichotomy): inf\n",
      "Objective function value at optimal point (Dichotomy): inf\n",
      "Number of iterations (Dichotomy): Method diverges\n"
     ]
    }
   ],
   "execution_count": 38
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## The Newton method to solve boundary value problems\n",
    "\n",
    "### Shooting method\n",
    "\n",
    "We consider the following non-linear ODE\n",
    "\n",
    "\\begin{equation}\n",
    "y''= f(x,y,y'),\\quad x\\in[a,b],\\quad y(a)=\\alpha, y(b)=\\beta.\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "To classically integrate such an ODE, we usually don’t have endpoints for $y$,\n",
    "but initial values for $y$ and $y'$. So, we cannot start at $x=a $ and integrate\n",
    "up to $x=b$. This is a _boundary value problem_.\n",
    "\n",
    "One approach is to approximate $y$ by somme finite difference and then arrive at\n",
    "a system for the discrete values $y(x_i)$ and finally solve large linear\n",
    "systems.\n",
    "\n",
    "Here, we will see how we can formulate the problem as a _shooting_ method, and\n",
    "use Newton method so solve it.\n",
    "\n",
    "The idea is to use a _guess_ for the initial value of $y'$. Let $s$ be a\n",
    "parameter for a fonction $y(\\;\\cdot\\;;s)$ solution of (1) such that\n",
    "$$y(a;s)=\\alpha,\\text{ and }y'(a;s)=s.$$\n",
    "\n",
    "There is no chance that $y(b;s)=y(b)=\\beta$ but we can adjust the value of $s$,\n",
    "refining the guess until it is (nearly equal to) the right value.\n",
    "\n",
    "This method is known as _shooting_ method in analogy to shooting a ball at a\n",
    "goal, determining the unknown correct velocity by throwing it too fast/too\n",
    "slow until it hits the goal exactly.\n",
    "\n",
    "#### In Practice\n",
    "\n",
    "For the parameter $s$, we integrate the following ODE:\n",
    "\n",
    "\\begin{equation}\n",
    "y''= f(x,y,y'),\\quad x\\in[a,b],\\quad y(a)=\\alpha, y'(a)=s.\n",
    "\\end{equation}\n",
    "\n",
    "We denote $y(\\;\\cdot\\;;s)$ solution of (2).\n",
    "\n",
    "Let us now define the _goal function_. Here, we want that $y(b;s)=\\beta$, hence,\n",
    "we define:\n",
    "$$g:s\\mapsto \\left.y(x;s)\\right|_{x=b}-\\beta$$\n",
    "\n",
    "We seek $s^*$ such that $g(s^*)=0$.\n",
    "\n",
    "Note that computing $g(s)$ involves the integration of an ODE, so each\n",
    "evaluation of $g$ is  expensive. Newton’s method seems then to be a good\n",
    "way due to its fast convergence.\n",
    "\n",
    "To be able to code a Newton’s method, we need to compute the derivative of $g$.\n",
    "For this purpose, let define\n",
    "$$z(x;s)=\\frac{\\partial y(x;s)}{\\partial s}.$$\n",
    "\n",
    "Then by differentiating (2) with respect to $s$, we get\n",
    "$$z''=\\frac{\\partial f}{\\partial y}z+\\frac{\\partial f}{\\partial y'}z',\\quad\n",
    "z(a;s)=0,\\text{ and }z'(a;s)=1.$$\n",
    "\n",
    "The derivative of $g$ can now be expressed in term of $z$:\n",
    "$$g'(z)=z(b;s).$$\n",
    "\n",
    "Putting this together, we can code the Newton's method:\n",
    "$$s_{n+1}=s_n-\\frac{g(s_n)}{g'(s_n)}.$$\n",
    "\n",
    "To sum up, a shooting method requires an ODE solver and a Newton solver.\n",
    "\n",
    "#### Example\n",
    "\n",
    "Apply this method to\n",
    "$$y''=2y^3-6y-2x^3,\\quad y(1)=2,y(2)=5/2,$$\n",
    "with standard library for integration, and your own Newton implementation.\n",
    "\n",
    "Note that you may want to express this with one order ODE. Moreover, it may be\n",
    "simpler to solve only one ODE for both $g$ and $g'$.\n",
    "\n",
    "With python, you can use `scipy.integrate.solve_ivp` function:\n",
    "https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp\n",
    "\n",
    "Plot the solution $y$.\n",
    "\n",
    "For numerical parameters, compute the solution up to a precision of $10^{-8}$\n",
    "and get the function on a grid of 1000 points.\n",
    "\n",
    "## Finite differences\n",
    "\n",
    "Here, we are going to use a different approach to solve the boundary value\n",
    "problem:\n",
    "\n",
    "\\begin{equation}\n",
    "y''= f(x,y,y'),\\quad x\\in[a,b],\\quad y(a)=\\alpha, y(b)=\\beta.\n",
    "\\end{equation}\n",
    "\n",
    "This problem can be solved by the following direct process:\n",
    "\n",
    "1. We discretize the domain choosing an integer $N$, grid points $\\{x_n\\}_{n=0,\\dots,N}$ and\n",
    "   we define the discrete solution $\\{y_n\\}_{n=0,\\dots,N}$.\n",
    "1. We discretize the ODE using derivative approximation with finite differences\n",
    "   in the interior of the domain.\n",
    "1. We inject the boundary conditions (here $y_0=\\alpha$ and $y_N=\\beta$) in the\n",
    "   discretized ODE.\n",
    "1. Solve the system of equation for the unknows $\\{y_n\\}_{n=1,\\dots,N-1}$.\n",
    "\n",
    "We use here a uniform grid :\n",
    "$$h:=(b-a)/N, \\quad \\forall n=0,\\dots, N\\quad x_n=hn.$$\n",
    "If we use a centered difference formula for $y''$ and $y'$, we obtain:\n",
    "$$\\forall n=1,\\dots,N-1,\\quad \\frac{y_{n+1}-2y_n+y_{n-1}}{h^2}=f\\left(x_n,y_n,\\frac{y_{n+1}-y_{n-1}}{2h}\\right).$$"
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "The result is a system of equations for $\\mathbf{y}=(y_1,\\dots,y_{N-1})$ :\n",
    "$$G(\\mathbf{y})=0,\\quad G:\\mathbb{R}^{N-1}\\to\\mathbb{R}^{N-1}.$$"
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "This system can be solved using Newton's method. Note that the Jacobian\n",
    "$\\partial G/\\partial \\mathbb{y}$ is _tridiagonal_.\n",
    "\n",
    "Of course, here, we are in the multidimensional context, so you will have to\n",
    "code a Newton algorithm well suited.\n",
    "\n",
    "#### Example\n",
    "\n",
    "Apply this method to\n",
    "$$y''=2y^3-6y-2x^3,\\quad y(1)=2,y(2)=5/2.$$\n",
    "\n",
    "Plot the solution $y$.\n",
    "\n",
    "For numerical parameters, compute the solution up to a precision of $10^{-8}$\n",
    "and get the function on a grid of 1000 points."
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "**Remark:** In the context of numerical optimal control, these two numerical\n",
    "methods are often called _indirect method_ (for the shooting method) and _direct\n",
    "method_ (for the finite difference method)."
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## Who’s the best?\n",
    "\n",
    "Compare the two methods playing with parameters (grid discretization, precision,\n",
    "initialization, etc.). Mesure the time computation."
   ]
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": ""
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "metadata": {},
   "cell_type": "code",
   "outputs": [],
   "execution_count": null,
   "source": "\n"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.21"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
