Plot phase portrait with MATLAB and Simulink
Plot phase portrait with MATLAB and Simulink
Nonlinear Systems
If a system includes one or more nonlinear devices, the system is called a nonlinear system. Every control system is essentially nonlinear. The characteristics of the nonlinear systems can not be described using linear differential equations.
Characteristics of nonlinear systems

The stablity of a nonlinear system depends on its initial conditions, the inputs to the system, and the system parameters

There may exist multiple equilibrium in a nonlinear system, in other words, there may have multiple solutions for $\dot{x} = 0$

The systems may be caught in the state of selfexcited oscillation (limit cycle) compared with linear systems

The outputs of the nonlinear system usually contain additional frequency components and my not contain the original input frequencies.
Analysis and Design Methods of Nonlinear Systems
 Phase plane (which will be discussed in this article)
 Describing function
 MATLAB (implemented in this article)
Phase plane anaylsis^{1}
Phase Plane Analysis is a graphical method for studying first and secondorder systems by

providing motion trajectories corresponding to various initial conditions.

then examine the qualitative features of the trajectories.

finally obtaining information regarding the stability and other motion patterns of the system.
Example 1
Q: Let’s start with a linear system:
$$\ddot{x} + x = 0$$
Method 1: Calculate by hands with phase plane analysis
First, find the eigenvalues of the characteristic equation: $$ \begin{aligned} &\lambda^{2}+1=0\\ &s_{1,2}=\pm i \end{aligned} $$
And we know that with such pole distribution, the phase portrait should look like:
Let’s see what we’ll get using MATLAB Simulink
Method 2: Running Simulink simulation
Steps to follow:
 Find out the input and output
 Describe the system with basic components, i.e. integrators, gain, etc.
 Implement the system
System Topology
In a second order system, we could cascade integrators in order to get $x, \dot{x}$ and $\ddot{x}$
Since the differential equation of the system is in the form of:
$$\ddot{x} + x = 0$$
we could find out the relationship between $\ddot{x}$ and $x$ by moving $x$ to the right of the equation.
$$\ddot{x} = x$$
Phase portrait
Example 2
Q: Find the phase portrait of this secondorder nonlinear system with such differential equation:
$$ \ddot{x}+0.5 \dot{x}+2 x+x^{2}=0 $$
Leaving $\ddot{x}$ on the left side of the equation, we can get the following equation $$ \ddot{x} = f(x, \dot{x}) = 0.5\dot{x}2xx^2 $$
Method 1: Calculate by hands with phase plane analysis
First, find the singularity points of the system, make
$$ \frac{d \dot{x}}{d x}=\frac{0.5 \dot{x}+2 x+x^{2}}{\dot{x}}=\frac{0}{0} $$
then we could find our two singularity points
$$ \left\{\begin{array}{l} x=2 \\ \dot{x}=0 \end{array}\right. $$ and $$ \left\{\begin{array}{l} x=0 \\ \dot{x}=0 \end{array}\right. $$
Therefore, the singularity points of the system lie at $(0, 0)$ and $(2, 0)$ on the phase plane.
Next, let’s linearize the differential equation $f(\dot{x}, x)$ and then solve for the roots for the equation (also serve as system poles) in order to get the phase portrait.
$$ \left.\frac{\partial f(\dot{x}, x)}{\partial x}\right_{x=0 \atop \dot{x}=0}=2 x+2=2 $$
$$ \left.\frac{\partial f(\dot{x}, x)}{\partial \dot{x}}\right_{x=0 \atop \dot{x}=0}=0.5 $$
The linearized system differential equation should look like:
$$ f(\dot{x}, x)=\ddot{x}+0.5 \dot{x}+2x = 0 $$
And the poles are:
$$ s_{1,2}=0.25\pm 1.39i$$
which means $(0,0)$ is a stable focus.
By the same token, we could get the poles for the other singularity point $(2, 0)$
$$ s_{1}=1.19$$ $$s_2 = 0.968i$$
Thus, $(2, 0)$ is a saddle point.
And we know that with such pole distribution, the phase portrait should look like:
Method 2: Running Simulink simulation
Method 3: MATLAB script^{2}

Step1. First create the file
vectfieldn.m
in your working directory% vectfieldn.m % vectfieldn vector field for system of 2 first order ODEs, % with arrows normalized to the same length % vectfield(func,y1val,y2val) plots the vector field for the system of % two first order ODEs given by func, using the grid of y1val and % y2 values given by the vectors y1val and y2val. func is either a % the name of an inline function of two variables, or a string % with the name of an mfile. % By default, t=0 is used in func. A t value can be specified as an % additional argument: vectfield(func,y1val,y2val,t) function vectfield(func,y1val,y2val,t) if nargin==3 t=0; end n1=length(y1val); n2=length(y2val); yp1=zeros(n2,n1); yp2=zeros(n2,n1); for i=1:n1 for j=1:n2 ypv = feval(func,t,[y1val(i);y2val(j)]); yp1(j,i) = ypv(1); yp2(j,i) = ypv(2); end end len=sqrt(yp1.^2+yp2.^2); quiver(y1val,y2val,yp1./len,yp2./len,.6,'r'); axis tight;

Step 2. Create another file named
main.m
and fill in the following code% main.m % init at (1, y20) where y20 are 2:0.3:2.7 could try other initial points f = @(t,y) [y(2);0.5*y(2)2*y(1)y(1)^2] vectfieldn(f,2:.5:8,2.5:.25:2.5) hold on for y20=2:0.3:2.7 [ts,ys] = ode45(f,[0,10],[1;y20]); plot(ys(:,1),ys(:,2)) end hold off

A good illustration on the phase plane analysis could be found at http://ele.aut.ac.ir/~abdollahi/Lec_1_N10.pdf ↩︎

Reference tutorial: http://www.javadtaghia.com/controlsystemmath/nonlinearcontrol/matlabphaseplot ↩︎