{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Solving equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Arrays (minimal)\n",
    "\n",
    "- Solvers make use of Python arrays\n",
    "- Arrays are collections of numbers, similar to mathematical vectors\n",
    "- Arrays are defined as comma-separated numbers contained in square brackets:\n",
    "    - `a = [1.1,-85,3.14,101325]`\n",
    "- Array elements are assigned and accessed using square brackets:\n",
    "    - Array elements are indexed (counted), starting at 0\n",
    "    - `a[0]` is `1.1`\n",
    "    - `a[1]` is `-85`\n",
    "    - `a[2]` is `3.14`\n",
    "    - `a[3]` is `101325`\n",
    "- Change elements:\n",
    "    - `a[1] = -8314.46`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1 -85 3.14 101325\n",
      "[1.1, -85, 3.14, 101325]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[8314.46, -85, 3.14, 101325]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = [1.1,-85,3.14,101325]\n",
    "print(a[0], a[1], a[2], a[3])\n",
    "print(a)\n",
    "\n",
    "a[0] = 8314.46\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2D Arrays (matrices)\n",
    "- A matrix is an array of arrays:\n",
    "    - `A = [[1,2,3], [4,5,6], [7,8,9]]`\n",
    "    $$A = \\begin{bmatrix}1 & 2 & 3 \\\\ 4 & 5 & 6 \\\\ 7 & 8 & 9\\end{bmatrix}$$\n",
    "- And accessed as $A_{i,j} = $ `A[i][j]` for row `i` and column `j`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.1, 2.2, 3.3], [4.4, 5.5, 6.6], [7.7, 8.8, 9.9]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "6.6"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [[1.1,2.2,3.3],[4.4,5.5,6.6],[7.7,8.8,9.9]]\n",
    "print(A)\n",
    "A[1][2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numpy arrays\n",
    "- The above are Python's built-in arrays.\n",
    "- It is generally better to use Numpy arrays for numerical work.\n",
    "- The syntax is basically the same, except we create arrays using the `np.array` function. \n",
    "    - **(More on this later.)**\n",
    "- Note, for 2D numpy arrays, we access and assign with `A[i,j]`, as opposed to `A[i][j]` for regular Python arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.1 2.2 3.3]\n",
      "6.6\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "b = np.array([1.1, 2.2, 3.3])\n",
    "A = np.array([[1.1, 2.2, 3.3],[4.4,5.5,6.6],[7.7,8.8,9.9]])\n",
    "\n",
    "print(b)\n",
    "print(A[1,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boilerplate\n",
    "- import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Solve Linear equations\n",
    "To solve a system $Ax=b$ for x, we do the following:\n",
    "* Define the matrix $A$ and the array $b$.\n",
    "* Use x = np.linalg.solve(A,b)\n",
    "\n",
    "**Example**\n",
    "\n",
    "Solve \n",
    "\n",
    "$$\\begin{bmatrix}\n",
    "-2 & 1 & 0 & 0 & 0 \\\\\n",
    "1 & -2 & 1 & 0 & 0 \\\\\n",
    "0 & 1 & -2 & 1 & 0 \\\\\n",
    "0 & 0 & 1 & -2 & 1 \\\\\n",
    "0 & 0 & 0 & 1 & -2\n",
    "\\end{bmatrix}\n",
    "\\begin{pmatrix}\n",
    "x_0 \\\\ x_1 \\\\ x_2 \\\\ x_3 \\\\ x_4\n",
    "\\end{pmatrix} = \n",
    "\\begin{pmatrix} 1 \\\\ 1 \\\\ 1 \\\\ 1 \\\\ 1\\end{pmatrix}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ax =  [1. 1. 1. 1. 1.]\n",
      "Ax =  [1. 1. 1. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([ [ -2, 1, 0, 0, 0 ],\n",
    "               [ 1, -2, 1, 0, 0 ],\n",
    "               [ 0, 1, -2, 1, 0 ],\n",
    "               [ 0, 0, 1, -2, 1 ],\n",
    "               [ 0, 0, 0, 1, -2 ] ])\n",
    "\n",
    "b = np.array([1, 1, 1, 1, 1])\n",
    "\n",
    "x = np.linalg.solve(A, b)\n",
    "\n",
    "#--------- verify\n",
    "\n",
    "print(\"Ax = \", np.dot(A,x))\n",
    "print(\"Ax = \", A@x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Here's a [link](http://docs.scipy.org/doc/numpy/reference/routines.linalg.html) to more linear algebra operations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Nonlinear equations\n",
    "### 1-D equations (1 equation in 1 unknown)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Recall\n",
    "* What are the key elements to solving nonlinear equations?\n",
    "* List the parts, the form, and things you need."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Exercise\n",
    "* Before learning the details about solving equations using Python, see if you can guess how to solve the following equation using Python's solver ```root```, based on what you need to solve equations:\n",
    "$$f(x) = x^2-5.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Solving a nonlinear equation involves first putting it in the form $f(x)=0$. We know the function $f$ and we want to find the value of $x$ that gives $f(x)=0$.\n",
    "\n",
    "To solve this, do the following:\n",
    "* Define the function f(x)\n",
    "* Set an initial guess for x.\n",
    "* We need to import a library: `from scipy.optimize import root`\n",
    "* Call the `root` function: `res = root(f, x0)`\n",
    "    * `f` is the name of the function you are solving, and `x0` is your initial guess.\n",
    "    * `res` stores a bunch of result information. Access the actual answer with `res.x`\n",
    "\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exercise\n",
    "\n",
    "Solve for both roots of $$f(x) = x^2-5.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.23606798]\n",
      "[-2.23606798]\n"
     ]
    }
   ],
   "source": [
    "def f(x):\n",
    "    return x*x - 5\n",
    "\n",
    "xguess = 1.0\n",
    "res = root(f, xguess)\n",
    "print(res.x)\n",
    "\n",
    "xguess = -1.0\n",
    "res = root(f, xguess)\n",
    "print(res.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Multiple equations in multiple unknowns\n",
    "* What do we do if we are solving multiple equations in multiple unknowns?\n",
    "* For example, suppose we have $h(y,z)=0$ and $g(y,z)=0$. Such as:\n",
    "$$h(y,z) = y+2z,$$\n",
    "$$g(y,z) = \\sin(y)/z.$$\n",
    "\n",
    "* Do you think Python's `root` function can handle this? \n",
    "* How will you go about specifying the problem?\n",
    "* What if there were three equations in three unknowns? Four and four? etc.?\n",
    "* If you were coding `root` how would you account for all these possibilities?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### General approach\n",
    "* Instead of $f(x)=0$, we solve\n",
    "* $\\vec{f}(\\vec{x})=\\vec{0}$.\n",
    "    * That is, $\\vec{x}$ is a vector of unknowns, and $\\vec{f}$ is a vector of functions.\n",
    "    * $\\vec{f}(\\vec{x})$ takes a vector $\\vec{x}$, and returns a vector $\\vec{f}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Details\n",
    "* Consider again the **two functions** for the **two unknowns**: $h(y,z) = 0$, and $g(y,z)=0$.\n",
    "$$h(y,z)=y+2z,$$ $$g(y,z) = \\sin(y)/z$$\n",
    "* For computer solvers, we write using a generic, standard notation. \n",
    "    * <font color=red>**We want all solvers to look like $f(x)=0.$**</font>\n",
    "* We would write these as $$f_0(x_0,x_1) = x_0+2x_1,$$ $$f_1(x_0,x_1) = \\sin(x_0)/x_1,$$\n",
    "which becomes: $$\\left[\\begin{matrix} f_0(x_0,x_1) \\\\ f_1(x_0,x_1)\\end{matrix}\\right] = \\left[\\begin{matrix}x_0+2x_1 \\\\ \\sin(x_0)/x_1\\end{matrix}\\right] = \\left[\\begin{matrix}0 \\\\ 0\\end{matrix}\\right],$$\n",
    "or $$\\vec{f}(\\vec{x}) = \\left[\\begin{matrix}x_0+2x_1 \\\\ \\sin(x_0)/x_1\\end{matrix}\\right] = \\vec{0},$$\n",
    "that is, $$\\vec{f}(\\vec{x})=\\vec{0}.$$\n",
    "\n",
    "* We solve the system just like before:\n",
    "    * Just define the vector function. It takes an array of unknowns $\\vec{x}$, and returns an array of function values $\\vec{f}$.\n",
    "    * We use the same solver, but we give an array of initial guesses $\\vec{x}_0$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Note** Here is a [link](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.root.html) to more information about the root function with extra arguments that can give more control over the solution process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exercise\n",
    "- Solve the following two equations in two unknowns.\n",
    "\n",
    "$$h(y,z)=y+2z,$$\n",
    "\n",
    "$$g(y,z) = \\sin(y)/z.$$\n",
    "\n",
    "- Use $y_0=1$, and $z_0=2$ as an initial guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = 9.424777960769402 y = -4.712388980384701\n"
     ]
    }
   ],
   "source": [
    "def hg(yz):\n",
    "    y = yz[0]\n",
    "    z = yz[1]\n",
    "\n",
    "    h = y + 2*z\n",
    "    g = np.sin(y)/z\n",
    "\n",
    "    return np.array([h,g])\n",
    "\n",
    "#------------\n",
    "\n",
    "xy_guess = np.array([1,2])\n",
    "\n",
    "res = root(hg, xy_guess)\n",
    "print(\"x =\", res.x[0], \"y =\", res.x[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Excercise\n",
    "\n",
    "Solve the following system of three equations in three unknowns:\n",
    "\n",
    "$$x^2 + y^2 = 1,$$\n",
    "$$xy + yz = -1.1,$$\n",
    "$$y^2 + z^2 = 2.$$\n",
    "\n",
    "A reasonable guess for all variables is $x=y=z=2.$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x,y,z =  [ 0.10056089 -0.99493091  1.00504353]\n"
     ]
    }
   ],
   "source": [
    "def F(xyz):\n",
    "    x = xyz[0]\n",
    "    y = xyz[1]\n",
    "    z = xyz[2]\n",
    "\n",
    "    eq0 = x*x + y*y - 1\n",
    "    eq1 = x*y + y*z + 1.1\n",
    "    eq2 = y*y + z*z - 2\n",
    "\n",
    "    return np.array([eq0, eq1, eq2])\n",
    "\n",
    "#-----------------\n",
    "\n",
    "xyz_guess = np.array([2,2,2])\n",
    "res = root(F, xyz_guess)\n",
    "print(\"x,y,z = \", res.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exercise\n",
    "\n",
    "The Colbrook equation is used to compute pressure drops across pipes, and therefore to compute pumping requirements and costs.\n",
    "\n",
    "The Colebrook equation is given by\n",
    "$$\\frac{1}{\\sqrt{f}} =-2\\log_{10}\\left(\\frac{\\epsilon/D}{3.7}+\\frac{2.51}{Re\\sqrt{f}}\\right),$$\n",
    "or,\n",
    "$$F(f, Re, \\epsilon/D) = \\frac{1}{\\sqrt{f}} + 2\\log_{10}\\left(\\frac{\\epsilon/D}{3.7}+\\frac{2.51}{Re\\sqrt{f}}\\right)=0.$$\n",
    "\n",
    "Solve for $f$ given $Re=20000$, and $\\epsilon/D=0.01$. Use 0.01 for a guess for $f$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f = [0.02794571]\n"
     ]
    }
   ],
   "source": [
    "Re     = 20000\n",
    "epsD   = 0.001\n",
    "fguess = 0.01\n",
    "\n",
    "#------------\n",
    "\n",
    "def F(f):\n",
    "    return 1.0/np.sqrt(f) + 2.0*np.log10(epsD/3.7 + 2.51/Re/np.sqrt(f))\n",
    "\n",
    "f = root(F, fguess).x\n",
    "print(\"f =\", f)\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "hide_input": false,
  "kernelspec": {
   "display_name": "base",
   "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.14.7"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}