Figure 2.18:

A limit cycle (thick dashed curve) and a trajectory (thin solid curve) approaching it.

Code for Figure 2.18

Text of the GNU GPL.

main.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
npts = 101;
time = linspace(0, 10, npts);
z0 = [0; 0.05];


[time, z] = ode15s (@rhs, time, z0);

% start at long time solution and compute the limit cycle; period is
% about 6.5; exact value doesn't matter

z0 = z(end,:)';
time = linspace(0, 6.5, npts);

[time, zcyc] = ode15s (@rhs, time, z0);

% save trajectory and limit cycle and annulus
theta = linspace(0, 2*pi, npts)';
anni = [cos(theta), sin(theta)];
anno = 2/sqrt(5)*anni;
traj = [z, zcyc, anni, anno];

figure()
plot(time,z)
figure()
plot(z(:,1),z(:,2), '-', zcyc(:,1),zcyc(:,2), 'o', anni(:,1), ...
     anni(:,2),  anno(:,1), anno(:,2))

save "cycle.dat" traj

rhs.m

1
2
3
4
5
function zdot = rhs(t, z)
  x = z(1);
  y = z(2);
  zdot = [x - y - x*(x^2+2*y^2);
          x + y - y*(x^2+y^2)];