drawing function phase space and draw the associated vector field

hello,

I'm trying to draw the following in phas-space, and to show the vector field, but I'm getting some very strange results. can anyone please explain to me how can I draw this function in phase-space and how to show the vector field?

I'm trying to plot elevation lines where H=constant, and to show the associated vector field.

my attempt:

m=1; R=1;Om=1;g=10;
x= linspace(-100,100,100000);
y=sqrt(m*(R^2)*(Om^2)*(sin(x).^2)-2*(m^2)*g*(R^3)*(cos(x)));

quiver(x,y);

thank you!

[David Goodmanson]

Hi Elinor,

could you describe what physical situation is represented by this eqn?

[Elinor Ginzburg]

hello, yes, I'm trying to draw the vector field and the elevetion line of the Hamiltonian describing Kapitza's pendulum for and

thank you for you attention!

[David Goodmanson]

since this represents a second order equation, are you looking to get a 2d contour plot in phase space? That would involve constructing the phase plane

th1 = -3*pi:0.1:3*pi;          % suggested limits
p1 = -p1_lim:spacing:p1_lim;   % for some appropriate limits and spacing
[th p] = meshgrid(th1,p1);

then calculating H and doing a contour plot with th,p, and H.

[Elinor Ginzburg]

thank you for the answer, unfortunately I still can't plot it' I've tried doing what you've suggested, but perhaps I misunderstood something? here's my code:

m=1; R=1; g=10; Om=g/R;
th1 = -3*pi:0.1:3*pi;          % suggested limits
p1 = -5:0.1:5;   % for some appropriate limits and spacing
H=sqrt(m*(R^2)*(Om^2)*(sin(th1).^2)-2*(m^2)*g*(R^3)*(cos(th1)));
[Th P] = meshgrid(th1,p1);

%x = linspace(-2*pi,2*pi);
%y = linspace(0,4*pi);
%[X,Y] = meshgrid(x,y);
%Z = sin(X)+cos(Y);
contour(real(Th),real(P),real(H))

I get the following error:

Error using contour (line 48)
Z must be at least a 2x2 matrix.

Error in pendulum_phas_space (line 11)
contour(real(Th),real(P),real(H))

But I have no idea how to make H a 2x2 matrix. how can I properly write H?

thank you very much!

[David Goodmanson]

Hi Elinor,

there are a couple of things here. The idea is not to calculate p as a function of theta in some manner, as you appear to be doing. Rather, both theta and p are independent variables. The idea is to set up a 2d plane with all possible pairs of values for p and theta, and then calculate H in the theta,p plane. Meshgrid allows you to do that, so meshgrid needs to be done before calculatiing H. And H is simply the expression you provided in the original question, except there is a mistake. The term involving Om^2 should have a + sign..

m = 2;
R = 5;
Om = 0;
g = 9.8;
thetavalues = -3*pi:.1:3*pi;
Lvalues = (-200:.2:200);
[theta,L] = meshgrid(thetavalues,Lvalues);
H = L.^2/(2*m*R^2) + (1/2)*m*R^2*Om^2*sin(theta).^2 - m*g*R*cos(theta);

Using p for the momentum is perfectly acceptable, but since the momentum in this problem is the angular momentum, I changed p to L.

Now you can do a contour plot with theta,L and H. The plot does not really provide enough contour lines, but if you take a look at 'help contour' you can see how to provide more contour lines with one more argument to the contour function. It helps to use 'grid on' just after the contour command, also 'colorbar'.

With Om = 0 you should see stable equilibrium points at theta = -2*pi, 0 and 2*pi, an ordinary pendulum at the low point of its swing. As Om is increased, eventually you should see more equilibrium points appear, at theta values corresponding to the pendulum at the top of its swing.

[Elinor Ginzburg]

thank you very much! I understand now how to use it. I'm still new to matlab, so sometimes it's hard to search the documentation for something specific. thank you again for your help!