Zymplectic

The Kapitza pendulum is a pendulum with a non-autonomous pivot whose motion depends explicitly on time. The motion of the pivot is typically vertical, although we will here consider the generalized Kapitza pendulum which oscillates in both the vertical and horizontal direction.

Derivation of the Hamiltonian

Define first the variables

l as the length of the pendulum from the pivot to the bob
m as the mass of the bob
θ is the angle between the vertical axis and the pendulum measured from top
g is the gravitational constant
a and b are the oscillation amplitudes in the first axis and second axis
ω₁ and ω₂ are the angular frequency of the oscillations
k₁ and k₂ are the phases of the oscillation

The coordinates of the pendulum bob are trivially expressed in terms of θ and the variables associated with the non-autonomous motion:

\(\displaystyle{\qquad x = l\sin \left( {\theta \left( t \right)} \right) + a\sin \left( {{\omega _1}t + {k_1}} \right)}\)

\(\displaystyle{\qquad y = l\cos \left( {\theta \left( t \right)} \right) + b\sin \left( {{\omega _2}t + {k_2}} \right)}\)

The velocities are obtained by differentiation:

\(\displaystyle{\qquad \dot x = l\cos \left( {\theta \left( t \right)} \right)\dot \theta \left( t \right) + a{\omega _1}\cos \left( {{\omega _1}t + {k_1}} \right)}\)

\(\displaystyle{\qquad \dot y = - l\sin \left( {\theta \left( t \right)} \right)\dot \theta \left( t \right) + b{\omega _2}\cos \left( {{\omega _2}t + {k_2}} \right)}\)

The kinetic energy T and the potential energy V are expressed by the velocity and height of the pendulum bob:

\(\displaystyle{ \qquad V = mgy }\)

\(\displaystyle{ \qquad T = \frac{1}{2}m ({\dot x}^2 +{\dot y}^2) }\)

Inserting positions and velocities yields:

\(\displaystyle{ \qquad V = mg\left( {l\cos \theta + b\sin \left( {{\omega _2}t + {k_2}} \right)} \right) }\)

\(\displaystyle{ \qquad T = \frac{1}{2}m\left( {{{\left( {a{\omega _1}\cos \left( {{\omega _1}t + {k_1}} \right)} \right)}^2} + {{\left( {b{\omega _2}\cos \left( {{\omega _2}t + {k_2}} \right)} \right)}^2} + {{\left( {l\dot \theta } \right)}^2}} \right) + ml\dot \theta \left( {a{\omega _1}\cos \theta cos\left( {{\omega _1}t + {k_1}} \right) - b{\omega _2}\sin \theta cos\left( {{\omega _2}t + {k_2}} \right)} \right) }\)

The Lagrangian L is now known L = T - V and utilized only to express the canonical momentum p_θ.

\(\displaystyle{ \qquad {p_\theta } = \frac{{\partial L}}{{\partial \dot \theta }} = lm\left( {a{\omega _1}\cos \theta \cos \left( {{\omega _1}t + {k_1}} \right) - b{\omega _2}\sin \theta \cos \left( {{\omega _2}t + {k_2}} \right) + l\dot \theta } \right) }\)

\(\displaystyle{ \qquad l\dot \theta = \frac{{{p_\theta }}}{{lm}} - a{\omega _1}\cos \theta \cos \left( {{\omega _1}t + {k_1}} \right) + b{\omega _2}\sin \theta \cos \left( {{\omega _2}t + {k_2}} \right) }\)

It follows from substitution of the angular velocity in the energies V and T that the energies are expressed only in terms of the canonical coordinates:

\(\displaystyle{ \qquad V\left( \theta \right) = mg\left( {l\cos \theta + b\sin \left( {{\omega _2}t + {k_2}} \right)} \right) }\)

\(\displaystyle{ \qquad T\left( {\theta ,{p_\theta }} \right) = \frac{1}{2}m\left( {{{\left( {a{\omega _1}\sin \theta \cos \left( {{\omega _1}t + {k_1}} \right)} \right)}^2} + {{\left( {b{\omega _2}\cos \theta \cos \left( {{\omega _2}t + {k_2}} \right)} \right)}^2}} \right) + abm{\omega _1}{\omega _2}\sin \theta cos\theta cos\left( {{\omega _1}t + {k_1}} \right)\cos \left( {{\omega _2}t + {k_2}} \right) + \frac{{{p_\theta }^2}}{{2{l^2}m}} }\)

The Hamiltonian H = T + V does not differentiate the contents of T and V, so the expressions may be re-arranged to become separable:

\(\displaystyle{ \qquad V\left( \theta \right) = \frac{1}{2}m\left( {{{\left( {a{\omega _1}\sin \theta \cos \left( {{\omega _1}t + {k_1}} \right)} \right)}^2} + {{\left( {b{\omega _2}\cos \theta \cos \left( {{\omega _2}t + {k_2}} \right)} \right)}^2}} \right) + abm{\omega _1}{\omega _2}\sin \theta cos\theta cos\left( {{\omega _1}t + {k_1}} \right)\cos \left( {{\omega _2}t + {k_2}} \right) + mg\left( {l\cos \theta + b\sin \left( {{\omega _2}t + {k_2}} \right)} \right) }\)

\(\displaystyle{ \qquad T\left( {{p_\theta }} \right) = \frac{{{p_\theta }^2}}{{2{l^2}m}} }\)

Finally we obtain the Hamiltonian gradients useful for simulating the system:

\(\displaystyle{ \qquad \frac{{\partial V}}{{\partial q}} = m\left( {ab{\omega _1}{\omega _2}\cos \left( {{\omega _1}t + {k_1}} \right)\cos \left( {{\omega _2}t + {k_2}} \right)\left( {2co{s^2}\left( \theta \right) - 1} \right) + \sin \theta \cos \theta \left( {{{\left( {a{\omega _1}\cos \left( {{\omega _1}t + {k_1}} \right)} \right)}^2} - {{\left( {b{\omega _2}\cos \left( {{\omega _2}t + {k_2}} \right)} \right)}^2}} \right) - gl\sin \theta } \right) }\)

\(\displaystyle{ \qquad \frac{{\partial T}}{{\partial p}} = \frac{{{p_\theta }}}{{{l^2}m}} }\)

Implementation

The system is implemented as any separable system. Zymplectic automatically detects the system as a non-autonomous system because the time "tau" is used in code. The following code simulates the motion of the pendulum bob with a circular motion of the pivot.

//Copy-paste into Zymplectic or load as file to run the simulation
#include "math.h"
const double m = 1.0; //pendulum bob mass
const double l = 2.0; //pendulum length
const double g = 9.8; //gravitational constant
const double a = 1.0, w1 = 1.0, k1 = 0.0; //oscillation amplitude, angular frequency and phase in the x-direction
const double b = 1.0, w2 = 1.0, k2 = 3.14159265358979323846/2; //... in the y-direction

double energy_V(Arg) {
	return m*g*(l*cos(q1) + b*sin(w2*tau + k2)) + 0.5*m*(pow(a*w1*sin(q1)*cos(w1*tau + k1),2) 
+ pow(b*w2*cos(q1)*cos(w2*tau + k2),2)) + a*b*m*w1*w2*sin(q1)*cos(q1)*cos(w1*tau + k1)*cos(w2*tau + k2);
}
void dV(Arg) {
	d1 = m*(a*b*w1*w2*cos(w1*tau + k1)*cos(w2*tau + k2)*(2*cos(q1)*cos(q1)-1.0) 
+ sin(q1)*cos(q1)*(pow(a*w1*cos(w1*tau + k1),2) + pow(b*w2*cos(w2*tau + k2),2)) - g*l*sin(q1));
}

double energy_T(Arg) {
	return p1*p1/(2.0*l*l*m);
}
void dT(Arg) {
	d1 = p1/(l*l*m);
}

void IterFunc(){
	double x = a*sin(w1*tau + k1) + l*sin(obj_q[0][0]);
	double y = b*sin(w2*tau + k2) + l*cos(obj_q[0][0]);
	draw.track_add(x,y); //obtains the trajectory of the pendulum bob
}

void DrawFunc(){
	double GLX[2],GLY[2];
	GLX[0] = a*sin(w1*tau + k1);
	GLY[0] = b*sin(w2*tau + k2);
	GLX[1] = GLX[0] + l*sin(obj_q[0][0]);
	GLY[1] = GLY[0] + l*cos(obj_q[0][0]);
	draw.points(GLX,GLY,2);
	draw.lines(GLX,GLY,2);
}

double _q_init[] = {0.5}; //pendulum angle (0 is vertical top)
double _p_init[] = {0.0}; //momentum

void main() { 
	Z.time_per_evaluation = 0.001;
	Z.track = 1;
	Z.delay = 2;
	Z.OnDraw = DrawFunc; //used for drawing the pendulum
	Z.OnIter = IterFunc; //used for obtaining the trajectory
	Z.xmin = -5;
	Z.xmax = 5;
	Z.ymin = -5;
	Z.ymax = 5;

	//--- required ---
	Z.T(energy_T,dT,0);
	Z.V(energy_V,dV,0);	
}

Qualitative graphical depictions of the Kapitza pendulum

The Kapitza pendulum with a circulating pivot is shown below:

The classical Kapitza pendulum is obtained by setting the oscillation amplitude a to 0. The pivot then moves only in the vertical direction, i.e. along the green line as shown in the following figure.

The oscillation frequency associated with w2 is large enough to create a stable trajectory of the pendulum bob which therefore never drops down. Balancing the pendulum in an inverted position is but one rich subject which require a comprehensive analysis and discussions omitted here.

Zymplectic project

Kapitza pendulum

Derivation of the Hamiltonian

Implementation

Qualitative graphical depictions of the Kapitza pendulum