Solving Differential Equations with Mathematica - Part I: Time Series

Download Mathematica notebook [right click / save as]

In this series, I will summarize some of the quickest methods to solve ordinary differential equations (ODEs hereafter) in Mathematica. To illustrate this, I will make use of the Lorenz dynamical system which is a set of three coupled ODEs. This dynamical system is notorious for its chaotic behavior and has given Lorenz the opportunity to be one of the first scientists who discovered chaos. In the process, I will also illustrate some of the techniques used to investigate chaotic behavior in dynamical systems; specifically, sensitivity to initial conditions, frequency domain transformations using the fast Fourier transform, and Poincaré sections.

The Lorenz equations are given bywhere x, y, and z are functions of time and sigma, rho, and beta are control parameters determined a priori. One could implement a fourth order Runga-Kutta method with adaptive time stepping to solve the above set of equations (and I would really recommend doing that). But in this article, we will use Mathematica which offers a super neat function called NDSolve[] that performs the numerical integration of ODEs. Without further adue, only a couple of lines of code are required in Mathematica to solve the above system of equations

(* Define control parameters *)
\[Sigma] = 3; \[Beta] = 1; \[Rho] = 10;
(* Define initial conditions for later use *)
x0 = 0; y0 = 1; z0 = 1;
(* Define interval of integration *)
Tend = 20 \[Pi];
(* Lump the initial conditions in one variable *)
initialConditions = {x[0] == x0, y[0] == y0, z[0] == z0};
(* Lump the Lorenz equations in one variable *)
LorenzEquations = {x'[t] == \[Sigma] (y[t] - x[t]),
y'[t] == \[Rho] x[t] - x[t] z[t] - y[t],
z'[t] == x[t] y[t] - \[Beta] z[t],
initialConditions};
(* Use NDSolve to integrate the Lorenz equations *)
s1 = NDSolve[LorenzEquations, {x[t], y[t], z[t]}, {t, 0, Tend}, MaxSteps -> \[Infinity]];
(* Plot the solution *)
Plot[Evaluate[x[t] /. s1], {t, 0, Tend}, PlotRange -> All]

Voila!

Note that \[Sigma] will automatically convert to the greek symbol sigma. The same applies for the rest. You can also generate the greek letters by pressing escape, typing a letter on the keyboard, and then pressing escape. For example, escape, s, escape will turn into sigma.

Going back to the previous code, the two important statements are the NDSolve[] and the Plot[Evaluate[]].

In the first one, we are solving the Lorenz equations for x[t], y[t], and z[t] from t = 0 to t = Tend with an infinite number of time steps (MaxSteps->Infinity).

As for the Plot[Evaluate[]], the "x[t] /. s1" means replace all x[t] with the data contained in s1, which holds the results of the numerical integration. One could have also chosen to plot y[t] or z[t].

For first or higher order ODEs, it is advisable to get rid of all derivatives by definig them as new variables. This will be helpful for phase space diagrams to be discussed in the next article. For example, if you have the following system (Ueda's oscillator)

it can be converted to


Download Mathematica notebook [right click / save as]