Pendulum - DAE formulation
> restart: libname := libname, "C:/MBSymba": with(MBSymba_r6):
> PDEtools[declare]([x(t),y(t),theta(t)],prime=t,quiet):
body definition
Define the body frame
> T1 := translate(x(t),y(t),0) . rotate('Z',theta(t)) ;
Define the body
> pendulum := make_BODY(T1,m,0,0,IG): show(pendulum);
Define acceleration due to gravity
> _gravity := make_VECTOR(ground,0,-g,0): show(_gravity);
pin joint
Define the pin joint A
> A := make_POINT(T1,0,L,0): show(A);
> Describe(spherical_JOINT);
# Create a spherical joint (3dof locked) between two points (kinematic
# constraint). When bodies are included, also the action-reaction
# forces/torques are considered. (mechanical constraint)
spherical_JOINT( P1::POINT, P2::POINT, B1::BODY := NULL,
B2::{BODY, mframe} := master_frame, suffix := "" )
> PDEtools[declare]([FxA(t),FyA(t)],prime=t,quiet):
FzA(t):=0; planar system
> revj_A, FA := spherical_JOINT( A,origin(ground),pendulum,"A"): show(revj_A); show(FA);
Define the damping torque
> make_VECTOR(ground,0,0,-c*diff(theta(t),t)):
DT := make_TORQUE(%, pendulum): show(DT);
Newton's equations
> eqnsN := newton_equations({ pendulum, FA, DT}): show(eqnsN);
Euler's equations with respect to the point A, in moving frame
> eqnsE := euler_equations({ pendulum, FA, DT},CoM(pendulum)): show(eqnsE);
> eqnsNE := [comp_Z(eqnsE), comp_X(eqnsN),comp_Y(eqnsN)]: 'eqnsNE'=Vector(%);
Lagrange's equations
> LF := kinetic_energy(pendulum) - gravitational_energy(pendulum) - sum( revj_A[expr][n] * revj_A[vars][n], n=1..nops(revj_A[vars]));
> lagrange(LF,theta(t),t)-generalized_force(DT,theta(t));
syntetically
> eqnsL1 := lagrange_equations({pendulum,revj_A,DT},[theta(t),x(t),y(t)]): 'eqnsL1'=Vector(eqnsL1);
In this case Lagrange's equations are equal to Newton's ones, indeed:
> simplify(leqns1 - dyna_eqns);
Alternatively, Lagrange's equations may be derived by considering reaction forces instead of constraint expressions
> eqnsL2 := lagrange_equations({pendulum,FA,DT},[theta(t),x(t),y(t)]): 'eqnsL2'=Vector(eqnsL2);
check
> simplify(eqnsL2 - eqnsL1);
DAE formulation and numerical solution
> eqns_DAE := eqnsNE union revj_A[expr][1..2]: 'eqns_DAE' = Vector(%);
> interface(displayprecision=2):
> params:= [L=1,IG=1,m=1,g=9.81,c=0.7];
> IC := {theta(0)=0,D(theta)(0)=3.5,
x(0)= 0, D(x)(0)=3.5,
y(0)=-1, D(y)(0)=0,
FxA(0)=0,FyA(0)=0};
> subs(t=0,IC,revj_A[expr]);
> eqns_DAEN:= convert(subs(params,eqns_DAE),set):
sys:= convert(subs(params,eqns_DAE),set) union IC:
solz:= dsolve(sys,numeric,range=0..5):
> solz(0); solz(5);
> npt:=200:
plots[odeplot](solz, [t,theta(t)], color=blue,labels=["t","theta"],numpoints=npt);
> plots[odeplot](solz, [x(t), y(t)], color=red, numpoints=npt,scaling=constrained);
> plots[odeplot](solz, [[t,FxA(t)],[t,FyA(t)]], color=[magenta,blue],labels=["t","forces"],legend=["FxA","FyA"],numpoints=npt);
>