Example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Playing with Runge-Kutta methods ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rkkit import *\n",
    "from sage.repl.ipython_kernel.interact import interact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us import somme predefined Runge-Kutta method descriptions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "for formula in (\"lobatto4.sage\",\"radau5.sage\",\"sdirk3.sage\",\"gauss4.sage\",\"radau2a.sage\",\"rk4.sage\",\"sdirk5.sage\"):\n",
    "    load(\"./methods/\"+formula)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choose one formula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#RK=Radau5\n",
    "#RK=Lobatto4\n",
    "#RK=SDIRK3\n",
    "#RK=Gauss4\n",
    "#RK=Radau2a\n",
    "#RK=RK4\n",
    "RK=SDIRK5()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SDIRKmethod, order 5\n"
     ]
    }
   ],
   "source": [
    "print(RK.Title)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and define the formula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "F=RKformula(RK)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, now let-us check different properties of the formula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.is_explicit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.is_A_stable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.is_L_stable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.is_stiffly_accurate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.is_algebraically_stable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.conserve_quadratic_invariants()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html>\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{-\\frac{28}{3} z^{4} - \\frac{32}{3} z^{3} + 128 z^{2} + 256 z - 1024}{z^{5} - 20 z^{4} + 160 z^{3} - 640 z^{2} + 1280 z - 1024}\\]</html>"
      ],
      "text/latex": [
       "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{-\\frac{28}{3} z^{4} - \\frac{32}{3} z^{3} + 128 z^{2} + 256 z - 1024}{z^{5} - 20 z^{4} + 160 z^{3} - 640 z^{2} + 1280 z - 1024}$$"
      ],
      "text/plain": [
       "(-28/3*z^4 - 32/3*z^3 + 128*z^2 + 256*z - 1024)/(z^5 - 20*z^4 + 160*z^3 - 640*z^2 + 1280*z - 1024)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%display latex\n",
    "F.stability_function()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html>\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left(4, 5\\right)\\right]\\]</html>"
      ],
      "text/latex": [
       "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left(4, 5\\right)\\right]$$"
      ],
      "text/plain": [
       "[(4, 5)]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.poles_of_stability_function()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the limit of the stability domain on $\\mathbb{R}^-$. It can be $-\\infty$, for example if the method is A-stable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html>\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\infty\\]</html>"
      ],
      "text/latex": [
       "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\infty$$"
      ],
      "text/plain": [
       "-Infinity"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.stability_on_real_negative_axis()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html>\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}4\\]</html>"
      ],
      "text/latex": [
       "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}4$$"
      ],
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.order_of_stability_function()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The true order (computed using rooted trees):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html>\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}4\\]</html>"
      ],
      "text/latex": [
       "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}4$$"
      ],
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.order()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let us plot the stability domain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4277acf8f6c14459aaf7ce45a90d243f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Interactive function <function P at 0x7fc7cf5539d0> with 2 widgets\n",
       "  WindowSize: FloatSlider(value=4.6, descri…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact\n",
    "def P(WindowSize=(0.1,10,0.5),Translate=(-50,50,5)):\n",
    "    RKplot(F,fill=True,ncurves=2,Enlarge=WindowSize,TranslateX=Translate).show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the order star (slow):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "867df88462a44499a0ef7fb5d3b13d18",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Interactive function <function PStar at 0x7fc7c64f4430> with 2 widgets\n",
       "  WindowSize: FloatSlider(value=4.51, d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact\n",
    "def PStar(WindowSize=(0.01,10,0.5),Translate=(-50,50,5)):\n",
    "    RKplot(F,fill=True,ncurves=1,type=\"star\",Enlarge=WindowSize,TranslateX=Translate).show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.4",
   "language": "sage",
   "name": "sagemath"
  },
  "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}