Solver_Explained.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "246235f0-e6b1-4c4c-879d-55010c4fe5ff",
   "metadata": {},
   "source": [
    "# Multi-Step Methods\n",
    "In these methods we need multiple previous points of the function we want to integrate. i. e. for $s\\text{-step}$ we need $f(t_{-s+1},y_{-s+1}) \\dots f(t_{0},y_{0})$ to solve for the equation $\\dot{y}(t)=f(t,y)$. \n",
    "\n",
    "## The Adams-Bashforth Method\n",
    "It is a multi-step method of solving first order ordinary differential equations. This method is as follows.\n",
    "\n",
    "\n",
    "First we denote that our *ODE* is given in the following form.\n",
    "\n",
    "$$ \\dot{y} = f(t,y) $$\n",
    "\n",
    "Assuming we can approximate the function on the *RHS* with interpolating it as a polynomial. Given that we have $s$ points the polynomial will be of the $(s-1)^{st}$ order. The procedure has been done in the following steps.\n",
    "\n",
    "$$p(t) \\approx f(t,y),\\quad p(t_{m})=f\\bigl(t_{m},y_{m}\\bigr)\\equiv f_{m} $$\n",
    "\n",
    "Then we can use the polynomial to get an approximation for the next step.\n",
    "\n",
    "$$ y_{n+1} = y_{n} + \\int_{t_{n}}^{t_{n+1}} f(t,y)\\mathrm{d}t\\approx y_{n} + \\int_{t_{n}}^{t_{n+1}} p(t)\\mathrm{d}t $$\n",
    "\n",
    "The procedure of obtaining the desired polynomial is written step-by-step:\n",
    "\n",
    "\n",
    "$$ p(t) = \\sum_{m=0}^{s} p_{m}f_{n-s+1+m} $$\n",
    "\n",
    "$$ p_{m}(t) = \\prod_{\\substack{l=0 \\\\ l\\neq m}}^{s-1} \\frac{t-t_{n-s+1+l}}{t_{n-s+1+m}-t_{n-s+1+l}} $$\n",
    "\n",
    "Then we will use it inthe integral (we integrate both sides of the ode and we get an approximation for $y_{n+1}$)\n",
    "\n",
    "$$ y_{n+1} = y_{n} + \\int_{t_{n}}^{t_{n+1}} f(t,y)\\mathrm{d}t\\approx y_{n} + \\int_{t_{n}}^{t_{n+1}} p(t)\\mathrm{d}t = y_{n} + \\int_{t_{n}}^{t_{n+1}} \\sum_{m=0}^{s-1}p_{m}(t)f_{n-s+1+m}\\mathrm{d}t= y_{n} + h\\times\\sum_{m=0}^{s-1} b_{m}f_{n-s+1+m} $$\n",
    "\n",
    "The integrals on the RHS of the equation can be reduced to some coefficients . \n",
    "\n",
    "$$ b_{m} = \\frac{1}{h}\\int_{t_{n}}^{t_{n+1}} p_{m}(t)\\mathrm{d}t = \\frac{1}{h}\\int_{t_{n}}^{t_{n+1}} \\prod_{\\substack{l=0 \\\\ l\\neq m}}^{s-1} \\frac{t-t_{n-s+1+l}}{t_{n-s+1+m}-t_{n-s+1+l}}\\mathrm{d}t $$\n",
    "\n",
    "\n",
    "$$ b_{n} = \\int_{h\\times n}^{h\\times(n+1)} \\prod_{\\substack{l=0 \\\\ l\\neq m}}^{s-1} \\frac{h\\times\\bigl(\\frac{t}{h}-(n-s+1+l)\\bigr)}{h\\times\\bigl((n-s+1+m)-(n-s+1+l)\\bigr)}\\mathrm{d}\\frac{t}{h} \\\\\n",
    "\\frac{t}{h}\\equiv \\tilde{t}\\Rightarrow b_{m} = \\int_{n}^{n+1} \\prod_{\\substack{l=0 \\\\ l\\neq m}}^{s-1} \\frac{\\tilde{t}-n+s-1-l}{m-l}\\mathrm{d}\\tilde{t} $$\n",
    "\n",
    "And the final form of the coefficients of Adams-Bashforth method is calculated as below:\n",
    "\n",
    "$$ \\tilde{t}\\equiv \\tau + n\\Rightarrow b_{m} = \\int_{0}^{1} \\prod_{\\substack{l=0 \\\\ l\\neq m}}^{s-1} \\frac{\\tau+s-1-l}{m-l}\\mathrm{d}\\tau $$\n",
    "\n",
    "So we get our method as:\n",
    "\n",
    "$$ y_{n+1} = y_{n} + h\\times\\sum_{m=0}^{s-1} b_{m}f_{n-s+1+m} $$\n",
    "\n",
    "## The Adams-Moulton Method\n",
    "It is a multi-step method of solving first order ordinary differential equations. This method is as follows. In this method we also take $f(t_{n+1},y_{n+1})$ into account.\n",
    "\n",
    "\n",
    "First we denote that our ode is given in the following form.\n",
    "\n",
    "$$ \\dot{y} = f(t,y) $$\n",
    "\n",
    "We can approximate the function on the *RHS* with interpolating it as a polynomial. Given that we have $s$ points the polynomial will be of the $(s-1)^{st}$ order. The procedure has been done in the following steps.\n",
    "\n",
    "$$q(t) \\approx f(t,y),\\quad q(t_{m})=f\\bigl(t_{m},y_{m}\\bigr)\\equiv f_{m} $$\n",
    "\n",
    "Then we can use the polynomial to get an approximation for the next step.\n",
    "\n",
    "$$ y_{n+1} = y_{n} + \\int_{t_{n}}^{t_{n+1}} f(t,y)\\mathrm{d}t\\approx y_{n} + \\int_{t_{n}}^{t_{n+1}} q(t)\\mathrm{d}t $$\n",
    "\n",
    "The procedure of obtaining the desired polynomial is written step-by-step:\n",
    "\n",
    "\n",
    "$$ q(t) = \\sum_{m=0}^{s} q_{m}f_{n-s+1+m} $$\n",
    "\n",
    "**NOTE** that this method generates an $s^{th}$ order polynomial, becuase it is using $s+1$ points to interpolate.\n",
    "\n",
    "$$ q_{m}(t) = \\prod_{\\substack{l=0 \\\\ l\\neq m}}^{s} \\frac{t-t_{n-s+1+l}}{t_{n-s+1+m}-t_{n-s+1+l}} $$\n",
    "\n",
    "Then we will use it inthe integral (we integrate both sides of the ode and we get an approximation for $y_{n+1}$)\n",
    "\n",
    "$$ y_{n+1} = y_{n} + \\int_{t_{n}}^{t_{n+1}} f(t,y)\\mathrm{d}t\\approx y_{n} + \\int_{t_{n}}^{t_{n+1}} q(t)\\mathrm{d}t = y_{n} + \\int_{t_{n}}^{t_{n+1}} \\sum_{m=0}^{s}q_{m}(t)f_{n-s+1+m}\\mathrm{d}t= y_{n} + h\\times\\sum_{m=0}^{s} c_{m}f_{n-s+1+m} $$\n",
    "\n",
    "The integrals on the RHS of the equation can be reduced to some coefficients . \n",
    "\n",
    "$$ c_{m} = \\frac{1}{h}\\int_{t_{n}}^{t_{n+1}} q_{m}(t)\\mathrm{d}t = \\frac{1}{h}\\int_{t_{n}}^{t_{n+1}} \\prod_{\\substack{l=0 \\\\ l\\neq m}}^{s} \\frac{t-t_{n-s+1+l}}{t_{n-s+1+m}-t_{n-s+1+l}}\\mathrm{d}t $$\n",
    "\n",
    "\n",
    "$$ c_{n} = \\int_{h\\times n}^{h\\times(n+1)} \\prod_{\\substack{l=0 \\\\ l\\neq m}}^{s} \\frac{h\\times\\bigl(\\frac{t}{h}-(n-s+1+l)\\bigr)}{h\\times\\bigl((n-s+1+m)-(n-s+1+l)\\bigr)}\\mathrm{d}\\frac{t}{h} \\\\\n",
    "\\frac{t}{h}\\equiv \\tilde{t}\\Rightarrow c_{m} = \\int_{n}^{n+1} \\prod_{\\substack{l=0 \\\\ l\\neq m}}^{s} \\frac{\\tilde{t}-n+s-1-l}{m-l}\\mathrm{d}\\tilde{t} $$\n",
    "\n",
    "And the final form of the coefficients of Adams-Moulton method is calculated as below:\n",
    "\n",
    "$$ \\tilde{t}\\equiv \\tau + n\\Rightarrow c_{m} = \\int_{0}^{1} \\prod_{\\substack{l=0 \\\\ l\\neq m}}^{s} \\frac{\\tau+s-1-l}{m-l}\\mathrm{d}\\tau $$\n",
    "\n",
    "So we get our method as:\n",
    "\n",
    "$$ y_{n+1} = y_{n} + h\\times\\sum_{m=0}^{s} c_{m}f_{n-s+1+m} $$\n",
    "\n",
    "## Adams-Bashforth-Moulton\n",
    "Obtaining the algebraic equation for $y_{n+1}$ in terms of the other components can be tricky for **Adams-Moulton** method, thus people often use a predictor-corrector method comprised of **Adams-Bashforth** method as the predictor and **Adams-Moulton** method of the same order as the predictor for correcting the answer obtained by the predictor. So the method is as follows:\n",
    "\n",
    "$$ \\tilde{y}_{n+1} = y_{n} + h\\times\\sum_{m=0}^{s-1} b_{m}f_{n-s+1+m} $$\n",
    "\n",
    "$$ y_{n+1} = y_{n} + h\\times\\Bigl(c_{s-1}f(t_{n+1},\\tilde{y}_{n+1})+\\sum_{m=0}^{s-2} c_{m}f_{n-s+2-m}\\Bigr) $$\n",
    "\n",
    "\n",
    "Where the coefficients are defined as:\n",
    "\n",
    "$$\\forall m \\in [0,\\dots,s-1],\\quad b_{m} = \\int_{0}^{1} \\prod_{\\substack{l=0 \\\\ l\\neq m}}^{s-1} \\frac{\\tau-n+s-1-l}{m-l}\\mathrm{d}\\tau $$\n",
    "\n",
    "and,\n",
    "\n",
    "$$\\forall m \\in [1,\\dots,s],\\quad c_{m} = \\int_{0}^{1} \\prod_{\\substack{l=1 \\\\ l\\neq m}}^{s} \\frac{\\tau+s-1-l}{m-l}\\mathrm{d}\\tau $$\n",
    "\n",
    "or\n",
    "\n",
    "$$\\forall m \\in [0,\\dots,s-1],\\quad c_{m} = \\int_{0}^{1} \\prod_{\\substack{l=0 \\\\ l\\neq m}}^{s-1} \\frac{\\tau+s-2-l}{m-l}\\mathrm{d}\\tau $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "a660d80c-800c-488f-9429-1b34fed6921c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "Please enter the degree of predictor-corrector method you want to be calculated:  6\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import sympy\n",
    "\n",
    "s = np.int64(input(\"Please enter the degree of predictor-corrector method you want to be calculated: \"))\n",
    "m = 1\n",
    "t = sympy.symbols('t')\n",
    "P = sympy.Function('P')(t, s, m)\n",
    "def P(t,s,m):\n",
    "    res = 1\n",
    "    for l in range(s):\n",
    "        if l!= m:\n",
    "            res*= (t+s-1-l)/(m-l)\n",
    "    return res\n",
    "Q = sympy.Function('Q')(t, s, m)\n",
    "def Q(t,s,m):\n",
    "    res = 1\n",
    "    # This loop could be from 1 to s or 0 to s-1, I chose the latter\n",
    "    for l in range(s):\n",
    "        if l!= m:\n",
    "            res*= (t+s-2-l)/(m-l)\n",
    "    return res\n",
    "\n",
    "# 'p's are the polynomial interpolations of the predictor\n",
    "p = sympy.symarray(('p'), shape=s)\n",
    "b = sympy.symarray(('b'),shape=s)\n",
    "# 'q's are the polynomial interpolations of the corrector\n",
    "q = sympy.symarray(('q'), shape=s)\n",
    "c = sympy.symarray(('c'),shape=s)\n",
    "bnum = np.zeros(s).reshape((1,s))\n",
    "bdenom = np.zeros(s).reshape((1,s))\n",
    "cnum = np.zeros(s).reshape((1,s))\n",
    "cdenom = np.zeros(s).reshape((1,s))\n",
    "\n",
    "for m in range(s):\n",
    "    p[m]=P(t,s,m)\n",
    "    b[m]=sympy.Rational(sympy.integrate(p[m],[t,0,1]).simplify())\n",
    "    bnum[0,m] = b[m].numerator\n",
    "    bdenom[0,m] = b[m].denominator\n",
    "    q[m]=Q(t,s,m)\n",
    "    c[m]=sympy.Rational(sympy.integrate(q[m],[t,0,1]).simplify())\n",
    "    cnum[0,m] = c[m].numerator\n",
    "    cdenom[0,m] = c[m].denominator\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "d5250aee-ac35-45d9-b058-d4499c28456c",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-t*(-t - 4)*(-t/2 - 3/2)*(-t/3 - 2/3)*(-t/4 - 1/4)/5,\n",
       "       -t*(-t - 3)*(-t/2 - 1)*(-t/3 - 1/3)*(t + 5)/4,\n",
       "       -t*(-t - 2)*(-t/2 - 1/2)*(t/2 + 5/2)*(t + 4)/3,\n",
       "       -t*(-t - 1)*(t/3 + 5/3)*(t/2 + 2)*(t + 3)/2,\n",
       "       -t*(t/4 + 5/4)*(t/3 + 4/3)*(t/2 + 3/2)*(t + 2),\n",
       "       (t/5 + 1)*(t/4 + 1)*(t/3 + 1)*(t/2 + 1)*(t + 1)], dtype=object)"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "5dcb060d-457d-424e-802c-70b38d565328",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-95/288, 959/480, -3649/720, 4991/720, -2641/480, 4277/1440],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "5badb4e9-309c-4579-809a-c09be131a358",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-t*(1/5 - t/5)*(-t - 3)*(-t/2 - 1)*(-t/3 - 1/3)/4,\n",
       "       -t*(1/4 - t/4)*(-t - 2)*(-t/2 - 1/2)*(t + 4)/3,\n",
       "       -t*(1/3 - t/3)*(-t - 1)*(t/2 + 2)*(t + 3)/2,\n",
       "       -t*(1/2 - t/2)*(t/3 + 4/3)*(t/2 + 3/2)*(t + 2),\n",
       "       (1 - t)*(t/4 + 1)*(t/3 + 1)*(t/2 + 1)*(t + 1),\n",
       "       t*(t/5 + 4/5)*(t/4 + 3/4)*(t/3 + 2/3)*(t/2 + 1/2)], dtype=object)"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "f30b290e-a345-4dfc-a523-ac0560a83fb9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3/160, -173/1440, 241/720, -133/240, 1427/1440, 95/288],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "07bdc39a-b122-4914-85e7-dcd140b216b0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  -95.,   959., -3649.,  4991., -2641.,  4277.]])"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bnum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "id": "e0cccfef-b89f-41a6-9000-a30596e67cf5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 288.,  480.,  720.,  720.,  480., 1440.]])"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bdenom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "f10688e7-498e-409f-b0f6-2a434e0c8aba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   3., -173.,  241., -133., 1427.,   95.]])"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cnum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "ce28c3c0-ff56-4ae4-865a-aa674ee87438",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 160., 1440.,  720.,  240., 1440.,  288.]])"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cdenom"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ffc8a09-666b-4c5f-bf56-ec9a43f4ab15",
   "metadata": {},
   "source": [
    "# Solving a system of equation via Adams-Bashforth and Adams-Bashforth-Moulton method\n",
    "\n",
    "## The Kuramoto model\n",
    "\n",
    "It is a special model of coupled phase oscillators. The model is defined in the following form:\n",
    "\n",
    "$$ \\dot{\\theta}_{i} = \\omega_{i} + \\frac{K}{N}\\sum_{j=0}^{N} A_{ij}\\sin\\bigl(\\theta_{j}-\\theta_{i}\\bigr) $$\n",
    "\n",
    "The dot indicates the derivative with respect to time, the $\\omega_{i}$ is the intrinsic frequency of the $i^{th}$ oscillator, $K$ is the coupling constant, $N$ is the number of oscillators and, $A$ is the adjacency matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "5b95b41b-7f7b-45f4-a100-a188f7ec3faf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "Please enter the number of oscillators:  12\n",
      "Please enter the frequency of the oscillators:  10\n",
      "Please enter the coupling constant:  10.\n",
      "Please enter the time step size:  0.5\n",
      "Please enter the integration method:  rk4\n"
     ]
    }
   ],
   "source": [
    "global ab; global abm\n",
    "ab = np.zeros(s).reshape((1,s))\n",
    "abm = np.zeros(s).reshape((1,s))\n",
    "for i in range(s):\n",
    "    ab[0,i] = 1.*np.int64(bnum[0,i])/np.int64(bdenom[0,i])\n",
    "    abm[0,i] = 1.*np.int64(cnum[0,i])/np.int64(cdenom[0,i])\n",
    "\n",
    "\n",
    "N = np.int64(input(\"Please enter the number of oscillators: \"))\n",
    "omega0 = np.float64(input(\"Please enter the frequency of the oscillators: \"))\n",
    "K = np.float64(input(\"Please enter the coupling constant: \"))\n",
    "dt = np.float64(input(\"Please enter the time step size: \"))\n",
    "\n",
    "global one_sixth\n",
    "one_sixth = 1./6\n",
    "def func(Theta,w0,N,k):\n",
    "    theta = Theta.reshape((1,N))\n",
    "    res = np.ones_like(theta)*w0\n",
    "    for i in range(N):\n",
    "        res[0,i] += k*np.sum(np.sin(theta[0,:]-theta[0,i]),axis = 0)\n",
    "    return res\n",
    "\n",
    "def rk4_step(Theta,w0,N,k,dt):\n",
    "    rk1 = func(Theta,w0,N,k)\n",
    "    rk2 = func(Theta + 0.5*dt*rk1,w0,N,k)\n",
    "    rk3 = func(Theta + 0.5*dt*rk2,w0,N,k)\n",
    "    rk4 = func(Theta + dt*rk3,w0,N,k)\n",
    "    res = Theta + one_sixth*(rk1+2*(rk2+rk3)+rk4)\n",
    "    return res, rk1\n",
    "\n",
    "def ab_step(dt,Theta_dot_History):\n",
    "    res = np.dot(ab,Theta_dot_History)*dt\n",
    "    return res\n",
    "\n",
    "def abm_step(dt,Theta_dot_History):\n",
    "    dummy = Theta_dot_History\n",
    "    for i in range(s-1):\n",
    "        dummy[i,:] = dummy[i+1,:]\n",
    "    dummy[s-1,:] = np.dot(ab,Theta_dot_History)*dt\n",
    "    res = np.dot(abm,dummy)*dt\n",
    "    return res\n",
    "    \n",
    "def initiate_kuramoto(N,omega0,K,dt):\n",
    "    # Initiating the phases as a splay with random fluctuation\n",
    "    Theta0 = np.arange(0., 2*np.pi, 2*np.pi/N).reshape((1,N))+np.random.uniform(0.,1.,N).reshape((1,N))*1.0e-010\n",
    "    Theta = np.zeros_like(Theta0)\n",
    "    Theta_dot = np.zeros_like(Theta0)\n",
    "    Theta_History = np.zeros((s,N))\n",
    "    Theta_dot_History = np.zeros((s,N))\n",
    "    #Creating the history\n",
    "    Theta_History[0,:] = Theta0[0,:]\n",
    "    k = 1.*N/K\n",
    "    for i in range(s-1):\n",
    "        Theta_History[i+1,:], Theta_dot_History[i,:] = rk4_step(Theta_History[i,:],omega0,N,k,dt)\n",
    "    Theta = Theta_History[s-1,:]\n",
    "    Theta_dot_History[s-1,:] = func(Theta,omega0,N,k)\n",
    "    Theta_dot = Theta_dot_History[s-1,:]\n",
    "    return Theta0,Theta,Theta_dot,Theta_History,Theta_dot_History,k\n",
    "\n",
    "def Kuramoto(Final_Time,N,omega0,K,dt,method):\n",
    "    Number_of_Steps = np.int64(Final_Time/dt)\n",
    "    Time = np.linspace(0.,Final_Time,Number_of_Steps+1)\n",
    "    Re_R = np.zeros(Number_of_Steps+s)\n",
    "    Im_R = np.zeros(Number_of_Steps+s)\n",
    "    rho = np.zeros(Number_of_Steps+s)\n",
    "    global k\n",
    "    Theta0,Theta,Theta_dot,Theta_History,Theta_dot_History,k = initiate_kuramoto(N,omega0,K,dt)\n",
    "    for i in range(s):\n",
    "        Re_R[i] = np.sum(np.cos(Theta_History[i,:]),axis=0)/N\n",
    "        Im_R[i] = np.sum(np.sin(Theta_History[i,:]),axis=0)/N\n",
    "    if method==\"rk4\":\n",
    "        for i in range(Number_of_Steps):\n",
    "            Theta, Theta_dot = rk4_step(Theta,omega0,N,k,dt)\n",
    "            for j in range(s-1):\n",
    "                Theta_History[j,:] = Theta_History[j+1,:]\n",
    "                Theta_dot_History[j,:] = Theta_dot_History[j+1,:]\n",
    "            Theta_History[s-1,:] = Theta\n",
    "            Theta_dot_History[s-1,:] = Theta       \n",
    "            Re_R[i+s] = np.sum(np.cos(Theta_History[s-1,:]),axis=0)/N\n",
    "            Im_R[i+s] = np.sum(np.sin(Theta_History[s-1,:]),axis=0)/N            \n",
    "    elif method==\"ab\":\n",
    "        for i in range(Number_of_Steps):\n",
    "            Theta = ab_step(dt,Theta_dot_History)\n",
    "            Theta_dot = func(Theta,omega0,N,k)\n",
    "            for j in range(s-1):\n",
    "                Theta_History[j,:] = Theta_History[j+1,:]\n",
    "                Theta_dot_History[j,:] = Theta_dot_History[j+1,:]\n",
    "            Theta_History[s-1,:] = Theta\n",
    "            Theta_dot_History[s-1,:] = Theta       \n",
    "            Re_R[i+s] = np.sum(np.cos(Theta_History[s-1,:]),axis=0)/N\n",
    "            Im_R[i+s] = np.sum(np.sin(Theta_History[s-1,:]),axis=0)/N                    \n",
    "    elif method==\"abm\":\n",
    "        for i in range(Number_of_Steps):\n",
    "            Theta = abm_step(dt,Theta_dot_History)\n",
    "            Theta_dot = func(Theta,omega0,N,k)\n",
    "            for j in range(s-1):\n",
    "                Theta_History[j,:] = Theta_History[j+1,:]\n",
    "                Theta_dot_History[j,:] = Theta_dot_History[j+1,:]\n",
    "            Theta_History[s-1,:] = Theta\n",
    "            Theta_dot_History[s-1,:] = Theta          \n",
    "            Re_R[i+s] = np.sum(np.cos(Theta_History[s-1,:]),axis=0)/N\n",
    "            Im_R[i+s] = np.sum(np.sin(Theta_History[s-1,:]),axis=0)/N            \n",
    "    else:\n",
    "        print(\"Invalid Method!!!\")\n",
    "    rho = np.sqrt(np.multiply(Re_R,Re_R)+np.multiply(Im_R,Im_R))\n",
    "    # phi = np.tan(np.divide(Im_R,Re_R))\n",
    "    return rho,Re_R,Im_R   \n",
    "\n",
    "Final_Time = 1000.\n",
    "method = input(\"Please enter the integration method: \")\n",
    "rho,Re_R,Im_R = Kuramoto(Final_Time,N,omega0,K,dt,method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "id": "08635050-3463-4871-ba7a-c8b4ca535e36",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.43605276e-11, 7.28203724e-10, 3.69263352e-08, ...,\n",
       "       4.60765993e-01, 4.08154081e-01, 1.12623650e-01])"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "ca3cb481-af9f-46e9-abce-07e0e4b203ca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(rho[80006:])/20000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "id": "c324b0b6-fa9d-4cc9-b0e7-eff1cb2cd303",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1081.4356319361536"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(rho)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fddd916-67ff-47d8-8e1a-82b1c7a0c2e9",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}