Dynamical Systems
Last revision: January 23, 2021Let us discuss how to use Octave to solve systems of differential equations, and to plot phase diagrams of dynamical systems.
Summary of the basic theory
(Important note: this section is merely a reminder, not a substitute for a real discussion of dynamical systems. For details please refer to a serious book, such as Steven Strogatz’s Nonlinear Dynamics and Chaos.)
Let us consider a system of $n$ differential equations, such as
Given an initial condition, $(y_1(0),y_2(0),\cdots,y_n(0))$, and under some basic regularity conditions, there exists a single set of curves $(y_1(t),y_2(t),\cdots,y_n(t))$ satisfying the previous system of equations. For concreteness, we will usually work with $n=2$.
A very simple dynamical system is given by
whose integral curves are simply circumferences around the origin. Notice that differentiating again the first equation and substituting into the other we obtain
which is the differential equation for the harmonic oscillator. A more interesting example is provided by the van der Pol’s oscillator,
Numerically, these systems of equations can by solved using standard algorithms, such as Euler’s method or RungeKutta’s method. Luckily for us, Octave has those methods programmed.
Dynamical sytems in Octave
Let’s play a bit with Octave. We can define the differential equation using a function handle. For example, the harmonic oscillator is written as
The @ symbol is shorthand for the function arguments: time $t$ and
the values of $y$. Now, we make use of the ode45
function, which
employs an algorithm which is equivalent to an order 4
RungeKutta. Function ode45
requires the following argument:
 Function handle.
 Values of t where we want the solution evaluated.
 Initial condition.
So, for example, we write
The syntax [t,y] on the lefthandside means that the output should be placed on variables t and y. Notice that, as promised, the first argument is the function handle df. The second argument is the set of values of $t$ where we need the function defined, from 0 to 1 in steps of 0.01, so [0:0.01:1]. And the last, we provide the initial values at the first point, $y_1(0)=1$ and $y_2(0)=0$ in the form $[1,0]$.
Now, we have the final curve in $y$, which is a matrix, as we can readily check by printing it on the screen. In order to plot it, we do the following
But in order to visualize a full dynamical system we should be able to plot many such curves. We will write a short program that does precisely that:

Chooses randomly $N=1000$ random initial points in a window, such as $[0,1]\times [0,1]$.

Obtains the integral curves for a short time.

Plots all the results together.
This is the final code:
And the result should look as follows
For van der Pol’s oscillator, use
in order to obtain
And you can try one more:
in the range $[1,1]\times [1,1]$ to obtain
It is convenient to use a full Octave script, such as
Thus, using these numerical techniques can help us understand the structure of the dynamical systems. For example, we can spot critical points and guess their nature. Of course, numerics is never a substitute for real analytical work, when it can be done. Thus, we can use Octave either to get some intuition about a new dynamical system, or to check what we have already found using analytical tools.