Test IVP

qualitative_bench.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import solve_ivp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. VDPOL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu = 1e3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(t, y):\n",
    "    return [\n",
    "        y[1],\n",
    "        mu * (1 - y[0] ** 2) * y[1] - y[0]\n",
    "    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def jac(t, y):\n",
    "    return [\n",
    "        [0, 1],\n",
    "        [-2 * mu * y[0] * y[1] - 1, mu * (1 - y[0] ** 2)]\n",
    "    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_span = [0, 2 * mu]\n",
    "y0 = [2, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res = solve_ivp(fun, t_span, y0, method='Radau')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(194,)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.t.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.status"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. ROBER"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(t, y):\n",
    "    return np.array([\n",
    "        -0.04 * y[0] + 1e4 * y[1] * y[2],\n",
    "        0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * y[1]**2,\n",
    "        3e7 * y[1]**2\n",
    "    ])\n",
    "\n",
    "y0 = np.array([1, 0, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def jac(t, y):\n",
    "    return np.array([\n",
    "        [-0.04, 1e4 * y[2], 1e4 * y[1]],\n",
    "        [0.04, -1e4 * y[2] - 6e7 * y[1], -1e4 * y[1]],         \n",
    "        [0, 6e7 * y[1], 0]\n",
    "    ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_span = [0, 1e11]\n",
    "y0 = [1, 0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res = solve_ivp(fun, t_span, y0, jac=jac, method='Radau')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(79,)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.t.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.status"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. OREGO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fun(t, y):\n",
    "    return [\n",
    "        77.27 * (y[1] + y[0] * (1 - 8.375e-6 * y[0] - y[1])),\n",
    "        (y[2] - (1 + y[0]) * y[1]) / 77.27,\n",
    "        0.161 * (y[0] - y[2])\n",
    "    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def jac(t, y):\n",
    "    J = np.array([\n",
    "        [1 - 2 * 8.375e-6 * y[0] - y[1], 1 - y[0], 0],\n",
    "        [-y[1], -(1 + y[0]), 1],\n",
    "        [1, 0, -1]\n",
    "    ])\n",
    "    J[0] *= 77.27\n",
    "    J[1] /= 77.27\n",
    "    J[2] *= 0.161\n",
    "    return J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_span = [0, 360]\n",
    "y0 = [1, 2, 3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res = solve_ivp(fun, t_span, y0, method='Radau')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(237,)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.t.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.status"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 4. HIRES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(t, y):\n",
    "    return np.array([\n",
    "        -1.71 * y[0] + 0.43 * y[1] + 8.32 * y[2] + 0.0007,\n",
    "        1.71 * y[0] - 8.75 * y[1],\n",
    "        -10.03 * y[2] + 0.43 * y[3] + 0.035 * y[4],\n",
    "        8.32 * y[1] + 1.71 * y[2] - 1.12 * y[3],\n",
    "        -1.745 * y[4] + 0.43 * y[5] + 0.43 * y[6],\n",
    "        -280 * y[5] * y[7] + 0.69 * y[3] + 1.71 * y[4] - 0.43 * y[5] + 0.69 * y[6],\n",
    "        280 * y[5] * y[7] - 1.81 * y[6],\n",
    "        -280 * y[5] * y[7] + 1.81 * y[6] \n",
    "    ])\n",
    "\n",
    "def jac(t, y):\n",
    "    return [\n",
    "        [-1.71, 0.43, 8.32, 0, 0, 0, 0, 0],\n",
    "        [1.71, -8.75, 0, 0, 0, 0, 0, 0],\n",
    "        [0, 0, -10.03, 0.43, 0.035, 0, 0, 0],\n",
    "        [0, 8.32, 1.71, -1.12, 0, 0, 0, 0],\n",
    "        [0, 0, 0, 0, -1.745, 0.43, 0.43, 0],\n",
    "        [0, 0, 0, 0.69, 1.71, -0.43 - 280 * y[7], 0.69, -280 * y[5]],\n",
    "        [0, 0, 0, 0, 0, 280 * y[7], -1.81, 280 * y[5]],\n",
    "        [0, 0, 0, 0, 0, -280 * y[7], 1.81, -280 * y[5]]\n",
    "    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_span = [0, 421.8122]\n",
    "y0 = np.array([1, 0, 0, 0, 0, 0, 0, 0.0057])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res = solve_ivp(fun, t_span, y0, method='Radau')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(52,)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.t.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.status"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. E5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = 7.89e-10\n",
    "B = 1.1e7\n",
    "C = 1.13e3\n",
    "M = 1e6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(t, y):  \n",
    "    return [\n",
    "        -A * y[0] - B * y[0] * y[2],\n",
    "        A * y[0] - M * C * y[1] * y[2],\n",
    "        A * y[0] - B * y[0] * y[2] - M * C * y[1] * y[2] + C * y[3],\n",
    "        B * y[0] * y[2] - C * y[3]\n",
    "    ]\n",
    "\n",
    "def jac(t, y):\n",
    "    return np.array([\n",
    "        [-A - B * y[2], 0, -B * y[0], 0],\n",
    "        [A, -M * C * y[2], -M * C * y[1], 0],\n",
    "        [A - B * y[2], -M * C * y[2], -B * y[0] - M * C * y[1], C],\n",
    "        [B * y[2], 0, B * y[0], -C]\n",
    "    ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_span = [0, 1e13]\n",
    "y0 = [0.00176, 0, 0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res = solve_ivp(fun, t_span, y0, method='Radau', jac=jac, atol=[1e-3, 1.7e-24, 1.7e-24, 1.7e-24])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(151,)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.t.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.status"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. CUSP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 32\n",
    "D = N ** 2 / 144"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(t, u):\n",
    "    y, a, b = u[:N], u[N: 2 * N], u[2 * N:]\n",
    "    y_ext = np.hstack((y[-1], y, y[0]))\n",
    "    a_ext = np.hstack((a[-1], a, a[0]))\n",
    "    b_ext = np.hstack((b[-1], b, b[0]))\n",
    "    u = (y - 0.7) * (y - 1.3)\n",
    "    v = u / (u + 0.1)    \n",
    "    y_dot = -1e4 * (y**3 + a * y + b) + D * (y_ext[:-2] - 2 * y + y_ext[2:])\n",
    "    a_dot = b + 0.07 * v + D * (a_ext[:-2] - 2 * a + a_ext[2:])\n",
    "    b_dot = (1 - a**2) * b - a - 0.4 * y + 0.035 * v + D * (b_ext[:-2] - 2 * b + b_ext[2:])\n",
    "    \n",
    "    return np.hstack((y_dot, a_dot, b_dot))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_span = [0, 1.1]\n",
    "i = np.arange(N) + 1\n",
    "y0 = np.zeros(N)\n",
    "a0 = -2 * np.cos(2 * np.pi * i / N)\n",
    "b0 = 2 * np.sin(2 * np.pi * i / N)\n",
    "y0 = np.hstack((y0, a0, b0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res = solve_ivp(fun, t_span, y0, method='Radau')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(76,)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.t.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.status"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. RINGMOD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "C = 1.6e-8\n",
    "Cs = 2e-12\n",
    "Cp = 1e-8\n",
    "Lh = 4.45\n",
    "Ls1 = 0.002\n",
    "Ls2 = 5e-4\n",
    "Ls3 = 5e-4\n",
    "gamma = 40.67286402e-9\n",
    "R = 25e3\n",
    "Rp = 50\n",
    "Rg1 = 36.3\n",
    "Rg2 = 17.3\n",
    "Rg3 = 17.3\n",
    "Ri = 50\n",
    "Rc = 600\n",
    "delta = 17.7493332"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def q(U):\n",
    "    with np.errstate(over='ignore'):\n",
    "        return gamma * np.expm1(delta * U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dq(U):\n",
    "    with np.errstate(over='ignore'):\n",
    "        return delta * gamma * np.exp(delta * U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(t, y):\n",
    "    Uin1 = 0.5 * np.sin(2000 * np.pi * t)\n",
    "    Uin2 = 2 * np.sin(20000 * np.pi * t)\n",
    "    \n",
    "    Ud1 = y[2] - y[4] - y[6] - Uin2\n",
    "    Ud2 = -y[3] + y[5] - y[6] - Uin2\n",
    "    Ud3 = y[3] + y[4] + y[6]  + Uin2\n",
    "    Ud4 = -y[2] - y[5] + y[6] + Uin2 \n",
    "    q1 = q(Ud1)\n",
    "    q2 = q(Ud2)\n",
    "    q3 = q(Ud3)\n",
    "    q4 = q(Ud4)\n",
    "    \n",
    "    return [\n",
    "        (y[7] - 0.5 * y[9] + 0.5 * y[10] + y[13] - y[0] / R) / C,\n",
    "        (y[8] - 0.5 * y[11] + 0.5 * y[12] + y[14] - y[1] / R) / C,\n",
    "        (y[9] - q1 + q4) / Cs,\n",
    "        (-y[10] + q2 - q3) / Cs,\n",
    "        (y[11] + q1 - q3) / Cs,\n",
    "        (-y[12] - q2 + q4) / Cs,\n",
    "        (-y[6] / Rp + q1 + q2 - q3 - q4) / Cp,\n",
    "        -y[0] / Lh,\n",
    "        -y[1] / Lh,\n",
    "        (0.5 * y[0] - y[2] - Rg2 * y[9]) / Ls2,\n",
    "        (-0.5 * y[0] + y[3] - Rg3 * y[10]) / Ls3,\n",
    "        (0.5 * y[1] - y[4] - Rg2 * y[11]) / Ls2,\n",
    "        (-0.5 * y[1] + y[5] - Rg3 * y[12]) / Ls3,\n",
    "        (-y[0] + Uin1 - (Ri + Rg1) * y[13]) / Ls1,\n",
    "        (-y[1] - (Rc + Rg1) * y[14]) / Ls1\n",
    "    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def jac(t, y):\n",
    "    Uin1 = 0.5 * np.sin(2000 * np.pi * t)\n",
    "    Uin2 = 2 * np.sin(20000 * np.pi * t)\n",
    "    \n",
    "    Ud1 = y[2] - y[4] - y[6] - Uin2\n",
    "    Ud2 = -y[3] + y[5] - y[6] - Uin2\n",
    "    Ud3 = y[3] + y[4] + y[6]  + Uin2\n",
    "    Ud4 = -y[2] - y[5] + y[6] + Uin2 \n",
    "    q1 = q(Ud1)\n",
    "    q2 = q(Ud2)\n",
    "    q3 = q(Ud3)\n",
    "    q4 = q(Ud4)\n",
    "    \n",
    "    dq1 = dq(Ud1)\n",
    "    dq2 = dq(Ud2)\n",
    "    dq3 = dq(Ud3)\n",
    "    dq4 = dq(Ud4)\n",
    "    \n",
    "    J = np.zeros((15, 15))\n",
    "    (y[7] - 0.5 * y[9] + 0.5 * y[10] + y[13] - y[0] / R) / C\n",
    "    J[0, 0] = -1 / R\n",
    "    J[0, 7] = 1\n",
    "    J[0, 9] = -0.5\n",
    "    J[0, 10] = 0.5\n",
    "    J[0, 13] = 1\n",
    "    J[0] /= C\n",
    "    \n",
    "    J[1, 1] = -1 / R\n",
    "    J[1, 8] = 1\n",
    "    J[1, 11] = -0.5\n",
    "    J[1, 12] = 0.5\n",
    "    J[1, 14] = 1\n",
    "    J[1] /= C\n",
    "    \n",
    "    J[2, 9] = 1\n",
    "    J[2, 2] = -dq1 - dq4\n",
    "    J[2, 4] = dq1\n",
    "    J[2, 5] = -dq4\n",
    "    J[2, 6] = dq1 + dq4\n",
    "    J[2] /= Cs\n",
    "    \n",
    "    J[3, 3] = -dq2 - dq3\n",
    "    J[3, 4] = -dq3\n",
    "    J[3, 5] = dq2\n",
    "    J[3, 6] = -dq2 - dq3\n",
    "    J[3, 10] = -1\n",
    "    J[3] /= Cs\n",
    "    \n",
    "    J[4, 2] = dq1\n",
    "    J[4, 3] = -dq3\n",
    "    J[4, 4] = -dq1 - dq3\n",
    "    J[4, 6] = -dq1 - dq3\n",
    "    J[4, 11] = 1\n",
    "    J[4] /= Cs\n",
    "    \n",
    "    J[5, 2] = -dq4\n",
    "    J[5, 3] = dq2\n",
    "    J[5, 5] = -dq2 - dq4\n",
    "    J[5, 6] = dq2 + dq4\n",
    "    J[5, 12] = -1\n",
    "    J[5] /= Cs\n",
    "        \n",
    "    J[6, 6] = -1 / Rp - dq1 - dq2 - dq3 - dq4\n",
    "    J[6, 2] = dq1 + dq4\n",
    "    J[6, 3] = -dq2 - dq3\n",
    "    J[6, 4] = -dq1 - dq3\n",
    "    J[6, 5] = dq2 + dq4\n",
    "    J[6] /= Cp\n",
    "    \n",
    "    J[7, 0] = -1 / Lh\n",
    "\n",
    "    J[8, 1] = -1 / Lh\n",
    "    \n",
    "    J[9, 0] = 0.5\n",
    "    J[9, 2] = -1\n",
    "    J[9, 9] = -Rg2\n",
    "    J[9] /= Ls2\n",
    "    \n",
    "    J[10, 0] = -0.5\n",
    "    J[10, 3] = 1\n",
    "    J[10, 10] = -Rg3\n",
    "    J[10] /= Ls3\n",
    "    \n",
    "    J[11, 1] = 0.5\n",
    "    J[11, 4] = -1\n",
    "    J[11, 11] = -Rg2\n",
    "    J[11] /= Ls2\n",
    "    \n",
    "    J[12, 1] = -0.5\n",
    "    J[12, 5] = 1\n",
    "    J[12, 12] = -Rg3\n",
    "    J[12] /= Ls3\n",
    "    \n",
    "    J[13, 0] = -1\n",
    "    J[13, 13] = -Ri - Rg1\n",
    "    J[13] /= Ls1\n",
    "    \n",
    "    J[14, 1] = -1\n",
    "    J[14, 14] = -Rc - Rg1\n",
    "    J[14] /= Ls1\n",
    "    \n",
    "    return J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_span = [0, 1e-3]\n",
    "y0 = np.zeros(15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/nmayorov/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:18: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/Users/nmayorov/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:19: RuntimeWarning: invalid value encountered in double_scalars\n",
      "/Users/nmayorov/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:21: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    }
   ],
   "source": [
    "res = solve_ivp(fun, t_span, y0, method='Radau')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(42382,)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.t.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.status"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 8. MEDAZKO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = 200\n",
    "k = 100\n",
    "v0 = 1\n",
    "c = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(t, y):\n",
    "    d = 1 / N\n",
    "    phi = 2 if t <= 5 else 0\n",
    "    j = 1 + np.arange(N)\n",
    "    alpha = 2 * (j * d  - 1) ** 3 / c**2\n",
    "    beta = (j * d - 1) ** 4 / c ** 2\n",
    "    \n",
    "    y = np.hstack((phi, 0, y, y[-2]))\n",
    "\n",
    "    j_2_p1 = 2 * j + 2\n",
    "    j_2_m3 = 2 * j - 2\n",
    "    j_2_m1 = 2 * j\n",
    "    j_2 = 2 * j + 1\n",
    "    \n",
    "    f = np.empty(2 * N)\n",
    "    f[::2] = alpha * (y[j_2_p1] - y[j_2_m3]) / (2 * d) + beta * (y[j_2_m3] - 2 * y[j_2_m1] + y[j_2_p1]) / d**2 - k * y[j_2_m1] * y[j_2]\n",
    "    f[1::2] = -k * y[j_2] * y[j_2_m1]\n",
    "    \n",
    "    return f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_span = [0, 20]\n",
    "y0 = np.zeros(2 * N)\n",
    "y0[1::2] = v0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res = solve_ivp(fun, t_span, y0, method='Radau')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(145,)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.t.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.status"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9. PLEI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m = 1 + np.arange(7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(t, u):\n",
    "    x = u[:7]\n",
    "    y = u[7:14]\n",
    "    x_dot = u[14:21]\n",
    "    y_dot = u[21:]\n",
    "    \n",
    "    dx = x[:, None] - x\n",
    "    dy = y[:, None] - y\n",
    "    denom = (dx ** 2 + dy ** 2) ** 1.5\n",
    "\n",
    "    with np.errstate(invalid='ignore'):\n",
    "        Fx = m * dx / denom\n",
    "        Fx = np.nan_to_num(Fx)\n",
    "        Fx = -np.sum(Fx, axis=1)\n",
    "        \n",
    "        Fy = m * dy / denom\n",
    "        Fy = np.nan_to_num(Fy)\n",
    "        Fy = -np.sum(Fy, axis=1)\n",
    "        \n",
    "    return np.hstack((x_dot, y_dot, Fx, Fy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t_span = [0, 3]\n",
    "x0 = np.array([3, 3, -1, -3, 2, -2, 2])\n",
    "y0 = np.array([3, -3, 2, 0, 0, -4, 4])\n",
    "x_dot0 = np.array([0, 0, 0, 0, 0, 1.75, -1.5])\n",
    "y_dot0 = np.array([0, 0, 0, -1.25, 1, 0, 0])\n",
    "u0 = np.hstack((x0, y0, x_dot0, y_dot0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res = solve_ivp(fun, [0, 3], u0, method='Radau')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(128,)"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.t.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.status"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}