>	restart:with(plots):with(DEtools):with(VariationalCalculus):

 

Here we set up the equation for the planar rotating chain in steps. First the catenary in polar coordinates, where the angle is measured from the negative y axis, as in spherical coordinates. 

>	f:=(r(phi)*cos(phi)+lambda)*sqrt(r(phi)^2+diff(r(phi),phi)^2);

 

(1)

 

These are the variables needed to do the EL substitutions so that Maple can differentiate with respect to a variable, and not a function. 

>	S1:={r(phi)=r,diff(r(phi),phi)=dr}; S2:={r=r(phi),dr=diff(r(phi),phi)};

 

 

	(2)

 

The Euler-Lagrange equation 

>	L1:=subs(S2,diff(subs(S1,f),r))-diff(subs(S2,diff(subs(S1,f),dr)),phi);

 

(3)

 

The Euler-Lagrange equation as calculated by Maple 

>	EulerLagrange(f,phi,[r(phi)]);

 

(4)

 

Put the EL eqn in the form of a standard differential equation 

>	DE1:=diff(r(phi),phi,phi)=simplify(solve(L1=0,diff(r(phi),phi,phi)));

 

(5)

 

We write the solution as a numerical procedure 

>	SDE1:=dsolve({DE1,r(0)=1,D(r)(0)=0},numeric,parameters=[lambda],output=listprocedure); SDE1(parameters=[lambda=-1.62]); SDE1(1);

 

 

 

	(6)

 

Now we can plot this as a curve in cartesian coordinates -- it looks just as we know it should. Note that we must choose lambda, the Lagrange multiplier, so that the endpoints land at 1, -1. Also, the initial conditions had to be chosen so that r = 1 and at phi = 0 because Maple won't solve this as a BVP. Might have to use a shooting method to get this done numerically. This restriction will give us trouble later. Luckily, it seems like the numerical procedure gives us sensible results for the next two variants, except for the unfortunate ICs necessary. 

>	evalf(Pi/2); plot({cos(k)*rhs(SDE1(k)[2]),sin(k)*rhs(SDE1(k)[2]),k=-Pi/2..Pi/2});

 

Now add a centrifugal potential 

>	g:=(r(phi)*cos(phi)+lambda+(1/2)*r(phi)^2*omega^2)*sqrt(r(phi)^2+diff(r(phi),phi)^2);

 

(7)

 

>	L2:=subs(S2,diff(subs(S1,g),r))-diff(subs(S2,diff(subs(S1,g),dr)),phi);

 

(8)

 

>	DE2:=diff(r(phi),phi,phi)=simplify(solve(L2=0,diff(r(phi),phi,phi)));

 

(9)

 

>	solve(DE2,lambda);

 

(10)

 

>	SDE2:=dsolve({DE2,r(0)=1,D(r)(0)=0},numeric,parameters=[lambda, omega],output=listprocedure); SDE2(parameters=[lambda=-15, omega=3]); SDE2(1);

 

 

 

	(11)

 

Once again, the plot seems realistic, but this symmetric solution should be unstable. Unfortunately, it is forced by the initial conditions. 

>	plot({cos(k)*rhs(SDE2(k)[2]),sin(k)*rhs(SDE2(k)[2]),k=-Pi/2..Pi/2},numpoints=300);

 

We can also specify different initial conditions and see we get realistic asymmetric solutions. 

>	SDE2a:=dsolve({DE2,r(Pi/2)=1,D(r)(Pi/2)=-0.32},numeric, parameters=[lambda, omega],output=listprocedure); SDE2a(parameters=[lambda=-22.4, omega=4.4]); SDE2a(1);

 

 

 

	(12)

 

Once again, the plot seems almost realistic, but we cannot easily adjust lambda to get the second endpoint to his the right spot. I conjecture that this is because the derivative and lambda have to be adjusted simultaneously. 

>	plot({cos(k)*rhs(SDE2a(k)[2]),sin(k)*rhs(SDE2a(k)[2]),k=-Pi/2..Pi/2},numpoints=300);

 

Here is the BVP version for the same thing 

>	SDE2b:=dsolve({subs({lambda=-10.8, omega=2.5},DE2),r(Pi/2)=1,r(-Pi/2)=1},numeric,method=bvp,output=listprocedure);

 

(13)

 

>	plot({cos(k)*rhs(SDE2b(k)[2]),sin(k)*rhs(SDE2b(k)[2]),k=-Pi/2..Pi/2},numpoints=300);

 

Now we let the plane rise up by an angle alpha. Gamma1 is the scaling factor for the centrifugal potential at plane angle alpha, gamma2 is the scaling factor for the gravitational force felt by the chain element. 

>	gamma1:=sqrt(sin(phi)^2+(cos(phi)^2)*sin(alpha)^2); gamma2:=cos(alpha);

 

 

	(14)

 

>	h:=(gamma2*r(phi)*cos(phi)+lambda+(1/2)*(gamma1*r(phi))^2*omega^2)*sqrt(r(phi)^2+diff(r(phi),phi)^2);

 

(15)

 

>	L3:=subs(S2,diff(subs(S1,h),r))-diff(subs(S2,diff(subs(S1,h),dr)),phi);

 

(16)

 

>	DE3:=diff(r(phi),phi,phi)=simplify(solve(L3=0,diff(r(phi),phi,phi)));

 

(17)

 

>	SDE3:=dsolve({DE3,r(Pi/2)=1,D(r)(Pi/2)=0.92},numeric,parameters=[lambda, omega, alpha],output=listprocedure); SDE3(parameters=[lambda=-11.9, omega=3, alpha=Pi/4]); SDE3(1);

 

 

 

	(18)

 

This time we were able to adjust lambda and the derivative at the endpoint simultaneously.  

>	plot({cos(k)*rhs(SDE3(k)[2]),sin(k)*rhs(SDE3(k)[2]),k=-Pi/2..Pi/2},numpoints=300);

 

Here is the BVP solution 

>	SDE3a:=dsolve({subs({alpha=Pi/8, lambda=-6, omega=2.75},DE3),r(Pi/2)=1,r(-Pi/2)=1},numeric,method=bvp,output=listprocedure);

 

(19)

 

>	plot({cos(k)*rhs(SDE3a(k)[2]),sin(k)*rhs(SDE3a(k)[2]),k=-Pi/2..Pi/2},numpoints=300);

 

>