home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / basic / library / trubasic / chaos / xpendsic.tru < prev    next >
Encoding:
Text File  |  1980-01-01  |  2.8 KB  |  94 lines
  1. !PROGRAM TITLE - "XPENDSIC"
  2. CLEAR
  3. PRINT"          ***PENDULUM - SENSITIVITY TO INITIAL CONDITIONS***"
  4. PRINT"THIS PROGRAM DISPLAYS THE 2-DIMENSIONAL PHASE DIAGRAM "
  5. PRINT"FOR TWO DIFFERENT SETS OF INTIAL CONDITIONS OF THE PENDULUM."
  6. PRINT"THE OBJECT IS TO ILLUSTRATE SENSITIVITY TO INITIAL CONDITIONS."
  7. !
  8. LIBRARY "SGLIB.TRC"
  9. DECLARE DEF ACCEL
  10. DIM A(1), B(1),XINT(2),VINT(2)
  11. !INPUT STATEMENTS
  12.  INPUT PROMPT"INPUT DRIVING FORCE STRENGTH:":g
  13.  INPUT PROMPT"INPUT DAMPING (IF NO DAMPING THEN INPUT 9999999):":q
  14. FOR I=1 TO 2 
  15.  INPUT PROMPT"INPUT: INITIAL ANGLE , ANGULAR VELOCITY:":XINT(I),VINT(I)
  16. NEXT I
  17.  INPUT PROMPT"INPUT: MINUMUM TIME , MAXIMUM TIME:":TMIN,TMAX
  18. CALL PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)   !SETS MISC AND GRAPH PARAMETERS 
  19. CALL SETXSCALE(XMIN,XMAX)  !FROM SGLIB
  20. CALL SETYSCALE(YMIN,YMAX)  !FROM SGLIB
  21. CALL SETTEXT("PENDULUM - 2-D PHASE DIAGRAM","ANGLE","ANGULAR VELOCITY")
  22. CALL RESERVELEGEND  !FROM SGLIB , SAVES SPACE FOR LEGENDS
  23.  
  24. DATA 0,0
  25. CALL DATAGRAPH(A,B,1,0,"RED")  !FROM SGLIB - PLOTS INITIAL POINT
  26. FOR J=1 TO 2
  27. LET T=0
  28. LET X=XINT(J)
  29. LET V=VINT(J)
  30. CALL GOTOCANVAS   !SETS SCREEN FOR GRAPH
  31. !
  32. !CALCULATION AND GRAPHNG BLOCK
  33. FOR I=1 TO 10000000
  34.   CALL RK4(X,V,TSTEP,XNEW,VNEW,T,W,G,Q)   !CALL RUNGE-KUTTA, STEP = TSTEP
  35.    LET TSHALF=TSTEP/2    ! SPLIT INTERVAL
  36.   CALL RK4(X,V,TSHALF,XNH,VNH,T,W,G,Q)  ! DO TWO HALF STEPS
  37.   CALL RK4(XNH,VNH,TSHALF,XN,VN,T+TSHALF,W,G,Q)
  38.   LET D1=ABS(XN-XNEW)
  39.   LET D2=ABS(VN-VNEW)
  40.   LET DELTA=MAX(D1,D2)
  41.   IF DELTA<EPS THEN
  42.     IF T>TMIN THEN
  43. !     IF ABS(X)>PI THEN LET X=X-2*PI*ABS(X)/X
  44.      IF J=1 THEN CALL GRAPHPOINT(X,V,3)
  45.      IF J=2 THEN CALL GRAPHPOINT(X,V,10)
  46.     END IF
  47.      LET X=XNEW
  48.      LET V=VNEW
  49.      LET T=T+TSTEP
  50.      LET TSTEP=TSTEP*.95*(EPS/DELTA)^.25
  51. !     IF ABS(X)>PI THEN LET X=X-2*PI*ABS(X)/X
  52.   ELSE
  53.      LET TSTEP=TSTEP*.95*(EPS/DELTA)^.2   !REDUCE STEP SIZE
  54.   END IF
  55.   IF T>TMAX THEN LET I=10000001
  56.   NEXT I
  57. NEXT J
  58.   CALL ADDLEGEND("G="&STR$(G)&"   Q="&STR$(Q),0,1,"WHITE")  
  59.   CALL DRAWLEGEND    !ADDS G AND Q VALUES TO LEGEND
  60. get key variable
  61. clear
  62. print"press <esc> key to finish"
  63. END 
  64. !
  65. SUB RK4(X,V,TSTEP,XNEW,VNEW,T,W,G,Q)   !RUNGE-KUTTA INTEGRATOR
  66.    DECLARE DEF ACCEL
  67.    LET XK1=TSTEP*V
  68.    LET VK1=TSTEP*ACCEL(X,V,T,W,G,Q)
  69.    LET XK2=TSTEP*(V+VK1/2)
  70.    LET VK2=TSTEP*ACCEL(X+XK1/2,V+VK1/2,T+TSTEP/2,W,G,Q)
  71.    LET XK3=TSTEP*(V+VK2/2)
  72.    LET VK3=TSTEP*ACCEL(X+XK2/2,V+VK2/2,T+TSTEP/2,W,G,Q)
  73.    LET XK4=TSTEP*(V+VK3)
  74.    LET VK4=TSTEP*ACCEL(X+XK3,V+VK3,T+TSTEP,W,G,Q)
  75.    LET VNEW=V+(VK1+2*VK2+2*VK3+VK4)/6
  76.    LET XNEW=X+(XK1+2*XK2+2*XK3+XK4)/6
  77. END SUB
  78. !
  79. DEF ACCEL(X,V,T,W,G,Q)
  80.    LET DAMP=1/Q
  81.    LET ACCEL=-SIN(X)-DAMP*V+G*COS(W*T)
  82. END DEF
  83. !
  84. SUB PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
  85.  LET W=0.66666666
  86.  LET EPS=1.0E-6
  87.  LET TSTEP=0.5
  88.  LET XMIN=-6
  89.  LET XMAX=6
  90.  LET YMIN=-4
  91.  LET YMAX=4
  92. END SUB
  93.  
  94.