{
"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
}