04a_Newton_Picard.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "![Logo TUBAF](https://tu-freiberg.de/sites/default/files/media/freiberger-alumni-netzwerk-6127/wbm_orig_rgb_0.jpg)\n",
    "\n",
    "Exercise material of the MSc-level course **Numerical Methods in Geotechnical Engineering**.\n",
    "Held at Technische Universität Bergakademie Freiberg.\n",
    "\n",
    "Comments to:\n",
    "\n",
    "*Prof. Dr. Thomas Nagel  \n",
    "Chair of Soil Mechanics and Foundation Engineering  \n",
    "Geotechnical Institute  \n",
    "Technische Universität Bergakademie Freiberg.*\n",
    "\n",
    "https://tu-freiberg.de/en/soilmechanics\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Iteration schemes and their convergence: Newton-Raphson and Picard"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example problem\n",
    "\n",
    "Given a function $f(x)$ we seek the solution $x^*$ such that $f(x^*) = 0$.\n",
    "\n",
    "Here, for illustrative purposes, we seek the solution to $f(x) = \\sin x$. Analytically, we find for this periodic function $x^*_n = n \\pi$ for $n \\in Z$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np #numerical methods\n",
    "import matplotlib.pyplot as plt #plotting\n",
    "\n",
    "from ipywidgets import widgets\n",
    "from ipywidgets import interact\n",
    "\n",
    "#Some plot settings\n",
    "import plot_functions.plot_settings\n",
    "%run plot_functions/iteration_schemes_plots.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define the function:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Newton-Raphson scheme\n",
    "\n",
    "Anticipating an iterative scheme, we depart from the Taylor series expansion around the given solution $x_i$ in order to find the next solution $x_{i+1}$, which we truncate at first order and for which we demand equality to zero, as we're seeking $f(x^*) = 0$:\n",
    "\n",
    "$$\n",
    "    f(x_{i+1}) \\approx f(x_i) + \\left. \\frac{\\partial f(x)}{\\partial x} \\right|_{x_i} (x_{i+1} - x_{i}) \\overset{!}{=} 0\n",
    "$$\n",
    "\n",
    "By rearrangement, we find the iterative update procecture\n",
    "\n",
    "$$\n",
    "    x_{i+i} = x_i - \\left( \\left. \\frac{\\partial f(x)}{\\partial x} \\right|_{x_i} \\right)^{-1} f(x_i)\n",
    "$$\n",
    "\n",
    "The iteration is started at a starting guess $x_0$ which, ideally, is close to the sought solution. We will investigate it's role in the sequel.\n",
    "\n",
    "We stop the iteration once we're sufficiently close to zero:\n",
    "\n",
    "$$\n",
    "    |f(x_i)| < \\varepsilon_\\text{tol} \\quad \\text{with } \\varepsilon_\\text{tol} \\ll 1 \\quad \\rightarrow \\quad x^* \\approx x_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the chosen problem we have:\n",
    "\n",
    "$$\n",
    "    f(x) = \\sin x \\quad \\rightarrow \\quad f'(x) = \\cos x\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f7bbb4bf6ee54370ae811597eeea131e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedFloatText(value=1.1, description='$x_0$', max=2.0), IntSlider(value=6, descriptio…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_interactive_newton()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tasks\n",
    "\n",
    "Increase the maximum iteration number one by one to see the iteration unfold. You can also play with the tolerance value.\n",
    "\n",
    "Now, keeping the tolerance fixed at $10^{-1}$ and the maximum iteration number at 10, answer the following questions.\n",
    "\n",
    "What happens -- and why -- if you \n",
    "\n",
    "* set the starting value to $x_0 = 2$?\n",
    "* set the starting value to $x_0 = 1.6$?\n",
    "* set the starting value to $x_0 = \\pi/2$ or multiples thereof?\n",
    "* increase the starting value from 1.15 to 1.2 in steps of 0.01?\n",
    "\n",
    "Note your observations in terms of the solution found, the rate of convergence, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Picard iterations\n",
    "\n",
    "You're already familiar with this type of iterations. Remember, for example, the iteration rule for the method of slices for the simplified Bishop's method:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dies führt auf die Formel nach DIN4084:\n",
    "\n",
    "\\begin{align*}\n",
    "    \\mu = \\frac{E_\\text{M,d}}{R_\\text{M,d}} = \\frac{\\sum \\limits_i \\left[ P_\\text{vi,d} + G_\\text{i,d} \\right] \\sin \\vartheta_i}{\\displaystyle \\sum \\limits_i  \\frac{(G_{i,\\text{d}} + P_{\\text{v},i,\\text{d}} - u_{i,\\text{d}}b_i) \\tan \\varphi_{i,\\text{d}} + c_{i,\\text{d}}b_i}{\\cos\\vartheta_i + \\mu \\tan \\varphi_{i,\\text{d}} \\sin \\vartheta_i}}\n",
    "\\end{align*}\n",
    "\n",
    "Vorgehen: Einsetzen von $\\mu_0 = 1$ oder einem anderen geschätzten Startwert. Dann wiederholtes Einsetzen und Ausrechnen von $\\mu$ bis\n",
    "\\begin{align*}\n",
    "    \\left| \\frac{\\mu_{i+1}-\\mu_i}{\\mu_i} \\right| < 0.03\n",
    "\\end{align*}\n",
    "\n",
    "Konvergenz wird i.d.R. innerhalb weniger Iterationen erzielt."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe, that we have an iteration update of the form\n",
    "\n",
    "$$\n",
    "    x_{i+1} = g(x_i)\n",
    "$$\n",
    "\n",
    "Therefore, we re-cast our example function $f(x) = \\sin x$ into the following form by using $f(x) = \\sin x + x - x$:\n",
    "\n",
    "$$\n",
    "    x = f(x) + x = g(x) = \\sin x + x\n",
    "$$\n",
    "\n",
    "In other words, finding $x^*$ for $f(x^*) = 0$ is now the fixed point problem of identifying $x^* = g(x^*)$.\n",
    "\n",
    "That's it -- we can start looking at how it works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a2d0e8f7b52d4a15b2ae2cf04c61d225",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedFloatText(value=0.7, description='$x_0$', max=3.0, min=-1.0), IntSlider(value=10,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_interactive_picard()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tasks\n",
    "\n",
    "* What do you observe for the identified solution as you vary the starting guess? Why is this?\n",
    "* What can you say about the convergence rate?\n",
    "\n",
    "Now let's repeat this by rearranging $f(x)$ slightly differently:\n",
    "\n",
    "$$\n",
    "    g(x) = x - \\sin(x)\n",
    "$$\n",
    "\n",
    "This definition is equivalent in terms of the roots ($\\sin x$ and $-\\sin x$ share the same roots)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ab69e3f754cd4900b8fd026aac26045e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedFloatText(value=0.7, description='$x_0$', max=3.0, min=-1.0), IntSlider(value=10,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_interactive_picard_inv()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The attractors change, however the principal problem remains the same as before ....\n",
    "\n",
    "\n",
    "## Additional material\n",
    "\n",
    "### Modified Newton scheme\n",
    "\n",
    "In the lecture we saw the modified Newton schemes and the initial stiffness schemes. Here we look at such a modification by keeping keeping the slope constant:\n",
    "\n",
    "$$\n",
    "    x_{i+i} = x_i - \\left( \\left. \\frac{\\partial f(x)}{\\partial x} \\right|_{x_0} \\right)^{-1} f(x_i)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "30b8438af415412d93ad2a04d69773b2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedFloatText(value=1.1, description='$x_0$', max=2.0), IntSlider(value=6, descriptio…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_interactive_newton_modified()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tasks\n",
    "\n",
    "Play with the starting values and observe the convergence rate. Especially, try $x_0 = [1.3, 1.1, 1.0, 0.8, 0.5, 0.2]$. What do you observe?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another function ..\n",
    "\n",
    "We now look at \n",
    "\n",
    "$$\n",
    "    f(x) = e^{-x} - x\n",
    "$$\n",
    "\n",
    "### Newton:\n",
    "\n",
    "The 'Jacobian' is given as $f'(x) = -e^{-x} - 1$\n",
    "\n",
    "### Picard:\n",
    "\n",
    "The fixed point equation is given ax $x = g(x) = e^{-x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "80ee6da02c8c481cb03609c8ef147bcc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedFloatText(value=1.0, description='$x_0$', max=1.5), IntSlider(value=10, descripti…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_interactive_picard_2()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "eeb46f6cdb764f369a036afc7b8bafcd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedFloatText(value=1.0, description='$x_0$', max=1.5), IntSlider(value=6, descriptio…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_interactive_newton_2()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tasks:\n",
    "\n",
    "* Do the methods converge to the same solution?\n",
    "* How quickly do they converge? Check how quickly the residuals decrease.\n",
    "* How would the modified Newton scheme perform in this simpler case? "
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}