{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This is a script to illustrate Laplace transforms in Sage, and using Laplace transforms to solve ODEs.\n",
    "\n",
    "#To Laplace transform an expression, e.g., t^2, declare \"s\" and \"t\" as symbolic variable and use the \"laplace\" command:\n",
    "var('s t')\n",
    "laplace(t^2,t,s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#More generally, to compute the Laplace transform of f(t) execute \"laplace(f(t),t,s)\". The second argument \"t\" indicates the independent\n",
    "#variable in f(t) and the third argument \"s\" indicates that the transform should have independent variable \"s\".\n",
    "\n",
    "#To inverse transform an expression use \"inverse_laplace()\", for example,\n",
    "inverse_laplace(2/s^3,s,t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Example: Solving a first-order ODE. Consider the ODE u'(t) = -2*u(t)+t with initial data u(0) = 3. To solve\n",
    "\n",
    "#Step 1: Define u(t) as a symbolic function and define the ODE\n",
    "var('s t');\n",
    "u = function('u')(t);\n",
    "ode = diff(u,t) == -2*u + t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Step 2: Laplace transform both sides of the ODE (which Sage does automatically):\n",
    "lapode = laplace(ode,t,s);\n",
    "print(lapode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Sage uses the notation \"laplace(u(t),t,s)\" for the Laplace transform of u(t). Let us replace this with the notation \"U\" (not necessary, just more\n",
    "#aesthetically pleasing) and also substitute in the initial data u(0) = 3.\n",
    "var('U');\n",
    "lapode2 = lapode.subs(u(0)==3,laplace(u(t),t,s)==U);\n",
    "print(lapode2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Sage may complain above  that the substitutions are deprecated and will not work #in a future release, so be aware. You \n",
    "#can always cut-and-paste to do the substitutions manually.\n",
    "\n",
    "#Step 3: We can now solve for U, which Sage returns in a list (enclosed in square brackets)\n",
    "Usol = solve(lapode2, U);\n",
    "print(Usol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Step 4: The last step is to inverse Laplace transform \"Usol\" to find the solution, which we call \"usol(t)\".\n",
    "sol = inverse_laplace(Usol[0],s,t);\n",
    "usol(t) = sol.rhs();\n",
    "print(\"The solution is \", usol(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#An easy alternative to the above steps is to use the desolve_laplace() command:\n",
    "usol(t) = desolve_laplace(ode,u,[0,3]);\n",
    "print(usol(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Example: A second-order ODE. Consider the second-order ODE\n",
    "ode = diff(u,t,2) + 4*diff(u,t) + 3*u == sin(t)\n",
    "#with u(0) = 1 and u'(0) = 2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Step 1: Laplace transform both sided of the ODE\n",
    "lapode = laplace(ode,t,s);\n",
    "print(lapode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Step 2: Substitute in the initial data and (optionally) use \"U\" to replace \"laplace(u(t),t,s)\". Again, Sage may complain about\n",
    "#the substitution.\n",
    "Du = diff(u,t); #Define derivative of u as a function\n",
    "lapode2 = lapode.subs(u(0)==1,Du(0)==2,laplace(u(t),t,s)==U);\n",
    "print(lapode2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Step 3: Solve for the Laplace transform U:\n",
    "Usol = solve(lapode2, U);\n",
    "print(Usol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Step 4: The last step is to inverse Laplace transform \"Usol\" to find the solution, which we call \"usol(t)\".\n",
    "sol = inverse_laplace(Usol[0],s,t);\n",
    "usol(t) = sol.rhs();\n",
    "print(\"The solution is \", usol(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Again, an easy alternative to the above steps is to use the desolve_laplace() command:\n",
    "usol(t) = desolve_laplace(ode, u, [0,1,2]);\n",
    "print(\"The solution really is \", usol(t));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The Heaviside and Dirac Delta Functions:  The Heaviside function in Sage is\n",
    "heaviside(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For example,\n",
    "plot(heaviside(t),(t,-2,2)).show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#HOWEVER, Sage does not (currently) have the internal ability to Laplace transform Heaviside or Dirac delta functions. To do\n",
    "#this we can call an external prorgram, e.g.,Scientific Python. Try\n",
    "laplace(heaviside(t),t,s,algorithm='sympy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The laplace transform returns the inverse transform as a function of s as well as conditions on s under #which the improper \n",
    "#integral defining the transform converges. We are interested in the transform (the first argument)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The Dirac delta function is\n",
    "dirac_delta(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Example: The mass in a damped spring-mass system with mass 1 kg, damping constant 2 newtons per meter per second, and \n",
    "#spring constant 5 newtons per meter is subjected to a constant force of 2 newtons for 0 < t < 4 seconds, at which point \n",
    "#the force drops to zero. At time t = 6 seconds the mass is subjected to a hammer blow with total impulse 8 newton-meters.\n",
    "#Find the motion of the mass if the mass starts at rest and at equilibrium, and plot this motion for 0 < t < 20 seconds.\n",
    "\n",
    "#The relevant ODE is\n",
    "ode = diff(u,t,2) + 2*diff(u,t) + 5*u == 2 - 2*heaviside(t-4) + 8*dirac_delta(t-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We can transform the left side using Sage's usual laplace commmand. This command recognizes the linearity of the ODE and breaks\n",
    "#the transform up into s^2*U + 2*s*U + 5*U, plus terms involving the initial conditions.\n",
    "lapode_left = laplace(ode.lhs(),t,s);\n",
    "print(lapode_left)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Then subsitute in the initial data, and use \"U\" for the transform\n",
    "Du = diff(u,t);\n",
    "lapode_left2 = lapode_left.subs(u(0)==0,Du(0)==0,laplace(u(t),t,s)==U);\n",
    "print(lapode_left2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Transform the right side of the ODE using a call to Scientific Python:\n",
    "lapode_right = laplace(ode.rhs(),t,s,algorithm='sympy');\n",
    "print(lapode_right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#And solve for U (using lapode_right[0] to pick off the transform as a function of s)\n",
    "sol = solve(lapode_left2==lapode_right[0],U);\n",
    "Usol = U.subs(sol[0]);\n",
    "print(Usol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Finally, inverse transform using the Scientific Python option (this can be a bit slow).\n",
    "usol(t) = inverse_laplace(Usol,s,t,algorithm=\"sympy\");\n",
    "print(usol(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A plot\n",
    "plot(usol(t), (t,0,20)).show()"
   ]
  },
  {
   "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
}