SymPy/Differential Equations

Introduction

This page focuses on using SymPy to find both the symbolic and the numeric solutions to differential equations obtained from electric circuits. The examples are available via a Google Drive at Basic Dsolve. If you have the Google Colaboratory app installed in your browser, you should be able to run the notebooks in your browser. If not, you can either install it or install a local program that can read notebook files such as Jupyter Notebooks, available with Anaconda. This page is based on SymPy/Simultaneous Equations, which is likely a good place to start before going through this page!

Very Basic Example

The example code below assumes you are either running code in a notebook of some kind or typing it into a console. If you want to follow along, the code is in the dsolve_1 notebook in the Basic Dsolve Google Drive folder.

Imagine you have the equation $$2\frac{dy(t)}{dt} + 3 y(t) = 4$$ with the initial condition $$y(0)=5$$, and you want to solve for $$y(t)$$. You can do this in SymPy as follows:

Initialization

You will need to import sympy to use symbolic calculations. Pratt Pundit will generally give that module the nickname sym:

import sympy as sym

You will likely want your notebook or script to try to use $$\LaTeX$$ to render the results, so add:

sym.init_printing(use_latex=True)

as well.

Declare Variables and Functions

You will need to let Python know what your variables and functions are. For the example equation above, the unknown variable is $$t$$ and the unknown function is $$y$$; you can declare them with the following lines of code:

t = sym.var('t')
y = sym.Function('y')

Note that the Function command starts with a capital F.

Define Equations

To define a derivative in SymPy, use sym.Derivative(FUN, VAR), where FUN is a function (with an argument, usually VAR) and VAR is the variable you are taking the derivative with respect to. For instance, sym.Derivative(y(t), t) represents $$\frac{dy(t)}{dt}$$. If you need to take a higher-order derivative, you can (a) append the independent variable to the argument list, (b) append the order of the derivative to the argument list, or (c) give a tuple with the independent argument and the order. For example, to get $$\frac{d^2y(t)}{dt^2}$$, you can use sym.Derivative(y(t), t, t), sym.Derivative(y(t), t, 2), or sym.Derivative(y(t), (t, 2)). See Defining Derivatives at SymPy for further examples.

When it comes to solving things later, there are two different ways to define an equation:

• If you have an expression that evaluates to other than zero, you can use sym.Eq(left, right) for that. Given that we are trying to solve $$2\frac{dy(t)}{dt} + 3 y(t) = 4$$, you can write that in code as:

eqn1 = sym.Eq(2 * sym.Derivative(y(t), t) + 3 * y(t), 4)

• If the equation is equal to zero, you do not explicitly have to include the "=0" part; you can simply define a variable to hold on to the expression. For instance, in this case, we can re-write the equation we are trying to solve as $$ax-d=0$$ and just define a variable to hold on to the left side:

eqn1 = 2 * sym.Derivative(y(t), t) + 3 * y(t) - 4

Solve Equations

There are two fundamentally different ways to solve one differential equation. If you give the solver the equation and the function to solve for, it will return an equation with the solution in it:

soln1 = sym.dsolve(eqn1, y(t))
soln1

will produce an equation that contains $$y(t)=\frac{C_1e^{-\frac{3t}{2}}}{3}+\frac{4}{3}$$. On the other hand, if you give the solver a list with an equation in it and a list with the function in it, it will return a list with the equation in it:

soln2 = sym.dsolve([eqn1], [y(t)])
soln2

will produce the list $$\left[ y(t)=\frac{C_1e^{-\frac{3t}{2}}}{3}+\frac{4}{3} \right]$$. This latter version is going to be much more useful if we are going to do anything else with the solution, so you should generally make sure to give dsolve two lists. Note that unlike solving linear algebra equations, if you give a list of differential equations, you must also give a list with the functions.

Initial Conditions

To incorporate initial conditions, you will give the dsolve command a dictionary of function values at particular times. For example, to solve our sample equation with $$y(0)=5$$, you will include the initial condition by adding ics={y(0):5} to the arguments.

soln3 = sym.solve([eqn1], [y(t)], ics = {y(0): 5})
soln3

will produce $$\left[ y(t)=\frac{4}{3} + \frac{11e^{-\frac{3t}{2}}}{3} \right]$$

Note that the initial condition does not have to be at time 0; if you know that $$y(6)=7$$, you can use that as well:

soln4 = sym.solve([eqn1], [y(t)], ics = {y(6): 7})
soln4

will produce $$\left[ y(t)=\frac{4}{3} + \frac{17e^9e^{-\frac{3t}{2}}}{3} \right]$$

Basic Symbolic Example

For some simple systems of differential equations, SymPy can actually solve with symbolic coefficients and undefined forcing functions. The result will be an unevaluated integral, but still - pretty impressive! Here is an example that solves $$a\,\frac{dy(t)}{dt} + b\,y(t)=x(t)$$ assuming $$y(t_0)=y_0$$. If you want to follow along, the code is in the dsolve_2 notebook in the Basic Dsolve Google Drive folder.

import sympy as sym
sym.init_printing(use_latex=True)
a, b, t, y0, t0 = sym.var('a, b, t, y0, t0')
x = sym.Function('x')
y = sym.Function('y')
eqn1 = sym.Eq(a*sym.Derivative(y(t), t) + b * y(t), x(t))
eqns = [eqn1]
soln = sym.dsolve(eqns, [y(t)], ics = {y(t0): y0})
soln

The solution is:

$$ \left[ y{\left(t \right)} = \left(y_{0} e^{\frac{b t_{0}}{a}} - \frac{\int x{\left(t_{0} \right)} e^{\frac{b t_{0}}{a}}\, dt_{0}}{a}\right) e^{- \frac{b t}{a}} + \frac{e^{- \frac{b t}{a}} \int x{\left(t \right)} e^{\frac{b t}{a}}\, dt}{a}\right] $$

Making Substitutions

Making substitutions into lists of equations is a little complicated - you need to go into each element in the list and make the substitutions there. Given that, there is a helper function in sym_helper called makesubseqnlist(EQN LIST, VARS TO SUB, VALS TO SUB) that will take a list of equations, a list of variables (or functions), and a list of values (or functions). The code for makesubseqnlist(EQN LIST, VARS TO SUB, VALS TO SUB) is

def makesubseqnlist(equations, vartosub, valtosub):
    sublist = list(zip(vartosub, valtosub))
    esub = []    
    for eqn in equations:
      esub.append(eqn.subs(sublist))
    return esub

and you can get by by adding the following to near the top of your notebook or script:

try:
  from sym_helper import *
except: # Get sym_helper.py from GitHub
  import requests
  url = 'https://tinyurl.com/mrgdemos/SymPy_Files/sym_helper.py'
  open('sym_helper.py', 'w').write(requests.get(url).text)
  from sym_helper import *

If we specifically want to set $$a=2$$, $$b=3$$, and $$x(t)=4\,\cos(5t)$$, and then solve that where $$y(0)=6$$, we can:

eqns2 = makesubseqnlist(eqns, [a, b, x(t)], [2, 3, 4*sym.cos(5*t)])
soln2 = sym.dsolve(eqns2, [y(t)], ics = {y(0): 6})

yields:

$$ \left[ y{\left(t \right)} = \frac{40 \sin{\left(5 t \right)}}{109} + \frac{12 \cos{\left(5 t \right)}}{109} + \frac{642 e^{- \frac{3 t}{2}}}{109}\right] $$

Making a Dictionary

For linear algebra, the results came in a dictionary that had each of the unknowns as keys and their solutions as variables. Since dsolve produces a list, versus a ductionary, it will be useful to convert that list to a dictionary for making substitutions later. sym_helper has a make_dict(LIST) command that will turn a list of equations into a dictionary. The code for make_dict is:

def makedict(sollist):
    d = {}
    for eqn in sollist:
        d[eqn.lhs] = eqn.rhs
    return d

and so the code you can add to your script is:

def makedict(sollist):
soln = makedict(solnlist)

and, as long as sym_helper has been imported, you can run it by writing code like:

solndict = makedict(soln2)

which yields:

$$ \left\{ y{\left(t \right)} : \frac{40 \sin{\left(5 t \right)}}{109} + \frac{12 \cos{\left(5 t \right)}}{109} + \frac{642 e^{- \frac{3 t}{2}}}{109}\right\} $$

Plotting Results

Plotting will require using lambdify to turn one or more functions into a lambda function and then setting up your plot. Here is the code for this example that does all those things (making a plot that runs from $$t$$=0 to 10 sec); note the additional required imports:

import matplotlib.pyplot as plt
import numpy as np
yfun = sym.lambdify(t, solndict[y(t)])
fig = plt.figure(num=1, clear=True)
ax = fig.add_subplot(1, 1, 1)
tn = np.linspace(0, 10, 1000)
ax.plot(tn, yfun(tn))

Lessons Learned

• dsolve is good for up to two coupled ODEs. Once there are three or more...not so much (even if all the coefficients are ints or floats). Once there are three or more (linear) coupled differential equations, consider going the Laplace route.
• If your equations are a combination of differential equations a non-differential equations, dsolve expects all the differential equations to come first.
• If you are using certain functions from SymPy that aren't exactly the same in NumPy, when you lambdify things you will need to provide the translation

Examples

• See SymPy/Differential Equations/RC Example for a complete example; this will reference IVP_RCDemo in the Google Drive folder.
• See SymPy/CAD for a list of SymPy examples from a free textbook. The ones most relevant to this page are from Chapter 5.
• There are several other files in the Basic Dsolve Google Drive folder:

• • Coupled ODEs (just to see how SymPy handles those!)
    • Coupled2
    • Coupled3
    • deq_coupled: symbolically solve two coupled first-order ODEs; substitute numbers later to make graph. Needs some editing - does not completely follow the general template.
    • deq_coupled_num: symbolically solve two coupled first-order ODEs with numbers substituted in before solving. Needs some editing - does not completely follow the general template.
  • First-order differential equations:
    • IVP_demo1st: symbolic solution to 1st order ODE with a constant forcing function; substitutes numbers later to make graph
    • IVP_demo1stfun: symbolic solution to 1st order ODE with a non-constant forcing function; substitutes numbers later to make graph (including forcing function on graph)
    • IVP_demo1st_pulse: symbolic solution to 1st order ODE with a non-constant forcing function comprised of step and ramp functions; substitutes numbers later to make graph (including forcing function on graph)
      • This one shows extra code needed to translate step and ramp functions into a lambdify function
  • Second-order differential equations:
    • IVP_demo2ndcon: symbolic solution to 2nd order ODE with a constant forcing function; substitutes numbers later to make graph. Needs some editing - does not completely follow the general template.
    • IVP_demo2ndfun: symbolic solution to 2nd order ODE with a non-constant forcing function; substitutes numbers later to make graph. Needs some editing - does not completely follow the general template.