C C 2ND ORDER NON-LINEAR DIFFERENTIAL EQUATIONS C SOLUTION BY RUNGE KUTTA 4TH ORDER METHOD C INTEGER IVECT (400), FIRST, LAST, BGN, END, FCN REAL VEL (2,50), DIST (2,50),AXES (4,30),Z(4) REAL XSCAL(21),LINE(110) DIMENSION GRID(4,30),TIME(50) COMMON/LIMITS/XMIN, XMAX, YMIN, YMAX DATA BLAN, DOT, STAR, PLUS/1H ,1H.,1H*,1H+/ DATA VP1,VP2,VP3/5H(F10.,5H5,5H ,5HVMAX)/ DATA DP1,DP2,DP3/5H(F10.,5H5,5H ,5HDMAX)/ DATA SECT,RESET,YES/1HN,1H ,1HY/ F(V,X,T)=+(D0+D1*T+D2*SIN(W*T))/AM-(B1*V+B2*V*V+C1*X+C2*X*X)/AM 300 CONTINUE REC=RESET DO 885 J=1,50 TIME(J)=0. DO 885 I=1,2 VEL(I,J)=0. 885 DIST(I,J)=0. C C******************************************** C TO INITIALIZE DISPLAY FILES INSERT THE FOLLOWING CALL CLEAR (.FALSE.) C C******************************************** C C.....K IS THE ADDRESS OF INPUT DEVICE C.....L IS THE ADDRESS OF OUTPUT DEVICE K=4 L=4 WRITE(L,800) WRITE(L,801) WRITE(L,802) 800 FORMAT(1X,44HSOLUTION OF 2ND ORDER DIFFERENTIAL EQUATIONS) 801 FORMAT(1X,44HBY RUNGE KUTTA 4TH ORDER METHOD. ENTER THE ) 802 FORMAT(1X,46HCONSTANTS BY TELETYPE, OR DELETE BY KEY RETURN) C ENTER THE X2D0T COEFFICIENT 112 WRITE (L,1) 01 FORMAT (1X,39H COEFFICIENT OF 2ND DERIVATIV, AM,F10.5) READ (K,2) AM IF(AM.LT.0.0001) GO TO 112 2 FORMAT (F10.5) C BM= +B1*V+B2*V*V WRITE (L,10) 10 FORMAT(1X,38H BM= B1*V +B2*V**2, ENTER B1,B2, F10.5) READ (K,2) B1 READ (K,2) B2 C C CM= +C1*X+C2*X*X C ENTRY OF DATA BY TELETYPE TERMINAL C EACH ENTRY IS TERMINATED BY CARRIAGE RETURN WRITE(L,18) 18 FORMAT(1X,38H CM= C1*X +C2*X**2, ENTER C1,C2, F10.5) READ (K,2) C1 READ (K,2) C2 C C DM=D0, + D1*T + D2*SIN(W*T) WRITE (L,8) 8 FORMAT (1X,46H DM=D0+D1*T+D2*SIN(W*T), ENTER D0,D1,D2, F10.5) READ (K,2) D0 READ (K,2) D1 READ (K,2) D2 WRITE(L,19) 19 FORMAT(1X,30HFREQUENCY, W IN RADIANS, F10.5) READ (K,2) W IF (REC.EQ.SECT) GO TO 112 REC=RESET 30 CONTINUE C C ENTER INITIAL CONDITIONS OF VELOCITY, DISPLACEMENT WRITE(L,35) 35 FORMAT (1X,28H INITIAL VELOCITY, , F10.5) READ (K,2) VT0 WRITE(L,40) 40 FORMAT (1X,31H INITIAL DISPLACEMENT , F10.5) READ (K,2) DT0 C C TOTAL TIME FOR SIMULATION AT WHICH PRINTOUT IS DESIRED C THE SCREEN DISPLAY IS SCALED IN BOTH AXES C C TIME INCREMENT AT WHICH PRINTOUT IS DESIRED OUTDEL 707 CONTINUE YMAX=0. YMIN=0. XMAX=0. XMIN=0. VMAX=0. VMIN=0. DMAX=0. DMIN=0. WRITE (L,50) 50 FORMAT (1X,45H TIME INCREMENT OUTDEL A MAX OF 50 INTERVALS) READ (K,2) OUTDEL IF(OUTDEL.LT.0.0001) GO TO 707 747 CONTINUE C C INTEGRATION INTERVAL, IN DECIMAL FRACTION OF OUTDEL ,F10.5 WRITE (L,55) 55 FORMAT (1X,35H INTEGRATION STEP DELT , F10.5 ) READ (K,2) DELT IF(DELT.LT.0.00001) GO TO 747 CYCLE=OUTDEL/DELT IF(CYCLE.LT.20.) GO TO 487 WRITE(L,488) 488 FORMAT(1X,28HINCREMENT DELT IS TOO SMALL.) GO TO 747 487 CONTINUE 727 CONTINUE FINTIM=50.*OUTDEL C IF(FINTIM.LT.0.0001) GO TO 707 C ANSWER NO RETURNS YOU TO DATA INPUT, YES PROCEEDS 15 FORMAT (1X,40H DO YOU WISH TO PROCEED? REPLY YES OR NO) 13 FORMAT (A1) WRITE(L,905)AM,VT0,DT0 905 FORMAT(1X,4H AM=,F10.5,4H VT0,F10.5,4H DT0,F10.5) WRITE(L,906)B1,B2 906 FORMAT(1X,4H B1=,F10.5,2X,4H B2=,F10.5) WRITE(L,907)C1,C2 907 FORMAT(1X,4H C1=,F10.5,2X,4H C2=,F10.5) WRITE(L,908)D0,D1,D2,W 908 FORMAT(1X,4H D0=,F10.5,1X,4HB1= ,F10.5,4H B2=,F10.5,2H W,F10.5) WRITE(L,909)DELT,OUTDEL,FINTIM 909 FORMAT(6H DELT ,F10.5,9H OUTDEL ,F10.5,9H FINTIM ,F10.5) WRITE(L,15) READ(K,13)REC IF(REC.EQ.SECT)GO TO 300 REC=RESET C C C AT THIS POINT THE PROGRAM ENTERS THE OUTER LOOP C WHICH CALLS ON THE RUNGE SUBROUTINE AND RETURNS C THE VEL. & DISPL VALUES T=0.0 VEL(1,1)=T VEL(2,1)=VT0 DIST(1,1)=T DIST(2,1)=DT0 TIME(1)=T V=VT0 X=DT0 C C........................................... C END=50 DO 900 I=2,50 C C.....L IS THE OUTPUT DEVICE I.E. TELETYPE =4,LINE PRINTER =6 IF(T.EQ.0.0) SLOPEV=F(V,X,T) TNET=0. H=DELT 150 TNET=TNET+H TM=T J=0 C C RESET VELOCITY AND DISPLACEMENT VM=V*H XM=X C C......4TH ORDER CONSTANTS ARE Z1,Z2,Z3,Z4 100 CONTINUE J=J+1 IF (J.LE.1) GO TO 105 C THE NEW ACCELERATION VALUES SLOPEV=F(V,X,T) C 105 CONTINUE Z(J)=((H**2.)*(SLOPEV))/2. IF (J-3) 106, 107, 108 C SET NEW VALUES OF T,X,V 106 T=TM+H/2. X=XM+VM/2.+Z(1)/4. V=(VM+Z(J))/H GO TO 100 C 107 T=TM+H X=XM+VM+Z(3) V=(VM+2.*Z(3))/H GO TO 100 C 108 V=VM/H+(Z(1)+2.*Z(2)+2.*Z(3)+Z(4))/(3.*H) X=XM+VM+(Z(1)+Z(2)+Z(3))/3. C C OUTPUT VALUES OF X,V, AND T IF(V.GT.VMAX)VMAX=V IF(V.LT.VMIN)VMIN=V IF(X.GT.DMAX)DMAX=X IF(X.LT.DMIN)DMIN=X IF(T.GT.XMAX)XMAX=T IF(T.LT.XMIN)XMIN=T IF(TNET.GE.OUTDEL) GO TO 115 GO TO 150 115 CONTINUE C.....NOW STORE VALUES OF V,X,T FOR PLOT ROUTINE VEL(1,I)=T VEL(2,I)=V DIST(1,I)=T DIST(2,I)=X TIME(I)=T 900 CONTINUE C C.....THIS TERMINATES THE MAIN LOOP #1 C C.....CHECK MAX Y VALUE AND SCALE CRT DISPLAY IF(VMAX.GE.DMAX)YMAX=VMAX IF(DMAX.GT.VMAX)YMAX=DMAX IF(DMIN.LT.VMIN)YMIN=DMIN IF(VMIN.LT.DMIN)YMIN=VMIN SCALX=(XMAX-XMIN)/1023. SCALY=(YMAX-YMIN)/1023. IF(ABS(YMAX).GT.ABS(YMIN))PERC=ABS(100./YMAX) IF(ABS(YMIN).GT.ABS(YMAX))PERC=ABS(100./YMIN) IF(ABS(VMAX).GT.ABS(VMIN))PERV=100./(ABS(VMAX)) IF(ABS(VMIN).GT.ABS(VMAX))PERV=100./(ABS(VMIN)) IF(ABS(DMAX).GT.ABS(DMIN))PERD=100./(ABS(DMAX)) IF(ABS(DMIN).GT.ABS(DMAX))PERD=100./(ABS(DMIN)) C.....SET UP GRID LINES HORIZONTAL AND VERTICAL GRID(1,1)=XMIN GRID(2,1)=YMAX GRID(3,1)=XMIN GRID(4,1)=YMIN GRID(1,2)=XMIN GRID(3,2)=XMAX GRID(2,2)=0. GRID(4,2)=0. C.....NOW THE DISPLAY IS TURNED ON AND HELD C 737 CONTINUE L=4 C******************************************** CALL VECTOR(3,VEL,1,50,IVECT,1,LAST,7,0,.FALSE.) CALL DISPLY(6,1,IVECT,FIRST,LAST) CALL VECTOR(3,DIST,1,50,IVECT,1,LAST,7,0,.FALSE.) CALL DISPLY(6,2,IVECT,FIRST,LAST) CALL VECTOR(1,GRID,1,2,IVECT,1,LAST,7,0,.FALSE.) CALL DISPLY(6,3,IVECT,FIRST,LAST) CALL TEXT(100,1000,7,1,VP1,VMAX) CALL TEXT(100,0950,7,1,DP1,DMAX) C C.....DISPLAY IS HELD UNTILREPLY TO TELETYPE INPUT C WRITE(L,941) 941 FORMAT(1X,47HCHANGE THE OUTPUT AND INCREMENT TIME? YES OR NO) READ(K,13)REC NEWT=0 IF(REC.EQ.YES)NEWT=1 REC=RESET WRITE(L,205) 205 FORMAT(1X,45H ENTER A NEW SET OF COEFF? REPLY YES OR NO. ) READ(K,13) REC CALL DISPLY(3,1,IVECT,FIRST,LAST) CALL DISPLY(3,2,IVECT,FIRST,LAST) CALL DISPLY(3,3,IVECT,FIRST,LAST) CALL TEXT(100,1000,0,1,VP1,VMAX) CALL TEXT(100,0950,0,1,DP1,DMAX) C******************************************** C IF(REC.EQ.YES)NEWT=-1 REC=RESET IF(NEWT)300,303,707 303 CONTINUE WRITE(L,210) 210 FORMAT(1X,45H DO YOU DESIRE A PERMANENT RECORD? YES OR NO. ) READ(K,13)REC IF(REC.EQ.SECT) STOP REC=RESET C C.....THE PLOT ROUTINE SCALES IN PERCENT AND PRINTS V,X,T) C WRITE(L,686)VMAX,VMIN,DMAX,DMIN 686 FORMAT(1X,4HVMAX,F10.4,2X,4HVMIN,F10.4,2X,4HDMAX,F10.4,2X,4HDMIN, ZF10.4) IF(PERV .GT.0.) GO TO 305 WRITE(L,578) 578 FORMAT(1H1) WRITE(L,310) 310 FORMAT(1X,41H THE VERTICAL SCALING IS ZERO OR NEGATIVE) 305 L=6 WRITE(L,905)AM,VT0,DT0 WRITE(L,906)B1,B2 WRITE(L,907)C1,C2 WRITE(L,908)D0,D1,D2,W WRITE(L,909)DELT,OUTDEL,FINTIM WRITE(L,686)VMAX,VMIN,DMAX,DMIN C....-L=6 IS THE LINE PRINTER XS=-110. DO 5 I=1,21 XS=XS+10. 5 XSCAL(I)=XS WRITE(L,14)(XSCAL(I),I=1,21) 14 FORMAT(1H0,22H TIME VEL DIST ,21F5.0) WRITE(L,21) 21 FORMAT(26X,21(1HT,4X)) DO51 I=1,50 DO 31 K=1,110 31 LINE(K)=BLAN C.....SCALE AND PLOT RESULTS VAL=DIST(2,I)*PERD IF(ABS(VAL)-101.)41,41,111 41 IN=50.+(VAL+1.)/2. VAK=VEL(2,I)*PERV IF(ABS(VAK)-101.)70,70,111 70 IM=50. + (VAK+1.)/2. LINE(50)= DOT LINE(IN)= STAR LINE(IM)= PLUS V=VEL(2,I) X=DIST(2,I) T=VEL(1,I) WRITE(L,80)T,V,X,(LINE(J),J=1,105) 80 FORMAT(F5.2,2F10.4,2X,105A1) 51 CONTINUE 111 WRITE(L,120) 120 FORMAT(16H END OF GRAPH ) GO TO 300 STOP END C C C===================================