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

Python part 2, solving ordinary differential equations


Numerical solution of differential equations.

The simplest autonomous ordinary differential equation of interest in science is the exponential decay equation

\[ \begin{gather*} dn/dt = r n \quad \text{ where $r < 0$.} \end{gather*} \]

Implement the following script

# Numerical integration of the exponential decay model
#       dn/dt = r n, where r < 0

from numpy import linspace, array
import scipy.integrate
from scipy import exp
from matplotlib.pyplot import \
        plot, figure, text, show, savefig, xlabel, ylabel

r = -0.05 # decay_rate per day

def vector_field(X, t):
    n = X[0]
    return array([ r * n ])

time_start = 0.
time_end = 10.
numsteps = 20

observation_times = linspace(time_start, time_end, numsteps)

X_initial = array([32.])

X = scipy.integrate.odeint(vector_field, X_initial, observation_times)

X_final_calc = X[-1, 0]
X_final_exct = X_initial[0]*exp(r*time_end)
X_final_error = abs( X_final_calc - X_final_exct )
error_str = "Final Error = %1.5g"%(X_final_error)

options = {'fontsize': 18}

plot(observation_times, X[:, 0],'x-')
xlabel('Time (days)',options)
ylabel('Mass (micrograms)',options)

Analysis of the dampled pendulum

A pendulum on a rod, hanging from a pivot with weak friction should move according to the following system of equations, according to Newton's laws.

\[ \begin{align*} \frac{d \theta}{dt} &= \omega, \\ \frac{d \omega}{dt} &= -\sin(\theta) - a \omega. \end{align*} \]

The angle is measured as different from the rest position. In our case, let's assume the friction coefficient is \(a=0.05\).