Why so complicated?
Why can’t the code be as simple as:
# Random code I wrote just to illustrate the issue - won't work from sympy import symbols, solve from sympy.integrals.transforms import laplace_transform, inverse_laplace_transform s = symbols('s') t = symbols('t', positive = True) y = Function('y')(t) L = Function('L')(s) # input y0 = [2, 1] deq = y.diff(t,2) + y.diff(t,1) - 2*y - 4 algeq = laplace_transform(eq, s,t, noconds = True) algsoln = solve(algeq, L) soln = inverse_laplace_transform(algsoln, s, t, noconds = True) print("Differential Equation: ", deq, "=0") print("Solution: y(t) = ", soln)
The trouble arises when you try to find the Laplace Transform of derivative of y, which laplace_trasform function cannot handle yet. (Try it yourself – print(laplace_transform(y.diff(t), t, s))
So we had to take a shortcut. Since it was supposed to be taught to students and it had to be as simple as possible (we are allowed to cheat a bit on code’s integrity), we settled for the following, which inputs the transformed equation instead of the differential equation itself. The friend who proposed it came up with this code (with minor touches from me for better readability):
from sympy import * from sympy.integrals.transforms import laplace_transform, inverse_laplace_transform s, L = symbols('s, L') t= symbols ('t', positive =True) y0 = 2 y10 = 1 # storing y(0) and y’(0) Ly2 = s**2*L-s*y0-y10 Ly1= s*L-y0 Ly = L algeq = Eq(Ly2 + Ly1 - 2*Ly, laplace_transform(4, t, s, noconds = True)) algsoln = solve(algeq, L) soln = inverse_laplace_transform(algsoln, s, t, noconds = True) print("Solution: y(t) = ", soln)
We decided to settle for this until we came to the third problem: . The program took on an average 5-10 minutes to give the output! The problem was with the inverse laplace transform. (Try tweaking the above code for this problem).