Better Code
Even now, we are inputting the transformed equation instead of giving the differential equation itself. I played around a bit to fix this. If laplace_transform cannot handle derivatives, let’s take things to our own hands:
The added function dLaplace finds the Laplace Transforms of differential equation up to 2nd order.
from sympy import symbols, solve
from sympy.integrals.transforms import laplace_transform, inverse_laplace_transform
s, y, y2, y1, L = symbols('s, y, y2, y1, L')
t = symbols('t', positive = True)
def dLaplace(expr):
ans = 0
for term in expr.args:
if y2 in term.free_symbols:
ans += term.subs(y2, s**2*L - s*y0-y10)
elif y1 in term.free_symbols:
ans += term.subs(y1, s*L - y0)
elif y in term.free_symbols:
ans += term.subs(y, L)
else:
trans_term = laplace_transform(term,t,s,noconds=True)
ans += term.subs(term, trans_term)
return ans
def InvLaplace(expr):
ans = 0
part_frac = expr.apart() # find partial fractions
for term in part_frac.args:
inv_term = inverse_laplace_transform(term,s,t,noconds=True)
ans = ans + inv_term
return(ans.simplify())
# input
y0 = 2
y10 = 1
deq = y2 + y1 - 2*y - 4
# solve
algeq = dLaplace(deq)
algsoln = solve(algeq, L)[0]
soln = InvLaplace(algsoln)
#print solution
print("Differential Equation: ", deq,"= 0\nSolution: y(t) = ", soln)
That settles it for us. Or does it? There is still something missing! Why are we using y2 and y1? Why not use differential itself?