KvN_NumOsc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Bad key text.latex.unicode in file /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle, line 112 ('text.latex.unicode : False # use \"ucs\" and \"inputenc\" LaTeX packages for handling')\n",
      "You probably need to get an updated matplotlibrc file from\n",
      "https://github.com/matplotlib/matplotlib/blob/v3.3.0/matplotlibrc.template\n",
      "or from the matplotlib source distribution\n",
      "\n",
      "Bad key savefig.frameon in file /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle, line 423 ('savefig.frameon : True')\n",
      "You probably need to get an updated matplotlibrc file from\n",
      "https://github.com/matplotlib/matplotlib/blob/v3.3.0/matplotlibrc.template\n",
      "or from the matplotlib source distribution\n",
      "\n",
      "Bad key pgf.debug in file /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle, line 444 ('pgf.debug           : False')\n",
      "You probably need to get an updated matplotlibrc file from\n",
      "https://github.com/matplotlib/matplotlib/blob/v3.3.0/matplotlibrc.template\n",
      "or from the matplotlib source distribution\n",
      "\n",
      "Bad key verbose.level in file /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle, line 475 ('verbose.level  : silent      # one of silent, helpful, debug, debug-annoying')\n",
      "You probably need to get an updated matplotlibrc file from\n",
      "https://github.com/matplotlib/matplotlib/blob/v3.3.0/matplotlibrc.template\n",
      "or from the matplotlib source distribution\n",
      "\n",
      "Bad key verbose.fileo in file /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle, line 476 ('verbose.fileo  : sys.stdout  # a log filename, sys.stdout or sys.stderr')\n",
      "You probably need to get an updated matplotlibrc file from\n",
      "https://github.com/matplotlib/matplotlib/blob/v3.3.0/matplotlibrc.template\n",
      "or from the matplotlib source distribution\n",
      "In /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n",
      "The text.latex.preview rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n",
      "In /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n",
      "The mathtext.fallback_to_cm rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n",
      "In /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: Support for setting the 'mathtext.fallback_to_cm' rcParam is deprecated since 3.3 and will be removed two minor releases later; use 'mathtext.fallback : 'cm' instead.\n",
      "In /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n",
      "The validate_bool_maybe_none function was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n",
      "In /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n",
      "The savefig.jpeg_quality rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n",
      "In /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n",
      "The keymap.all_axes rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n",
      "In /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n",
      "The animation.avconv_path rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n",
      "In /home/arurz/.local/lib/python3.6/site-packages/matplotlib/mpl-data/stylelib/_classic_test.mplstyle: \n",
      "The animation.avconv_args rcparam was deprecated in Matplotlib 3.3 and will be removed two minor releases later.\n",
      "┌ Info: Precompiling DifferentialEquations [0c46a032-eb83-5123-abaf-570d42b7fbaa]\n",
      "└ @ Base loading.jl:1273\n"
     ]
    }
   ],
   "source": [
    "using PyCall\n",
    "using PyPlot\n",
    "using DifferentialEquations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spy = pyimport(\"scipy.interpolate\")\n",
    "np = pyimport(\"numpy\")\n",
    "mp = pyimport(\"mpmath\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Δ = 1/2\n",
    "ω = 2.5\n",
    "\n",
    ";"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function k(t)\n",
    "    return Δ + cos(ω*t)\n",
    "end\n",
    "\n",
    "function k1(t,κ)\n",
    "    return 2*κ^2/cosh(κ*t)^2\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function osc!(du,u,p,t)\n",
    " du[1] = u[2]\n",
    " du[2] = (1/u[1]^3)-k(t)*u[1]\n",
    " du[3] = (1/u[1]^2)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u0 = [1.0;0.0;0.0]\n",
    "tspan = (0.0,20)\n",
    "prob = ODEProblem(osc!,u0,tspan)\n",
    "sol = solve(prob,Vern7(),dt=0.00000000001)\n",
    ";"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(np.min(sol.t), np.max(sol.t), 200)  \n",
    "\n",
    "spl1 = spy.make_interp_spline(sol.t, sol'[:,1], k=3)\n",
    "spl2 = spy.make_interp_spline(sol.t, sol'[:,2], k=3)\n",
    "spl3 = spy.make_interp_spline(sol.t, sol'[:,3], k=3)\n",
    "\n",
    "ps1 = spl1(t)\n",
    "ps2 = spl2(t)\n",
    "ps3 = spl3(t)\n",
    "\n",
    ";"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rc(\"text\", usetex=true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,(ax1,ax2,ax3)=plt.subplots(3,1,figsize=(5,5),sharex=true)\n",
    "plt.subplots_adjust(hspace = 0.1)\n",
    "\n",
    "ax1.plot(t,ps1,\"b-\", linewidth=2,label=L\"\\omega = 2.5\\qquad \\Delta=1/2\")\n",
    "ax2.plot(t,ps2,\"b-\", linewidth=2)\n",
    "ax3.plot(t,ps3,\"b-\", linewidth=2)\n",
    "\n",
    "ax1.set_aspect(\"auto\")\n",
    "ax2.set_aspect(\"auto\")\n",
    "ax3.set_aspect(\"auto\")\n",
    "\n",
    "ax1.tick_params(direction=\"out\",length=5,width=2,labelsize=14)\n",
    "ax2.tick_params(direction=\"out\",length=5,width=2,labelsize=14)\n",
    "ax3.tick_params(direction=\"out\",length=5,width=2,labelsize=14)\n",
    "\n",
    "ax1.grid(linewidth=1)\n",
    "ax2.grid(linewidth=1)\n",
    "ax3.grid(linewidth=1)\n",
    "\n",
    "ax1.set_ylabel(L\"\\rho(t)\",fontsize=16)\n",
    "ax2.set_ylabel(L\"\\dot{\\rho}(t)\",fontsize=16)\n",
    "ax3.set_ylabel(L\"\\omega_{\\rho}(t)\",fontsize=16)\n",
    "\n",
    "ax3.set_xlabel(L\"t\\ [s]\",fontsize=16,labelpad=0)\n",
    "xlim(0,20)\n",
    "\n",
    "ax1.legend(bbox_to_anchor=(0.3,0.3,0.5,1),loc=9,ncol=1, mode=\"expand\", borderaxespad=0,fontsize=11)\n",
    "\n",
    "show()\n",
    "\n",
    "savefig(\"osckvn.png\", transparent = \"true\", dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function β1(t,ρ0)\n",
    "    return convert(Float64,(ρ0/sol'[t,1])*mp.cos(sol'[t,3])+ρ0*sol'[t,2]*mp.sin(sol'[t,3]))\n",
    "end\n",
    "\n",
    "function β2(t,ρ0)\n",
    "    return convert(Float64,ρ0*sol'[t,1]*mp.sin(sol'[t,3]))\n",
    "end\n",
    "\n",
    "function β3(t,ρ0,dρ0)\n",
    "    return convert(Float64,(1/(ρ0*sol'[t,1])+dρ0*sol'[t,2])*mp.sin(sol'[t,3])+((dρ0/sol'[t,1])-(sol'[t,2]/ρ0))*mp.cos(sol'[t,3]))\n",
    "end\n",
    "\n",
    "function β4(t,ρ0,dρ0)\n",
    "    return convert(Float64,(sol'[t,1]/ρ0)mp.cos(sol'[t,3])-dρ0*sol'[t,1]*mp.cos(sol'[t,3]))\n",
    "end\n",
    "\n",
    ";"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function qc(q,p,t)\n",
    "    return q*β1(t,1)-p*β2(t,1)\n",
    "end\n",
    "    \n",
    "function pc(q,p,t)\n",
    "    return q*β3(t,1,0)+p*β4(t,1,0)\n",
    "end\n",
    "\n",
    ";"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qv = zeros(Float64,length(sol.t))\n",
    "for i in 1:length(sol.t)\n",
    "    qv[i] = qc(2,2,i)\n",
    "end\n",
    "pv = zeros(Float64,length(sol.t))\n",
    "for i in 1:length(sol.t)\n",
    "    pv[i] = pc(2,2,i)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = np.linspace(np.min(sol.t), np.max(sol.t), 200)  \n",
    "\n",
    "sp1 = spy.make_interp_spline(sol.t, qv, k=3)\n",
    "sp2 = spy.make_interp_spline(sol.t, pv, k=3)\n",
    "\n",
    "ps1 = sp1(t)\n",
    "ps2 = sp2(t)\n",
    "\n",
    ";"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figps = figure(figsize=(5,5))\n",
    "axps = gca(projection=\"3d\")\n",
    "plot3D(ps1,ps2,t1,\"b-\", linewidth=2)\n",
    "axps.view_init(30,-30)\n",
    "xlabel(L\"q\",fontsize=16)\n",
    "ylabel(L\"p\",fontsize=16)\n",
    "zlabel(L\"t\\ [s]\",fontsize=16)\n",
    "axps.tick_params(axis=\"both\",direction=\"out\",length=10,width=3,grid_linewidth=3,labelsize=14)\n",
    "axps.set_zticks(range(0,stop=20,step=5))\n",
    "axps.set_zticklabels([L\"0\",L\"5\",L\"10\",L\"15\",L\"20\"])\n",
    "show()\n",
    "tight_layout()\n",
    "\n",
    "savefig(\"oscckvn.png\", transparent = \"true\", dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.4.2",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.3.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}