Solving Linear ODE using Python

The Complete Program

The following program solves $n^{th}$ order differential equation:

Problem: $y''+y' - 2y = 4, y(0)=2, y'(0) = 1$

```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):
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()	# 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()

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))
```

This site uses Akismet to reduce spam. Learn how your comment data is processed.