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][2]:=0:

y[0][2]:=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][2]:=0: #the x position of the end point of each beam

y[n][2]:=0: #the x position of the end point of each beam

x[n][1]:=0: #the y position of the center point of each beam

y[n][1]:=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][2]:=l[i]*cos(q[i](t))+x[(i-1)][2]: #x[i][2],y[i][2] is the positon of the end of beam i

y[i][2]:=l[i]*sin(q[i](t))+y[(i-1)][2]:

x[i][1]:=x[i][2]-0.5*l[i]*cos(q[i](t)): #x[i][1],y[i][1] is the center of beam i

y[i][1]:=y[i][2]-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][1],t):

vy[i]:=diff(y[i][1],t):

ax[i]:=diff(x[i][1],t$2):

ay[i]:=diff(y[i][1],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[1];

Lk[2];

Ll[1];

Ll[2];

L[1];

L[2];

ics:={theta[1](0)=Pi/6,theta[2](0)=Pi/3,D(theta[1])(0)=0,D(theta[2])(0)=0};

> #below is the problem, the code of solving the equations doesn't work.

eqs:={Ll[1],Ll[2],Lk[1],Lk[2],L[1],L[2]};#equations

ics:={theta[1](0)=Pi/6,theta[2](0)=Pi/3,D(theta[1])(0)=0,D(theta[2])(0)=0};#initial conditions

#ic:=theta[1](0)=0,theta[2](0)=0,theta[3](0)=0,theta[4](0)=0,

#D(theta[1])(0)=0,D(theta[2])(0)=0,D(theta[3])(0)=0,D(theta[3])(0)=0,D(theta[4])(0)=0;

#theta[1](0)=0,theta[2](0)=0,theta[3](0)=0,theta[4](0)=0, #D(theta[1])(0)=0,D(theta[2])(0)=0,D(theta[3])(0)=0,D(theta[4])(0)=0;

#loes:=dsolve({Ll[1],Ll[2],Lk[1],Lk[2],L[1],L[2],theta[1](0)=0,theta[2](0)=0,D(theta[1])(0)=0,D(theta[2])(0)=#0},{theta[1](t),theta[2](t),fx[0],fy[0],fx[1],fy[1]},numeric);

#sol:=dsolve(eqs union ics, numeric);

loes:=dsolve({Ll[1],Ll[2],Lk[1],Lk[2],L[1],L[2],theta[1](0)=.5,theta[2](0)=1,D(theta[1])(0)=0,D(theta[2])(0)=0},{theta[1](t),theta[2](t),fx[1],fx[0],fy[1],fy[0]},type=numeric,output=listprocedure);#solving code that"suppose to work"

#Diff(theta[1](t),t),Diff(theta[2](t),t),Diff(theta[3](t),t),Diff(theta[4](t),t),

Error, (in dsolve/numeric/process_input) invalid argument: {theta[2](t), theta[1](t), fy[1], fx[1], fx[0], fy[0]}

> #dsolve[interactive](L[4]);

**> **

This post generated using the online HTML conversion tool

Download the original worksheet