{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This worksheet illustrates using the Laplace transform to solve linear system of ODEs (and also just using Sage's\n",
    "#desolve_system() command).\n",
    "\n",
    "#Example System: Consider the linear constant-coefficient nonhomogeneous linear system of ODEs\n",
    "var('t')\n",
    "x1 = function('x1')(t);\n",
    "x2 = function('x2')(t);\n",
    "de1 = diff(x1,t) + 5*x1 - 6*x2 == -2*sin(t) + 10*cos(t) - 6*exp(t)\n",
    "de2 = diff(x2,t) + 3*x1 - x2 == 6*cos(t)\n",
    "#for functions x1(t) and x2(t), with initial data x1(0) = 2, x2(0) = 1.\n",
    "\n",
    "#NOTE: The laplace command (and desolve) seems to Be more successful when all dependent variables appear on the left side of \n",
    "#the ODEs, all independent variables on the right as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#To solve, Laplace transform both sides of both equations\n",
    "var('s');\n",
    "de1lap = laplace(de1,t,s)\n",
    "de2lap = laplace(de2,t,s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Substitute in the initial conditions, and (for convenience) let X1 = laplace(x1(t),t,s) and X2 = laplace(x2(t),t,s):\n",
    "var('X1 X2');\n",
    "de1lap2 = de1lap.substitute(x1(0)==2,x2(0)==1,laplace(x1(t),t,s)==X1,laplace(x2(t),t,s) == X2)\n",
    "de2lap2 = de2lap.substitute(x1(0)==2,x2(0)==1,laplace(x1(t),t,s)==X1,laplace(x2(t),t,s) == X2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Solve for the transforms X1 and X2\n",
    "sol = solve([de1lap2,de2lap2],[X1,X2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Assign transforms to variables X1sol and X2sol:\n",
    "X1sol = X1.substitute(sol[0]); print(X1sol)\n",
    "X2sol = X2.substitute(sol[0]); print(X2sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Inverse transform to find the solutions\n",
    "x1sol(t) = inverse_laplace(X1sol,s,t)\n",
    "x2sol(t) = inverse_laplace(X2sol,s,t)\n",
    "print(\"x1(t) = \",x1sol(t))\n",
    "print(\"x2(t) = \",x2sol(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Solution via dsolve: The solution can be obtained with Sage's dsolve_system() command:\n",
    "sol = desolve_system([de1, de2], [x1,x2], ics=[0,2,1]);\n",
    "print(\"Solution via desolve_system() is \",sol)"
   ]
  },
  {
   "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
}