31.3: Ordinary Differential Equations

Ordinary Differential Equations (ODEs) are yet another for of linear dynamical systems and are a scientific model used in a wide range of problems of the basic form:

\[\dot{x} = Ax \nonumber \]

These are equations such that the is the instantaneous rate of change in \(x\) (i.e. \(\dot{x}\) is the derivative of \(x\)) is dependent on \(x\). Many systems can be modeled with these types of equations.

Here is a quick video that introduces the concepts of Differential Equations. The following is a good review of general ODEs.

from IPython.display import YouTubeVideo
YouTubeVideo("8QeCQn7uxnE",width=640,height=360, cc_load_policy=True)

Now consider an ODE as a system of linear equations:

\[\dot{x_t} = A x_t \nonumber \]

Based on the current \(x\) vector at time \(t\) and the matrix \(A\), we can calculate the derivative at \(\dot{x}\) at time \(t\). Once we know the derivative, we can increment the time to by some small amount \(dt\) and calculate a new value of \(x\) as follows:

\[x_{t+1} = x_t + \dot{x_t}dt \nonumber \]

Then we can do the exact sequence of calculations again for \(t+2\). The following function has the transition matrix (\(A\)), the starting state vector (\(x_0\)) and a number of time steps (\(N\)) and uses the above equations to calculate each state and return all of the \(x\) statues:

The following code generates a trajectory of points starting from x_0, applying the matrix \(A\) to get \(x_1\) and then applying \(A\) again to see how the system progresses from the start state.

%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
import sympy as sym
sym.init_printing()

def traj(A, x, n):
    dt = 0.01
    x_all = np.matrix(np.zeros((len(x),n)))   # Store all points on the trajectory
    for i in range(n):  
        x_dot = A*x         # First we transform x into the derrivative
        x = x + x_dot*dt    # Then we estimate x based on the previous value and a small increment of time.
        x_all[:,i] = x[:,0] 
    return x_all

For example the following code uses the matrix \(A= \begin{bmatrix}1 & 1 \\ 1 & -2\end{bmatrix}\) and the starting point (0,0) over 50 timesteps to get a graph:

A = np.matrix([[1,1],[1,-2]])
x0 = np.matrix([[1],[1]])

x_all = traj(A, x0, 50)
plt.scatter(np.asarray(x_all[0,:]),np.asarray(x_all[1,:]))

plt.scatter(list(x0[0,:]),list(x0[1,:])) #Plot the start point as a refernce

Do This

Let \(A= \begin{bmatrix}2 & 3 \\ 4 & -2\end{bmatrix}\)

Write a loop over the points \((−1.5,−1), (−1,2)(−1,2)\) and plot the results of the traj function:

A = np.matrix([[2,3],[4,-2]])
x0 = np.matrix([[1.5, -1.5, -1, 1, 2],[1, -1, 2, -2, -2]])

# Put your code here

Do This

Let \(A= \begin{bmatrix}6 & -1 \\ 1 & 4\end{bmatrix}\)

Write a loop over the points \((−1.5,−1), (−1,2)(−1,2)\) and plot the results of the traj function:

# Put your code here

Do This

Let \(A= \begin{bmatrix}5 & 2 \\ -4 & 1\end{bmatrix}\)

Write a loop over the points \((−1.5,−1), (−1,2)(−1,2)\) and plot the results of the traj function:

# Put your code here