Question: Why my MAPLE program gives different result everytime I run it?

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;

> #Here you should set the initial condition of the system,and give the value of the characteristic variables
beams(1,2,3,1/6*Pi,0):
beams(2,2,3,1/3*Pi,0):
beams(3,2,3,-1/6*Pi,0):
beams(4,2,3,-1/3*Pi,0):
Force(4,100,100);

> la:=proc(nn)
global x,y,q,vy,vx,T,lg,lm,lp,L,Ll,Lk,ax,ay,dsys,dsys2,ini,dsn,ss,plt:
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
ini:={seq(icq[i],i=1..n),seq(icqt[i],i=1..n)}:
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 n 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:
#because diff(T,x)=Diff(T,omega)*Diff(omega,x)
Ll[k]:=fx[k-1]-fx[k]=m[k]*ax[k]:      #newton equations
Lk[k]:=fy[k-1]-fy[k]-m[k]*g=m[k]*ay[k]:
end do:
#rearrange the equation system,
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(q[r](t)=theta[r](t),Lk[s]):
Ll[s]:=subs(q[r](t)=theta[r](t),Ll[s]):
L[s]:=subs(omega[r](t)=diff(theta[r](t),t),L[s]):
L[s]:=combine(L[s],trig):
Ll[s]:=combine(Ll[s],trig):
Lk[s]:=combine(Lk[s],trig):
#creat intermedia variables for plotting
plt[s]:=[t,theta[s](t)];
end do:
end do:
#solve the ODE system
dsys:={seq(Lk[i],i=1..n),seq(Ll[i],i=1..n)}:
#get the explicit form of the forces,
ss:=solve(dsys,{seq(fx[i],i=0..n-1),seq(fy[i],i=0..n-1)}):
#replace fx and fy with explicit form, then we get 2n equations(lagrangian equations with angle theta) with 2n variables(theta and D(theta,t))
dsys2:=eval({seq(L[i],i=1..n)},ss):
#solve the ODE system numericlly
dsn := dsolve({op(dsys2)} union ini,numeric,implicit=true):
#plot the angle/ time graph
plots[odeplot](dsn,[seq(plt[i],i=1..n)],0..20):
end proc;

> la(4);

200+              DDDDDDDDDD   DDDDD                                          
   +          DDDD          DDD     DDDDDD                                    
   +   D*****2*C******4*        6        8DDDDD   10       12       14        
  -*************************************************************************-
   +                               C**********BBB    DD                      
-200                                          CCC*****B**BBBBB                
   +                                                  CCC*CCCC***********BBB  
   +                                                      DD             CCC  
-400                                                       DD                
   +                                                        DD                
-600                                                          DD              
   +                                                            DDD          
-800                                                               DDD        
   +                                                                  DDDD    
   +                                                                     DDD  
                                                                              
                                                                              
                                                                              
                      AAAAAAAAAAAAAAAA   theta[1]                            
                      BBBBBBBBBBBBBBBB   theta[2]                            
                      CCCCCCCCCCCCCCCC   theta[3]                            
                      DDDDDDDDDDDDDDDD   theta[4]                            
                                                                              

>

This post generated using the online HTML conversion tool
Download the original worksheet