A DAE formulation of the pendulum problem

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));


> 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);


> 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,

> 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);


Leave a Reply

Your email address will not be published. Required fields are marked *


You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>