On Oct. 25, Juan Pablo Carbajal wrote:
-------------------------
> I am trying to deduce the equations of motions of a kinematic chain
> given the Potential and Kinetic energy. I was wondering if any of the
> chapters deals with the Lagrangian or Hamiltonian formalism of
> classical mechanics.
------------------------------
Hi Juan,
None of my chapters cover lagrangian and hamiltonian
mechanics methods, as yet.
Some misc. code is in share/contrib/rand, small specialized
programs for mechanics analysis, but advanced and not
very useful unless you have the book by Richard Rand,
"Topics in Nonlinear Dynamics with Computer Algebra",
Gordon and Breach Sci. Publishers, 1994.
I have some code based on his lagrangian approach on
pages38-39, and I ignore the possibility of
constraint equations and incorporating them with
lagrange multipliers in what follows.
lag1.mac is a batch file which defines a one at a time
eqn. of motion constructor le(q), and also loads
file lag2.mac, which defines leN(qL,mylag) which
wants a list qL of generalized coordinates and
the lagrangian (ke - pe) as a function of q1,q2,...
and expects the corresponding generalized
velocities in the lagrangian to be q1d, q2d, etc.
Here is output from batch file lag1.mac, using
display2d:false so I can get the output into this
mailing list message coherently, but the results
look better it you run the batch file using
display2d:true (the default).
=====================
(%i1) batch("lag1.mac")$
read and interpret file: #pC:/work5/lag1.mac
(%i2) kill(all)
(%i1) display2d : false
(%i2) le(q):=expand(diff(diff(lag,diff(q,t)),t)-diff(lag,q)) = 0
(%i3) load("lag2.mac")
" leN (qL,Lag) " = " leN (qL,Lag) "
(%i4) " use le(q) for deriving one equation at a time "
(%i5) depends(x,t)
(%i6) xd:diff(x,t)
(%i7) ke:m*xd^2/2
(%i8) pe:V(x)
(%i9) lag:ke-pe
(%i10) eqx:le(x)
(%o10) 'diff(V(x),x,1)+m*'diff(x,t,2) = 0
(%i11) " case simple harmonic motion about x = 0 "
(%i12) V(x):=k*x^2/2
(%i13) eqx1:ev(eqx,diff)
(%o13) m*'diff(x,t,2)+k*x = 0
(%i14) " case vertical motion in grav. field, x axis up "
(%i15) V(x):=m*g*x
(%i16) eqx1:(ev(eqx,diff),expand(%%/m))
(%o16) 'diff(x,t,2)+g = 0
(%i17) "-------------------------------------------------"
(%i18) " case projectile motion in (x,y) plane, y axis up "
(%i19) depends(y,t)
(%i20) yd:diff(y,t)
(%i21) ke:m*yd^2/2+ke
(%i22) pe:V2(x,y)
(%i23) lag:ke-pe
(%i24) eqx:le(x)
(%o24) 'diff(V2(x,y),x,1)+m*'diff(x,t,2) = 0
(%i25) eqy:le(y)
(%o25) m*'diff(y,t,2)+'diff(V2(x,y),y,1) = 0
(%i26) " case no friction "
(%i27) V2(x,y):=m*g*y
(%i28) eqx1:(ev(eqx,diff),expand(%%/m))
(%o28) 'diff(x,t,2) = 0
(%i29) eqy1:(ev(eqy,diff),expand(%%/m))
(%o29) 'diff(y,t,2)+g = 0
(%i30) " ---------------------------------------------"
(%i31) " use leN(qL,Lag) for a list of equations of motion "
(%i32) " corresponding to list qL "
(%i33) kill(x,y)
(%i34) qL:[x,y]
(%i35) " set mass m = 1 here "
(%i36) mylag:-g*y+yd^2/2+xd^2/2
(%i37) leN(qL,mylag)
qdL = [xd,yd]
qLdiff = ['diff(x,t,1),'diff(y,t,1)]
qLR = [xd = 'diff(x,t,1),yd = 'diff(y,t,1)]
_Lag2% = ('diff(y,t,1))^2/2-g*y+('diff(x,t,1))^2/2
(%o37) ['diff(x,t,2) = 0,'diff(y,t,2)+g = 0]
===================================
Here is file lag1.mac:
===================
/* file lag1.mac batch file */
/* Richard Rand: Lagrangian, no constraint eqns, method
in Topics in Nonlinear Dynamics with Computer
Algebra,1994, p.38,39 */
kill(all)$
display2d:false$
/* generate each degree of freedom eqn. of motion
one at a time with: */
le(q) := expand (diff (diff (lag, diff (q,t)),t) - diff (lag,q)) = 0$
/* file lag2.mac has code for a list of gen. coordinates */
load ("lag2.mac")$
/* examples of le(q) use */
" use le(q) for deriving one equation at a time "$
depends(x,t)$
xd : diff(x,t)$
ke : m*xd^2/2$
pe : V(x)$
lag : ke - pe$
eqx : le (x);
" case simple harmonic motion about x = 0 "$
V(x) := k*x^2/2$
eqx1 : ev (eqx,diff);
" case vertical motion in grav. field, x axis up "$
V(x) := m*g*x$
eqx1 : (ev (eqx,diff), expand (%%/m));
"-------------------------------------------------"$
" case 2 degrees of freedom "$
" case projectile motion in (x,y) plane, y axis up "$
depends (y,t)$
yd : diff(y,t)$
ke : ke + m*yd^2/2$
pe : V2(x,y)$
lag : ke - pe$
eqx : le(x);
eqy : le(y);
" case no friction "$
V2(x,y) := m*g*y$
eqx1 : (ev(eqx,diff),expand (%%/m));
eqy1 : (ev(eqy,diff),expand (%%/m));
" ---------------------------------------------"$
" use leN(qL,Lag) for a list of equations of motion "$
" corresponding to list qL "$
kill(x,y)$
qL : [x,y]$
" set mass m = 1 here "$
mylag : xd^2/2 + yd^2/2 - g*y$
leN (qL, mylag);
===========================
and here is file lag2.mac:
=========================
/* file lag2.mac */
_le% (_q%,_lag%) := expand (diff (diff (_lag%, diff (_q%,t)),t) - diff
(_lag%,_q%)) = 0$
leN (_qL%,_Lag%) :=
block ([qdL,qLdiff,_Lag2% ],
qdL : map ('lambda ([x],concat (x,d)), _qL%),
display (qdL),
depends (_qL%,t),
qLdiff : diff (_qL%,t),
display (qLdiff),
qLR : map ("=",qdL, qLdiff),
display (qLR),
_Lag2% : subst (qLR, _Lag%),
display (_Lag2%),
map ('lambda ([x],_le% (x,_Lag2%)),_qL%))$
display (" leN (qL,Lag) ")$
===========================
I will put these two mechanics files on my webpage
tomorrow.
Best Wishes,
Ted Woollett
http://www.csulb.edu/~woollett
-------------- next part --------------
A non-text attachment was scrubbed...
Name: lag1.mac
Type: image/x-macpaint
Size: 1519 bytes
Desc: not available
URL: <http://www.math.utexas.edu/pipermail/maxima/attachments/20101025/8537be45/attachment.bin>;
-------------- next part --------------
A non-text attachment was scrubbed...
Name: lag2.mac
Type: image/x-macpaint
Size: 574 bytes
Desc: not available
URL: <http://www.math.utexas.edu/pipermail/maxima/attachments/20101025/8537be45/attachment-0001.bin>;