# Question:a multibody dynamic problem(solving an 2nd order ODE system)

## Question:a multibody dynamic problem(solving an 2nd order ODE system)

Maple
several beams are connected by gemel one after another, the 1st is fixed at one end on the roof, and the last is driven by a force(fx,fy vector summation) at the end point. i am using lagrange method for the virtual displacement of angle theta[i] , and momentum balance in 2 direction (x,and y)to build the model. in the last step, i got an ODE system of 2nd order, but i cant code in maple to solve it. i think it is solveable, (6 equations, 6 unknowns). Kindly let me know how to write the code to solve it! below is the code.

> restart:
interface(warnlevel=0):
with( plots ): with( plottools ):
with( DEtools ): with( PDEtools ):
g:=9.81:
x:=0:
y:=0:
#initialization function of each beam,in this function, we should give 3 variables, n=the number of this #beam, la=the length of this beam, ma=the  mass of this beam
beams:=proc(n,la,ma)
global l,m,x,y:
l[n]:=la:           #length of  beam n
m[n]:=ma:           #mass of beam n
#initialize the positions of each beam
x[n]:=0:        #the x position of the end point of each beam
y[n]:=0:        #the x position of the end point of each beam
x[n]:=0:        #the y position of the center point of each beam
y[n]:=0:        #the y position of the center point of each beam
end proc;
#define the forces that applied on the end of  beam n
Force:=proc(n,fm,fn)
global fx,fy:
fx[n]:=fm:
fy[n]:=fn:
end proc:

> beams(1,1,1):
beams(2,2,1):
beams(3,3,1):
beams(4,4,1);

> Force(2,80,100);

> la:=proc(nn)
global x,y,q,vy,vx,T,lg,lm,lp,L,Ll,Lk,ax,ay:
n:=nn: #transfer the function varibale to local function
# generate the position of the beams, nn is the number of beams
# x,y is the position vector, q is the angle vector
for i from 1 to n do:
x[i]:=l[i]*cos(q[i](t))+x[(i-1)]:     #x[i],y[i] is the positon of the end of beam i
y[i]:=l[i]*sin(q[i](t))+y[(i-1)]:
x[i]:=x[i]-0.5*l[i]*cos(q[i](t)):     #x[i],y[i] is the center of beam i
y[i]:=y[i]-0.5*l[i]*sin(q[i](t)):
end do:
#generate the velocity of  the center of the beams
for i from 1 to n do:
vx[i]:=diff(x[i],t):
vy[i]:=diff(y[i],t):
ax[i]:=diff(x[i],t\$2):
ay[i]:=diff(y[i],t\$2):
end do:
#generate the kinetic energy of  the system
T:=0:
for i from 1 to n do:
T:=T+0.5*m[i]*(vx[i]^2+vy[i]^2)+1/24*m[i]*l[i]^2*diff(q[i](t),t)^2:
end do:
# substitue the angle and the angle velocity  with time independent variables, for later diffrenciation
for i from 1 to 4 do:
T:=subs(diff(q[i](t),t)=omega[i],T);   # angle velocity(q(t) is angle function)
T:=subs(q[i](t)=theta[i],T):           #angle (q(t) is angle function))
end do;
#generate the lagrange equations
#lg,lm,lp,L,Ll,Lk,those variables is the intermedia variables for Lagrange equation, explained detailly below
for i from 1 to n do:
lg[i]:=diff(T,omega[i]):   #derivation of "T" with respect to angle velocity

#substitute the angle and angle velocity in each lg[i] with time function, for later derivation over time
for j from 1 to n do:
lg[i]:=subs({theta[j]=theta[j](t),omega[j]=omega[j](t)},lg[i]):
end do;
end do:
#---------------------------------------------------------
for p from 1 to n do:
lp[p]:=diff(lg[p],t)  :#derivation of "lg" term  over time,make lp the first term in lagrange eq.
end do:
#--------------------------------------------------------
for o from 1 to n do:
lm[o]:=diff(T,theta[o]):  #derivation of "T" over angle, which is the 2nd term in the eq.
for k from 1 to n do:
#substitute the angle and angle velocity in each lm[o] with time function, for later derivation over time
lm[o]:=subs({theta[k]=theta[k](t),omega[k]=omega[k](t)},lm[o]):
end do:
end do:
#---------------------------------------------------------
for k from 1 to n do:
#L is the lagrange eq. with the virtual displacement  of angle
#Ll is the lagrange eq. with the virtual displacement of x
#Lk is the lagrange eq. with the virtual displacement of y
#unlike L,Ll and Lk are written directly,no deduction by programming
L[k]:=-lp[k]+lm[k]+(fy[k]+0.5*m[k]*g)*cos(theta[k](t))*l[k]-fx[k]*sin(theta[k](t))*l[k]=0:
#L[k]:=-lp[k]+lm[k]+(fy[k]+0.5*m[k]*g)*cos(theta[k](t))*l[k]-fx[k]*sin(theta[k](t))*l[k]-(fy[k]+0.5*m[k]*g)*theta[k](t)*sin(theta[k]##(t))*l[k]-fx[k]*cos(theta[k](t))*l[k]*theta[k](t)=0:
#because diff(T,x)=Diff(T,omega)*Diff(omega,x)
#Ll[k]:=lp[k]*(1/l[k])*cos(theta[k](t))*omega[k](t)/sin(theta[k](t))^2-lm[k]/(l[k]*sin(theta[k](t)))-fx[k]:
#because diff(T,y)=Diff(T,omega)*Diff(omega,y)
#Lk[k]:=lp[k]*(1/l[k])*sin(theta[k](t))*omega[k](t)/cos(theta[k](t))^2+lm[k]/(l[k]*cos(theta[k](t)))-fy[k]-m[#k]*g:
Ll[k]:=fx[k-1]-fx[k]=m[k]*ax[k]:
Lk[k]:=fy[k-1]-fy[k]=m[k]*ay[k]:
end do:
for s from 1 to n do:
for r from 1 to n do:
L[s]:=subs(diff(omega[r](t),t)=diff(theta[r](t),t,t),L[s]):
#Lk[s]:=subs(diff(omega[r](t),t)=diff(theta[r](t),t,t),Lk[s]):
Lk[s]:=subs(q[r](t)=theta[r](t),Lk[s]):
Ll[s]:=subs(q[r](t)=theta[r](t),Ll[s]):
#Ll[s]:=subs(diff(omega[r](t),t)=diff(theta[r](t),t,t),Ll[s]):
L[s]:=subs(omega[r](t)=diff(theta[r](t),t),L[s]):
#Lk[s]:=subs(omega[r](t)=diff(theta[r](t),t),Lk[s]):
#Ll[s]:=subs(omega[r](t)=diff(theta[r](t),t),Ll[s]):
L[s]:=combine(L[s],trig):
Ll[s]:=combine(Ll[s],trig):
Lk[s]:=combine(Lk[s],trig):
end do:
end do:

end proc;

> la(2);
Lk;
Lk;
Ll;
Ll;
L;
L;
ics:={theta(0)=Pi/6,theta(0)=Pi/3,D(theta)(0)=0,D(theta)(0)=0};

> #below is the problem, the code of solving the equations doesn't work.
eqs:={Ll,Ll,Lk,Lk,L,L};#equations
ics:={theta(0)=Pi/6,theta(0)=Pi/3,D(theta)(0)=0,D(theta)(0)=0};#initial conditions
#ic:=theta(0)=0,theta(0)=0,theta(0)=0,theta(0)=0,
#D(theta)(0)=0,D(theta)(0)=0,D(theta)(0)=0,D(theta)(0)=0,D(theta)(0)=0;
#theta(0)=0,theta(0)=0,theta(0)=0,theta(0)=0, #D(theta)(0)=0,D(theta)(0)=0,D(theta)(0)=0,D(theta)(0)=0;
#loes:=dsolve({Ll,Ll,Lk,Lk,L,L,theta(0)=0,theta(0)=0,D(theta)(0)=0,D(theta)(0)=#0},{theta(t),theta(t),fx,fy,fx,fy},numeric);
#sol:=dsolve(eqs union ics, numeric);
loes:=dsolve({Ll,Ll,Lk,Lk,L,L,theta(0)=.5,theta(0)=1,D(theta)(0)=0,D(theta)(0)=0},{theta(t),theta(t),fx,fx,fy,fy},type=numeric,output=listprocedure);#solving code that"suppose to work"
#Diff(theta(t),t),Diff(theta(t),t),Diff(theta(t),t),Diff(theta(t),t),

Error, (in dsolve/numeric/process_input) invalid argument: {theta(t), theta(t), fy, fx, fx, fy}

> #dsolve[interactive](L);

>

This post generated using the online HTML conversion tool 