home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / basic / library / trubasic / chaos / xpendlya.tru < prev    next >
Encoding:
Text File  |  1980-01-01  |  4.2 KB  |  169 lines
  1. !PROGRAM TITLE "XPENDLYA"
  2. LIBRARY "SGLIB.TRC"
  3. DIM A(1),B(1)
  4. CLEAR
  5. PRINT"                      ***PENDULUM - LYAPUNOV EXPONENTS***"
  6. PRINT
  7. PRINT"THIS PROGRAM CALCULATES THE 3 LYAPUNOV EXPONENTS FOR THE DRIVEN PENDULUM"
  8. PRINT"USING THE ALGORITHM OF WOLF ET AL.  THE EXPONENTS ARE SMOOTHED OVER "
  9. PRINT"MANY DRIVE CYCLES (ORBITS)."
  10. !
  11. !N=# OF NONLINEAR EQUTNS., NN= TOTAL # OF EQUATIONS
  12. LET N=3
  13. LET NN=12
  14. DECLARE DEF YPRIM
  15. !
  16. DIM Y(12), ZNORM(3), GSC(3), CUM(3), YNEW(12)
  17. !
  18. !INITIAL CONDITIONS FOR NONLINEAR SYSTEM
  19. LET Y(1)=.5
  20. LET Y(2)=0
  21. LET Y(3)=0
  22. !
  23. !INITIAL CONDITIONS FOR LINEAR SYSTEM (ORTHONORMAL FRAME)
  24. FOR  I= N+1 TO NN
  25.     LET Y(I)=0.0
  26. NEXT I
  27. FOR I=1 TO n
  28.     LET Y((N+1)*I) = 1.0
  29.     LET CUM(I)=0.0
  30. NEXT I
  31. !
  32. INPUT PROMPT"INTEGRATION STEPSIZE (TRY 0.5):":TSTEP
  33. INPUT PROMPT"NUMBER OF ORBITS:":NUMORBITS
  34. INPUT PROMPT"DRIVING FORCE:":G
  35. INPUT PROMPT"DRIVE FREQUENCY:":W
  36. INPUT PROMPT"DAMPING FACTOR:":Q
  37. INPUT PROMPT"LOG BASE (1)NATURAL (2)INF0-2 CHOOSE 1 OR 2:":BASE
  38. !
  39. !SET UP GRAPHICS
  40. IF BASE=1 THEN CALL SETYSCALE(-.7,.3)
  41. IF BASE=2 THEN CALL SETYSCALE(-1,.5)
  42. CALL SETXSCALE(0,NUMORBITS)
  43. CALL SETTEXT("PENDULUM - LYAPUNOV EXPONENTS","# DRIVE CYCLES","LYAP. EXPS.")
  44. CALL SETAXES(0)
  45. CALL RESERVELEGEND
  46. DATA 0,0
  47. CALL DATAGRAPH(A,B,1,0,"WHITE")
  48. CALL GOTOCANVAS
  49. ! 
  50. DO WHILE Y(3)<2*PI*NUMORBITS
  51.    !INITIALIZE INTEGRATOR
  52.    CALL  RKK4(Y(),NN,TSTEP,Q,W,G,YNEW())
  53.    FOR K=1 TO 12
  54.        LET Y(K)=YNEW(K)
  55.    NEXT K
  56.    ! 
  57.    !CONSTRUCT A NEW ORTHONORMAL BASIS BY GRAM-SCHMIDT METHOD
  58.    !NORMALIZE FIRST VECTOR
  59.    LET ZNORM(1)=0.0
  60.    FOR J = 1 TO N
  61.        LET ZNORM(1) = ZNORM(1) +Y(N*J+1)^2
  62.    NEXT J
  63.    LET ZNORM(1)=(ZNORM(1))^.5
  64.    FOR J=1 TO N
  65.        LET Y(N*J+1)=Y(N*J+1)/ZNORM(1)
  66.    NEXT J
  67.    ! 
  68.    !GENERATE THE NEW ORTHONORMAL SET OF VECTORS
  69.    FOR J=2 TO N
  70.        !  GENERATE J-1 COEFFICIENTS
  71.        FOR K = 1 TO (J-1)
  72.            LET GSC(K)=0.0
  73.            FOR L = 1 TO N
  74.                LET GSC(K) = GSC(K) + Y(N*L+J)*Y(N*L+K)
  75.            NEXT L
  76.        NEXT K
  77.        !
  78.        ! CONSTRUCT A NEW VECTOR
  79.        FOR K = 1 TO N
  80.            FOR L= 1 TO J-1
  81.                LET Y(N*K+J) = Y(N*K+J) -GSC(L)*Y(N*K+L)
  82.            NEXT L
  83.        NEXT K
  84.        !
  85.        ! CALCULATE THE VECTOR'S NORM
  86.        LET ZNORM(J) = 0.0
  87.        FOR K= 1 TO N
  88.            LET ZNORM(J) = ZNORM(J)+Y(N*K+J)^2
  89.        NEXT K
  90.        LET ZNORM(J) = (ZNORM(J))^.5
  91.        !
  92.        ! NORMALIZE THE NEW VECTOR
  93.        FOR K = 1 TO N
  94.            LET Y(N*K+J) = Y(N*K+J)/ZNORM(J)
  95.        NEXT K
  96.    NEXT J
  97.    !
  98.    ! UPDATE RUNNING VECTOR MAGNITUDES
  99.    FOR K = 1 TO N
  100.        LET CUM(K) = CUM(K) +LOG(ZNORM(K))
  101.        IF BASE = 2 THEN LET CUM(K) = CUM(K)/LOG(2)
  102.    NEXT K
  103.    !
  104.    !NORMALIZE EXPONENT AND PLOT EXPONENTS
  105.    IF Y(3)>0 THEN
  106.       LET T=Y(3)/W
  107.       FOR K = 1 TO N
  108.           LET CUMKT=CUM(K)/T
  109.           CALL GRAPHPOINT(T*W/(2*PI),CUMKT,1)
  110.       NEXT K
  111.    END IF
  112.  
  113. LOOP
  114. !
  115. CALL ADDLEGEND("G="&STR$(G)&"  Q="&STR$(Q)&"  W="&STR$(W),0,1,"WHITE")
  116. CALL DRAWLEGEND
  117. END
  118. !
  119. SUB RKK4(Y(),NN,TSTEP,Q,W,G,YNEW())
  120.     DIM Y1(12), Y2(12), Y3(12), Y4(12), YY1(12), YY2(12), YY3(12)
  121.     DECLARE DEF YPRIM
  122.     FOR K= 1 TO NN
  123.         LET Y1(K)=TSTEP*YPRIM(Y(),K,Q,W,G)
  124.     NEXT K
  125.     FOR K= 1 TO NN
  126.         LET YY1(K)=Y(K)+Y1(K)/2
  127.     NEXT K
  128.     FOR K=1 TO NN
  129.         LET Y2(K)=TSTEP*YPRIM(YY1(),K,Q,W,G)
  130.     NEXT K
  131.     FOR K = 1 TO NN
  132.         LET YY2(K)=Y(K)+Y2(K)/2
  133.     NEXT K
  134.     FOR K = 1 TO NN
  135.         LET Y3(K)=TSTEP*YPRIM(YY2(),K,Q,W,G)
  136.     NEXT K
  137.     FOR K = 1 TO NN
  138.         LET YY3(K)=Y(K)+Y3(K)
  139.     NEXT K
  140.     FOR K= 1 TO NN
  141.         LET Y4(K)=TSTEP*YPRIM(YY3(),K,Q,W,G)
  142.     NEXT K
  143.     FOR K = 1 TO NN
  144.         LET YNEW(K) = Y(K)+(Y1(K) +2*Y2(K)+2*Y3(K)+Y4(K))/6
  145.     NEXT K
  146. END SUB
  147.  
  148.  
  149.  
  150. !DEFINE DERIVATIVES FUNCTIONS
  151. DEF YPRIM(Y(),K,Q,W,G)
  152.     IF K=1 THEN LET YPRIM=-Y(1)/Q-SIN(Y(2))+G*COS(Y(3))
  153.     IF K=2 THEN LET YPRIM=Y(1)
  154.     IF K=3 THEN LET YPRIM=W
  155.     IF K>3 THEN
  156.        IF K<7 THEN
  157.           LET I = K-4
  158.           LET YPRIM=-Y(4+I)/Q-Y(7+I)*COS(Y(2))+G*Y(10+I)*COS(Y(3))
  159.        END IF
  160.     END IF
  161.     IF K>6 THEN
  162.        IF K<10 THEN
  163.           LET I = K-7
  164.           LET YPRIM=Y(4+I)
  165.        END IF
  166.     END IF
  167.     IF K>9  THEN LET YPRIM = 0
  168. END DEF
  169.