Eberly College of Science Mathematics Department, Center for Interdisciplinary Mathematics
REU banner

Python part 3, analyzing ordinary differential equations

First, a reminder that the material here builds on the previous material.

An interlude in Sympy for symbolic computations

Like Maxima(free), Maple(good), Mathematica(bad), python can also do symbolic mathematical calculations, thanks to the sympy module. To use it, we have to import both the functions and a set of symbols. We use permissive import so we don't have to type "sympy." and the beginning of everything.

>>> from sympy import *
>>> from sympy.abc import *
>>> eq = 3*x**2 + 3 - 5*sin(2*x)
>>> eq
>>> deq = eq.diff(x)
>>> deq
>>> ideq = Integral(deq,x)
>>> ideq
>>> ideq.doit()

So sympy is a free way to do all your calculus homework if you learn how to use it. (but user beware -- sympy is still in development and has plenty of bugs.I'm one of it's developers, so I know ;) If you use the right fonts, everything will come out looking very much like LaTex equations do. In FACT, sympy will automatically convert your equations into LaTex markup for you!

>>> J = Integral( exp(x)/(1+x**2), (x,-1,2))
>>> J
>>> pprint(J)
>>> print latex(J)

The function pprint does "pretty printing" of mathematical expressions to make them easier to read on the screen. Unfortunately, it still messes things up a little sometimes. In general, such type setting is un-reliable, which is why we still use LaTex typesetting.

Here's another example.

>>> A = Matrix([[1,2],[3,0]])
>>> A
>>> pprint(A)
>>> A**-1

Multiplication of sympy matrices works the way you might expect matrix multiplication to work. If we multiply a matrix by its inverse, we expect to get the identity matrix, and we do.

>>> A * (A**-1)
>>> (A**-1) * A

Part A. Watt's governor

Diagram of Watt's centrifugal flying-ball governor

Diagram of Watt's centrifugal flying-ball governor


Watt's governor was the first automatic control system used at the beginning of the industrial revolution for steam engines. It opened up whole new fields of both engineering and dynamical system analysis. In particular, James Maxwell showed that these governors could become unstable. Here are there equations, which we'll study a little.

\[ \begin{align*} \frac{dx}{dt} =& y \\ \frac{dy}{dt} =& z^{2} \sin{\left (x \right )} \cos{\left (x \right )} - \sin{\left (x \right )} - f y \\ \frac{dz}{dt} =& a \left(\cos{\left (x \right )}-b\right) \\ \text{parameter values} \quad & a = 1, \quad b = 1/2 \\ \text{initial condition} \quad & x(0) = 1.0472, \quad y(0) = 0., \quad z(0) = 1.386 \end{align*} \] Here, the state \((x, y, z) \in [0,\pi] \times (-\infty, \infty) \times [0, \infty)\). where

The parameters \((a, b, f) \in [0,\infty) \times [0,a] \times [0,\infty)\). where


  1. I. A. Vyshnegradskii and J. C. Maxwell independently showed that if the friction \(f\) is large enough, then the governor settles on a stable control. Plot timeseries for x,y,z when \(f = 4/5\) over 100 seconds.

  2. But if f is small enough, then the governor "hunts" for the best level without ever settling! Plot timeseries for x,y,z when \(f = 1/2\), and describe how the behavior differs from the previous case.

  3. Enter the vector field system into a column vector in sympy. Use sympy to solve for the steady-states of the system. Then, numerically estimate the steady-states for \((a,b,f) = (1,1/2,4/5)\) and \((a,b,f) = (1,1/2,1/2)\).

  4. Create a new function called "jacobian" to calculate the Jacob's derivative matrix of a vector field

    jacobian = lambda es, xs : Matrix([[ e.diff(x) for x in xs] for e in es])

    Use your new function to calculate the Jacobian of vector field for Watt's governor.

  5. Using the results of your calculates so far, create a new python the determines the sign of the largest eigenvalue of the Jacobian when evaluated at the steady-state solution. Numerically estimate the steady-states eigenvalues at steady-state for parameter sets \((a,b,f) = (1,1/2,4/5)\) and \((a,b,f) = (1,1/2,1/2)\). In which case is there an exponentially-growing part to solutions near the steady-state? At what value of \(f\) does this eigenvalue switch from having positive real part instead of a negative real part?

Part B. Behavior Change in an epidemic

A simple extension of classic epidemic theory is a case where resistance to disease is behavioral, and individuals who become sick tend to abandon their protective behaviors. This kind of phenomena can be described with the system \[ \begin{align} \frac{dS}{dt} &= \mu - \beta S I + f \gamma I + a R - v S - \mu S, \\ \frac{dI}{dt} &= \beta (S+\sigma R) I - \gamma I - \mu I, \\ \frac{dR}{dt} &= - \sigma \beta R I + (1-f) \gamma I - a R + v S - \mu R. \end{align} \] where \(S\) is the fraction of susceptibles, \(I\) is the fraction of infecteds, and \(R\) is the fraction of recovered people (\(1=S+I+R\)). Calculate all steady-states when \(f=9/10\), \(\gamma=61\), \(a=1/6\), \(\mu =1/78\), \(\beta = 570\), \(\nu = 1/3\), and determine their local stability.

Part C. Autocatalytic reactions

Once upon a time, it was believed that chemical reactions would always got to equilibrium straight-away, but some visionary chemists eventually showed that, in fact, chemical reactions can oscillate. The simpliest mathematical model of this phenomena is given by the bi-molecular mass-action reaction equations \[ \begin{align*} \frac{dp}{dt} =& a - p + p^2 q \\ \frac{dq}{dt} =& b - p^2 q \end{align*} \] where \(a\) and \(b\) are positive constants.

  1. Determine the steady-states of this system

  2. Find values of \(a\) and \(b\) for which the system is unstable and oscillates. Show with a time-series plot.

  3. Find different values of \(a\) and \(b\) for which the system is stable and converges to equilibrium. Show with a time-series plot.