Dynamical Systems¶

Many physical systems are explained by an ordinary differential equation (ODE) and it is often needed to solve for a solution of the differential equation. This solution will explain the trajectory behaviour and characteristics of the system. Some types of ODE can be certainly solved analytically such as linear systems. However, a numerical solution can provide an approximate solution to a general equation. This section, we will show built-in commands in MATLAB used for solving differential equations.

An ordinary differential equation (ODE) has the form:

(1) $\dot{x} = f(x,t), \quad x(0) = x_0$

where $x(\cdot) \in \mathbf{R}^n$ , $\dot{x} = \frac{dx}{dt}$ and $f:\mathbf{R}^n \rightarrow \mathbf{R}^n$ . An ODE is an equation that contains only ordinary derivatives of the dependent variable. Examples of ODE are

(2) $\ddot{y} + 2y \dot{y} = y^2$

(3) $\ddot{y}+ 3t^2 \dot{y} = -2\sin y$

(4) $\dot{y} + 2ty = 0$

(5) $\ddot{y} + 4 y = \cos t$

(6) $\frac{d^n y}{d t^n} + a_1 \frac{d^{n-1} y}{d t^{n-1}} + a_2 \frac{d^{n-2} y}{d t^{n-2}} + \cdots + a_n y = 0$

If an equation contains partial derivatives with respect to two or more independent variables, then it is called a partial differential equation (PDE). For example, a heat equation

$\frac{\partial u(x,t)}{\partial t} = a \frac{\partial^2 u(x,t)}{\partial^2 x^2}$

where $u(x,t)$ is the temperature at a given location $x$ , and at time $t$ . Finding a numerical solution of a PDE is an advanced topic. In this document, we will cover how to solve ODE only.

Short review on dynamical systems¶

The order of a differential equation is the order of its highest derivative. The systems (2), (3), and (5) are second-order, system (4) is first-order, and system (6) is $n$ -order.

Formulating to a first-order system¶

Any $n$ -order ODEs can be arranged to the form (1) with $n$ components in variable $x$ , which is referred to as state variable. Consider the system (2). If we define $x_1 = y, x_2 = \dot{y}$ , then

(7) $\ddot{y} + 2y \dot{y} = y^2 \quad \Longleftrightarrow \quad \left [ \begin{array}{l} \dot{x}_1 \\ \dot{x}_2 \end{array} \right ] = \left [ \begin{array}{l} x_2 \\ -2x_1 x_2 + x_1^2 \end{array} \right ].$

Similarly for system (3), we define $x_1 = y, x_2 = \dot{y}$ . Then

(8) $\ddot{y}+ 3t^2 \dot{y} = -2\sin y \quad \Longleftrightarrow \quad \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = \begin{bmatrix} x_2 \\ -3t^2 x_2 - 2 \sin x_1 \end{bmatrix}.$

The system (4) is already a first-order system with $f(x,t) = -2ty$ . For system (5), we define $x_1 = y, x_2 = \dot{y}$ . Then

(9) $\ddot{y} + 4 y = \cos t \quad \Longleftrightarrow \quad \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = \begin{bmatrix} x_2 \\ -4 x_1 + \cos t \end{bmatrix}.$

Last example is the system (6). Define

$x_1 = y ,\quad x_2 = \frac{dy}{dt} ,\quad \cdots ,\quad x_n = \frac{d^{n-1} y }{dt^{n-1}}$

then, we can transform (?) into

(10) $\frac{d^n y}{d t^n} + a_1 \frac{d^{n-1} y}{d t^{n-1}} + a_2 \frac{d^{n-2} y}{d t^{n-2}} + \cdots + a_n y = 0 \Longleftrightarrow \frac{d}{dt} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_{n-1} \\ x_n \end{bmatrix} = \begin{bmatrix} x_2 \\ x_3 \\ \vdots \\ x_n \\ -a_1 x_n - a_2 x_{n-1} - \cdots - a_n x_1 \end{bmatrix}.$

Linear systems¶

The system of the form (1) is said to be linear or the ODE (1) is a linear equation if $f(x,t)$ is a linear function in $x$ . A linear ODE of order $n$ can also be written as

(11) $\frac{d^n y}{d t^n} + a_1(t) \frac{d^{n-1} y}{d t^{n-1}} + a_2(t) \frac{d^{n-2} y}{d t^{n-2}} + \cdots + a_n(t) y = g(t).$

For example, the systems (4), (9), and (10) are linear systems, while the systems (7) and (8) are called nonlinear systems.

Homogenous linear equations¶

A linear ODE of the form (?) is called homogeneous linear equation if $g(t)=0$ (identically zero). For example, the linear systems (4) and (10) are homogenous equations, while the linear system (5) is nonhomogenous.

Autonomous systems¶

The system of the form (1) is said to be autonomous or time-invariant if $f(x,t)$ is not explicit dependent on time $t$ . For example, the systems (7) and (10) are autonomous systems, while the systems (8), (4) and (9) are not.

Linear systems and examples on chemical processes¶

Any linear autonomous system can be expressed as

(12) $\dot{x} = Ax$

where $x \in \mathbf{R}^n$ and $A \in \mathbf{R}^{n \times n}$ .

Example: $n$ -order linear ODE¶

The ODE described in (10) can be arranged into (12) with

$\frac{d}{dt} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_{n-1} \\ x_n \end{bmatrix} = \begin{bmatrix} x_2 \\ x_3 \\ \vdots \\ x_n \\ -a_1 x_n - a_2 x_{n-1} - \cdots - a_n x_1 \end{bmatrix} \quad \Longrightarrow \quad A = \begin{bmatrix} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & & & \\ 0 & 0 & 0 & & 1 \\ -a_n & -a_{n-1} & -a_{n-2} & \ldots & -a_1 \end{bmatrix}$

Example: Liquid level systems¶

Suppose $Q$ is a steady state liquid flow to a tank and $q_i$ is a deviation of inflow rate. Define $h$ is the deviation of the liquid level in the tank. The differential equation of this system is

$\frac{dh }{dt} + \frac{1}{RC} h = \frac{1}{C} q_i$

where $R,C$ are the resistance and capacitance of the tank (parameters of the system). When the inflow rate $q_i$ is zero, this system is autonomous and can be described by the form (12) with $A = -1/RC$ .

Example: Electrically heated stirred tank¶

Consider a stirred-bank heating system with a mass inflow rate $w_i$ and an inlet temperature $T_i$ . Let $T$ be the temperature of the tank contents and $T_e$ be the temperature of the heating element. The dynamics of the these temperatures are

(13) $\begin{eqnarray*} mC \dot{T} & = & wC (T_i - T) + H_e (T_e - T) \\ m_e C_e \dot{T_e} & =& Q - H_e (T_e-T) \end{eqnarray*}$

where $mC, m_e C_e$ are the thermal capacitances of the tank contents and the heating element, respectively. $Q$ is an input variable, $h_eA_e$ is the heat transfer coefficient.

Suppose the input variables $w,Q,T_i$ are specified to be zero. Then the system (13) can be formulated to (12) as

$\begin{bmatrix} \dot{T} \\ \dot{T_e} \end{bmatrix} = \begin{bmatrix} -H_e/mC & H_e/mC \\ H_e/m_eC_e & -H_e/m_e C_e \end{bmatrix} \begin{bmatrix} T \\ T_e \end{bmatrix}$

Solution to an autonomous linear system¶

Given an initial condition $x(0)$ . We take a Laplace transform of (12) (elementwise) and we obtain

$sX(s) - X(0) = A X(s) \quad \Longrightarrow \quad (sI-A) X(s) = X(0) \quad \Longrightarrow \quad X(s) = (sI-A)^{-1} X(0)$

We define the matrix exponential of $A$ as

$e^{tA} = \mathcal{L}^{-1} (sI-A)^{-1}$

Then the solution to (12) is

$x(t) = e^{tA} x(0)$

Comments

The matrix exponential $e^{tA}$ is an $n \times n$ matrix and is a function of $t$ .

It is a generalization of $e^{at}$ to an $n$ -dimensional system.

$e^{tA}$ is a power series of the matrix $A$

$e^{tA} = I+At + \frac{A^2 t^2}{2!} + \frac{A^3 t^3}{3!} + \ldots$

Example¶

Suppose a linear system is given by

(14) $\dot{x} = \begin{bmatrix} -1 & 1 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 2 \end{bmatrix}$

with initial condition $x(0) = \begin{bmatrix} 1 & -1 & 2 \end{bmatrix}^T$ . We simulate a time response of $x(t)$ from $t = 0$ to $t=10$ by:

>> A=[-1 1 0;0 -1 0;0 0 -2];
>> t=linspace(0,10,300);
>> x=zeros(3,300);
>> matA = zeros(3,3,300);
>> x0 = [1 -1 2]';
>> for k=1:300, matA(:,:,k) = expm(A*t(k)); x(:,k) = matA(:,:,k)*x0;end;
>> plot(t,x)

We simulate 300 time points, so we must have 300 matrices of $e^{At}$ (which is stored in matA). We create a 3-D array matA and for each matA(:,:,k), it is a 3 x 3 matrix of $e^{At_k}$ .

ODE solvers for nonlinear differential equations¶

This section provides MATLAB commands used to solve an ODE with the form

$\dot{x} = f(x,t)$

for $t_0 \leq t \leq t_f$ with a given initial condition $x(t_0)$ .

An example of a problem to solve is $\dot{x} = -x^3$ with $x(0)= 1$ . We will explain steps to solve ODE as follows.

1. Create a user-defined function¶

The description of $f(x,t)$ must be written as a user-defined function (in a separate M-file) or as an anonymous function:

function [dx] = ode1(t,x)

dx = -x^3;

If an anonymous function is used, we can define it in the command window as:

>> ode2 = @(t,x) x^3

Note that though $f(x,t)$ is not an explicit function of $t$ , the ODE function must be defined as a function of two variables; time and state respectively.

2. Select an ODE solvers¶

There are many numerical methods to solve ODEs. Typically, a numerical method starts by dividing the time interval into small time pieces. The solution starts at the known point $x(t_0)$ and the value of $x$ at each time step is calculated using an integration method according to the selected solver. A short list of ODE solvers is given by:

ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb

The solvers that contain s or t in its name are used to solve stiff problems. We will discuss stiff and nonstiff problems later. For now, we can start by ode45 solver, which can be a good try for most problems.

3. Solve the ODE¶

An ODE solver command has the form:

[T,X] = solver(odefun,tspan,x0,options)

Description

solver is the name of the solver such as ode45, ode23.
odefun is the function descrption of $f(x,t)$ obtained from step 1.
tspan is a vector that specifies the interval of the solution
x0 is the initial value of $x$ (the first value of $x$ at the first point of the interval.
[T,X] is the solution of the ODE where T is the time variable.
options is the structure of optional parameters that change the default integration properties. The options can be set using the odeset function.

For example, we simulate the system $\dot{x} = -x^3$ with $x(0)= 1$ from t = 0 to t=10 by:

>> [t,x]=ode45(@ode1,[1 10],1);
>> plot(t,x);

If we refer to the ODE description written as anonymous function, then we type:

>> [t,x]=ode45(ode2,[1 10],1);
>> plot(t,x);

(without the @ sign.)

We can also solve the ODE in this example analytically. The closed-form solution to $\dot{x} = -x^3$ is given by

$x(t) = \sqrt{\frac{0.5}{0.5/x(t_0)^2+t - t_0}}$

Hence, we can verify if the two methods produce the same results:

>> x0=1
>> y=sqrt(0.5./(0.5/x0^2+t-t(1)));
>> plot(t,y,t,x)

From these steps, a response from the system (14) can be certainly solved from ode45 as:

>> A=[-1 1 0;0 -1 0;0 0 -2];
>> lin1 = @(t,x) A*x
>> [t,x]=ode45(lin1,[0 10],[1 -1 2]');
>> plot(t,x)

The result must be the same as using the expm command obtained previously.

For the rest of this section, we will provide more examples.

Time-varying systems¶

Given a time-varying system satisfying an ODE

$t^3\dddot{y} -3t^2 \ddot{y} + 6t \dot{y} - 6y = t^4 \log t ,\quad t >0$

The system can be expressed in a state-space form as

$\begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \\ \dot{x}_3 \end{bmatrix} = \begin{bmatrix} x_2 \\ x_3 \\ \frac{3}{t} x_3 - \frac{6}{t^2} x_2 + \frac{6}{t^3} x_1 + t^4 \log t \end{bmatrix}$

We write the system description in ode3.m as

function dx = ode3(t,x)

dx = zeros(3,1);
dx(1) = x(2);
dx(2) = x(3);
dx(3) = (3/t)*x(3)-(6/t^2)*x(2)+(6/t^3)*x(1) + t^4*log(t);

In this example, the number of state variables is 3, so ode3 must return a column vector of size 3.

We simulate the system from t = 1 to t = 5, with $x(1) = \begin{bmatrix} 1 & 0 & -1 \end{bmatrix}^T$ by

>> [t,x] = ode45(@ode3,[1 5],[1 0 -1]');
>> plot(t,x);

Stiff differential equations¶

Stiffness is ODEs refers to a wide range in the time scales of the components in the vector solution. For example, the following equation may be considered stiff

$\ddot{y} + 1001 \dot{y} + 1000 y = 0$

The closed-form solution to this system is

$y(t) = C_1 e^{-t} + C_2 e^{-1000 t}$

So the response due to $e^{-1000 t}$ decays very quickly compared to $e^{-t}$ . A numerical solver would need a small step size to compute the rapid changes due to the fast mode. Hence, some numerical procedures that are quite satisfactory in general will perform poorly on stiff equations.

MATLAB also provides stiff solvers such as ode15s, ode23s, ode23t, ode23tb. Consult MATLAB help for detail of each solver.

As an example of a stiff system, we consider the van der Pol equations in relaxation oscillation

$\ddot{y} - \mu (1-y^2) \dot{y} + y = 0$

for $\mu = 1000$ . To simulate this system, create a function vdp containing the equations:

function dy = vdp(t,y)
global mu

dy = zeros(2,1);    % a column vector
dy(1) = y(2);
dy(2) = 1000*(1 - y(1)^2)*y(2) - y(1);

We solve on a time interval of [0 3000] with initial condition vector [2 0] at time 0:

>> global mu
>> mu = 1000;
>> [t,x] = ode15s(@vdp,[0 3000],[2 0]);
>> plot(t,x(:,1),'-o')

We have solved using a stiff solver ode15s. Then we plot the first component of vector x which is the variable y. If we plot the second component of x, we will see a rapid change in the response:

>> plot(t,x(:,2));

In this example, if we were try to solve using a nonstiff solver such as ode45, the solver would take a long time to complete the task.

Phase portrait¶

A second-order time-invariant system is represented by

$\begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = \begin{bmatrix} f_1(x_1,x_2) \\ f_2(x_1,x_2) \end{bmatrix}$

We are interested in a plot of solution $x_1(t)$ versus $x_2(t)$ for all $t \geq 0$ . This curve is called a trajectory or orbit of the system. The family of all trajectories or solution curves (which started by different initial points) is called phase portrait. In control systems, it is used to illustrate qualitative behaviour of a 2-dimensional nonlinear system.

Consider a pendulum system described by

$\ddot{\theta} + c \dot{\theta} + k \sin \theta = 0$

where $c,k$ are system parameters. Both are positive numbers. So we write the ODE in pendulum file as:

function dx = pendulum(t,x)
global c k

dx = zeros(2,1);
dx(1) = x(2);
dx(2) = -c*x(2)-k*sin(x(1));

Suppose $c = 0, k=1$ and we simulate the response from t = 0 to t=20 as

>> global c k
>> c = 0; k = 1;
>> [t,x] = ode45(@pendulum,[0 20],[1 0]);
>> plot(x(:,1),x(:,2))

From the plot, the trajectory starts from the initial condition x=(1,0) and then converges to zero. We can make a several curves using different initial conditions by performing a loop:

c = 0; k = 1;
intval = linspace(-1,1,10);

figure;hold on;
for ii=1:10,
for jj=1:10,
        [t,x] = ode45(@pendulum,[0 20],[intval(ii) intval(jj)]);
        plot(x(:,1),x(:,2));
end
end

This produces a phase portrait as shown in the figure (left). The trajectories are in a circle shape. Now let $c = 1, k =1$ (make c nonzero) and repeat the same process, we have the following plot (right). If $c > 0$ , it is known that the trajectories eventually converge to the origin (which is the equilibrium point of this system.)

Exercises¶

Solve the ODE

$\begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = \begin{bmatrix} x_1 + x_2 -x_1(x_1^2+x_2^2) \\ -2x_1 + x_2 - x_2(x_1^2+x_2^2) \end{bmatrix}$

using the initial point $x = [-2 \;\; 3]^T$ .

Solve the ODE from t = 0 to t = 10

$\ddot{y}(t) + \frac{1}{9}y(t) = 0$

with $y(0) = 1$ .

Solve the ODE from t = 0 to t = 10

$\dddot{y} + 2\ddot{y} - \dot{y} - 2y = 1-4 t^3$

with $y(0) = 1, \dot{y}(0) = 0, \ddot{y}(0) =0$ .

Solve the Duffing equation

$\ddot{y} + \omega_0^2 y + \beta y^3 = 0$

where usually $|\beta |$ is small. Plot a phase portrait for both $\beta > 0$ and $\beta < 0$ .

Plot phase portraits of the Van der Pol equation using $\mu = 0.2, 1.0, 5$ . Use different initial points in the interval [-3,3].
Plot phase portraits of the Hopt bifurcation equation

$\begin{eqnarray*} \dot{x}_1 &=& x_1 [ \mu + (x_1^2 + x_2^2)-(x_1^2+x_2^2)^2] - x_2 \\ \dot{x}_2 &=& x_2 [ \mu + (x_1^2 + x_2^2)- (x_1^2+x_2^2)] +x_1 \end{eqnarray*}$

for $\mu > 0$ and $\mu < 0$ .

Compare phase portraits of the three equations

$\begin{eqnarray*} \ddot{y}+2\dot{y} + 5 y &=& 0 \\ \ddot{y}+3\dot{y} + 2y & =& 0 \\ \ddot{y}-3\dot{y}+2 y &=& 0 \end{eqnarray*}$

How are the plots related to the poles of the auxillary equations ?

Navigation

Table Of Contents