Using Matlab for Autonomous Systems

clearvars; clf

You need to download the file vectorfield.m and put it in the directory where your mlx files are.

Example: Consider the following autonomous system:

Define the right hand side function f

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

Note that you must use @(t,y) ..., despite the fact that the functions do not depend on t.

Plot the vector field

for

between -1 and 3 and

between -2 and 2 using steps of size .5:

vectorfield(f,-1:.5:3,-2:.5:2); xlabel('y_1'); ylabel('y_2')

Plot trajectories for many initial points

We can solve the ODE with initial condition

for t going from 0 to 4 and plot the trajectory using

[ts,ys] = ode45(f,[0,4],[a1;a2]);

plot(ys(:,1),ys(:,2))

We can solve this initial value problem t going from 0 to 4 and plot the trajectory using

[ts,ys] = ode45(f,[0,-4],[a1;a2]);

plot(ys(:,1),ys(:,2))

The ode45 commands may generate a number of warning messages since the solutions does not exist for t in the whole interval [-4,4]). We can turn off these warnings with the command

warning('off','MATLAB:ode45:IntegrationTolNotMet')

Now we plot the trajectories going through the points

with

and

for a1=-1:.5:3

for a2=-2:.5:2

[ts,ys] = ode45(f,[0,4],[a1;a2]);

plot(ys(:,1),ys(:,2),'c'); hold on

[ts,ys] = ode45(f,[0,-4],[a1;a2]);

plot(ys(:,1),ys(:,2),'c')

end

vectorfield(f,-1:.5:3,-2:.5:2); % plot vector field

hold off; axis([-1 3 -2 2]); xlabel('y_1'); ylabel('y_2')

Investigation of stationary points

Define symbolic variables

syms t y1 y2

Define symbolic expression F by using f with the vector [y1;y2]:

F = f(t,[y1;y2])

F =

Find the stationary points

Find

where

using solve:

[y1s,y2s] = solve(F,[y1;y2]) % find all solutions (y1,y2) of F=[0;0]

y1s =

y2s =

m = length(y1s) % number of solutions found

m =      2

We have found two stationary points: (0,1) and (2,-1)

Find the Jacobian matrix:

J = jacobian(F,[y1;y2])

J =

Investigate the stationary points

Result:

stationary point (0,1): eigenvalues 1,2, nodal source
stationary point (2,-1): eigenvalues -2,1, saddle

for k=1:m
  stationary_point = [y1s(k),y2s(k)]

Evaluate the Jacobian at the stationary point:

A = subs(J,{y1,y2},{y1s(k),y2s(k)})

[eigvect,eigvals] = eig(A)

figure

phaseportrait(A);

title(sprintf('linearized problem for stationary point (%g,%g)',stationary_point))

end

stationary_point =

A =

eigvect =

eigvals =

stationary_point =

A =

eigvect =

eigvals =

Function for drawing phase portrait

function phaseportrait(A)
syms y(t) [2 1]
syms C1 C2
sol = dsolve(diff(y)==A*y);             % solve ODE 
ysol = subs(y(t),sol);
f = @(t,y) A*y;                         % define function f(t,y) for vectorfield
vectorfield(f,-5:5,-5:5); hold on
for c1=[-1 0 1]
  for c2=[-1 0 1]
    ys = subs(ysol,{C1,C2},{c1,c2});    % use C1=c1, C2=c2
    fplot(ys(1),ys(2),[-5 5]);          % plot trajectory for t=-5...5
  end
end
hold off; axis equal; axis([-5 5 -5 5]); xlabel('y_1'); ylabel('y_2')
end