home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / basic / library / trubasic / chaos / xphase3d.tru < prev    next >
Encoding:
Text File  |  1980-01-01  |  2.8 KB  |  94 lines
  1. !PROGRAM PEND/PHASE3D
  2. LIBRARY "3dlib.trc"
  3. CLEAR
  4. PRINT"                         ***PENDULUM - 3D PHASE DIAGRAM***"
  5. PRINT
  6. PRINT"THIS PROGRAM SHOWS THE COMPLETE 3D PHASE SPACE OF THE PENDULUM.  THE"
  7. PRINT"TIME COORD IS 3/2 OF THE CYCLIC PHI COORD, AND THEREFORE THE BOUNDARY"
  8. PRINT"CONDITIONS ON TIME ARE ALSO CYCLIC."
  9. !
  10. declare def accel
  11.  
  12. input prompt"Input driving force strength: ":g
  13. input prompt"Input damping :":q
  14. input prompt"Input initial position: ": xint
  15. input prompt"Input initial angular velocity:":v
  16. input Prompt"Input min. and max. time:":tmin,tmax
  17. input prompt"Input 'cyclic' time step as fraction of phi: ":phistep
  18.  
  19. let w=0.6666667
  20. let eps=1.0e-8
  21. let tstep=0.001
  22. let t=0
  23. let x=xint
  24.  
  25.  
  26. Call Parawindow(0,3*pi,-pi,pi,-3,3,s$)
  27. Set mode "hires"
  28. Call Ticks3(1,1,1,s$)
  29. CALL PLOTTEXT3(0,2,5,"PENDULUM - 3D PHASE DIAGRAM",S$)
  30. CALL PLOTTEXT3(0,5,2,"G="&STR$(G)&" Q=2",S$)
  31.    let phistep=2*pi*phistep
  32.  for i=1 to 1000000
  33.      call rk4(x,v,tstep,xnew,vnew,t,w,g,q)! Take a step of size tstep
  34.      let tshalf=tstep/2
  35.      call rk4(x,v,tshalf,xnh,vnh,t,w,g,q)!Take two half steps
  36.      call rk4(xnh,vnh,tshalf,xn,vn,t+tshalf,w,g,q)
  37.      let d1=abs(xn-xnew)
  38.      let d2=abs(vn-vnew)
  39.      let delta=max(d1,d2)
  40.      if delta<eps then
  41.        if t>tmin then
  42.           let tnew=t+tstep
  43.           let w1=mod(-w*t,phistep)  !Check for time axis step
  44.           let w2=mod(w*tnew,phistep)
  45.           if w1<w*tstep then
  46.           if w2<w*tstep then
  47.             let ts=w1/w
  48.             call rk4(x,v,ts,xp,vp,t,w,g,q)
  49.             if abs(xp)>pi then let xp=xp-2*pi*abs(xp)/xp
  50.             let period=2*pi/w
  51.             let tyme=mod(t+ts,period)
  52.             Call PlotOff3(tyme,xp,vp,s$)
  53.            end if
  54.           end if
  55.        end if
  56.      let x=xnew
  57.      let v=vnew 
  58.      let t=t+tstep              !Expand step size
  59.      let tstep=tstep*.95*(eps/delta)^.25
  60.         if abs(x)>pi then      !bring theta back into range
  61.           let x=x-2*pi*abs(x)/x
  62.         end if
  63.      else       !else reduce step size
  64.        let tstep=tstep*.95*(eps/delta)^.2
  65.      end if
  66.      if t>tmax then let i=1000001
  67. next i
  68. call plottext3(10,0,0,"TIME",s$)
  69. call plottext3(0,4,0,"ANGLE",s$)
  70. call plottext3(0,-1,4,"ANG. VEL.",s$)
  71. end
  72. !
  73. !
  74. sub rk4(x,v,tstep,xnew,vnew,t,w,g,q)
  75. declare def accel
  76.     let xk1=tstep*v
  77.     let vk1=tstep*accel(x,v,t,w,g,q)
  78.     let xk2=tstep*(v+vk1/2)
  79.     let vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
  80.     let xk3=tstep*(v+vk2/2)
  81.     let vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
  82.     let xk4=tstep*(v+vk3)
  83.     let vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
  84.     let vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
  85.     let xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
  86. end sub
  87. !
  88. !
  89. def accel(x,v,t,w,g,q)
  90.  let accel=-sin(x)-(1/q)*v+g*cos(w*t)
  91. end def
  92.  
  93.             
  94.