cBook‎ > ‎Non-linear Control‎ > ‎

Matlab Phase plot

posted Sep 13, 2015, 7:38 PM by Javad Taghia

Using Matlab for Higher Order ODEs and Systems of ODEs

(Continuation of Using Matlab for First Order ODEs)

Contents

Numerical Solution
Converting problems to first order systems
Plotting the solution
Finding numerical values at given t values
Making phase plane plots
Vector fields for autonomous problems
Plotting the vector field
Plotting the vector field together with solution curves
Symbolic Solution
Example with 2nd order problem
Plotting the solution
Finding numerical values at given t values
Example with first order system
Plotting the solution
Finding numerical values at given t values
Making phase plane plots

Numerical solution

Example problem: The angle y of an undamped pendulum with a driving force sin(5 t) satisfies the differential equation

y'' = -sin(y) + sin(5 t)

and the initial conditions

y(0) = 1
y'(0) = 0.

If your problem is of order 2 or higher: rewrite your problem as a first order system. Let y1=y and y2=y', this gives the first order system

y1' = y2,
y2' = -sin(y1) + sin(5 t)

with the initial conditions

y1(0) = 1
y2(0) = 0.

Define an @-function f for the right hand side of the first order system: Note that f must be specified as a column vector using [...;...] (not a row vector using [...,...]). In the definition of f, use y(1) for y1, use y(2) for y2. The definition off should have the form

f = @(t,y) [expression for y1' ; expression for y2' ];

You have to use f = @(t,y)... even if t does not occur in your formula. For our example the first component of f is y2, the second component is -sin(y1) + sin(5 t) and we define

f = @(t,y) [y(2); -sin(y(1))+sin(5*t)]

To plot the numerical solution: To obtain plots of all the components of the solution for t going from t0 to t1 use ode45(f,[t0,t1],[y10;y20]) where y10y20 are the initial values for y1y2 at the starting point t0. In our example we get a plot for t going from 0 to 20 with the initial values [1;0] using

ode45(f,[0,20],[1;0])

This shows the two functions y1(t)=y(t) (blue) and y2(t)=y'(t) (green).

The circles mark the values which were actually computed (the points are chosen by Matlab to optimize accuracy and efficiency). You can obtain a vector ts and a matrix ys with the coordinates of these points using [ts,ys] = ode45(f,[t0,t1],[y10;y20]). You can then plot the solution curves using plot(ts,ys) (this is a way to obtain a plot without the circles). Note that each row of the matrix ys contains 2 entries corresponding to the two components of the solution at that time:

[ts,ys] = ode45(f,[0,20],[1;0]); % find ts, ys, but don't showplot(ts,ys) % make plot of y1 and y2 vs. t[ts,ys] % show table with 3 columns for t, y1, y2

You can obtain the vector of y1 values in the first column of ys by using ys(:,1), therefore plot(ts,ys(:,1)) plots only y1(t).

To obtain numerical values at specific t values: You can specify a vector tv of t values and use [ts,ys] = ode45(f,tv,[y10;y20]). The first element of the vector tv is the initial t value; the vector tv must have at least 3 elements. E.g., to obtain the solution with the initial values [1;0] at t = 0, 0.5, ..., 10 and display the results as a table with 3 columns ty1y2, use

[ts,ys]=ode45(f,0:0.5:10,[1;0]);[ts,ys]

To plot trajectories in the phase plane: To see the points with coordinates ( y1(t), y2(t) ) in the y1y2 plane for t going from 0 to 20 type

options=odeset('OutputFcn','odephas2');ode45(f,[0,20],[1;0],options)

This shows the points while they are being computed (the plotting can be stopped with the stop button). To first compute the numerical values and then plot the curve without the circles, type

[ts,ys] = ode45(f,[0,20],[1;0]); % find ts, ys, but don't showplot(ys(:,1),ys(:,2)) % make plot of y2 vs. y1

Vector fields for autonomous systems of two first order ODEs

If the right hand side function f(ty) does not depend on t, the problem is called autonomous. In this case the behavior of the differential equation can be visualized by plotting the vector f(ty) at each point y = (y1,y2) in the y1,y2 plane (the so-called phase plane).

First save the files vectfield.m and vectfieldn.m into the same directory where your m-files are.

To plot the vector field for y1 going from a1 to b1 with a spacing of d1 and y2 going from a2 to b2 with a spacing of d2 use vectfield(f,a1:d1:b1,a2:d2:b2). The command vectfieldn works in the same way, but produces arrows which all have the same length. This makes it easier to see the direction of the vector field.

Example: The undamped pendulum problem without a driving force has the differential equation y'' = -sin(y). This gives the first order system

y1' = y2,
y2' = -sin(y1)

Here we define

f = @(t,y) [y(2);-sin(y(1))]

We can plot the vector field and 10 trajectories with starting points (0,0), (0,0.3), ..., (0,2.7) in the phase plane as follows:

vectfield(f,-2:.5:8,-2.5:.25:2.5)
hold on
for y20=0:0.3:2.7
[ts,ys] = ode45(f,[0,10],[0;y20]);
plot(ys(:,1),ys(:,2))
end
hold off

Symbolic solution

Use the dsolve command. Specify all differential equations as strings, using Dy for y'(t), D2y for y''(t) etc. .

For an initial value problem specify the initial conditions in the form 'y(t0)=y0''Dy(t0)=y1' etc. . The last argument of dsolve is the name of the independent variable, e.g., 't'.

Example: For the differential equation

y'' = -y + sin(5 t)

use

sol = dsolve('D2y = -y + sin(5*t)','t')

In this case, the answer can be simplified by typing

s = simple(sol)

which gives the general solution -1/24*sin(5*t)+C1*cos(t)+C2*sin(t) with two constants C1C2.

To solve the ODE with initial conditions y(0) = 1, y'(0) = 0 use

sol = dsolve('D2y = -y + sin(5*t)','y(0)=1','Dy(0)=0','t')

Then s = simple(sol) gives the solution -1/24*sin(5*t)+5/24*sin(t)+cos(t).

To plot the solution curve use ezplot:

ezplot(sol, [0,20])

Example for system of ODEs: For the system

y1' = y2,
y2' = -y1 + sin(5 t)

with the initial conditions

y1(0) = 1
y2(0) = 0

type

sol = dsolve('Dy1=y2','Dy2=-y1+sin(5*t)','y1(0)=1','y2(0)=0','t')

which gives the somewhat strange response

sol =y1: [1x1 sym] y2: [1x1 sym]

This means that the two components of the solution can be accessed as sol.y1 and sol.y2 : Typing

sol.y1

gives -2/3*sin(t)*cos(t)^4+1/2*sin(t)*cos(t)^2+1/6*sin(t)+cos(t). This can be simplified by typing

s1 = simple(sol.y1)

which gives -1/24*sin(5*t)+5/24*sin(t)+cos(t). For the second component y2 of the solution we proceed in the same way: Typing

s2 = simple(sol.y2)

gives -sin(t)-5/24*cos(5*t)+5/24*cos(t).

To plot the solution curves use ezplot:

ezplot(sol.y1, [0,20])hold onezplot(sol.y2, [0,20])hold off

To obtain numerical values use double(subs(sol.y1,'t',tval)) where tval is a number or a vector of numbers. E.g., the following commands generate a table with three columns ty1y2 for t=0, 0.5, ..., 10:

tval = (0:0.5:10)'; % column vector with t-valuesyval = double(subs([sol.y1,sol.y2],'t',tval)); % 2 columns with y1,y2 [tval, yval] % display 3 columns together

To plot the solution in the phase plane use a similar approach with more t values:

tval = (0:0.1:10)'; % column vector with t-valuesyval = double(subs([sol.y1,sol.y2],'t',tval)); % 2 columns with y1,y2 plot(yval(:,1),yval(:,2)) % plot col.2 of yval vs. col.1 of yval

clear all; close all; clc


Problem setup

An undamped pendulum with no driving force is described by

We reduce this to standard matlab form of a system of first order ODEs by letting and . This leads to:

The phase portrait is a plot of a vector field which qualitatively shows how the solutions to these equations will go from a given starting point. here is our definition of the differential equations:

f = @(t,Y) [Y(2); -sin(Y(1))];


To generate the phase portrait, we need to compute the derivatives  and  at  on a grid over the range of values for  and  we are interested in. We will plot the derivatives as a vector at each (y1, y2) which will show us the initial direction from each point. We will examine the solutions over the range -2 < y1 < 8, and -2 < y2 < 2 for y2, and create a grid of 20x20 points.

y1 = linspace(-2,8,20);
y2 = linspace(-2,2,20);

% creates two matrices one for all the x-values on the grid, and one for
% all the y-values on the grid. Note that x and y are matrices of the same
% size and shape, in this case 20 rows and 20 columns
[x,y] = meshgrid(y1,y2);
size(x)
size(y)

ans =

20    20

ans =

20    20



computing the vector field

we need to generate two matrices, u and v. u will contain the value of y1' at each x and y position, while v will contain the value of y2' at each x and y position of our grid.

we preallocate the arrays so they have the right size and shape

u = zeros(size(x));
v = zeros(size(x));

% we can use a single loop over each element to compute the derivatives at
% each point (y1, y2)
t=0; % we want the derivatives at each point at t=0, i.e. the starting time
for i = 1:numel(x)
Yprime = f(t,[x(i); y(i)]);
u(i) = Yprime(1);
v(i) = Yprime(2);
end


now we use the quiver command to plot our vector field

quiver(x,y,u,v,'r'); figure(gcf)
xlabel('y_1')
ylabel('y_2')
axis tight equal;


Plotting solutions on the vector field

Let's plot a few solutions on the vector field. We will consider the solutions where y1(0)=0, and values of y2(0) = [0 0.5 1 1.5 2 2.5], in otherwords we start the pendulum at an angle of zero, with some angular velocity.

hold on
for y20 = [0 0.5 1 1.5 2 2.5]
[ts,ys] = ode45(f,[0,50],[0;y20]);
plot(ys(:,1),ys(:,2))
plot(ys(1,1),ys(1,2),'bo') % starting point
plot(ys(end,1),ys(end,2),'ks') % ending point
end
hold off


what do these figures mean? For starting points near the origin, and small velocities, the pendulum goes into a stable limit cycle. For others, the trajectory appears to fly off into y1 space. Recall that y1 is an angle that has values from  to . The y1 data in this case is not wrapped around to be in this range.

'done'

% categories: ODEs
% tags: math
% post_id = 799; %delete this line to force new post;

ans =

done