Solving Linear ODE using Python

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)

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!

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.