lagrangian (mech3),peter.vlasschaert@gmail.com, 25/07/2017
double pendulum 'two degrees of freedom q1,q2 '
---------------------------------------------------------------------------------------------
=> l1,theta1 = q1(ref vertical= y axis),'deviation'),m1=mass 1
=> l2,theta2 = q2(ref vertical= y axis),'deviation'),m2=mass 2
connection(origin) : fourth quadrant.
(l1,theta1,m1,l2,theta2,m2)
---------------------------------------------------------------------------------------------
ref : pE = 0  'pos x axis'
---------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------
rem: pE = potential Energy.
qi ,i=1,2 ('see below')

kill(all)$
depends(lag,[q1,q2,q1_t,q2_t],[q1,q2,q1_t,q2_t],t);
/*  ' (=prime),prevent evalution of the equation */
p1:'diff(lag,q1)-'diff(diff(lag,q1_t),t)=0;
p2:'diff(lag,q2)-'diff(diff(lag,q2_t),t)=0;

potential Energy (=pE),kinetic Energy(=kE)

given :

assumption :variables 'physics'

depends([theta1,theta2],t);
/*depends(theta1,t);*/

assume(l1>0,l2>0,m1>0,m2>0,g>0);

pE1:-m1*g*l1*cos(theta1)-m2*l2*g*cos(theta2)-m2*l1*g*cos(theta1);
kE1a:1/2*(m1+m2)*l1^2*(diff(theta1,t))^2+1/2*m2*l2^2*(diff(theta2,t))^2;
kE2a:m2*l1.l2*diff(theta1,t)*diff(theta2,t)*cos(theta1-theta2);
kE1:kE1a+kE2a;

conversion : q,q_t

pE:ev(pE1,theta1=q1,theta2=q2);
kE:ev(kE1,diff(theta1,t)=q1_t,diff(theta2,t)=q2_t,theta1=q1,theta2=q2);

definition : Lagrangian =kinetic Energy -potential Energy

p3:lag=kE-pE;

no evalution of differentation operators

p4:'diff(lag,q1)-'diff(diff(lag,q1_t),t)=0;
p41:ev(p4,p3);
p5:'diff(lag,q2)-'diff(diff(lag,q2_t),t)=0;
p51:ev(p5,p3);

 evalution of differentation operators  : 'diff ' option with 'command :ev'

p6:ev(p41,diff);

p7:ev(p51,diff);

definition : acceleration (velocity=qi_t,position=qi,i=1,2)

p61:diff(q1_t,t)=diff(q1,t,2);
p71:diff(q2_t,t)=diff(q2,t,2);

simplify:radcan,trigreduce

p8:radcan(ev(p6,p61,p71));
p9:radcan(ev(p7,p61,p71));

two differential equations p10,p11 for 'q1(t),q2(t)'

1e) velocity

p91:q1_t=diff(q1,t);
p92:q2_t=diff(q2,t);

2e)  ode's

p10:ev(trigreduce(p8),p91,p92);
p11:ev(trigreduce(p9),p91,p92);

rem: 2 ode's,with use ( command :factor)

p12:factor(p10);
p13:factor(p11);
p14:part(p13,1,1,2)=rhs(p13);