lagrangian (mech8),peter.vlasschaert@gmail.com, 02/08/2017
q=general coordinates   (qi ,i=1,2)  here x,y
q_t= time derivative from general coordinates.(qi_t,i=1,2) here x_t,y_t

Compound Atwood Machine:

ref : ceiling ,eP=0

! positive direction is from ceiling to the floor. 
compound => First Pulley  1            :l=l_1 ,x
                      => Second  Pulley 2      :l=l_2 ,y
l_1 = Abitrary constant
l_2 = Abitrary constant
--------------    =  (reason choose ,other ref :x+fixedvalue (= center  Pulley 1 to the ceiling))
=>                 eP=0 (--------------------)

assume(m_1>0,m_2>0,m_3>0,g>0);
depends(lag,[x,y,x_t,y_t],[x,y,x_t,y_t],t);

position m_i   of  x_i  (i=1,2,3)

x_1:x;
x_2:l_1-x+y;
x_3:l_1-x+l_2-y;

velocity x_i =>  vi_t

a) m1 < m2+m3
b) m3 > m2
=>   follows from the physics vi_t (i=1,2,3)

assume(m1<m2+m3)$
assume(m3>m2)$

you don't need  %i7 gives the same result.

v1_t :-x_t;
v2_t:-y_t+x_t;
v3_t:+y_t+x_t;

eK:1/2*m_1*(v1_t)^2+1/2*m_2*(v2_t)^2+1/2*m_3*(v3_t)^2;
eP:-(m_1*g*x_1+m_2*g*x_2+m_3*g*x_3);

below the ceiling :reason minus sign in the front.

definition : Lagrangian(=lag) =kinetic Energy -potential Energy=eK-eP

lag:eK-eP;

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

p1:double(lag,q,q_t):=diff(lag,q)-diff(diff(lag,q_t),t);

Euler-Lagrange equation(x,x_t)

p2:double(lag,x,x_t)=0;
p21:ev(p2,x_t=diff(x,t),y_t=diff(y,t));
p211:factor(ev(part(p21,1),diff(x,t,2)=0,diff(y,t,2)=0));
p212:ratcoef(part(p21,1),diff(x,t,2));
p213:ratcoef(part(p21,1),diff(y,t,2));

p214:(p212*'diff(x,t,2)+p213*'diff(y,t,2)+p211-p211=-p211);

Euler-Lagrange equation(y,y_t)

p3:double(lag,y,y_t)=0;
p31:ev(p3,x_t=diff(x,t),y_t=diff(y,t));
p311:factor(ev(part(p31,1),diff(x,t,2)=0,diff(y,t,2)=0));
p312:ratcoef(part(p31,1),diff(x,t,2));
p313:ratcoef(part(p31,1),diff(y,t,2));

p314:p312*'diff(x,t,2)+p313*'diff(y,t,2)+p311-p311=-p311;

How to solve for diff(x,t,2) or diff(y,t,2) from 'p214 and p314'

1e) solve for x=x(t)

p4:part(solve([p314], [diff(y,t,2)]),1);

p41:ev(p214,p4);

p42:part(solve(p41,diff(x,t,2)),1);

p43:ode2(p42,x,t);

1e) solve for y=y(t)

p5:part(solve([p314], [diff(x,t,2)]),1);
p51:ev(p214,p5);
p52:part(solve(p51,diff(y,t,2)),1);
p53:ode2(p52,y,t)$
p54:ev(p53,%k1 = %k3,%k2=%k4);

find v_i => mi (i=1,2,3), v=velocity

p61:v_1=ev(v1_t,x_t=diff(x,t),y_t=diff(y,t));
p62:v_2=ev(v2_t,x_t=diff(x,t),y_t=diff(y,t));
p63:v_3=ev(v3_t,x_t=diff(x,t),y_t=diff(y,t));

p611:ev(ev(p61,p43,p54),diff);
p622:ev(ev(p62,p43,p54),diff);
p633:ev(ev(p63,p43,p54),diff);

find a_i => mi (i=1,2,3), a=acceleration

p711:a_1=diff(part(p611,2),t);
p722:a_2=diff(part(p622,2),t);
p733:a_3=diff(part(p633,2),t);

check solution x=x(t),y=y(t)

check x=x(t),y=y(t) in 'p214'

p8:ratsimp(ev(ev(p214,p43,p54),diff));
p81:is(p8);

check x=x(t),y=y(t) in 'p314'

p9:ratsimp(ev(ev(p314,p43,p54),diff));
p91:is(p9);