This is a program discribing the dynamic behavior of a N pendelum system. When I tried to comparae the result of different initial conditions, I found that every time i run it, the graph will be different.

I know that double pendelum system is a typical chaotic system, is this the reason why i have different result for each time?

Here I attach the file, but plz rename it to ***.mws file before run, b/c i wrote it in maple classic worksheet.

> restart:
interface(warnlevel=0):
with(SolveTools):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,q=the initial angle of the beam,and qt=the initial angle velocity of the beam
beams:=proc(n,la,ma,q,qt)
global l,m,x,y,theta,icq,icqt:
icq[n]:=theta[n](0)=q:
icqt[n]:=D(theta[n])(0)=qt:
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;

i got the solution of a ODE system (dsn), the results are time dependent variable theta[i](t),and its deriavtives diff(theta[i](t),t). and i have one set of expressions, which are explicit form of theta[i](t),and diff(theta[i](t),t), such as ss[1]:=fx[1]=sin(theta[1](t))*3+cos(theta[1](t))*diff(theta[2](t),t). I want to put the solution i got into the explicit expression of fx[i], and then polt them, how can I do it? i use eval(ss[1],dsn)but it does not work. here is the code: > restart: interface(warnlevel=0): with(SolveTools):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:

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:

i have this im plicit ODE system in modeling a multibody dynamic system. theta[1](t), theta[2](t), are the angle displacement of the beams, and they are time depentent, in these 6 equations, i have the term of cos(theta[1](t)), sin(theta[2](t)),cos(theta[2](t)), sin(theta[1](t)), multipling (diff(theta[1](t),t))^2, and diff(theta[1](t),t$2), which make them implicit equations. there are 6 unknown variables fx1,fx0,fy1,fy0(which are the forces at the end of the beams),and theta[1](t), theta[2](t). i gave the initial condition as ics:={theta[1](0) = 1/6*Pi, theta[2](0) = 1/3*Pi, D(theta[1])(0) = 0, D(theta[2])(0) = 0} Kindly let me know how to write the code to solve it! Thanks. Still below is the code: >