Visualisation of systems of first order ODEs

This page is intended to be a supplement to the lectures on systems of first order ODEs for the UNC course MATH383: First Course in Differential Equations. The general system of first order ODEs

$\mathbf{x}^\prime = \mathbf{P}(t) \mathbf{x} + \mathbf{g}(t)$ ,

with an initial condition
$\mathbf{x}(t_0) = \mathbf{x}_0$ ,

defines the general initial value problem for a system of first order ODEs. Here we will look at the case of homogeneous equations with constant coefficients, given by

$\mathbf{x}^\prime = \mathbf{A} \mathbf{x}$ ,

where the entries of the matrix $\mathbf{A}$ are constants. We will visualize solutions to these initial value problems with a Python script. For an introduction to using Python, have a look at Numerical Approximations to First Order ODEs.

Running the script

Download the script and rename it to “systems.py”. When you run the script the first time, you should get an animated plot that ends on a frame that looks like this:

On the plot:

The title gives the initial value problem being solved.
The black line shows the solution from a highly precise numerical approximation. This is precise enough that it should be visually indistinguishable from the analytic solution.
The red line shows an analytic solution that the user provides.

The buttons above the plot allow you to pan, zoom in and out, and save the figure.

Modifying the problem

The code between

 #BEGIN CODE TO BE EDITED#

and

 #END CODE TO BE EDITED#

controls the initial value problem, the analytic solution, and how the solutions are displayed.

System matrix

The lines

 P = [1.0, 0.0,
      0.0, -1.0]

describe the system matrix, which is currently

$P = \left( \begin{array}{cc} 1 & 0 \\ 0 & -1\end{array} \right).$

Let’s change our matrix to

$P = \left( \begin{array}{cc} 0.1 & 1 \\ -1 & 0.1\end{array} \right)$ ,

by changing the lines to

 P = [0.1, 1.0,
      -1.0, 0.1]

Initial condition

The line

x_init = [0.1, 1.0]

describes the initial condition, which is currently

$\mathbf{x}_0 = \mathbf{x}(0) = \begin{pmatrix} 0.1 \\ 1.0 \end{pmatrix}.$

Let’s change our initial condition to

$\mathbf{x}_0 = \mathbf{x}(0) = \begin{pmatrix} 0.1 \\ 0.1 \end{pmatrix}$

by changing the line to

x_init = [0.1, 0.1]

End time

Currently, the line

t_n = 1.0

sets the visualization to go from t=0.0 to t=1.0. Lets change the end time to t=15.0 by changing the line to

t_n = 15.0

Analytic solution

The lines

x0 = x_init[0]*exp(t)
x1 = x_init[1]*exp(-t)

describe the analytic solution that we put in. In this case, if we call our initial condition

$\mathbf{x}_0 = \begin{pmatrix} \alpha_1 \\ \alpha_2 \end{pmatrix}$

these lines describe the analytic solution as

$\mathbf{x}(t) = \begin{pmatrix} \alpha_1 e^t \\ \alpha_2 e^{-t} \end{pmatrix}$

After the changes we’ve just made to the problem, this solution is no longer valid. Once we’ve covered systems with complex eigenvalues, we’ll be able to calculate the analytic solution and place it here. For now, we will remove the plotting of the analytic solution by changing the line

do_plot_analytic = True

do_plot_analytic = False

Plot after changes

After making these changes and running the script again, the plot should look like this:

Modifying other aspects

The line

t_anim = 2.5

sets the number of real-world seconds to spend animating the trajectory. If we don’t want any animation at all, and just want to display the final time, we can set it to

t_anim = None

Feel free to change anything else in the code to your heart’s content.