First download the file vectfield.m .
Here we consider the following example of an autonomous system:
y1' = 4y1 + y2 + y1 y2
y2' = y1 + 4 y2 + y22
First define the right hand side function
f of the differential equation as
f = @(t,y) [4*y(1)+y(2)+y(1)*y(2) ; y(1)+4*y(2)+y(2)^2]
Note that you must use
@(t,y) ..., despite the fact that the
functions do not depend on t.
To plot the vector field for y1 between -6 and 5 and y2 between -6 and 1 type
We now want to plot the trajectory through a point (a1,a2), say (a1,a2)=(1,-3):
vectfield(f,-6:.5:5,-6:.5:1)
a1=1; a2=-3
To plot a trajectory in the phase
plane starting at a point (a1, a2) at
time t=0 for increasing values of t going from 0 to
4 type
[ts,ys] = ode45(f,[0,4],[a1;a2]); plot(ys(:,1),ys(:,2))
To plot the trajectory in the phase plane starting at the point
(a1, a2) at time t=0 for
decreasing values of t going from 0 to -4 type
[ts,ys] = ode45(f,[0,-4],[a1;a2]); plot(ys(:,1),ys(:,2))
To get an idea of the behavior of the ODE, we plot the direction
field and several trajectories together: We can choose e.g. the 63
starting points (a1, a2) with
a1=-7,-5,-3,-1,1,3,5 and a2=-7,-6,-5,-4,-3,-2,-1,0,1. From each
of those points we plot the trajectory for t going from 0 to 4 and
for t going from 0 to -4:
vectfield(f,-6:.5:5,-6:.5:1)
hold on
for a1=-7:2:5
for a2=-7:1
[ts,ys] = ode45(f,[0,4],[a1;a2]);
plot(ys(:,1),ys(:,2))
[ts,ys] = ode45(f,[0,-4],[a1;a2]);
plot(ys(:,1),ys(:,2))
end
end
hold off
These commands will generate a number of warning messages (since some solutions don't exist for t in the whole interval [-4,4]) which you can ignore. You can use the command
warning('off','MATLAB:ode45:IntegrationTolNotMet')
at the beginning of your m-file to suppress these warnings. Then
your code should run without any warnings or errors.
From the picture it appears that there are three critical points:
Recall that we defined the right hand side function
f of the differential equation as
f = @(t,y) [4*y(1)+y(2)+y(1)*y(2) ; y(1)+4*y(2)+y(2)^2]
We define symbolic variables
syms t y1 y2
and define F by using f with the vector
[y1;y2]
F = f(t,[y1;y2])
which gives the symbolic vector
[ 4*y1+y2+y1*y2] [ y1+4*y2+y2^2]
To find critical points where both components
F(1) and F(2) of F are zero we use
solve:
[y1s,y2s] = solve(F(1),F(2),y1,y2)
This returns
y1s = [ 0]
[ -5]
[ 3]
y2s = [ 0]
[ -5]
[ -3]
Therefore we have found three critical points: (0,0), (-5,-5) and (3,-3). Your version of Matlab may give the answer in a different order.
To find the Jacobian matrix we type
jac = jacobian(F,[y1;y2])
This returns the matrix
[ 4+y2, 1+y1] [ 1, 4+2*y2]
As an example we will study the second critical point with coordinates
y1s(2)=-5, y2s(2)=-5.
To evaluate the Jacobian at the critical point we
substitute y1s(2) and y2s(2) for y1 and
y2:
A = subs(jac,{y1,y2},{y1s(2),y2s(2)})
which gives the matrix
[ -1, -4] [ 1, -6]
To find the eigenvectors and eigenvalues of this matrix type
[eigvect,eigval] = eig(A)
which gives
eigvect = [ 1, 4]
[ 1, 1]
eigval = [ -5, 0]
[ 0, -2]
This means that the eigenvalues are -5 and -2. Using the table for critical points we see that this critical
point is a stable node. The eigenvector for the first eigenvalue -5
is given by the first column of the matrix eigvect, it has
components 1, 1. The eigenvector for the second eigenvalue -5 is given by the
second column of the matrix eigvect, it has components 4, 1.
Using the eigenvectors we can sketch by hand the behavior of the trajectories
close to the critical point:
Compare this picture with the trajectories near (-5,-5) in the picture above.
You should try to analyze the other two critical two points by yourself in the same way.