home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / basic / library / trubasic / chaos / xbifurca.tru < prev    next >
Encoding:
Text File  |  1980-01-01  |  4.7 KB  |  141 lines
  1. !PROGRAM TITLE "XBIFURCA"
  2. CLEAR
  3. PRINT"              ***PENDULUM - BIFURCATION DIAGRAM***"
  4. PRINT"THIS PROGRAM DISPLAYS THE BIFURCATION DIAGRAM FOR THE PENDULUM."
  5. PRINT"FOR EACH VALUE OF FORCING AMPLITUDE, G, THE SYSTEM COMES TO A "
  6. PRINT"STEADY STATE (AFTER MIN.TIME) AND THEN THE ANGULAR VELOCITY AT"
  7. PRINT"AT THE BEGINNING OF EACH FORCING CYCLE IS DISPLAYED FOR A NUMBER OF"
  8. PRINT"FURTHER CYCLES (GOVERNED BY MAX. TIME).  THE DATA CAN BE SAVED TO A FILE"
  9. PRINT
  10.  
  11. DIM XINT(10), VINT(10)
  12. LIBRARY "SGLIB.TRC"
  13. !
  14. DECLARE DEF accel
  15. DIM A(1),B(1)
  16. INPUT prompt"Input LOWEST DRIVING FORCE STRENGTH: ":GMIN
  17. INPUT PROMPT"INPUT HIGHEST DRIVING FORCE STRENGTH:":GMAX
  18. INPUT PROMPT"INPUT G STEPSIZE:":DELTAG
  19. INPUT PROMPT"INPUT NUMBER OF SETS OF INITIAL CONDITIONS:":NUMSETS
  20. FOR I=1 TO NUMSETS
  21.     INPUT PROMPT"INPUT INITIAL ANGLE:":XINT(i)
  22.     INPUT PROMPT"INPUT INITIAL ANGULAR VELOCITY:":VINT(I)
  23. NEXT I
  24. INPUT prompt"Input damping (If no damping then input 9999999):":q
  25. INPUT Prompt"Input min. and max. time:":tmin,tmax
  26. INPUT prompt"Input phase angle/(2*pi) - USE ZERO IF WANT BEGINNING OF CYCLE:":PHI
  27. PRINT
  28. PRINT"SINCE THE RUNTIME IS VERY LONG THE NEXT SET OF INPUTS GIVE AN OPTION"
  29. PRINT"TO SAVE THE DATA TO A FILE"
  30. INPUT PROMPT"SAVE TO A FILE? YES(1), NO(2):":SV
  31. IF SV=1 THEN
  32.    PRINT "A REASONABLE 8 CHARACTER FILE NAME (USE A NUMBER) MIGHT INCLUDE"
  33.    PRINT"1)FIRST 2 DIGITS FOR Q VALUE"
  34.    PRINT"2)NEXT 3 DIGITS FOR LOWEST G VALUE"
  35.    PRINT"3)LAST 3 DIGITS FOR HIGHES G VALUE"
  36.    PRINT" EXAMPLE  20145150"
  37.    INPUT PROMPT"FILE NAME =>":FILENAME
  38.    INPUT PROMPT"DATA FILE DRIVE (A/B/C/D):":B$
  39.    LET NAME$=STR$(FILENAME)
  40. END IF
  41. CLEAR
  42. CALL PARAMS(W,EPS,TSTEP)
  43. CALL SETAXES(0)
  44. CALL SETXSCALE(GMIN,GMAX)
  45. CALL SETYSCALE(-1,3)
  46. CALL SETTEXT("PENDULUM BIFURCATION DIAGRAM","FORCING-G","ANGULAR VELOCITY")
  47. CALL RESERVELEGEND
  48. !
  49. DATA O,O
  50. CALL DATAGRAPH(A,B,1,O,"WHITE")
  51. !
  52. IF SV=1 THEN                      !OPENS A FILE AND SPECIFIES CHARACTERISTICS
  53.    OPEN #1:NAME B$&":"&NAME$,ORGANIZATION RECORD, CREATE NEWOLD
  54.    ASK #1:FILESIZE LENGTH
  55.    IF LENGTH=0 THEN SET#1: RECSIZE 10
  56.    SET #1: POINTER END
  57. END IF
  58. !
  59. FOR II=1 TO NUMSETS               !LOOPS FOR ALL INITIAL CONDITIONS
  60.     LET T=0
  61.     LET XP =XINT(ii)
  62.     LET VP= VINT(II)
  63.     FOR G=GMIN TO GMAX STEP DELTAG     !LOOPS FOR ALL G VALUES
  64.         LET t=0
  65.         LET x=xP
  66.         LET v=vP
  67.         CALL GOTOCANVAS
  68.         !
  69.         !CALCULATION AND GRAPHING BLOCK
  70.         LET phi=phi*2*pi
  71.         FOR i=1 to 1000000
  72.             CALL rk4(x,v,tstep,xnew,vnew,t,w,g,q)     ! Take a step of size tstep
  73.             LET tshalf=tstep/2
  74.             CALL rk4(x,v,tshalf,xnh,vnh,t,w,g,q)      !Take two half steps
  75.  
  76.             CALL rk4(xnh,vnh,tshalf,xn,vn,t+tshalf,w,g,q)
  77.             LET d1=abs(xn-xnew)
  78.             LET d2=abs(vn-vnew)
  79.             LET delta=max(d1,d2)
  80.             IF delta<eps then
  81.                IF t>tmin then
  82.                   LET tnew=t+tstep
  83.                   LET w1=mod(phi-w*t,2*pi)  !Check for Poincare section
  84.                   LET w2=mod(w*tnew-phi,2*pi)
  85.                   IF w1<w*tstep then
  86.                      IF w2<w*tstep then
  87.                         LET ts=w1/w
  88.                         CALL rk4(x,v,ts,xp,vp,t,w,g,q)     !CALCULATES POINT AT SECTION
  89.                         IF abs(xp)>pi then LET xp=xp-2*pi*abs(xp)/xp
  90.                         CALL GRAPHPOINT(G,VP,1)
  91.                         IF SV=1 THEN WRITE #1:G,VP
  92.                      END IF
  93.                   END IF
  94.                END IF
  95.                LET x=xnew
  96.                LET v=vnew
  97.                LET t=t+tstep      !Expand step size
  98.                LET tstep=tstep*.95*(eps/delta)^.25
  99.                IF abs(x)>pi then  !bring theta back into range
  100.                   LET x=x-2*pi*abs(x)/x
  101.                END IF
  102.             ELSE                  !else reduce step size
  103.                LET tstep=tstep*.95*(eps/delta)^.2
  104.             END IF
  105.             IF t>tmax then LET i=1000001
  106.         NEXT i
  107.     NEXT G
  108. NEXT II
  109. LET G$=STR$(G)
  110. LET Q$=STR$(Q)
  111. CALL ADDLEGEND(" Q="&STR$(Q)&"   PHI="&STR$(PHI),0,1,"WHITE")
  112. CALL DRAWLEGEND
  113. get key variable
  114. clear
  115. print"press <esc> key to finish"
  116. END
  117. !
  118. !
  119. SUB rk4(x,v,tstep,xnew,vnew,t,w,g,q)
  120.     DECLARE DEF accel
  121.     LET xk1=tstep*v
  122.     LET vk1=tstep*accel(x,v,t,w,g,q)
  123.     LET xk2=tstep*(v+vk1/2)
  124.     LET vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
  125.     LET xk3=tstep*(v+vk2/2)
  126.     LET vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
  127.     LET xk4=tstep*(v+vk3)
  128.     LET vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
  129.     LET vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
  130.     LET xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
  131. END SUB
  132. DEF accel(x,v,t,w,g,q)
  133.     LET accel=-sin(x)-(1/q)*v+g*cos(w*t)
  134. END DEF
  135. !
  136. SUB PARAMS(W,EPS,TSTEP)
  137.     LET W=0.66666666
  138.     LET EPS=1.0E-6
  139.     LET TSTEP=0.5
  140. END SUB
  141.