lagrangian (mech10),peter.vlasschaert@gmail.com, 03/08/2017
example : constraint lagrangian,curvilinear (use cylindrical coordinates)
particle of mass(= m ) moves on fixed hemisphere
r        = distance from the center to the particle              (r)
θ        = place where particle leaves hemisphere           (theta)
constraint
=>    r=a  => f(r,θ)= r-a    'one constraint' =>  find :λ
q=general coordinates   (qi ,i=1,2)  here r,θ
q_t= time derivative from general coordinates.(qi_t,i=1,2) here r_t,θ_t
given : see 'cylindrical  coordinates'  z=0 (2D)  , ρ = r,φ = θ
EK = 1/2*m*((diff(r,t))^2+r^2*(diff(θ,t))^2)

assume(r>0,a>0,m>0,g>0);
depends(lag,[r,θ,r_t,θ_t],[r,θ,r_t,θ_t],t);

constraint : 'moves on fixed hemisphere' f(r,θ) = r-a

f1:r=a;

f:part(f1,1)-part(f1,2);

height : particle at θ => h = r*cos(θ)

eK:1/2*m*((r_t)^2+r^2*(θ_t)^2);
eP:m*g*r* cos(θ);

lag:eK-eP;

=>   " Euler-Lagrange equation (qi,i=1,2)  "

p1:constraint(lag,q,q_t):=diff(lag,q)-diff(diff(lag,q_t),t)+λ*diff(f,q);

Euler-Lagrange equation(r,r_t)

p2:constraint(lag,r,r_t)=0;
p21:ev(p2,r_t=diff(r,t),θ_t=diff(θ,t));

Euler-Lagrange equation(θ,θ_t)

p3:constraint(lag,θ,θ_t)=0;
p31:ev(p3,r_t=diff(r,t),θ_t=diff(θ,t));

How to solve for   'λ'

p41:p21;
p42:p31;
p43:f1;

p5:diff(p43,t);
p6:diff(p5,t);

p51:ev(p41,p5,p6,p43);
p52:ev(p42,p5,p6,p43);

p61:part(solve(p52,diff(θ,t,2)),1);

p62:integrate(p61, θ);

? problem with first part integral , 'p62'
 =>    diff(θ,t,2) => trick
          rem :  diff(θ,t,2) = θ_t*d(θ_t)/dθ

p63:diff(θ,t,2) *dθ= θ_t*d(θ_t);
p64:p61*dθ;
p65:ev(p64,p63);
p66:integrate(part(p61,2), θ);
p67:ratsubst((diff(θ,t))^2/2,part(p65,1),part(p65,1));
p68:p67=p66+%k1;
;

find  %k1 :  begin position (t=0) then (θ = 0 )of particle with (mass=m)
         =>  diff(θ,t)  = 0 

p7:ev(p68,θ=0,diff(θ,t)=0);
p8:part(solve(p7,%k1),1);
p9:ev(p68,p8);

find :λ from     =>  'p9,p51'

p11:part(solve(p9,(diff(θ,t))^2),1);
p12:ev(p51,p11);
p13:part(solve(p12,λ),1);

what angle θ : particle falls from the hemisphere.
=> find angle θ when λ=0

p14:factor(ev(p13,λ=0));
p15:part(solve(p14,θ),1);

rem : m*g <> 0  only (3*cos(θ)-2) = 0