Python (SymPy) Programming for solving Ordinary Differential Equations

This article demonstrates how to solve ordinary differential equations with SymPy. Please refer to the previous article for introductory stuffs of SymPy:

pianofisica.hatenablog.com

Indefinite integral
Ordinary Differential Equations
Linear ODE of 2nd order: Free Fall
One-dimensional harmonic oscillator
Airy equation

Indefinite integral

Let it start with reviewing an indefinite integral. The indefinite integral of a function $y=f(x)$

$\qquad\displaystyle{F(x)=\int dx \, f(x) +C}$

where $C$ is the constant of integral, is given as follows. As an example, let $\ f(x)=\exp(-ax)$ whose integral shall be obtained by running the following code in Python

import sympy as sp
sp.var('a, x')
sp.integrate(sp.exp(-a*x), x)

Let I be a function returning its indefinite integral of the input:

def I(f):
    return sp.integrate(f, x)

The example above is also found by applying this function I:

I(sp.exp(-a*x))

resulting in

$\qquad\displaystyle{\int dx \,e^{-ax} =\begin{cases} \ -\frac{e^{-ax}}{a}+C &\qquad (a\ne0)\\ \ x +C &\qquad (a=0)\end{cases}}$

For another example,

sp.var('n')
I(x**n)

one finds the indefinite integral of $f(x)=x^n$ as

$\qquad\displaystyle{\int dx \, x^n =\begin{cases} \ \frac{x^{n+1}}{n+1}+C &\qquad (n\neq-1)\\ \ \log x +C &\qquad (n=-1)\end{cases}}$

For $f(x)=\cos(nx)$

I(sp.cos(n*x))

one sees

$\qquad\displaystyle{\int dx \,\cos(nx) =\begin{cases} \ \frac{\sin(nx)}{n}+C &\qquad (n\neq0)\\ \ x+C &\qquad (n=0)\end{cases}}$

and so on.
　

Ordinary Differential Equations

Linear ODE of 1st order

As a simpler example, consider the following equation

$\qquad\displaystyle{\frac{dy}{dx}=-ay}$

where $y$ is an unknown function of $x$ and $a$ is a constant.

This ODE is easily solved by a method of separation of variables:

$\qquad\displaystyle{\int dx\, \frac{dy}{dx}\,\frac{1}{y}=-a\int dx}$

implies

$\qquad\displaystyle{\log y=-ax + c \quad (c:\text{constant of integral})}$

Taking logarithms on both sides, one may find

$\qquad\displaystyle{y=C\exp\left(-ax\right) \quad (C:\text{constant of integral})}$
　

dsolve: general solution

SymPy has a function to solve ODE "dsolve" which is run as follows:

import sympy as sp
sp.var('a, x')

y = sp.Function('y')(x)
eq = sp.Eq( sp.diff(y, x),  -a*y )

sp.dsolve(eq)

gives

$\qquad\displaystyle{y(x)=Ce^{-ax} \quad (C:\text{積分定数})}$

dsolve: solution with boundary condition

In order to impose (a) boundary condition(s), for example, $y(0)=1$ putting

sp.dsolve(eq, ics={y.subs(x,0):1})

shows

$\qquad\displaystyle{y(x)=e^{-ax}}$
　

Linear ODE of 2nd order: Free Fall

One-dimensional free fall is described by the following ODE (equation of motion)

$\qquad\displaystyle{m\frac{d^2x}{dt^2}=-mg}$

where $x \ [{\rm m}]$ coordinates the position of the mass $m \ [{\rm kg}]$ along the vertical direction (positive upwards) and $g \ [{\rm m/s^2}]$ is gravitational acceleration.

This equation is solved by Python running

import sympy as sp
sp.var('g, t')

x = sp.Function('x')(t)
eq1 = sp.Eq( sp.diff(x,t,2),  -g )

sp.dsolve(eq1)

resulting in

$\qquad\displaystyle{x(t)=-\frac12 g t^2+C_2 t+C_1}$

where $C_2$ is initial speed and $C_1$ is initial position.

Once imposing the initial conditions, e.g. $\ x(0)=0,\ \frac{dx}{dt}(0)=v$

sp.var('v')
sp.dsolve(eq1, ics={x.subs(t,0):0, sp.diff(x,t,1).subs(t,0):v})

gives

$\qquad\displaystyle{x(t)=-\frac12 g t^2+v t}$

One-dimensional harmonic oscillator

The position of mass $m$ shall be labeled by $x$ . The EoM is given by

$\qquad\displaystyle{m\frac{d^2x}{dt^2}=F_x \qquad F_x=-kx}$

which is implemented as

import sympy as sp
sp.var('m,k, t')

x = sp.Function('x')(t)
F = -k*x
eq2 = sp.Eq( m*sp.diff(x,t,2),  F )

sp.dsolve(eq2)

showing general solution

$\qquad\displaystyle{x(t)=C_{1} e^{-t \sqrt{-\frac{k}{m}}}+C_{2} e^{t \sqrt{-\frac{k}{m}}}}$

where $C_1, \ C_2$ are constants of integral.

Since both $k \ {\rm [ N/m]}$ and $m\ {\rm [ m]}$ are positive constants, introducing

$\qquad\displaystyle{i\omega=\sqrt{-\frac{k}{m}}}$

and applying the Euler formula $\ e^{i\theta}=\cos\theta+i\sin\theta$ , one may find the general solution is also written as

$\qquad\displaystyle{x(t)=A \cos (\omega t)+B\sin(\omega t)}$

Or, more explicitly, considering

$\qquad\displaystyle{\frac{d^2x}{dt^2}=-x}$

import sympy as sp
sp.init_printing()
sp.var('t')
x = sp.Function('x')(t)
eq3 = sp.Eq( sp.diff(x,t,2),  -x)
sp.dsolve(eq3)

one gets

$\qquad\displaystyle{x(t)=C_1 \sin (t)+C_2\cos(t)}$

Airy equation

Consider the Airy equation

$\qquad\displaystyle{\frac{d^2y}{dx^2}-xy=0}$

running the following code on Python

import sympy as sp
sp.init_printing()
sp.var('x')
y = sp.Function('y')(x)
eq4 = sp.Eq( sp.diff(y,x,2), x*y)
sp.dsolve(eq4)

one may obtain

$\qquad\displaystyle{y(x)=C_{1}{\rm Ai}(x)+C_{2}{\rm Bi}(x)}$

which is by definition given by linear combination of ${\rm Ai}$ and ${\rm Bi}$ .
　

　

Keywords: Differential equation, Python