Why Settle for Less?
Now, Laplace transform of nth derivative is given by:
![\mathscr{L}[y^{(n)}]=s^n\mathscr{L}[y] - s^{n-1} y(0) - \cdots - s y^{(n-2)}(0) - y^{(n-1)}(0)](https://s0.wp.com/latex.php?latex=%5Cmathscr%7BL%7D%5By%5E%7B%28n%29%7D%5D%3Ds%5En%5Cmathscr%7BL%7D%5By%5D+-+s%5E%7Bn-1%7D+y%280%29+-+%5Ccdots+-+s+y%5E%7B%28n-2%29%7D%280%29+-+y%5E%7B%28n-1%29%7D%280%29&bg=ffffff&fg=6a6c6e&s=0&c=20201002)
I decided to use this equation to generate the nth derivative, using for loop.
Now the function dLaplace is capable of finding transform of nth derivative, so, why stop at second order? I expanded it further, to the form it is now.
Now, since we are making it more general purpose, why not add a line or two to check if the input problem is solvable? A few more optimizations can be spotted there (like the use of ‘powers’ list).
from sympy import symbols, solve, Function, ode_order, classify_ode
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)
def dLaplace(expr, ics):
ans = 0
powers = [i for i in range(len(ics),-1,-1)]
for term in expr.args:
trans_term = 0
if y in term.atoms(Function): # Calculates L[y^(n)]
for i in powers:
if y.diff(t,i) in term.args or y.diff(t,i) == term:
trans_term += s**i * L
for j in range(0, i):
trans_term -= s**(i - j - 1) *ics[j]
powers.remove(i)
trans_term = term.subs(y.diff(t, i), trans_term)
else:
trans_term = laplace_transform(term,t,s,noconds=True)
ans += trans_term
return ans.simplify()
def InvLaplace(expr):
ans = 0
part_frac = expr.apart()
for term in part_frac.args:
inv_term = inverse_laplace_transform(term,s,t,noconds=True)
ans = ans + inv_term
return ans.simplify()
def dsolve_laplace(eq, ics):
if 'nth_linear_constant_coeff_variation_of_parameters' in classify_ode(eq):
if ode_order(eq,y) > len(ics):
raise Exception('Not enough initial conditions')
elif ode_order(eq, y) < len(ics):
raise Exception('More initial conditions provided')
else:
algeq = dLaplace(eq, y0)
algsoln = solve(algeq, L)[0]
soln = InvLaplace(algsoln)
return soln
else:
raise Exception('Can solve only linear differential equations.')
# input
y0 = [2, 1]
deq = y.diff(t,2) + y.diff(t,1) - 2*y - 4
print("Differential Equation: ", deq, "=0")
print("Solution: y(t) = ", dsolve_laplace(deq, y0))
There are chances for mistakes, if any, please do point it out in comments.
Leaving it here hoping that someone will find it useful sometime in future!
Mathematicos 