Solving Linear ODE using Python

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?

Leave a Reply

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

WordPress.com Logo

You are commenting using your WordPress.com 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.