{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A script to illustrate how to solve systems of ODEs numerically, and to draw a direction field for a pair of autonomous ODEs in Sage.\n",
    "\n",
    "#Consider a pair of ODEs x1' = f(t,x1,x2), x2' = g(t,x1,x2) where f(t,x1,x2) = x1-x2^2, g(t,x1,x2) = x1*x2+x1 (these happen \n",
    "#to be autonomous, but this isn't necessary for the numerical solver \"dsolve_system_rk4()\" below.)\n",
    "def f(t,x1,x2):\n",
    "    return(x1-x2^2)\n",
    "def g(t,x1,x2):\n",
    "    return(x1*x2+x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Numerical Solution: We can solve the system numerically with initial data x1(0) = 4, x2(0) = 1 on the interval 0 <= t <= 2\n",
    "#with the commands\n",
    "var('t x1 x2');\n",
    "sol = desolve_system_rk4([f(t,x1,x2),g(t,x1,x2)],[x1,x2],ics=[0,4,1],ivar=t,end_points=2,step=0.05)\n",
    "#The last argument \"step\" is the stepsize (can be omitted, defaults to 0.1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The output \"sol\" is a list of the form (t,x1,x2) at specified times. To isolate the values of x1 and plot, execute\n",
    "x1vals = [ [i,j] for i,j,k in sol]\n",
    "plt1 = list_plot(x1vals,plotjoined=True, color='red')\n",
    "show(plt1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Or to also plot x2 and overlay the graph of both x1(t) and x2(t)\n",
    "x2vals = [ [i,k] for i,j,k in sol]\n",
    "plt2 = list_plot(x2vals,plotjoined=True, color='blue')\n",
    "show(plt1+plt2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A parametric plot of the solution trajectory can be obtained with\n",
    "x1x2vals = [ [j,k] for i,j,k in sol]\n",
    "list_plot(x1x2vals,plotjoined=True, color='blue').show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Direction Fields: If the ODEs are autonomous, we can sketch a direction field on a range a < x1 < b, c < x2 < d, with or \n",
    "#without solution trajectories, by using the supplied command \"draw_auto_field.sage\".\n",
    "load('draw_auto_field.sage')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define f and g so they do not depend on t.\n",
    "def f(x1,x2):\n",
    "    return(x1-x2^2)\n",
    "def g(x1,x2):\n",
    "    return(x1*x2+x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Construct a direction field on the range -5 <= t <= 5, -5 <= u <= 5.\n",
    "draw_auto_field(f, g, [-5, 5, -5, 5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Construct a direction field with solutions satisfying x1(0) = 4, x2(0) = 1, and x1(0) = 3, x2(0) = -2. Use stepsize 0.01 \n",
    "#(last argument) to sketch trajectories.\n",
    "draw_auto_field(f, g, [-5, 5, -5, 5], [[4,1],[3,-2]],0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}