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]
					trans_term = term.subs(y.diff(t, i), trans_term)
			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')
			algeq = dLaplace(eq, y0)
			algsoln = solve(algeq, L)[0]
			soln = InvLaplace(algsoln)
			return soln
		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))

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

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