{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"toc": true
},
"source": [
"\n",
"# Numerical Routines: SciPy and NumPy\n",
"\n",
"SciPy is a Python library of mathematical routines. Many of the SciPy\n",
"routines are Python \"wrappers\", that is, Python routines that provide a\n",
"Python interface for numerical libraries and routines originally written\n",
"in Fortran, C, or C++. Thus, SciPy lets you take advantage of the\n",
"decades of work that has gone into creating and optimizing numerical\n",
"routines for science and engineering. Because the Fortran, C, or C++\n",
"code that Python accesses is compiled, these routines typically run very\n",
"fast. Therefore, there is no real downside---no speed penalty---for\n",
"using Python in these cases.\n",
"\n",
"single: SciPy\n",
"\n",
"We have already encountered one of SciPy's routines,\n",
"`scipy.optimize.leastsq`, for fitting nonlinear functions to\n",
"experimental data, which was introduced in the the chapter on `chap8`.\n",
"Here we will provide a further introduction to a number of other SciPy\n",
"packages, in particular those on special functions, numerical\n",
"integration, including routines for numerically solving ordinary\n",
"differential equations (ODEs), discrete Fourier transforms, linear\n",
"algebra, and solving non-linear equations. Our introduction to these\n",
"capabilities does not include extensive background on the numerical\n",
"methods employed; that is a topic for another text. Here we simply\n",
"introduce the SciPy routines for performing some of the more frequently\n",
"required numerical tasks.\n",
"\n",
"One final note: SciPy makes extensive use of NumPy arrays, so NumPy\n",
"should always be imported with SciPy\n",
"\n",
"single: SciPy; special functions\n",
"\n",
"Special functions\n",
"-----------------\n",
"\n",
"SciPy provides a plethora of special functions, including Bessel\n",
"functions (and routines for finding their zeros, derivatives, and\n",
"integrals), error functions, the gamma function, Legendre, Laguerre, and\n",
"Hermite polynomials (and other polynomial functions), Mathieu functions,\n",
"many statistical functions, and a number of other functions. Most are\n",
"contained in the `scipy.special` library, and each has its own special\n",
"arguments and syntax, depending on the vagaries of the particular\n",
"function. We demonstrate a number of them in the code below that\n",
"produces a plot of the different functions called. For more information,\n",
"you should consult the SciPy web site on the `scipy.special` library.\n",
"\n",
"> \n",
"> Plots of a few selected special functions\n",
"> \n",
"\n",
"``` python\n",
"import numpy as np\n",
"import scipy.special\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# create a figure window\n",
"fig = plt.figure(1, figsize=(9,8))\n",
"\n",
"# create arrays for a few Bessel functions and plot them\n",
"x = np.linspace(0, 20, 256)\n",
"j0 = scipy.special.jn(0, x)\n",
"j1 = scipy.special.jn(1, x)\n",
"y0 = scipy.special.yn(0, x)\n",
"y1 = scipy.special.yn(1, x)\n",
"ax1 = fig.add_subplot(321)\n",
"ax1.plot(x,j0, x,j1, x,y0, x,y1)\n",
"ax1.axhline(color=\"grey\", ls=\"--\", zorder=-1)\n",
"ax1.set_ylim(-1,1)\n",
"ax1.text(0.5, 0.95,'Bessel', ha='center', va='top',\n",
" transform = ax1.transAxes)\n",
"\n",
"# gamma function\n",
"x = np.linspace(-3.5, 6., 3601)\n",
"g = scipy.special.gamma(x)\n",
"g = np.ma.masked_outside(g, -100, 400)\n",
"ax2 = fig.add_subplot(322)\n",
"ax2.plot(x,g)\n",
"ax2.set_xlim(-3.5, 6)\n",
"ax2.axhline(color=\"grey\", ls=\"--\", zorder=-1)\n",
"ax2.axvline(color=\"grey\", ls=\"--\", zorder=-1)\n",
"ax2.set_ylim(-20, 100)\n",
"ax2.text(0.5, 0.95,'Gamma', ha='center', va='top',\n",
" transform = ax2.transAxes)\n",
"\n",
"# error function\n",
"x = np.linspace(0, 2.5, 256)\n",
"ef = scipy.special.erf(x)\n",
"ax3 = fig.add_subplot(323)\n",
"ax3.plot(x,ef)\n",
"ax3.set_ylim(0,1.1)\n",
"ax3.text(0.5, 0.95,'Error', ha='center', va='top',\n",
" transform = ax3.transAxes)\n",
"\n",
"# Airy function\n",
"x = np.linspace(-15, 4, 256)\n",
"ai, aip, bi, bip = scipy.special.airy(x)\n",
"ax4 = fig.add_subplot(324)\n",
"ax4.plot(x,ai, x,bi)\n",
"ax4.axhline(color=\"grey\", ls=\"--\", zorder=-1)\n",
"ax4.axvline(color=\"grey\", ls=\"--\", zorder=-1)\n",
"ax4.set_xlim(-15,4)\n",
"ax4.set_ylim(-0.5,0.6)\n",
"ax4.text(0.5, 0.95,'Airy', ha='center', va='top',\n",
" transform = ax4.transAxes)\n",
"\n",
"# Legendre polynomials\n",
"x = np.linspace(-1, 1, 256)\n",
"lp0 = np.polyval(scipy.special.legendre(0),x)\n",
"lp1 = np.polyval(scipy.special.legendre(1),x)\n",
"lp2 = np.polyval(scipy.special.legendre(2),x)\n",
"lp3 = np.polyval(scipy.special.legendre(3),x)\n",
"ax5 = fig.add_subplot(325)\n",
"ax5.plot(x,lp0, x,lp1, x,lp2, x,lp3)\n",
"ax5.axhline(color=\"grey\", ls=\"--\", zorder=-1)\n",
"ax5.axvline(color=\"grey\", ls=\"--\", zorder=-1)\n",
"ax5.set_ylim(-1,1.1)\n",
"ax5.text(0.5, 0.9,'Legendre', ha='center', va='top',\n",
" transform = ax5.transAxes)\n",
"\n",
"# Laguerre polynomials\n",
"x = np.linspace(-5, 8, 256)\n",
"lg0 = np.polyval(scipy.special.laguerre(0),x)\n",
"lg1 = np.polyval(scipy.special.laguerre(1),x)\n",
"lg2 = np.polyval(scipy.special.laguerre(2),x)\n",
"lg3 = np.polyval(scipy.special.laguerre(3),x)\n",
"ax6 = fig.add_subplot(326)\n",
"ax6.plot(x,lg0, x,lg1, x,lg2, x,lg3)\n",
"ax6.axhline(color=\"grey\", ls=\"--\", zorder=-1)\n",
"ax6.axvline(color=\"grey\", ls=\"--\", zorder=-1)\n",
"ax6.set_xlim(-5,8)\n",
"ax6.set_ylim(-5,10)\n",
"ax6.text(0.5, 0.9,'Laguerre', ha='center', va='top',\n",
" transform = ax6.transAxes)\n",
"\n",
"plt.show()\n",
"```\n",
"\n",
"The arguments of the different functions depend, of course, on the\n",
"nature of the particular function. For example, the first argument of\n",
"the two types of Bessel functions called in lines 10-13 is the so-called\n",
"*order* of the Bessel function, and the second argument is the\n",
"independent variable. The Gamma and Error functions take one argument\n",
"each and produce one output. The Airy function takes only one input\n",
"argument, but returns four outputs, which correspond the two Airy\n",
"functions, normally designated $\\mathrm{Ai}(x)$ and $\\mathrm{Bi}(x)$,\n",
"and their derivatives $\\mathrm{Ai}^\\prime(x)$ and\n",
"$\\mathrm{Bi}^\\prime(x)$. The plot shows only $\\mathrm{Ai}(x)$ and\n",
"$\\mathrm{Bi}(x)$.\n",
"\n",
"The polynomial functions shown have a special syntax that uses NumPy's\n",
"`polyval` function for generating polynomials. If `p` is a list or array\n",
"of `N` numbers and `x` is an array, then\n",
"\n",
"``` python\n",
"polyval(p, x) = p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x +\n",
" p[N-1]\n",
"```\n",
"\n",
"For example, if `p = [2.0, 5.0, 1.0]`, `polyval(p, x)` generates the\n",
"following quadratic polynomial: $2x^2 + 5x +1$.\n",
"\n",
"SciPy's `special.legendre(n)` and `special.laguerre(n)` functions output\n",
"the coefficients `p` needed in `polyval` to produce the\n",
"$n^\\mathrm{th}$-order Legendre and Laguerre polynomials, respectively.\n",
"The `scipy.special` library has functions that specify many other\n",
"polynomial functions in this same way.\n",
"\n",
"single: SciPy; numerical integration single: numerical integration;\n",
"single integrals\n",
"\n",
"Numerical integration\n",
"-------------------------\n",
"\n",
"When a function cannot be integrated analytically, or is very difficult\n",
"to integrate analytically, one generally turns to numerical integration\n",
"methods. SciPy has a number of routines for performing numerical\n",
"integration. Most of them are found in the same `scipy.integrate`\n",
"library. We list them here for reference.\n",
"\n",
"+---------------+-------------------------------------------------------+\n",
"| **Function** | **Description** |\n",
"+===============+=======================================================+\n",
"| `quad` | single integration |\n",
"+---------------+-------------------------------------------------------+\n",
"| `dblquad` | double integration |\n",
"+---------------+-------------------------------------------------------+\n",
"| `tplquad` | triple integration |\n",
"+---------------+-------------------------------------------------------+\n",
"| `nquad` | $n$-fold multiple integration |\n",
"+---------------+-------------------------------------------------------+\n",
"| `fixed_quad` | Gaussian quadrature, order n |\n",
"+---------------+-------------------------------------------------------+\n",
"| `quadrature` | Gaussian quadrature to tolerance |\n",
"+---------------+-------------------------------------------------------+\n",
"| `romberg` | Romberg integration |\n",
"+---------------+-------------------------------------------------------+\n",
"+---------------+-------------------------------------------------------+\n",
"| `trapz` | trapezoidal rule |\n",
"+---------------+-------------------------------------------------------+\n",
"| `cumtrapz` | trapezoidal rule to cumulatively compute integral |\n",
"+---------------+-------------------------------------------------------+\n",
"| `simps` | Simpson's rule |\n",
"+---------------+-------------------------------------------------------+\n",
"| `romb` | Romberg integration |\n",
"+---------------+-------------------------------------------------------+\n",
"+---------------+-------------------------------------------------------+\n",
"| `polyint` | Analytical polynomial integration (NumPy) |\n",
"+---------------+-------------------------------------------------------+\n",
"| `poly1d` | Helper function for `polyint` (NumPy) |\n",
"+---------------+-------------------------------------------------------+\n",
"\n",
"### Single integrals\n",
"\n",
"The function `quad` is the workhorse of SciPy's integration functions.\n",
"Numerical integration is sometimes called *quadrature*, hence the name.\n",
"It is normally the default choice for performing single integrals of a\n",
"function $f(x)$ over a given fixed range from $a$ to $b$\n",
"\n",
"$$\\int_a^b f(x)\\, dx$$\n",
"\n",
"The general form of `quad` is `scipy.integrate.quad(f, a, b)`, where `f`\n",
"is the name of the function to be integrated and `a` and `b` are the\n",
"lower and upper limits, respectively. The routine uses *adaptive\n",
"quadrature* methods to numerically evaluate integrals, meaning it\n",
"successively refines the subintervals (makes them smaller) until a\n",
"desired level of numerical precision is achieved. For the `quad`\n",
"routine, this is about $10^{-8}$, although it usually does even better.\n",
"\n",
"As an example, let's integrate a Gaussian function over the range from 0\n",
"to 1\n",
"\n",
"$$\\int_0^1 e^{-x^2} dx$$\n",
"\n",
"We first need to define the function $f(x)=e^{-x^2}$, which we do using\n",
"a lambda expression, and then we call the function `quad` to perform the\n",
"integration.\n",
"\n",
"``` ipython\n",
"In [1]: import scipy.integrate\n",
"\n",
"In [2]: f = lambda x : exp(-x**2)\n",
"\n",
"In [3]: scipy.integrate.quad(f, 0, 1)\n",
"Out[3]: (0.7468241328124271, 8.291413475940725e-15)\n",
"```\n",
"\n",
"The function call `scipy.integrate.quad(f, 0, 1)` returns two numbers.\n",
"The first is `0.7468...`, which is the value of the integral, and the\n",
"second is `8.29...e-15`, which is an estimate of the absolute error in\n",
"the value of the integral, which we see is quite small compared to\n",
"`0.7468`.\n",
"\n",
"Because `quad` requires a function *name* as its first argument, we\n",
"can't simply use the expression `exp(-x**2)`. On the other hand, we\n",
"could use the usual `def` statement to create a normal function, and\n",
"then use the name of that function in `quad`. However, it's simpler here\n",
"to use a lambda expression. In fact, we can just put the lambda\n",
"expression directly into the first argument, as illustrated here\n",
"\n",
"``` ipython\n",
"In [4]: scipy.integrate.quad(lambda x : exp(-x**2), 0, 1)\n",
"Out[4]: (0.7468241328124271, 8.291413475940725e-15)\n",
"```\n",
"\n",
"That works too! Thus we see a `lambda` expression used as an *anonymous\n",
"function*, a function with no name, as promised in the section `lambda`.\n",
"\n",
"Note\n",
"\n",
"The `quad` function accepts positive and negative infinity as limits.\n",
"\n",
"``` ipython\n",
"In [5]: scipy.integrate.quad(lambda x : exp(-x**2), 0, inf)\n",
"Out[5]: (0.8862269254527579, 7.101318390472462e-09)\n",
"\n",
"In [6]: scipy.integrate.quad(lambda x : exp(-x**2), -inf, 1)\n",
"Out[6]: (1.6330510582651852, 3.669607414547701e-11)\n",
"```\n",
"\n",
"The `quad` function handles infinite limits just fine. The absolute\n",
"errors are somewhat larger but still well within acceptable bounds for\n",
"practical work.\n",
"\n",
"The `quad` function can integrate standard predefined NumPy functions of\n",
"a single variable, like `exp`, `sin`, and `cos`.\n",
"\n",
"``` ipython\n",
"In [7]: scipy.integrate.quad(exp, 0, 1)\n",
"Out[7]: (1.7182818284590453, 1.9076760487502457e-14)\n",
"\n",
"In [8]: scipy.integrate.quad(sin, -0.5, 0.5)\n",
"Out[8]: (0.0, 2.707864644566304e-15)\n",
"\n",
"In [9]: scipy.integrate.quad(cos, -0.5, 0.5)\n",
"Out[9]: (0.9588510772084061, 1.0645385431034061e-14)\n",
"```\n",
"\n",
"Let's integrate the first order Bessel function of the first kind,\n",
"usually denoted $J_1(x)$, over the interval from 0 to 5. Here is how we\n",
"do it, using `scipy.special.jn(v,x)` where `v` is the (real) order of\n",
"the Bessel function:\n",
"\n",
"``` ipython\n",
"In [10]: import scipy.special\n",
"\n",
"In [11]: scipy.integrate.quad(lambda x: scipy.special.jn(1,x),0,5)\n",
"Out[11]: (1.177596771314338, 1.8083362065765924e-14)\n",
"```\n",
"\n",
"Because the SciPy function `scipy.special.jn(v, x)` is a function of two\n",
"variables, `v` and `x`, we cannot use the function name\n",
"`scipy.special.jn` in `quad`. So we use a `lambda` expression, which is\n",
"a function of only one variable, `x`, because we have set the `v`\n",
"argument equal to 1.\n",
"\n",
"single: numerical integration; integrals of polynomials\n",
"\n",
"#### Integrating polynomials\n",
"\n",
"Working in concert with the NumPy `poly1d`, the NumPy function `polyint`\n",
"takes the $n^\\mathrm{th}$ antiderivative of a polynomial and can be used\n",
"to evaluate definite integrals. The function `poly1d` essentially does\n",
"the same thing as `polyval` that we encountered in the section\n",
"`specFunc`, but with a different syntax. Suppose we want to make the\n",
"polynomial function $p(x) = 2x^2 + 5x +1$. Then we write\n",
"\n",
"``` ipython\n",
"In [12]: p = np.poly1d([2, 5, 1])\n",
"\n",
"In [13]: p\n",
"Out[13]: poly1d([2, 5, 1])\n",
"```\n",
"\n",
"The polynomial $p(x) = 2x^2 + 5x +1$ is evaluated using the syntax\n",
"`p(x)`. Below, we evaluate the polynomial at three different values of\n",
"`x`.\n",
"\n",
"``` ipython\n",
"In [14]: p(1), p(2), p(3.5)\n",
"Out[14]: (8, 19, 43.0)\n",
"```\n",
"\n",
"Thus `polyval` allows us to define the function $p(x) = 2x^2 + 5x +1$.\n",
"Now the antiderivative of $p(x) = 2x^2 + 5x +1$ is\n",
"$P(x) = \\frac{2}{3}x^3 + \\frac{5}{2}x^2 +x+C$ where $C$ is the\n",
"integration constant. The NumPy function `polyint`, which takes the\n",
"$n^\\mathrm{th}$ antiderivative of a polynomial, works as follows\n",
"\n",
"``` ipython\n",
"In [15]: P = polyint(p)\n",
"\n",
"In [16]: P\n",
"Out[16]: poly1d([ 0.66666667, 2.5 , 1. , 0. ])\n",
"```\n",
"\n",
"When `polyint` has a single input, `p` is this case, `polyint` returns\n",
"the coefficients of the antiderivative with the integration constant set\n",
"to zero, as `Out[16]` illustrates. It is then an easy matter to\n",
"determine any definite integral of the polynomial $p(x) = 2x^2 + 5x +1$\n",
"since\n",
"\n",
"$$q \\equiv \\int_a^b p(x)\\, dx = P(b) - P(a) \\;.$$\n",
"\n",
"For example, if $a=1$ and $b=5$,\n",
"\n",
"``` ipython\n",
"In [17]: q=P(5)-P(1)\n",
"\n",
"In [18]: q\n",
"Out[18]: 146.66666666666666\n",
"```\n",
"\n",
"or\n",
"\n",
"$$\\int_1^5 \\left(2x^2 + 5x +1\\right)\\, dx = 146\\tfrac{2}{3} \\;.$$\n",
"\n",
"single: numerical integration; double integrals\n",
"\n",
"### Double integrals\n",
"\n",
"The `scipy.integrate` function `dblquad` can be used to numerically\n",
"evaluate double integrals of the form\n",
"\n",
"$$\\int_{y=a}^{y=b} dy \\int_{x=g(y)}^{x=h(y)} dx\\,f(x,y)$$\n",
"\n",
"The general form of `dblquad` is\n",
"\n",
"``` ipython\n",
"scipy.integrate.dblquad(func, a, b, gfun, hfun)\n",
"```\n",
"\n",
"where `func` if the name of the function to be integrated, `a` and `b`\n",
"are the lower and upper limits of the `x` variable, respectively, and\n",
"`gfun` and `hfun` are the *names* of the functions that define the lower\n",
"and upper limits of the `y` variable.\n",
"\n",
"As an example, let's perform the double integral\n",
"\n",
"$$\\int_0^{1/2} dy \\int_0^{\\sqrt{1-4y^2}} 16xy\\, dx$$\n",
"\n",
"We define the functions f, g, and h, using\n",
"lambda expressions. Note that even if g,\n",
"and h are constants, as they may be in\n",
"many cases, they must be defined as functions, as we have done here for\n",
"the lower limit.\n",
"\n",
"``` ipython\n",
"In [19]: f = lambda x, y : 16*x*y\n",
"\n",
"In [20]: g = lambda x : 0\n",
"\n",
"In [21]: h = lambda y : sqrt(1-4*y**2)\n",
"\n",
"In [22]: scipy.integrate.dblquad(f, 0, 0.5, g, h)\n",
"Out[22]: (0.5, 5.551115123125783e-15)\n",
"```\n",
"\n",
"Once again, there are two outputs: the first is the value of the\n",
"integral and the second is its absolute uncertainty.\n",
"\n",
"Of course, the lower limit can also be a function of $y$, as we\n",
"demonstrate here by performing the integral\n",
"\n",
"$$\\int_0^{1/2} dy \\int_{1-2y}^{\\sqrt{1-4y^2}} 16xy\\, dx$$\n",
"\n",
"The code for this is given by\n",
"\n",
"``` ipython\n",
"In [23]: g = lambda y : 1-2*y\n",
"\n",
"In [24]: scipy.integrate.dblquad(f, 0, 0.5, g, h)\n",
"Out[24]: (0.33333333333333326, 3.700743415417188e-15)\n",
"```\n",
"\n",
"#### Other integration routines\n",
"\n",
"In addition to the routines described above, `scipy.integrate` has a\n",
"number of other integration routines, including `nquad`, which performs\n",
"$n$-fold multiple integration, as well as other routines that implement\n",
"other integration algorithms. You will find, however, that `quad` and\n",
"`dblquad` meet most of your needs for numerical integration.\n",
"\n",
"single: ODEs; numerical solutions single: SciPy; ODEs\n",
"\n",
"Solving ODEs\n",
"------------\n",
"\n",
"The `scipy.integrate` library has two powerful powerful routines, `ode`\n",
"and `odeint`, for numerically solving systems of coupled first order\n",
"ordinary differential equations (ODEs). While `ode` is more versatile,\n",
"`odeint` (ODE integrator) has a simpler Python interface works very well\n",
"for most problems. It can handle both stiff and non-stiff problems. Here\n",
"we provide an introduction to `odeint`.\n",
"\n",
"A typical problem is to solve a second or higher order ODE for a given\n",
"set of initial conditions. Here we illustrate using `odeint` to solve\n",
"the equation for a driven damped pendulum. The equation of motion for\n",
"the angle $\\theta$ that the pendulum makes with the vertical is given by\n",
"\n",
"$$\\frac{d^2\\theta}{dt^2} = -\\frac{1}{Q} \\frac{d\\theta}{dt} +\n",
" \\sin\\theta + d \\cos\\Omega t$$\n",
"\n",
"where $t$ is time, $Q$ is the quality factor, $d$ is the forcing\n",
"amplitude, and $\\Omega$ is the driving frequency of the forcing. Reduced\n",
"variables have been used such that the natural (angular) frequency of\n",
"oscillation is 1. The ODE is nonlinear owing to the $\\sin\\theta$ term.\n",
"Of course, it's precisely because there are no general methods for\n",
"solving nonlinear ODEs that one employs numerical techniques, so it\n",
"seems appropriate that we illustrate the method with a nonlinear ODE.\n",
"\n",
"The first step is always to transform any $n^\\mathrm{th}$-order ODE into\n",
"a system of $n$ first order ODEs of the form:\n",
"\n",
"$$\\begin{aligned}\n",
"\\frac{dy_1}{dt} &= f_1(t, y_1, ..., y_n) \\\\\n",
"\\frac{dy_2}{dt} &= f_2(t, y_1, ..., y_n) \\\\\n",
" \\vdots\\quad &= \\quad\\vdots \\\\\n",
"\\frac{dy_n}{dt} &= f_n(t, y_1, ..., y_n) \\;.\n",
"\\end{aligned}$$\n",
"\n",
"We also need $n$ initial conditions, one for each variable $y_i$. Here\n",
"we have a second order ODE so we will have two coupled ODEs and two\n",
"initial conditions.\n",
"\n",
"We start by transforming our second order ODE into two coupled first\n",
"order ODEs. The transformation is easily accomplished by defining a new\n",
"variable $\\omega \\equiv d\\theta/dt$. With this definition, we can\n",
"rewrite our second order ODE as two coupled first order ODEs:\n",
"\n",
"$$\\begin{aligned}\n",
"\\frac{d\\theta}{dt} &= \\omega \\\\\n",
"\\frac{d\\omega}{dt} &= -\\frac{1}{Q}\\,\\omega + \\sin\\theta \n",
" + d \\cos\\Omega t \\;.\n",
"\\end{aligned}$$\n",
"\n",
"In this case the functions on the right hand side of the equations are\n",
"\n",
"$$\\begin{aligned}\n",
"f_1(t, \\theta, \\omega) &= \\omega \\\\\n",
"f_2(t, \\theta, \\omega) &= -\\frac{1}{Q}\\,\\omega + \\sin\\theta \n",
" + d \\cos\\Omega t \\;.\n",
"\\end{aligned}$$\n",
"\n",
"Note that there are no explicit derivatives on the right hand side of\n",
"the functions $f_i$; they are all functions of $t$ and the various\n",
"$y_i$, in this case $\\theta$ and $\\omega$.\n",
"\n",
"The initial conditions specify the values of $\\theta$ and $\\omega$ at\n",
"$t=0$.\n",
"\n",
"SciPy's ODE solver `scipy.integrate.odeint` has three required arguments\n",
"and many optional keyword arguments, of which we only need one, `args`,\n",
"for this example. So in this case, `odeint` has the form\n",
"\n",
"``` ipython\n",
"odeint(func, y0, t, args=())\n",
"```\n",
"\n",
"The first argument `func` is the name of a Python function that returns\n",
"a list of values of the $n$ functions $f_i(t, y_1, ..., y_n)$ at a given\n",
"time $t$. The second argument `y0` is an array (or list) of the values\n",
"of the initial conditions of $y_1, ..., y_n)$. The third argument is the\n",
"array of times at which you want `odeint` to return the values of\n",
"$y_1, ..., y_n)$. The keyword argument `args` is a tuple that is used to\n",
"pass parameters (besides `y0` and `t`) that are needed to evaluate\n",
"`func`. Our example should make all of this clear.\n",
"\n",
"After having written the $n^\\mathrm{th}$-order ODE as a system of $n$\n",
"first-order ODEs, the next task is to write the function `func`. The\n",
"function `func` should have three arguments: (1) the list (or array) of\n",
"current `y` values, the current time `t`, and a list of any other\n",
"parameters `params` needed to evaluate `func`. The function `func`\n",
"returns the values of the derivatives $dy_i/dt = f_i(t, y_1, ..., y_n)$\n",
"in a list (or array). Lines 5-11 illustrate how to write `func` for our\n",
"example of a driven damped pendulum. Here we name the function simply\n",
"`f`, which is the name that appears in the call to `odeint` in line 33\n",
"below.\n",
"\n",
"The only other tasks remaining are to define the parameters needed in\n",
"the function, bundling them into a list (see line 22 below), and to\n",
"define the initial conditions, and bundling them into another list (see\n",
"line 25 below). After defining the time array in lines 28-30, the only\n",
"remaining task is to call `odeint` with the appropriate arguments and a\n",
"variable, `psoln` in this case to store output. The output `psoln` is an\n",
"$n$ element array where each element is itself an array corresponding\n",
"the the values of $y_i$ for each time in the time `t` array that was an\n",
"argument of `odeint`. For this example, the first element `psoln[:,0]`\n",
"is the $y_0$ or `theta` array, and the second element `psoln[:,1]` is\n",
"the $y_1$ or `omega` array. The remainder of the code simply plots out\n",
"the results in different formats. The resulting plots are shown in the\n",
"figure `fig:odePend` after the code.\n",
"\n",
"``` python\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy.integrate import odeint\n",
"\n",
"def f(y, t, params):\n",
" theta, omega = y # unpack current values of y\n",
" Q, d, Omega = params # unpack parameters\n",
" derivs = [omega, # list of dy/dt=f functions\n",
" -omega/Q + np.sin(theta) + d*np.cos(Omega*t)]\n",
" return derivs\n",
"\n",
"# Parameters\n",
"Q = 2.0 # quality factor (inverse damping)\n",
"d = 1.5 # forcing amplitude\n",
"Omega = 0.65 # drive frequency\n",
"\n",
"# Initial values\n",
"theta0 = 0.0 # initial angular displacement\n",
"omega0 = 0.0 # initial angular velocity\n",
"\n",
"# Bundle parameters for ODE solver\n",
"params = [Q, d, Omega]\n",
"\n",
"# Bundle initial conditions for ODE solver\n",
"y0 = [theta0, omega0]\n",
"\n",
"# Make time array for solution\n",
"tStop = 200.\n",
"tInc = 0.05\n",
"t = np.arange(0., tStop, tInc)\n",
"\n",
"# Call the ODE solver\n",
"psoln = odeint(f, y0, t, args=(params,))\n",
"\n",
"# Plot results\n",
"fig = plt.figure(1, figsize=(8,8))\n",
"\n",
"# Plot theta as a function of time\n",
"ax1 = fig.add_subplot(311)\n",
"ax1.plot(t, psoln[:,0])\n",
"ax1.set_xlabel('time')\n",
"ax1.set_ylabel('theta')\n",
"\n",
"# Plot omega as a function of time\n",
"ax2 = fig.add_subplot(312)\n",
"ax2.plot(t, psoln[:,1])\n",
"ax2.set_xlabel('time')\n",
"ax2.set_ylabel('omega')\n",
"\n",
"# Plot omega vs theta\n",
"ax3 = fig.add_subplot(313)\n",
"twopi = 2.0*np.pi\n",
"ax3.plot(psoln[:,0]%twopi, psoln[:,1], '.', ms=1)\n",
"ax3.set_xlabel('theta')\n",
"ax3.set_ylabel('omega')\n",
"ax3.set_xlim(0., twopi)\n",
"\n",
"plt.tight_layout()\n",
"plt.show()\n",
"```\n",
"\n",
"\n",
"Pendulum trajectory\n",
"\n",
"\n",
"The plots above reveal that for the particular set of input parameters\n",
"chosen, `Q = 2.0`, `d = 1.5`, and `Omega = 0.65`, the pendulum\n",
"trajectories are chaotic. Weaker forcing (smaller $d$) leads to what is\n",
"perhaps the more familiar behavior of sinusoidal oscillations with a\n",
"fixed frequency which, at long times, is equal to the driving frequency.\n",
"\n",
"single: discrete Fourier transforms single: SciPy; discrete Fourier\n",
"transforms see: fast Fourier transforms; discrete Fourier transforms\n",
"see: FFTs; discrete Fourier transforms\n",
"\n",
"Discrete (fast) Fourier transforms\n",
"----------------------------------\n",
"\n",
"The SciPy library has a number of routines for performing discrete\n",
"Fourier transforms. Before delving into them, we provide a brief review\n",
"of Fourier transforms and discrete Fourier transforms.\n",
"\n",
"### Continuous and discrete Fourier transforms\n",
"\n",
"The Fourier transform of a function $g(t)$ is given by\n",
"\n",
"$$G(f) = \\int_{-\\infty}^\\infty g(t)\\, e^{-i\\, 2\\pi f t}\\, dt \\;,$$\n",
"\n",
"where $f$ is the Fourier transform variable; if $t$ is time, then $f$ is\n",
"frequency. The inverse transform is given by\n",
"\n",
"$$g(t) = \\int_{-\\infty}^\\infty G(f)\\, e^{i\\, 2\\pi ft}\\, df$$\n",
"\n",
"Here we define the Fourier transform in terms of the frequency $f$\n",
"rather than the angular frequency $\\omega = 2\\pi f$.\n",
"\n",
"The conventional Fourier transform is defined for continuous functions,\n",
"or at least for functions that are dense and thus have an infinite\n",
"number of data points. When doing numerical analysis, however, you work\n",
"with *discrete* data sets, that is, data sets defined for a finite\n",
"number of points. The discrete Fourier transform (DFT) is defined for a\n",
"function $g_n$ consisting of a set of $N$ discrete data points. Those\n",
"$N$ data points must be defined at *equally-spaced* times\n",
"$t_n=n\\Delta t$ where $\\Delta t$ is the time between successive data\n",
"points and $n$ runs from 0 to $N-1$. The discrete Fourier transform\n",
"(DFT) of $g_n$ is defined as\n",
"\n",
"$$G_l = \\sum_{n=0}^{N-1} g_n\\, e^{-i\\,(2\\pi/N)\\,ln}$$\n",
"\n",
"where $l$ runs from 0 to $N-1$. The inverse discrete Fourier transform\n",
"(iDFT) is defined as\n",
"\n",
"$$g_n = \\frac{1}{N} \\sum_{l=0}^{N-1} G_l\\, e^{i\\,(2\\pi/N)\\,ln} \\;.$$\n",
"\n",
"The DFT is usually implemented on computers using the well-known Fast\n",
"Fourier Transform (FFT) algorithm, generally credited to Cooley and\n",
"Tukey who developed it at AT&T Bell Laboratories during the 1960s. But\n",
"their algorithm is essentially one of many independent rediscoveries of\n",
"the basic algorithm dating back to Gauss who described it as early as\n",
"1805.\n",
"\n",
"### The SciPy FFT library\n",
"\n",
"The SciPy library `scipy.fftpack` has routines that implement a\n",
"souped-up version of the FFT algorithm along with many ancillary\n",
"routines that support working with DFTs. The basic FFT routine in\n",
"`scipy.fftpack` is appropriately named `fft`. The program below\n",
"illustrates its use, along with the plots that follow.\n",
"\n",
"``` python\n",
"import numpy as np\n",
"from scipy import fftpack\n",
"import matplotlib.pyplot as plt\n",
"\n",
"width = 2.0\n",
"freq = 0.5\n",
"\n",
"t = np.linspace(-10, 10, 101) # linearly space time array\n",
"g = np.exp(-np.abs(t)/width) * np.sin(2.0*np.pi*freq*t)\n",
"\n",
"dt = t[1]-t[0] # increment between times in time array\n",
"\n",
"G = fftpack.fft(g) # FFT of g\n",
"f = fftpack.fftfreq(g.size, d=dt) # frequenies f[i] of g[i]\n",
"f = fftpack.fftshift(f) # shift frequencies from min to max\n",
"G = fftpack.fftshift(G) # shift G order to coorespond to f\n",
"\n",
"fig = plt.figure(1, figsize=(8,6), frameon=False)\n",
"ax1 = fig.add_subplot(211)\n",
"ax1.plot(t, g)\n",
"ax1.set_xlabel('t')\n",
"ax1.set_ylabel('g(t)')\n",
"\n",
"ax2 = fig.add_subplot(212)\n",
"ax2.plot(f, np.real(G), color='dodgerblue', label='real part')\n",
"ax2.plot(f, np.imag(G), color='coral', label='imaginary part')\n",
"ax2.legend()\n",
"ax2.set_xlabel('f')\n",
"ax2.set_ylabel('G(f)')\n",
"\n",
"plt.show()\n",
"```\n",
"\n",
"\n",
"Function g(t) and its DFT G(f).\n",
"\n",
"\n",
"The DFT has real and imaginary parts, both of which are plotted in the\n",
"figure.\n",
"\n",
"The `fft` function returns the $N$ Fourier components of $G_n$ starting\n",
"with the zero-frequency component $G_0$ and progressing to the maximum\n",
"positive frequency component $G_{(N/2)-1}$ (or $G_{(N-1)/2}$ if $N$ is\n",
"odd). From there, `fft` returns the maximum *negative* component\n",
"$G_{N/2}$ (or $G_{(N-1)/2}$ if $N$ is odd) and continues upward in\n",
"frequency until it reaches the minimum negative frequency component\n",
"$G_{N-1}$. This is the standard way that DFTs are ordered by most\n",
"numerical DFT packages. The `scipy.fftpack` function `fftfreq` creates\n",
"the array of frequencies in this non-intuitive order such that `f[n]` in\n",
"the above routine is the correct frequency for the Fourier component\n",
"`G[n]`. The arguments of `fftfreq` are the size of the the orignal array\n",
"`g` and the keyword argument `d` that is the spacing between the\n",
"(equally spaced) elements of the time array (`d=1` if left unspecified).\n",
"The package `scipy.fftpack` provides the convenience function `fftshift`\n",
"that reorders the frequency array so that the zero-frequency occurs at\n",
"the middle of the array, that is, so the frequencies proceed\n",
"monotonically from smallest (most negative) to largest (most positive).\n",
"Applying `fftshift` to both `f` and `G` puts the frequencies `f` in\n",
"ascending order and shifts `G` so that the frequency of `G[n]` is given\n",
"by the shifted `f[n]`.\n",
"\n",
"The `scipy.fftpack` module also contains routines for performing\n",
"2-dimensional and $n$-dimensional DFTs, named `fft2` and `fftn`,\n",
"respectively, using the FFT algorithm.\n",
"\n",
"As for most FFT routines, the `scipy.fftpack` FFT routines are most\n",
"efficient if $N$ is a power of 2. Nevertheless, the FFT routines are\n",
"able to handle data sets where $N$ is not a power of 2.\n",
"\n",
"`scipy.fftpack` also supplies an inverse DFT function `ifft`. It is\n",
"written to act on the *unshifted* FFT so take care! Note also that\n",
"`ifft` returns a *complex* array. Because of machine roundoff error, the\n",
"imaginary part of the function returned by `ifft` will, in general, be\n",
"very near zero but not exactly zero even when the original function is a\n",
"purely real function.\n",
"\n",
"single: linear algebra single: SciPy; linear algebra\n",
"\n",
"Linear algebra\n",
"--------------\n",
"\n",
"Python's mathematical libraries, NumPy and SciPy, have extensive tools\n",
"for numerically solving problems in linear algebra. Here we focus on two\n",
"problems that arise commonly in scientific and engineering settings: (1)\n",
"solving a system of linear equations and (2) eigenvalue problems. In\n",
"addition, we also show how to perform a number of other basic\n",
"computations, such as finding the determinant of a matrix, matrix\n",
"inversion, and $LU$ decomposition. The SciPy package for linear algebra\n",
"is called `scipy.linalg`.\n",
"\n",
"### Basic computations in linear algebra\n",
"\n",
"SciPy has a number of routines for performing basic operations with\n",
"matrices. The determinant of a matrix is computed using the\n",
"`scipy.linalg.det` function:\n",
"\n",
"``` ipython\n",
"In [1]: import scipy.linalg\n",
"In [2]: a = array([[-2, 3], [4, 5]])\n",
"In [3]: a\n",
"Out[4]: array([[-2, 3],\n",
" [ 4, 5]])\n",
"\n",
"In [5]: scipy.linalg.det(a)\n",
"Out[5]: -22.0\n",
"```\n",
"\n",
"The inverse of a matrix is computed using the `scipy.linalg.inv`\n",
"function, while the product of two matrices is calculated using the\n",
"NumPy `dot` function:\n",
"\n",
"``` ipython\n",
"In [6]: b = scipy.linalg.inv(a)\n",
"\n",
"In [6]: b\n",
"Out[6]: array([[-0.22727273, 0.13636364],\n",
" [ 0.18181818, 0.09090909]])\n",
"\n",
"In [7]: dot(a,b)\n",
"Out[7]: array([[ 1., 0.],\n",
" [ 0., 1.]])\n",
"```\n",
"\n",
"single: linear algebra; solving systems of equations single: SciPy;\n",
"solving systems of equations\n",
"\n",
"### Solving systems of linear equations\n",
"\n",
"Solving systems of equations is nearly as simple as constructing a\n",
"coefficient matrix and a column vector. Suppose you have the following\n",
"system of linear equations to solve:\n",
"\n",
"$$\\begin{aligned}\n",
"2x_1 + 4x_2 + 6x_3 &= 4\\\\\n",
" x_1 - 3x_2 - 9x_3 &= -11\\\\\n",
"8x_1 + 5x_2 - 7x_3 &= 1\\\\\n",
"\\end{aligned}$$\n",
"\n",
"The first task is to recast this set of equations as a matrix equation\n",
"of the form $\\mathsf{A}\\, \\mathbf{x} = \\mathbf{b}$. In this case, we\n",
"have:\n",
"\n",
"$$\\begin{aligned}\n",
"\\mathsf{A} = \\left(\\begin{array}{ccc}2 & 4 & 6 \\\\\n",
" 1 & -3 & -9 \\\\\n",
" 8 & 5 & -7 \\end{array}\\right)\n",
" \\;,\\quad\n",
"\\mathbf{x} = \\left(\\begin{array}{c}x_1 \\\\x_2 \\\\x_3\\end{array}\\right)\n",
" \\;,\\quad\n",
"\\mathbf{b} = \\left(\\begin{array}{c}4 \\\\-11 \\\\1\\end{array}\\right) \\;.\n",
"\\end{aligned}$$\n",
"\n",
"Next we construct the array $\\mathsf{A}$ and vector $\\mathbf{b}$ as\n",
"NumPy arrays:\n",
"\n",
"``` ipython\n",
"In [8]: A = array([[2, 4, 6], [1, -3, -9], [8, 5, -7]])\n",
"In [9]: b = array([4, -11, 2])\n",
"```\n",
"\n",
"Finally we use the SciPy function `scipy.linalg.solve` to find $x_1$,\n",
"$x_2$, and $x_3$.\n",
"\n",
"``` ipython\n",
"In [10]: scipy.linalg.solve(A,b)\n",
"Out[10]: array([ -8.91304348, 10.2173913 , -3.17391304])\n",
"```\n",
"\n",
"which gives the results: $x_1=-8.91304348$, $x_2= 10.2173913$, and\n",
"$x_3= -3.17391304$. Of course, you can get the same answer by noting\n",
"that $\\mathbf{x}=\\mathsf{A}^{-1}\\mathbf{b}$. Following this approach, we\n",
"can use the scipy.linalg.inv introduced\n",
"in the previous section:\n",
"\n",
"``` ipython\n",
"Ainv = scipy.linalg.inv(A)\n",
"\n",
"In [10]: dot(Ainv, b)\n",
"Out[10]: array([ -8.91304348, 10.2173913 , -3.17391304])\n",
"```\n",
"\n",
"which is the same answer we obtained using `scipy.linalg.solve`. Using\n",
"`scipy.linalg.solve` is numerically more stable and a faster than using\n",
"$\\mathbf{x}=\\mathsf{A}^{-1}\\mathbf{b}$, so it is the preferred method\n",
"for solving systems of equations.\n",
"\n",
"You might wonder what happens if the system of equations are not all\n",
"linearly independent. For example if the matrix $\\mathsf{A}$ is given by\n",
"\n",
"$$\\begin{aligned}\n",
"\\mathsf{A} = \\left(\\begin{array}{ccc}2 & 4 & 6 \\\\\n",
" 1 & -3 & -9 \\\\\n",
" 1 & 2 & 3 \\end{array}\\right)\n",
"\\end{aligned}$$\n",
"\n",
"where the third row is a multiple of the first row. Let's try it out and\n",
"see what happens. First we change the bottom row of the matrix\n",
"$\\mathsf{A}$ and then try to solve the system as we did before.\n",
"\n",
"``` ipython\n",
"In [11]: A[2] = array([1, 2, 3])\n",
"\n",
"In [12]: A\n",
"Out[12]: array([[ 2, 4, 6],\n",
" [ 1, -3, -9],\n",
" [ 1, 2, 3]])\n",
"\n",
"In [13]: scipy.linalg.solve(A,b)\n",
"LinAlgError: Singular matrix\n",
"\n",
"In [14]: Ainv = scipy.linalg.inv(A)\n",
"LinAlgError: Singular matrix\n",
"```\n",
"\n",
"Whether we use `scipy.linalg.solve` or `scipy.linalg.inv`, SciPy raises\n",
"an error because the matrix is singular.\n",
"\n",
"single: linear algebra; eigenvalue problems single: eigenvalue problems\n",
"single: anonymous functions; lambda expressions single: lambda\n",
"expressions\n",
"\n",
"### Eigenvalue problems\n",
"\n",
"One of the most common problems in science and engineering is the\n",
"eigenvalue problem, which in matrix form is written as\n",
"\n",
"$$\\mathsf{A}\\mathbf{x} = \\lambda \\mathbf{x}$$\n",
"\n",
"where $\\mathsf{A}$ is a square matrix, $\\mathbf{x}$ is a column vector,\n",
"and $\\lambda$ is a scalar (number). Given the matrix $\\mathsf{A}$, the\n",
"problem is to find the set of eigenvectors $\\mathbf{x}$ and their\n",
"corresponding eigenvalues $\\lambda$ that solve this equation.\n",
"\n",
"We can solve eigenvalue equations like this using `scipy.linalg.eig`.\n",
"the outputs of this function is an array whose entries are the\n",
"eigenvalues and a matrix whose rows are the eigenvectors. Let's return\n",
"to the matrix we were using previously and find its eigenvalues and\n",
"eigenvectors.\n",
"\n",
"``` ipython\n",
"A = array([[2, 4, 6],[1, -3, -9],[8, 5, -7]])\n",
"\n",
"In [15]: A\n",
"Out[15]: array([[ 2, 4, 6],\n",
" [ 1, -3, -9],\n",
" [ 8, 5, -7]])\n",
"\n",
"In [16]: lam, evec = scipy.linalg.eig(A)\n",
"\n",
"In [17]: lam\n",
"Out[17]: array([ 2.40995356+0.j, -8.03416016+0.j,\n",
" -2.37579340+0.j])\n",
"\n",
"In [18]: evec\n",
"Out[18]: array([[-0.77167559, -0.52633654, 0.57513303],\n",
" [ 0.50360249, 0.76565448, -0.80920669],\n",
" [-0.38846018, 0.36978786, 0.12002724]])\n",
"```\n",
"\n",
"The first eigenvalue and its corresponding eigenvector are given by\n",
"\n",
"``` ipython\n",
"In [19]: lam[0]\n",
"Out[19]: (2.4099535647625494+0j)\n",
"\n",
"In [20]: evec[:,0]\n",
"Out[20]: array([-0.77167559, 0.50360249, -0.38846018])\n",
"```\n",
"\n",
"We can check that they satisfy the\n",
"$\\mathsf{A}\\mathbf{x} = \\lambda \\mathbf{x}:$\n",
"\n",
"``` ipython\n",
"In [21]: dot(A,evec[:,0])\n",
"Out[21]: array([-1.85970234, 1.21365861, -0.93617101])\n",
"\n",
"In [22]: lam[0]*evec[:,0]\n",
"Out[22]: array([-1.85970234+0.j, 1.21365861+0.j, \n",
" -0.93617101+0.j])\n",
"```\n",
"\n",
"Thus we see by direct substitution that the left and right sides of\n",
"$\\mathsf{A}\\mathbf{x} = \\lambda \\mathbf{x}:$ are equal. In general, the\n",
"eigenvalues can be complex, so their values are reported as complex\n",
"numbers.\n",
"\n",
"single: linear algebra; generalized eigenvalue problem\n",
"\n",
"#### Generalized eigenvalue problem\n",
"\n",
"The `scipy.linalg.eig` function can also solve the *generalized*\n",
"eigenvalue problem\n",
"\n",
"$$\\mathsf{A}\\mathbf{x} = \\lambda \\mathsf{B}\\mathbf{x}$$\n",
"\n",
"where $\\mathsf{B}$ is a square matrix with the same size as\n",
"$\\mathsf{A}$. Suppose, for example, that we have\n",
"\n",
"``` ipython\n",
"In [22]: A = array([[2, 4, 6], [1, -3, -9], [8, 5, -7]])\n",
"Out[22]: B = array([[5, 9, 1], [-3, 1, 6], [4, 2, 8]])\n",
"```\n",
"\n",
"Then we can solve the generalized eigenvalue problem by entering\n",
"$\\mathsf{B}$ as the optional second argument to `scipy.linalg.eig`\n",
"\n",
"``` ipython\n",
"In [23]: lam, evec = scipy.linalg.eig(A,B)\n",
"```\n",
"\n",
"The solutions are returned in the same fashion as before, as an array\n",
"`lam` whose entries are the eigenvalues and a matrix `evac` whose rows\n",
"are the eigenvectors.\n",
"\n",
"``` ipython\n",
"In [24]: lam\n",
"Out[24]: array([-1.36087907+0.j, 0.83252442+0.j,\n",
" -0.10099858+0.j])\n",
"\n",
"In [25]: evec\n",
"Out[25]: array([[-0.0419907 , -1. , 0.93037493],\n",
" [-0.43028153, 0.17751302, -1. ],\n",
" [ 1. , -0.29852465, 0.4226201 ]])\n",
"```\n",
"\n",
"single: linear algebra; Hermitian and banded matrices\n",
"\n",
"#### Hermitian and banded matrices\n",
"\n",
"SciPy has a specialized routine for solving eigenvalue problems for\n",
"Hermitian (or real symmetric) matrices. The routine for hermitian\n",
"matrices is `scipy.linalg.eigh`. It is more efficient (faster and uses\n",
"less memory) than `scipy.linalg.eig`. The basic syntax of the two\n",
"routines is the same, although some of the *optional* arguments are\n",
"different. Both routines can solve generalized as well as standard\n",
"eigenvalue problems.\n",
"\n",
"SciPy also has a specialized routine `scipy.linalg.eig_banded` for\n",
"solving eigenvalue problems for real symmetric or complex hermitian\n",
"banded matrices.\n",
"\n",
"single: non-linear equations single: non-linear equations; solving\n",
"single: SciPy; non-linear equations see: solving non-linear equations;\n",
"non-linear equations see: roots of equations; non-linear equations\n",
"\n",
"Solving non-linear equations\n",
"----------------------------\n",
"\n",
"SciPy has many different routines for numerically solving non-linear\n",
"equations or systems of non-linear equations. Here we will introduce\n",
"only a few of these routines, the ones that are relatively simple and\n",
"appropriate for the most common types of nonlinear equations.\n",
"\n",
"### Single equations of a single variable\n",
"\n",
"Solving a single nonlinear equation is enormously simpler than solving a\n",
"system of nonlinear equations, so that is where we start. A word of\n",
"caution: solving non-linear equations can be a tricky business so it is\n",
"important that you have a good sense of the behavior of the function you\n",
"are trying to solve. The best way to do this is to plot the function\n",
"over the domain of interest before trying to find the solutions. This\n",
"will greatly assist you in finding the solutions you seek and avoiding\n",
"spurious solutions.\n",
"\n",
"We begin with a concrete example. Suppose we want to find the solutions\n",
"to the equation\n",
"\n",
"$$\\tan x=\\sqrt{(8/x)^2-1}$$\n",
"\n",
"Plots of $\\tan x$ and $\\sqrt{(8/x)^2-1}$ *vs* $x$ are shown in the top\n",
"plot in the figure `fig-subplotDemo`, albeit with $x$ replaced by\n",
"$\\theta$. The solutions to this equation are those $x$ values where the\n",
"two curves $\\tan x$ and $\\sqrt{(8/x)^2-1}$ cross each other. The first\n",
"step towards obtaining a numerical solution is to rewrite the equation\n",
"to be solved in the form $f(x)=0$. Doing so, the above equation becomes\n",
"\n",
"$$\\tan x - \\sqrt{(8/x)^2-1} = 0$$\n",
"\n",
"Obviously the two equations above have the same solutions for $x$.\n",
"Parenthetically we mention that the problem of finding the solutions to\n",
"equations of the form $f(x)=0$ is often referred to as *finding the\n",
"roots* of $f(x)$.\n",
"\n",
"Next, we plot $f(x)$ over the domain of interest, in this case from\n",
"$x=0$ to 8. We are only interested in positive solutions and for $x>8$,\n",
"the equation has no real solutions as the argument of the square root\n",
"becomes negative. The solutions, the points where $f(x)=0$ are indicated\n",
"by green circles; there are three of them. Another notable feature of\n",
"the function is that it diverges to $\\pm\\infty$ at\n",
"$x = \\{0, \\pi/2, 3\\pi/2, 5\\pi/2\\}$.\n",
"\n",
"\n",
"Roots of a nonlinear function\n",
"\n",
"\n",
"single: non-linear equations; Brent method\n",
"\n",
"#### Brent method\n",
"\n",
"One of the workhorses for finding solutions to a single variable\n",
"nonlinear equation is the method of Brent, discussed in many texts on\n",
"numerical methods. SciPy's implementation of the Brent algorithm is the\n",
"function `scipy.optimize.brentq(f, a, b)`, which has three required\n",
"arguments. The first `f` is the name of the user-defined function to be\n",
"solved. The next two, `a` and `b` are the $x$ values that bracket the\n",
"solution you are looking for. You should choose `a` and `b` so that\n",
"there is only one solutions in the interval between `a` and `b`. Brent's\n",
"method also requires that `f(a)` and `f(b)` have opposite signs; an\n",
"error message is returned if they do not. Thus to find the three\n",
"solutions to $\\tan x - \\sqrt{(8/x)^2-1} = 0$, we need to run\n",
"`scipy.optimize.brentq(f, a, b)` three times using three different\n",
"values of `a` and `b` that bracket each of the three solutions. The\n",
"program below illustrates the how to use `scipy.optimize.brentq`\n",
"\n",
"``` python\n",
"import numpy as np\n",
"import scipy.optimize\n",
"import matplotlib.pyplot as plt\n",
"\n",
"def tdl(x):\n",
" y = 8./x\n",
" return np.tan(x) - np.sqrt(y*y-1.0)\n",
"\n",
"# Find true roots\n",
"\n",
"rx1 = scipy.optimize.brentq(tdl, 0.5, 0.49*np.pi)\n",
"rx2 = scipy.optimize.brentq(tdl, 0.51*np.pi, 1.49*np.pi)\n",
"rx3 = scipy.optimize.brentq(tdl, 1.51*np.pi, 2.49*np.pi)\n",
"rx = np.array([rx1, rx2, rx3])\n",
"ry = np.zeros(3)\n",
"# print using a list comprehension\n",
"print('\\nTrue roots:')\n",
"print('\\n'.join('f({0:0.5f}) = {1:0.2e}'.format(x, tdl(x)) for x in rx)) \n",
"\n",
"# Find false roots\n",
"\n",
"rx1f = scipy.optimize.brentq(tdl, 0.49*np.pi, 0.51*np.pi)\n",
"rx2f = scipy.optimize.brentq(tdl, 1.49*np.pi, 1.51*np.pi)\n",
"rx3f = scipy.optimize.brentq(tdl, 2.49*np.pi, 2.51*np.pi)\n",
"rxf = np.array([rx1f, rx2f, rx3f])\n",
"# print using a list comprehension\n",
"print('\\nFalse roots:')\n",
"print('\\n'.join('f({0:0.5f}) = {1:0.2e}'.format(x, tdl(x)) for x in rxf)) \n",
"\n",
"# Plot function and various roots\n",
"\n",
"x = np.linspace(0.7, 8, 128)\n",
"y = tdl(x)\n",
"# Create masked array for plotting\n",
"ymask = np.ma.masked_where(np.abs(y)>20., y)\n",
"\n",
"plt.figure(figsize=(6, 4))\n",
"plt.plot(x, ymask)\n",
"plt.axhline(color='black')\n",
"plt.axvline(x=np.pi/2., color=\"gray\", linestyle='--', zorder=-1)\n",
"plt.axvline(x=3.*np.pi/2., color=\"gray\", linestyle='--', zorder=-1)\n",
"plt.axvline(x=5.*np.pi/2., color=\"gray\", linestyle='--', zorder=-1)\n",
"plt.xlabel(r'$x$')\n",
"plt.ylabel(r'$\\tan x - \\sqrt{(8/x)^2-1}$')\n",
"plt.ylim(-8, 8)\n",
"\n",
"plt.plot(rx, ry, 'og', ms=5, label='true roots')\n",
"\n",
"plt.plot(rxf, ry, 'xr', ms=5, label='false roots')\n",
"plt.legend(numpoints=1, fontsize='small', loc = 'upper right',\n",
" bbox_to_anchor = (0.92, 0.97))\n",
"plt.tight_layout()\n",
"plt.show()\n",
"```\n",
"\n",
"Running this code generates the following output:\n",
"\n",
"``` ipython\n",
"In [1]: run rootbrentq.py\n",
"\n",
"True roots:\n",
"f(1.39547) = -6.39e-14\n",
"f(4.16483) = -7.95e-14\n",
"f(6.83067) = -1.11e-15\n",
"\n",
"False roots:\n",
"f(1.57080) = -1.61e+12\n",
"f(4.71239) = -1.56e+12\n",
"f(7.85398) = 1.16e+12\n",
"```\n",
"\n",
"The Brent method finds the three true roots of the equation quickly and\n",
"accurately when you provide values for the brackets `a` and `b` that are\n",
"valid. However, like many numerical methods for finding roots, the Brent\n",
"method can produce spurious roots as it does in the above example when\n",
"`a` and `b` bracket singularities like those at\n",
"$x = \\{\\pi/2, 3\\pi/2, 5\\pi/2\\}$. Here we evaluated the function at the\n",
"purported roots found by `brentq` to verify that the values of $x$ found\n",
"were indeed roots. For the true roots, the values of the function were\n",
"very near zero, to within an acceptable roundoff error of less than\n",
"$10^{-13}$. For the false roots, exceedingly large numbers on the order\n",
"of $10^{12}$ were obtained, indicating a possible problem with these\n",
"roots. These results, together with the plots, allow you to\n",
"unambiguously identify the true solutions to this nonlinear function.\n",
"\n",
"The `brentq` function has a number of optional keyword arguments that\n",
"you may find useful. One keyword argument causes `brentq` to return not\n",
"only the solution but the value of the function evaluated at the\n",
"solution. Other arguments allow you to specify a tolerance to which the\n",
"solution is found as well as a few other parameters possibly of\n",
"interest. Most of the time, you can leave the keyword arguments at their\n",
"default values. See the `brentq` entry online on the SciPy web site for\n",
"more information.\n",
"\n",
"single: non-linear equations; Newton-Raphson method\n",
"\n",
"#### Other methods for solving equations of a single variable\n",
"\n",
"SciPy provides a number of other methods for solving nonlinear equations\n",
"of a single variable. It has an implementation of the Newton-Raphson\n",
"method called `scipy.optimize.newton`. It's the racecar of such methods;\n",
"its super fast but less stable that the Brent method. To fully realize\n",
"its speed, you need to specify not only the function to be solved, but\n",
"also its first derivative, which is often more trouble than its worth.\n",
"You can also specify its second derivative, which may further speed up\n",
"finding the solution. If you do not specify the first or second\n",
"derivatives, the method uses the secant method, which is usually slower\n",
"than the Brent method.\n",
"\n",
"single: non-linear equations; Ridder method single: non-linear\n",
"equations; Bisection method\n",
"\n",
"Other methods, including the Ridder (`scipy.optimize.ridder`) and\n",
"bisection (`scipy.optimize.bisect`), are also available, although the\n",
"Brent method is generally superior. SciPy let's you use your favorite.\n",
"\n",
"single: non-linear equations; systems of nonlinear equations\n",
"\n",
"### Solving systems of nonlinear equations\n",
"\n",
"Solving systems of nonlinear equations is not for the faint of heart. It\n",
"is a difficult problem that lacks any general purpose solutions.\n",
"Nevertheless, SciPy provides quite an assortment of numerical solvers\n",
"for nonlinear systems of equations. However, because of the complexity\n",
"and subtleties of this class of problems, we do not discuss their use\n",
"here.\n",
"\n",
"Exercises\n",
"---------\n",
"\n",
"1. Use NumPy's `polyval` function together with SciPy to plot the\n",
" following functions:\n",
" 1. The first four Chebyshev polynomials of first kind. Plot these\n",
" over the interval from -1 to +1.\n",
" 2. The first four Hermite polynomials *multiplied* by $e^{-x^2/2}$.\n",
" Plot these on the interval from -5 to +5. These are the first\n",
" four wave functions of the quantum mechanical simple harmonic\n",
" oscillator."
]
}
],
"metadata": {
"jupytext": {
"cell_metadata_filter": "all",
"formats": "ipynb,py:percent",
"notebook_metadata_filter": "all,-language_info,-toc,-latex_envs"
},
"kernelspec": {
"display_name": "Python 3",
"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.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
![\"\"](\"attachment:specFuncPlots.jpg\")
![\"\"](\"attachment:odePend.png\")
![\"\"](\"attachment:fftSimple.png\")
![\"\"](\"attachment:rootsBrent.png\")