Figure 2.15:

Contours of an energy function V(x_{1},x_{2}) or H(x_{1},x_{2}).

Code for Figure 2.15

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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
% contour plot of V(x,y)= 10y^2 - kappa*cos(pi/2*x) + 0.2*(x-2)^2
%
% jbr, 9/3/2011
%


n = 201;
x = linspace(-5, 5, n);
y = linspace(-5, 5, n);
kappa = 2;

xx = kron(x, ones(n,1));
yy = kron(y', ones(1,n));

xx2 = (xx-2).*(xx-2);
yy2 = yy.*yy;

z = 10*yy2 - kappa*cos(xx*pi/2) + 0.2*xx2;

val = [-1,  0,  1, 2, 3, 5, 8, 10];
[c, lev] = contour(x,y,z,val);
axis([-1.75, 3, -1, 1])
xlabel ("$q$");
% avoid rotating y-axis label
text (-5.75, 0, "$p$")

% write contour data for gnuplot
col = 1;
myfile = fopen('energy.dat', 'w');
while col <= size(c, 2)
  coln = col + c(2,col);
  data = c(:,col+1:coln);
% print this contour data and move to next contour line
  fprintf(myfile, '%e %e \n', data)
  fprintf(myfile, '\n\n')
  col = coln + 1;
end
fclose(myfile);

% now solve the gradient system to locate the arrows on the contour
% plot; start on the V=3 contour
w0 = [1.60000; 0.36666; 3];
npts = 101;
time = linspace(0, 0.5, npts);


[time, w] = ode15s (@rhs, time, w0);

hold on
plot(w(:,1),w(:,2),'-o')

% choose rows of w matrix corresponding to energy contours in val(1:5)
for i = 1:5
  [Vval, Vrows(i)]  = min (abs(w(:,3)-val(i)));
end
Vrows = fliplr(Vrows);
diff = w(Vrows+1,1:2) - w(Vrows,1:2);
scale = 0.15;
arrowloc = [w(Vrows,1:2), w(Vrows,1:2)+scale*diff./kron(vecnorm(diff,2,2),ones(1,2))]';
disp ("write out some good gnuplot arrow commands -- gradient system")
for i = 1: length(Vrows)
  txt = sprintf ("set arrow %i from %g, %g to %g, %g filled size graph 0.025, 15 lw 1\n", ...
          i+5, arrowloc(1:4,i))
end

rhs.m


1
2
3
4
5
6
7
function wdot = rhs(t, w)
  x = w(1);
  y = w(2);
  a = -pi*sin(pi/2*x) -0.4*(x-2);
  b = -20*y;
  c = -(a^2+b^2);
  wdot = [a; b; c];