Differential Equations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vector Fields and Euler's Method\n",
    "\n",
    "by Mike Hansen (tested and updated by William Stein, and later by Dan Drake)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "80d7ba3bee494584a697ec194d78476d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "SW50ZXJhY3RpdmUgZnVuY3Rpb24gPGZ1bmN0aW9uIF8gYXQgMHg3ZjRhZmYzMjVlZDg+IHdpdGggMTAgd2lkZ2V0cwogIGY6IEV2YWxUZXh0KHZhbHVlPXUneScsIGRlc2NyaXB0aW9uPXUnZifigKY=\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x,y = var('x,y')\n",
    "from sage.ext.fast_eval import fast_float\n",
    "@interact\n",
    "def _(f = input_box(default=y), g=input_box(default=-x*y+x^3-x),\n",
    "      xmin=input_box(default=-1), xmax=input_box(default=1),\n",
    "      ymin=input_box(default=-1), ymax=input_box(default=1),\n",
    "      start_x=input_box(default=0.5), start_y=input_box(default=0.5),\n",
    "      step_size=(0.01,(0.001, 0.2)), steps=(1000, (0, 2000)) ):\n",
    "    ff = fast_float(f, 'x', 'y')\n",
    "    gg = fast_float(g, 'x', 'y')\n",
    "    steps = int(steps)\n",
    "\n",
    "    points = [ (start_x, start_y) ]\n",
    "    for i in range(steps):\n",
    "        xx, yy = points[-1]\n",
    "        try:\n",
    "            points.append( (xx + step_size * ff(xx,yy), yy + step_size * gg(xx,yy)) )\n",
    "        except (ValueError, ArithmeticError, TypeError):\n",
    "            break\n",
    "\n",
    "    starting_point = point(points[0], pointsize=50)\n",
    "    solution = line(points)\n",
    "    vector_field = plot_vector_field( (f,g), (x,xmin,xmax), (y,ymin,ymax) )\n",
    "\n",
    "    result = vector_field + starting_point + solution\n",
    "\n",
    "    pretty_print(html(r\"$\\displaystyle\\frac{dx}{dt} = %s$  $ \\displaystyle\\frac{dy}{dt} = %s$\" % (latex(f),latex(g))))\n",
    "    result.show(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vector Field with Runga-Kutta-Fehlberg\n",
    "\n",
    "by Harald Schilly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8f9f8a80b9c54fb89f1469801d909124",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "SW50ZXJhY3RpdmUgZnVuY3Rpb24gPGZ1bmN0aW9uIF8gYXQgMHg3ZjRhZmYzMjU1ZjA+IHdpdGggMTIgd2lkZ2V0cwogIGZpbjogRXZhbFRleHQodmFsdWU9dScteCp5XjMgLSAxLzEyKig0KnnigKY=\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "var('x y')\n",
    "@interact\n",
    "def _(fin = input_box(default=y+exp(x/10)-1/3*((x-1/2)^2+y^3)*x-x*y^3), gin=input_box(default=x^3-x+1/100*exp(y*x^2+x*y^2)-0.7*x),\n",
    "      xmin=input_box(default=-1), xmax=input_box(default=1.8),\n",
    "      ymin=input_box(default=-1.3), ymax=input_box(default=1.5),\n",
    "      x_start=(-1,(-2.0,2.0)), y_start=(0,(-2.0,2.0)), error=(0.5,(0.0,1.0)),\n",
    "      t_length=(23,(0, 100)) , num_of_points = (1500,(5,2000)),\n",
    "      algorithm = selector([\n",
    "         (\"rkf45\" , \"runga-kutta-felhberg (4,5)\"),\n",
    "         (\"rk2\" , \"embedded runga-kutta (2,3)\"),\n",
    "         (\"rk4\" , \"4th order classical runga-kutta\"),\n",
    "         (\"rk8pd\" , 'runga-kutta prince-dormand (8,9)'),\n",
    "         (\"rk2imp\" , \"implicit 2nd order runga-kutta at gaussian points\"),\n",
    "         (\"rk4imp\" , \"implicit 4th order runga-kutta at gaussian points\"),\n",
    "         (\"bsimp\" , \"implicit burlisch-stoer (requires jacobian)\"),\n",
    "         (\"gear1\" , \"M=1 implicit gear\"),\n",
    "         (\"gear2\" , \"M=2 implicit gear\")\n",
    "      ])):\n",
    "    f(x,y)=fin\n",
    "    g(x,y)=gin\n",
    "\n",
    "    ff = f._fast_float_(*f.args())\n",
    "    gg = g._fast_float_(*g.args())\n",
    "\n",
    "    #solve\n",
    "    path = []\n",
    "    err = error\n",
    "    xerr = 0\n",
    "    for yerr in [-err, 0, +err]:\n",
    "      T=ode_solver()\n",
    "      T.algorithm=algorithm\n",
    "      T.function = lambda t, yp: [ff(yp[0],yp[1]), gg(yp[0],yp[1])]\n",
    "      T.jacobian = lambda t, yp: [[diff(fun,dval)(yp[0],yp[1]) for dval in [x,y]] for fun in [f,g]]\n",
    "      T.ode_solve(y_0=[x_start + xerr, y_start + yerr],t_span=[0,t_length],num_points=num_of_points)\n",
    "      path.append(line([p[1] for p in T.solution]))\n",
    "\n",
    "    #plot\n",
    "    vector_field = plot_vector_field( (f,g), (x,xmin,xmax), (y,ymin,ymax) )\n",
    "    starting_point = point([x_start, y_start], pointsize=50)\n",
    "    show(vector_field + starting_point + sum(path), aspect_ratio=1, figsize=[8,9])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Euler-Maruyama method and geometric Brownian motion (a common simple model of the stock market)\n",
    "\n",
    "by Marshall Hampton"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f7fcc33e99784e4b9cafbfa1896bd892",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "TWFudWFsIGludGVyYWN0aXZlIGZ1bmN0aW9uIDxmdW5jdGlvbiBFdWxlck1hcnV5YW1hRXhhbXBsZSBhdCAweDdmNGFmZDhlMzQxMD4gd2l0aCA1IHdpZGdldHMKICBtdTogU2VsZWN0aW9uU2zigKY=\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def EulerMaruyama(xstart, ystart, xfinish, nsteps, f1, f2): \n",
    "    sol = [ystart] \n",
    "    xvals = [xstart] \n",
    "    h = N((xfinish-xstart)/nsteps) \n",
    "    for step in range(nsteps): \n",
    "        sol.append(sol[-1] + h*f1(sol[-1]) + h^(.5)*f2(sol[-1])*normalvariate(0,1)) \n",
    "        xvals.append(xvals[-1] + h) \n",
    "    return zip(xvals,sol)\n",
    "    \n",
    "out = Graphics()\n",
    "save(out,'temp')\n",
    "@interact\n",
    "def EulerMaruyamaExample(mu = slider(srange(0,10,.1),default=2.0),\n",
    "                        sigma = slider(srange(0,10,.1),default=0.5),\n",
    "                        plots_at_a_time = slider(range(1,100),default=10), \n",
    "                        number_of_steps = slider(range(1,1000),default=100), \n",
    "                        clear_plot = checkbox(True), \n",
    "                        auto_update=False):\n",
    "    html('<center>Example of the Euler-Maruyama method applied to<br>the stochastic differential equation for geometric Brownian motion</center>')\n",
    "    html('<center>$S = S_0 + \\int_0^t \\mu S dt + \\int_0^t \\sigma S dW$</center>')\n",
    "    emplot = list_plot(EulerMaruyama(0,1,1,number_of_steps,lambda x: mu*x,lambda x:sigma*x),plotjoined=True)\n",
    "    for i in range(1,plots_at_a_time):\n",
    "        emplot = emplot + list_plot(EulerMaruyama(0,1,1,100,lambda x: mu*x,lambda x:sigma*x),plotjoined=True)\n",
    "    if clear_plot:\n",
    "        out = emplot\n",
    "        save(out,'temp')\n",
    "    else:\n",
    "        out = load('temp')\n",
    "        out = out + emplot\n",
    "        save(out,'temp')\n",
    "    show(out, figsize = [8,5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Autonomous equations and stable/unstable fixed points\n",
    "\n",
    "by Marshall Hampton"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "def RK4_1d(f, double t_start, double y_start, double t_end, int steps, double y_upper = 10**6, double y_lower = -10**6):\n",
    "    '''\n",
    "    Fourth-order scalar Runge-Kutta solver with fixed time steps. f must be a function of t,y, \n",
    "    where y is just a scalar variable.\n",
    "    '''\n",
    "    cdef double step_size = (t_end - t_start)/steps\n",
    "    cdef double t_current = t_start\n",
    "    cdef double y_current = y_start\n",
    "    cdef list answer_table = []\n",
    "    cdef int j\n",
    "    answer_table.append([t_current,y_current])\n",
    "    for j in range(0,steps):\n",
    "        k1=f(t_current, y_current)\n",
    "        k2=f(t_current+step_size/2, y_current + k1*step_size/2)\n",
    "        k3=f(t_current+step_size/2, y_current + k2*step_size/2)\n",
    "        k4=f(t_current+step_size, y_current + k3*step_size)\n",
    "        t_current += step_size\n",
    "        y_current = y_current + (step_size/6)*(k1+2*k2+2*k3+k4)\n",
    "        if y_current > y_upper or y_current < y_lower: \n",
    "            j = steps\n",
    "        answer_table.append([t_current,y_current])\n",
    "    return answer_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cd0435a2c6424b3fb67d49466b236c95",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "SW50ZXJhY3RpdmUgZnVuY3Rpb24gPGZ1bmN0aW9uIGF1dG9ub21vdXNfcGxvdCBhdCAweDdmNGFmZjdiMjIzMD4gd2l0aCAyIHdpZGdldHMKICBwb2x5OiBFdmFsVGV4dCh2YWx1ZT11Jyh4IC3igKY=\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sage.rings.polynomial.real_roots import real_roots\n",
    "var('x')\n",
    "@interact\n",
    "def autonomous_plot(poly=input_box(x*(x-1)*(x-2),label='polynomial'), t_end = slider(1,10,step_size = .1)): \n",
    "    var('x')   \n",
    "    y = polygen(ZZ)\n",
    "    ypoly = sage_eval(repr(poly).replace('x','y'),locals=locals())\n",
    "    rr = real_roots(ypoly, max_diameter = 1/100)\n",
    "    eps = 0.2\n",
    "    delvals = .04\n",
    "    minr = min([x[0][0] for x in rr])\n",
    "    maxr = max([x[0][1] for x in rr])\n",
    "    svals = [(minr-eps)*t+(1-t)*(maxr+eps) for t in srange(0,1+delvals,delvals)]\n",
    "    def polyf(t,xi):\n",
    "        return poly(x=xi)\n",
    "    paths = [RK4_1d(polyf,0.0,q,t_end,100.0) for q in svals]    \n",
    "    miny = max(minr-eps,min([min([q1[1] for q1 in q]) for q in paths]))\n",
    "    maxy = min(maxr+eps,max([max([q1[1] for q1 in q]) for q in paths]))\n",
    "    solpaths = sum([line(q) for q in paths])\n",
    "    fixedpoints = sum([line([[0,(q[0][0]+q[0][1])/2],[t_end,(q[0][0]+q[0][1])/2]], rgbcolor = (1,0,0)) for q in rr])\n",
    "    var('t')\n",
    "    html(\"Autonomous differential equation $x' = p(x)$\")\n",
    "    show(solpaths+fixedpoints, ymin = miny, ymax = maxy, xmin = 0, xmax = t_end, figsize = [6,4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heat equation using Fourier series\n",
    "\n",
    "by Pablo Angulo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d241f0a0b00d40afb6c8476efe7155cf",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "SW50ZXJhY3RpdmUgZnVuY3Rpb24gPGZ1bmN0aW9uIF8gYXQgMHg3ZjRhZmY3YjIzOTg+IHdpdGggMSB3aWRnZXQKICB0OiBTZWxlY3Rpb25TbGlkZXIoZGVzY3JpcHRpb249dSd0Jywgb3B0aW/igKY=\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "var('x')\n",
    "x0  = 0\n",
    "k=1\n",
    "f   = x*exp(-x^2)\n",
    "p   = plot(f,0,2*pi, thickness=2)\n",
    "c   = 1/pi\n",
    "orden=10\n",
    "alpha=[(n,c*numerical_integral(f*sin(x*n/2),0,2*pi)[0] ) for n in range(1,orden)]\n",
    "\n",
    "@interact\n",
    "def _(t = (0.1*j for j in (0..10)) ):\n",
    "    ft = sum( a*sin(x*n/2)*exp(-k*(n/2)^2*t) for n,a in alpha)\n",
    "    pt = plot(ft, 0, 2*pi, color='green', thickness=2)\n",
    "    show( p + pt, ymin = -.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 8.5.beta1",
   "language": "",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}