3. Solving an equation (or systems of equations).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook I show you \n",
    "- how to solve equations of one variable using the package Roots, \n",
    "- how to solve systems of equations using NLSolve.\n",
    "\n",
    "Of course NLSolve can also solve univariate equations but it is a bit more complicated to use than Roots. Therefore, I show you Roots first which might be all you need."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving univariate equations using Roots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual I start with a ridiculously simple equation, namely $x+3=0.$ The way this equation is presented to Roots is that we define a function $f(x)=x+3$ and then we want to know where this function is zero, i.e. intersects the x-axis. Roots is made to give the answer to exactly such problems.\n",
    "\n",
    "To start with we install the Roots package (in case you have not done so already in the past)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mCloning cache of Roots from https://github.com/JuliaMath/Roots.jl.git\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling Roots v0.4.1\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mBuilding SpecialFunctions\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage database updated\n",
      "\u001b[39m"
     ]
    }
   ],
   "source": [
    "Pkg.add(\"Roots\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we tell julia to use Roots, then we define the function $f$ and finally tell Roots to solve our equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-3.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Roots\n",
    "f(x)=x+3\n",
    "fzero(f,0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The command to solve the equation is fzero (note that the command is still \"fzero\" if your function is not called $f$). The first argument you put into this command is the function and the second argument is a starting point for the algorithm, i.e. a rough guess of what the solution could be.\n",
    "\n",
    "If your function is well behaved, e.g. monotone, then the starting value will not matter much. However, some equations have several solutions. The roots package will only give you one solution. Which solution this is depends on the starting value. Let me illustrate this with the equation $x^2=0.25$ which has the two solutions 0.5 and -0.5. I will write this first as a function $g(x)=x^2-0.25$ and then tell Roots to solve this equation with different starting values. As you will see, the solution I get differs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5000000000000001"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g(x)=x^2-0.25\n",
    "fzero(g,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.5"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fzero(g,-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes you know that an equation has several solutions but you are only interested in one of them. For example, your variable is a probability and therefore the solution you are interested in is the one in $[0,1]$. You can tell Roots to look only for solutions in this interval, i.e. instead of a starting value you can give Roots a lower and an upper bound. However, Roots will only accept this if the $f(upper bound)$ and $f(lower bound)$ have opposite signs, i.e. one has to be positive and the other negative.\n",
    "\n",
    "I use the same example as above: $x^2=0.25$ which has the two solutions 0.5 and -0.5. Say, we are only interested in solutions between 0 and 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fzero(g,[0.0,1.0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, I use a vector [lower bound, upper bound] as second argument instead of the starting value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ...but I really have to know all solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, so you suspect your equation to have several solutions and you want to know all of them. Of course, you can try to solve with a lot of different starting values and see whether your result changes but you still could miss out on a solution. In general finding all solutions of a (system of) equation(s) is a hard problem but in case of one equation of one variable there is an easy way out (at least if you can restrict your variable to an interval): the solution is to write your equation as a function (as above) and the plot this function! At every point at which the function intersects the x-axis your equation has a solution. The plot gives you immediately an overview how many solutions you have and where (approximately) they are. \n",
    "\n",
    "For example, we could plot the $g$ function above to see that the equation $x^2=0.25$ has two solutions 8and that they are -0.5 and 0.5). I do this here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<script src=\"/home/christoph/.julia/v0.6/Plots/src/backends/../../deps/plotly-latest.min.js\"></script>    <div id=\"84356822-a29f-4770-8271-7817e981a1fd\" style=\"width:600px;height:400px;\"></div>\n",
       "    <script>\n",
       "    PLOT = document.getElementById('84356822-a29f-4770-8271-7817e981a1fd');\n",
       "    Plotly.plot(PLOT, [{\"showlegend\":true,\"mode\":\"lines\",\"xaxis\":\"x\",\"colorbar\":{\"title\":\"\"},\"line\":{\"color\":\"rgba(0, 154, 250, 1.000)\",\"shape\":\"linear\",\"dash\":\"solid\",\"width\":1},\"y\":[3.711032095141974,3.672254932726204,3.6336685127526884,3.5952728352214263,3.2492851269548884,2.919605272941246,2.606233273180499,2.3091691276726483,1.9737023223879113,1.6617939020608075,1.416458800375875,1.1879643802382422,0.9588382884766735,0.7495842829497379,0.5742756023340709,0.41585662319809324,0.264775108269963,0.1331067953122433,0.004831842283663679,-0.0973838091208333,-0.1663515607439204,-0.2147388310040407,-0.2402335967050226,-0.24990257019278828,-0.24241122306608404,-0.21610799023332192,-0.17066850581744425,-0.10619286281008558,-0.012748316992100522,0.1039662765275377,0.23662093515814903,0.3903420664616095,0.5725161634443228,0.777465591526814,0.967644704281261,1.1739625384061412,1.4477329038452074,1.745559069976096,2.047459551183071,2.3706170717921715,2.93264790881421,3.5492381784562914,3.6489730006783123],\"type\":\"scatter\",\"name\":\"y1\",\"yaxis\":\"y\",\"x\":[-1.9902341809802118,-1.9804683619604238,-1.9707025429406357,-1.9609367239208475,-1.8706376257722628,-1.780338527623678,-1.6900394294750933,-1.5997403313265088,-1.4912083430520067,-1.3826763547775045,-1.2909139399572207,-1.1991515251369371,-1.0994718225023656,-0.9997921198677943,-0.9078962508646409,-0.8160003818614874,-0.717478298117764,-0.6189562143740406,-0.5048087185099557,-0.3906612226458709,-0.2892203991008926,-0.18777957555591424,-0.09882511469751709,-0.009870653839119958,0.08711358639107895,0.18409782662127785,0.2816584708162631,0.3792191150112484,0.48708488275443274,0.5949506504976171,0.6975822067384955,0.8002137629793739,0.9069267685123881,1.0136397740454022,1.1034693943563914,1.1932990146673805,1.3029707993064186,1.4126425839454564,1.5157372962301452,1.618832008514834,1.7839977322895368,1.9491634560642397,1.9745817280321198]}], {\"showlegend\":true,\"xaxis\":{\"gridwidth\":0.5,\"tickvals\":[-1.5,-1.0,-0.5,0.0,0.5,1.0,1.5],\"visible\":true,\"ticks\":\"inside\",\"tickmode\":\"array\",\"domain\":[0.02137649460484106,0.9934383202099737],\"linecolor\":\"rgba(0, 0, 0, 1.000)\",\"showgrid\":true,\"title\":\"\",\"mirror\":false,\"tickangle\":0,\"showline\":true,\"gridcolor\":\"rgba(0, 0, 0, 0.100)\",\"titlefont\":{\"color\":\"rgba(0, 0, 0, 1.000)\",\"family\":\"sans-serif\",\"size\":15},\"tickcolor\":\"rgb(0, 0, 0)\",\"ticktext\":[\"-1.5\",\"-1.0\",\"-0.5\",\"0.0\",\"0.5\",\"1.0\",\"1.5\"],\"zeroline\":false,\"type\":\"-\",\"tickfont\":{\"color\":\"rgba(0, 0, 0, 1.000)\",\"family\":\"sans-serif\",\"size\":11},\"zerolinecolor\":\"rgba(0, 0, 0, 1.000)\",\"anchor\":\"y\"},\"paper_bgcolor\":\"rgba(255, 255, 255, 1.000)\",\"annotations\":[],\"height\":400,\"margin\":{\"l\":0,\"b\":20,\"r\":0,\"t\":20},\"plot_bgcolor\":\"rgba(255, 255, 255, 1.000)\",\"yaxis\":{\"gridwidth\":0.5,\"tickvals\":[0.0,1.0,2.0,3.0],\"visible\":true,\"ticks\":\"inside\",\"tickmode\":\"array\",\"domain\":[0.03762029746281716,0.9901574803149606],\"linecolor\":\"rgba(0, 0, 0, 1.000)\",\"showgrid\":true,\"title\":\"\",\"mirror\":false,\"tickangle\":0,\"showline\":true,\"gridcolor\":\"rgba(0, 0, 0, 0.100)\",\"titlefont\":{\"color\":\"rgba(0, 0, 0, 1.000)\",\"family\":\"sans-serif\",\"size\":15},\"tickcolor\":\"rgb(0, 0, 0)\",\"ticktext\":[\"0\",\"1\",\"2\",\"3\"],\"zeroline\":false,\"type\":\"-\",\"tickfont\":{\"color\":\"rgba(0, 0, 0, 1.000)\",\"family\":\"sans-serif\",\"size\":11},\"zerolinecolor\":\"rgba(0, 0, 0, 1.000)\",\"anchor\":\"x\"},\"legend\":{\"bordercolor\":\"rgba(0, 0, 0, 1.000)\",\"bgcolor\":\"rgba(255, 255, 255, 1.000)\",\"font\":{\"color\":\"rgba(0, 0, 0, 1.000)\",\"family\":\"sans-serif\",\"size\":11},\"y\":1.0,\"x\":1.0},\"width\":600});\n",
       "    </script>\n"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots\n",
    "plot(g,-2,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The only problem with the plotting solution is that you have to determine a domain for the plot (e.g. I plotted for x between -2 and 2 above). In principal, there could be solutions out side this interval. In many economic applications that is not much of a problem because the variables have natural lower and upper bounds that you can use, e.g. a consumer can buy not less than zero units and not more units than his budget divided by the price of a good."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving systems of equations using NLsolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving systems of equations requires a somewhat more involved syntax. I will take an example system consisting of 2 equations in two variables:\n",
    "$$y-z=5$$\n",
    "$$y+z=3$$\n",
    "which has the solution $y=4$ and $z=-1$. The first step is to rearrange this such that we can write this system as $f(x)=0$ where $x$ is a vector with two components (the first corresponding to $y$ and the second to $z$):\n",
    "$$f(x)=[x[1]-x[2]-5,x[1]+x[2]-3]$$\n",
    "This is similar to multivariate functions that we looked at before with the slight difference that the function $f$ now returns a vector where the first component corresponds to the first and the second component to the second equation. Unfortunately, NLsolve wants us not to enter this $f$ function as it is into julia but wants us to use a somewhat different form. I input it below and explain then what it does: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "f! (generic function with 1 method)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function f!(x,f)\n",
    "    f[1] = x[1]-x[2]-5\n",
    "    f[2] = x[1]+x[2]-3\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function $f!$ is a bit different from what we have seen so far: It takes as an argument a vector x that has two elements and a vector f that has also two elements. It then overwrites the first element of the f vector with the value x[1]-x[2]-5 that correspons to our first equation. Then it overwrites the second element of the f vector with the value x[1]+x[2]-3 that correspons to our second equation. The function does not return anything! It over writes f in place instead. It is a convention in julia to let the name of such an overwriting function end with a \"!\" (read \"f bang\").\n",
    "\n",
    "NLsolve wants us to write our system of equations in this form and can solve it then. But first we have to install the NLsolve package in case you have not done so earlier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mCloning cache of CommonSubexpressions from https://github.com/rdeits/CommonSubexpressions.jl.git\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mUpdating cache of DiffBase...\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mCloning cache of DiffResults from https://github.com/JuliaDiff/DiffResults.jl.git\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mCloning cache of DiffRules from https://github.com/JuliaDiff/DiffRules.jl.git\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mUpdating cache of ForwardDiff...\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mUpdating cache of LineSearches...\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mCloning cache of NLSolversBase from https://github.com/JuliaNLSolvers/NLSolversBase.jl.git\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mCloning cache of NLsolve from https://github.com/JuliaNLSolvers/NLsolve.jl.git\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mUpdating cache of NaNMath...\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mCloning cache of Parameters from https://github.com/mauro3/Parameters.jl.git\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mCloning cache of StaticArrays from https://github.com/JuliaArrays/StaticArrays.jl.git\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling Calculus v0.2.2\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling CommonSubexpressions v0.0.1\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling DiffBase v0.3.2\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling DiffResults v0.0.1\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling DiffRules v0.0.1\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling ForwardDiff v0.6.0\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling LineSearches v3.2.0\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling NLSolversBase v3.1.0\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling NLsolve v0.12.1\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling NaNMath v0.2.6\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling Parameters v0.8.0\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mInstalling StaticArrays v0.6.6\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mBuilding SpecialFunctions\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage database updated\n",
      "\u001b[39m"
     ]
    }
   ],
   "source": [
    "Pkg.add(\"NLsolve\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we tell julia to use NLsolve and tell NLsolve to solve our equation system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Results of Nonlinear Solver Algorithm\n",
       " * Algorithm: Trust-region with dogleg and autoscaling\n",
       " * Starting Point: [0.0, 0.0]\n",
       " * Zero: [4.0, -1.0]\n",
       " * Inf-norm of residuals: 0.000000\n",
       " * Iterations: 3\n",
       " * Convergence: true\n",
       "   * |x - x'| < 0.0e+00: false\n",
       "   * |f(x)| < 1.0e-08: true\n",
       " * Function Calls (f): 4\n",
       " * Jacobian Calls (df/dx): 4"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using NLsolve\n",
    "nlsolve(f!,[0.0,0.0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first argument for NLsolve is the \"bang function\" representing our equation system. The second argument is a starting point for the algorithm, i.e. an initial approximate guess of the solution. The solution of our problem is given in the output field \"Zero\"."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.7.0",
   "language": "julia",
   "name": "julia-0.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}