*APP PHYS01 SUBROUTINE INITAL C... C... TRANSIENT RESPONSE OF A RESONANT RCL CIRCUIT C... C... THE EQUATIONS WHICH DEFINE THE VOLTAGE OF THE CIRCUIT AS A FUNC- C... TION OF TIME ARE C... C... L*C*D2V/DT2 + R*C*DV/DT + V = 0 (1) C... C... V(0) = V0 (2) C... C... DV(0)/DT = 0 (3) C... WHERE C... C... L INDUCTANCE = 0.5 HENRIES C... C... C CAPACITANCE = 2.0E-06 FARADS C... C... R RESISTANCE = 100 OHMS C... C... V0 INITIAL VOLTAGE = 10 VOLTS C... C... THE NUMERICAL SOLUTION TO EQUATIONS (1) TO (3) IS PLOTTED AND C... AND COMPARED WITH THE ANALYTICAL SOLUTION IN SUBROUTINE PRINT. C... COMMON/T/T,NFIN,NRUN/Y/V1,V2/F/DV1DT,DV2DT/PARM/L,C,R,V0 REAL L C... C... THE CIRCUIT PARAMETERS ARE SET IN A BLOCK DATA ROUTINE AND PASSED C... THROUGH LABELLED COMMON/PARM/ C... C... INITIAL CONDITIONS C... V(0) = V0, DV(0)/DT = 0 V1=V0 V2=0. RETURN END SUBROUTINE DERV COMMON/T/T,NFIN,NRUN/Y/V1,V2/F/DV1DT,DV2DT/PARM/L,C,R,V0 REAL L C... C... CIRCUIT DIFFERENTIAL EQUATION DV2DT=-1./(L*C)*(R*C*V2+V1) DV1DT=V2 RETURN END BLOCK DATA COMMON/T/T,NFIN,NRUN/Y/V1,V2/F/DV1DT,DV2DT/PARM/L,C,R,V0 REAL L C... C... THE CIRCUIT PARAMETERS ARE SET IN THIS ROUTINE DATA L, C, R, V0/ 1 0.5, 2.0E-06, 100., 10./ END SUBROUTINE PRINT(NI,NO) COMMON/T/T,NFIN,NRUN/Y/V1,V2/F/DV1DT,DV2DT/PARM/L,C,R,V0 REAL L C... C... DIMENSION THE ARRAYS FOR THE SOLUTION TO BE PLOTTED DIMENSION TPLOT(51), V1PLOT(51), V2PLOT(51) C... C... INITIALIZE A COUNTER FOR THE PLOTTED SOLUTION DATA IP/0/ C... C... PRINT A HEADING IF(IP.EQ.0)WRITE(NO,2)L,C,R 2 FORMAT(23H RESONANT CIRCUIT, L = ,E9.2,3X,4HC = ,E9.2, 1 3X,4HR = ,E9.2,//,21X,1HT,11X,4HV(T),5X,10HV(T)(ANAL)) C... C... EVALUATE THE ANALYTICAL SOLUTION 1 AS=1./(L*C)-R**2/(4.*L**2) A=SQRT(ABS(AS)) IF(AS)3,4,5 C... C... UNDERDAMPED SOLUTION 5 VA=V0/(A*SQRT(C*L))*EXP(-R*T/(2.*L))*COS(A*T-ATAN(R/(2.*L*A))) GO TO 6 C... C... CRITICALLY DAMPED SOLUTION 4 VA=V0*EXP(-R*T/(2.*L))*(1.+R*T/(2.*L)) GO TO 6 C... C... OVERDAMPED SOLUTION 3 VA=V0*EXP(-R*T/(2.*L))*((0.5+R/(4.*L*A))*EXP(A*T) 1 + (0.5-R/(4.*L*A))*EXP(-A*T)) C... C... PRINT THE NUMERICAL, ANALYTICAL SOLUTIONS 6 WRITE(NO,7)T,V1,VA 7 FORMAT(7X,E15.2,2E15.3) C... C... STORE THE NUMERICAL SOLUTION FOR SUBSEQUENT PLOTTING IP=IP+1 TPLOT(IP)=T V1PLOT(IP)=V1 V2PLOT(IP)=V2 C... C... PLOT THE NUMERICAL SOLUTION AS V(T) VS T IF(IP.LT.51)RETURN CALL TPLOTS(1,IP,TPLOT,V1PLOT) WRITE(NO,9)L,C,R,A 9 FORMAT(/,16X,16HRESONANT CIRCUIT, 1 //,16X,4HL = ,E9.2,3X,4HC = ,E9.2, 2 3X,4HR = ,E9.2,3X,4HA = ,E9.2) C... C... PLOT THE NUMERICAL SOLUTION AS DV(T)/DT VS V(T) CALL TPLOTS(1,IP,V1PLOT,V2PLOT) WRITE(NO,9)L,C,R,A C... C... ****************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND MUST BE MODIFIED FOR OTHER COMPUTERS OR CONVERTED TO COMMENTS C... C... CALL QIKSAX(3,3) C... CALL QIKPLT(TPLOT,V1PLOT,IP,3H*T*,6H*V(T)*,2H**,1) C... CALL QIKSAX(3,3) C... CALL QIKPLT(V1PLOT,V2PLOT,IP,6H*V(T)*,10H*DV(T)/DT*,2H**,1) C... C... ****************************************************************** C... IP=0 RETURN END CARNAHAN, ET AL, APPLIED NUMERICAL METHODS, EX. 6.3, PP. 367-380 0. 0.025 0.0005 2 100 4 1 ABS 0.01 END OF RUNS *APP PHYS02 SUBROUTINE INITAL C... C... SOLUTION OF THE STEADY STATE AND TRANSIENT SINE-GORDON EQUATIONS C... BY QUASILINEARIZATION AND THE METHOD OF LINES C... C... THE SINE-GORDON EQUATION HAS SEVERAL APPLICATIONS IN PHYSICS, C... INCLUDING THE MODELING OF A JOSEPHSON JUNCTION. THE STARTING C... POINT FOR THE ANALYSIS IS THE PARTIAL DIFFERENTIAL EQUATION C... (PDE) (1): C... C... B**PHI + PHI - PHI - A*PHI = SIN(PHI) (1) C... XXT XX TT T C... C... WHERE C... C... PHI TRANSVERSE MAGNETIC FLUX C... C... X NORMALIZED DISTANCE ALONG THE JUNCTION C... C... T TIME C... C... A,B CONSTANTS WHICH INTRODUCE DISSIPATION INTO THE MODEL C... C... (1) SCOTT, ALWYN C., FLORA Y. F. CHU AND STANLEY A. REIBLE, C... MAGNETIC FLUX PROPAGATION ON A JOSEPHSON TRANSMISSION LINE, C... MATHEMATICS RESEARCH CENTER REPORT 1615, NIVERSITY OF C... WISCONSIN, MADISON, WI C... C... A LETTER SUBSCRIPT DENOTES A PARTIAL DERIVATIVE, E.G., PHI IS C... T C... THE PARTIAL DERIVATIVE OF PHI WITH RESPECT TO T, PHI IS THE C... XXT C... THIRD ORDER PARTIAL DERIVATIVE WITH RESPECT TO X (TWICE) AND T C... (ONCE). C... C... FOR THE NONDISSIPATIVE CASE, A = B = 0, AND EQUATION (1) REDUCES C... TO THE SINE-GORDON EQUATION C... C... PHI - PHI = SIN(PHI) (2) C... XX TT C... C... IF WE NOW CONSIDER THE SPECIAL STEADY STATE CASE OF EQUATION (2) C... C... PHI = SIN(PHI) (3) C... XX C... C... WE START THE ANALYSIS BY CONSIDERING THE SOLUTION TO EQUATION (3) C... SUBJECT TO THE BOUNDARY CONDITIONS C... C... PHI(0) = 1, PHI (1) = 0 (4)(5) C... X C... C... EQUATION (3) IS A SECOND ORDER, BOUNDARY VALUE, NONLINEAR ORDINARY C... DIFFERENTIAL EQUATION. C... C... IF THE DERIVATIVE IN EQUATION (2) IS REPLACED BY THE SECOND ORDER C... CENTERED APPROXIMATION C... C... PHI = (PHI(I+1) - 2*PHI(I) + PHI(I-1))/(DX**2) + O(DX**2) C... XX C... C... WHERE DX IS THE GRID SPACING IN X, AND I IS THE GRID INDEX, EQUA- C... TION (3) CAN BE APPROXIMATED BY A SET OF NONLINEAR EQUATIONS C... WRITTEN FOR THE GRID POINTS I = 0, 1, 2,..., N C... C... 1 - (2*PHI(1) + (DX**2)*SIN(PHI(1))) + PHI(2) = 0 C... C... PHI(1) - (2*PHI(2) + (DX**2)*SIN(PHI(2))) + PHI(3) = 0 C... . . C... . . (6) C... . . C... C... PHI(N-2) - (2*PHI(N-1) + (DX**2)*SIN(PHI(N-1))) + PHI(N) = 0 C... C... 2*PHI(N-1) - (2*PHI(N) + (DX**2)*SIN(PHI(N))) = 0 C... C... NOTE THAT BOUNDARY CONDITIONS (4) AND (5) HAVE BEEN INCLUDED IN C... THE FIRST AND LAST EQUATIONS FOR GRID POINTS I = 1 AND N. FOR C... BOUNDARY CONDITION (5), THE SECOND ORDER CENTERED APPROXIMATION C... C... PHI (1) = (PHI(N+1) - PHI(N-1))/(2*DX) = 0 C... X C... C... OR C... C... PHI(N+1) = PHI(N-1) C... C... WAS USED. C... C... (6) IS A SET OF NONLINEAR, TRIDIAGONAL TRANSCENDENTAL EQUATIONS. C... WE WILL COMPUTE A SOLUTION BY NEWTON'S METHOD, WHICH IS ALSO C... CALLED QUASILINEARIZATION. IF THE NONLINEAR SIN TERM IS EX- C... PANDED IN A TAYLOR SERIES C... C... K+1 K C... SIN(PHI(I) ) = SIN(PHI(I) ) + C... C... K K+1 K C... + COS(PHI(I) )*(PHI(I) - PHI(I) ) C... C... K+1 C... = A(I) + B(I)*PHI(I) + C(I) C... C... WHERE C... C... K ITERATION COUNTER C... C... K C... A(I) SIN(PHI(I) ) C... C... K C... B(I) COS(PHI(I) ) C... C... K K C... C(I) -COS(PHI(I) )*PHI(I) C... C... THEN EQUATIONS (6) CAN BE WRITTEN IN APPROXIMATE LINEAR FORM AS C... (WITH THE ITERATION COUNTER K DROPPED) C... C... - (2 + (DX**2)*B(1))*PHI(1) + PHI(2) = C... C... - 1 + (DX**2)*(A(1) + C(1)) C... C... PHI(1) - (2 + (DX**2)*B(2))*PHI(2) + PHI(3) = C... C... (DX**2)*(A(2) + C(2)) C... . . C... . . (7) C... . . C... C... PHI(N-2) - (2 + (DX**2)*B(N-1))*PHI(N-1) + PHI(N) = C... C... (DX**2)*(A(N-1) + C(N-1)) C... C... 2*PHI(N-1) - (2 + (DX**2)*B(N))*PHI(N) = C... C... (DX**2)*(A(N) + C(N)) C... C... NOTE THAT A, B AND C IN EACH EQUATION ARE FUNCTIONS OF I, AND C... THEREFORE MUST BE EVALUATED FOR EACH EQUATION. C... C... (7) IS A SET OF LINEAR, TRIDIAGONAL ALGEBRAIC EQUATIONS WHICH C... CAN BE SOLVED BY A BANDED SOLVER. IT WILL BE NECESSARY TO ASSUME C... A TRIAL SOLUTION FOR ITERATION K = 0. THEN (7) CAN BE SOLVED C... FOR THE K = 1 SOLUTION, AND THIS SOLUTION CAN BE USED IN PLACE OF C... THE ASSUMED K = 0 SOLUTION TO COMPUTE THE K = 2 SOLUTION, ETC., C... UNTIL CONVERGENCE IS ACHIEVED. C... C... TWO SOLUTIONS TO EQUATION (3) ARE COMPUTED BY THIS METHOD FOR C... N = 10 AND N = 20. SOLUTIONS ARE ALSO COMPUTED BY IMSL ROUTINE C... ZSPOW WHICH CAN ACCEPT THE NONLINEAR EQUATIONS (EQUATIONS (6)) C... DIRECTLY. C... C... ANOTHER APPROACH TO THE SOLUTION OF EQUATION (3) IS ALSO PRO- C... GRAMMED. IF A TIME DERIVATIVE IS ADDED TO EQUATION (3) TO FORM C... AN UNSTEADY STATE OR TRANSIENT FORM OF THE SINE-GORDON EQUATION, C... C... PHI = PHI - SIN(PHI) (8) C... T XX C... C... EQUATION (8) CAN THEN BE INTEGRATED UNTIL PHI ---+ 0, SO THAT C... T C... A SOLUTION TO EQUATION (3) RESULTS. THIS APPROACH IS CALLED C... THE METHOD OF FALSE TRANSIENTS. C... C... ONCE THIS METHOD IS UNDERSTOOD, IT CAN BE EXTENDED TO THE SOLU- C... TION OF EQUATION (1) WITH ANY COMBINATION OF DERIVATIVES IN X C... AND T. C... C... THE CODING FOR THESE VARIOUS METHODS FOLLOWS. THE SOLUTIONS C... TO EQUATION (3) CAN THEN BE COMPARED, AND THE TRANSIENT SOLUTION C... TO EQUATION (8) CAN BE OBSERVED. C... C... ACCESS THE INPUT/OUTPUT UNIT NUMBERS COMMON/IO/ NI, NO C... C... DSS/2 COMMON AREA COMMON/T/ T, NSTOP, NORUN 1 /Y/ PHIU(21) 2 /F/ PHIT(21) 3 /S/ PHIX(21),PHIXX(21) 4 /C/ N, IP C... C... COMMON AREA FOR ACCESS IN FUNCTIONS AND SUBROUTINES USED BY C... LEQT1B COMMON/C1/ PHI0(20) C... C... COMMON/PAR/ CONTAINS THE PARAMETERS USED IN SUBROUTINE FCN C... CALLED BY ZSPOW COMMON/PAR/ DXS C... C... DIMENSION THE ARRAYS REQUIRED BY SUBROUTINE LEQT1B DIMENSION ABSM(20,3), BRHS(20), XL(40), X(20), II(20) C... C... DIMENSION THE ARRAYS REQUIRED BY SUBROUTINE ZSPOW DIMENSION PHI(20), F(20), WK(750), PAR(1) C... C... EXTERNAL THE ROUTINE CALLED BY ZSPOW TO EVALUATE THE FUNCTIONS EXTERNAL FCN C... C... FOR THE FIRST RUN (NORUN = 1), SOLUTIONS TO EQUATIONS (3) ARE C... COMPUTED BY IMSL ROUTINES LEQT1F AND ZSPOW IF(NORUN.GT.1)GO TO 103 C... C... STEP THROUGH TWO CASES DO 7 NCASE=1,2 C... C... FOR THE FIRST CASE, N = 10 IF(NCASE.EQ.1)N=10 C... C... FOR THE SECOND CASE, N = 20 IF(NCASE.EQ.2)N=20 C... C... ***************************************************************** C... C... COMPUTE A SOLUTION TO EQUATION (3) BY IMSL ROUTINE LEQT1B C... BASED ON QUASILINEARIZATION C... C... GRID SPACING AND ITS SQUARE DX=1./FLOAT(N) DXS=DX*DX C... C... DEFINE THE INITIAL TRIAL SOLUTION, SPATIAL GRID AND INDEX DO 6 I=1,N PHI0(I)=0.5 X(I)=DX*FLOAT(I) II(I)=I 6 CONTINUE C... C... STEP THROUGH A SERIES OF ITERATIONS DO 8 ITER=1,20 C... C... SET UP THE COEFFICIENT MATRIX IN BAND STORAGE MODE (STORED IN C... ABSM) - - - C... C... LOWER DIAGONAL DO 1 I=1,N IF(I.EQ.1)ABSM(1,1)=0. IF((I.GT.1).AND.(I.LT.N))ABSM(I,1)=1.0 IF(I.EQ.N)ABSM(N,1)=2.0 1 CONTINUE C... C... MAIN DIAGONAL DO 2 I=1,N ABSM(I,2)=-(2.0+DXS*B(I)) 2 CONTINUE C... C... UPPER DIAGONAL DO 3 I=1,N IF(I.NE.N)ABSM(I,3)=1.0 IF(I.EQ.N)ABSM(N,3)=0. 3 CONTINUE C... C... RIGHT HAND SIDE VECTOR DO 4 I=1,N IF(I.EQ.1)BRHS(1)=-1.0+DXS*(A(I)+C(I)) IF(I.NE.1)BRHS(I)= DXS*(A(I)+C(I)) 4 CONTINUE C... C... INPUT PARAMETERS FOR SUBROUTINE LEQT1B NLC=1 NUC=1 IA=20 M=1 IB=20 IJOB=0 C... C... SOLVE THE LINEAR ALGEBRAIC EQUATIONS BY IMSL SUBROUTINE LEQT1B CALL LEQT1B(ABSM,N,NLC,NUC,IA,BRHS,M,IB,IJOB,XL,IER) C... C... CHECK FOR CONVERGENCE DO 9 I=1,N IF(ABS(BRHS(I)-PHI0(I)).GT.0.001)GO TO 10 9 CONTINUE C... C... CONVERGENCE HAS BEEN ACHIEVED, SO PRINT THE SOLUTION VECTOR WRITE(NO,5)N,ITER,IER,(II(I),X(I),BRHS(I),I=1,N) 5 FORMAT(1H ,//, 1 ' N = ',I5,3X,' ITER = ',I5,3X,' IER = ',I5,//, 2 9X,'I',9X,'X',9X,'PHI(X)',/,(I10,F10.3,F15.4)) GO TO 12 C... C... CONVERGENCE HAS NOT BEEN ACHIEVED, SO CONTINUE THE ITERATIONS 10 DO 11 I=1,N PHI0(I)=BRHS(I) 11 CONTINUE C... C... GO ON TO THE NEXT ITERATION 8 CONTINUE C... C... ***************************************************************** C... C... COMPUTE A SOLUTION TO EQUATION (3) BY IMSL ROUTINE ZSPOW C... C... DEFINE THE PARAMETERS FOR ZSPOW 12 NSIG=4 ITMAX=200 C... C... SET THE INITIAL ESTIMATE OF THE SOLUTION DO 100 I=1,N PHI(I)=0.5 100 CONTINUE C... C... COMPUTE THE SOLUTION TO EQUATION (6) CALL ZSPOW(FCN,NSIG,N,ITMAX,PAR,PHI,FNORM,WK,IER) C... C... PRINT THE SOLUTION WRITE(NO,101)N,FNORM,IER,(II(I),X(I),PHI(I),I=1,N) 101 FORMAT(1H ,//, 1 ' N = ',I5,3X,' FNORM = ',E11.3,3X,' IER = ',I5,//, 2 9X,'I',9X,'X',9X,'PHI(X)',/,(I10,F10.3,F15.4)) C... C... GO ON TO THE NEXT CASE 7 CONTINUE C... C... ***************************************************************** C... C... COMPUTE A SOLUTION TO EQUATION (8) BY THE METHOD OF FALSE C... TRANSIENTS C... C... INITIAL CONDITION 103 N=21 DO 102 I=1,N PHIU(I)=0.0 102 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END FUNCTION A(I) C... C... FUNCTION A COMPUTES THE FIRST TERM IN THE LINEAR APPROXIMATION C... OF SIN(PHI) AS DISCUSSED PREVIOUSLY C... COMMON/C1/ PHI0(20) A=SIN(PHI0(I)) RETURN END FUNCTION B(I) C... C... FUNCTION B COMPUTES THE SECOND TERM IN THE LINEAR APPROXIMATION C... OF SIN(PHI) AS DISCUSSED PREVIOUSLY C... COMMON/C1/ PHI0(20) B=COS(PHI0(I)) RETURN END FUNCTION C(I) C... C... FUNCTION C COMPUTES THE THIRD TERM IN THE LINEAR APPROXIMATION C... OF SIN(PHI) AS DISCUSSED PREVIOUSLY C... COMMON/C1/ PHI0(20) C=-COS(PHI0(I))*PHI0(I) RETURN END SUBROUTINE FCN(PHI,F,N,PAR) C... C... SUBROUTINE IMPLEMENTS THE NONLINEAR EQUATIONS C... C... COMMON/PAR/ CONTAINS THE PARAMETERS USED IN SUBROUTINE FCN C... CALLED BY ZSPOW COMMON/PAR/ DXS C... C... DIMENSION THE ARRAYS USED BY FCN DIMENSION PHI(N), F(N), PAR(1) C... C... IMPLEMENT THE NONLINEAR EQUATIONS DO 1 I=1,N IF(I.EQ.1)THEN F(1)=(PHI(2)-2.0*PHI(1)+1.0)/DXS-SIN(PHI(1)) ELSE IF(I.EQ.N)THEN F(N)=(2.0*PHI(N-1)-2.0*PHI(N))/DXS-SIN(PHI(N)) ELSE F(I)=(PHI(I+1)-2.0*PHI(I)+PHI(I-1))/DXS-SIN(PHI(I)) END IF 1 CONTINUE RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN C... C... FOR THE FIRST RUN (NORUN = 1), PHI IN EQUATION (8) IS COMPUTED C... BY EXPLICIT PROGRAMMING OF SECOND-ORDER APPROXIMATIONS IF(NORUN.EQ.1)CALL DERV1 C... C... FOR THE SECOND RUN (NORUN = 2), PHI IN EQUATION (8) IS COMPUTED C... BY PREPROGRAMMED SECOND-ORDER APPROXIMATIONS IF(NORUN.EQ.2)CALL DERV2 C... C... FOR THE THIRD RUN (NORUN = 3), PHI IN EQUATION (8) IS COMPUTED C... BY PREPROGRAMMED FOURTH-ORDER APPROXIMATIONS IF(NORUN.EQ.3)CALL DERV3 RETURN END SUBROUTINE DERV1 COMMON/T/ T, NSTOP, NORUN 1 /Y/ PHIU(21) 2 /F/ PHIT(21) 3 /S/ PHIX(21),PHIXX(21) 4 /C/ N, IP COMMON/PAR/ DXS C... C... BOUNDARY CONDITION PHIU(1)=1.0 C... C... EQUATION (8) DO 1 I=1,N IF(I.EQ.1)PHIT(1)=0. IF((I.GT.1).AND.(I.LT.N))PHIT(I)= 1 (PHIU(I+1)-2.0*PHIU(I)+PHIU(I-1))/DXS-SIN(PHIU(I)) IF(I.EQ.N)PHIT(N)=2.0*(PHIU(N-1)-PHIU(N))/DXS-SIN(PHIU(N)) 1 CONTINUE RETURN END SUBROUTINE DERV2 COMMON/T/ T, NSTOP, NORUN 1 /Y/ PHIU(21) 2 /F/ PHIT(21) 3 /S/ PHIX(21),PHIXX(21) 4 /C/ N, IP C... C... BOUNDARY CONDITIONS PHIU(1)=1.0 PHIX(N)=0.0 C... C... PHI BY SECOND-ORDER APPROXIMATIONS C... XX CALL DSS042(0.,1.,N,PHIU,PHIX,PHIXX,1,2) C... C... EQUATION (8) DO 1 I=1,N PHIT(I)=PHIXX(I)-SIN(PHIU(I)) 1 CONTINUE RETURN END SUBROUTINE DERV3 COMMON/T/ T, NSTOP, NORUN 1 /Y/ PHIU(21) 2 /F/ PHIT(21) 3 /S/ PHIX(21),PHIXX(21) 4 /C/ N, IP C... C... BOUNDARY CONDITIONS PHIU(1)=1.0 PHIX(N)=0.0 C... C... PHI BY FOURTH-ORDER APPROXIMATIONS C... XX CALL DSS044(0.,1.,N,PHIU,PHIX,PHIXX,1,2) C... C... EQUATION (8) DO 1 I=1,N PHIT(I)=PHIXX(I)-SIN(PHIU(I)) 1 CONTINUE RETURN END SUBROUTINE PRINT(NI,NO) C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ PHIU(21) 2 /F/ PHIT(21) 3 /S/ PHIX(21),PHIXX(21) 4 /C/ N, IP C... C... DIMENSION THE ARRAYS FOR PLOTTING DIMENSION XP(21), PHIP(6,21) C... C... PRINT A HERADING FOR THE SOLUTION IP=IP+1 IF(IP.EQ.1)WRITE(NO,1) 1 FORMAT(1H ,//,9X,'T', 1 4X,'PHI(0,T)',2X,'PHI(1/4,T)',2X,'PHI(1/2,T)', 2 2X,'PHI(3/4,T)',4X,'PHI(1,T)') C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,2)T,(PHIU(I),I=1,N,5) 2 FORMAT(F10.2,5F12.4) C... C... STORE THE SOLUTION FOR PLOTTING DX=1.0/FLOAT(N-1) DO 3 I=1,N XP(I)=DX*FLOAT(I-1) PHIP(IP,I)=PHIU(I) 3 CONTINUE C... C... PLOT THE SOLUTION IF(IP.LT.6)RETURN CALL TPLOTS(IP,N,XP,PHIP) C... RETURN END STEADY STATE AND TRANSIENT SINE-GORDON EQUATION 0. 2.5 0.5 2110000 1 1 REL 0.001 STEADY STATE AND TRANSIENT SINE-GORDON EQUATION 0. 2.5 0.5 2110000 1 1 REL 0.001 STEADY STATE AND TRANSIENT SINE-GORDON EQUATION 0. 2.5 0.5 2110000 1 1 REL 0.001 END OF RUNS *APP PHYS03 SUBROUTINE INITAL C... C... 1 + 1 SINE-GORDON EQUATION C... C... THE 1 + 1 SINE-GORDON EQUATION IS DISCUSSED IN DETAIL IN ROGERS, C... C. AND W. F. SHADWICK, BACKLUND TRANSFORMATIONS AND THEIR APPLI- C... CATIONS, ACADEMIC PRESS, NY, 1982. THE SINE-GORDON EQUATON HAS C... MANY PHYSICAL APPLICATIONS SUCH AS THE ANALYSIS OF ULTRASHORT C... OPTICAL PULSE PROPAGATION PHENOMENA. C... C... CONSIDER THE FOLLOWING ONE-DIMENSIONAL SINE-GORDON EQUATION C... C... U = SIN(U) (1) C... XT C... C... WITH THE BOUNDARY AND INITIAL CONDITIONS C... C... U(0,T) = 0, T GE 0, U (X,0) = SIN(PI*X), 0 LE X LE 1 (2)(3) C... X C... C... THE FOLLOWING CODE PRODUCES A NUMERICAL SOLUTION TO EQUATIONS C... (1) TO (3) WITH PLOTTING OF THE SOLUTION. C... C... COMMON AREA COMMON/T/ T, NSTOP, NORUN 1 /Y/ UX(0:20) 2 /F/UXT(0:20) 3 /C/ U(0:20), X(0:20), DX, N C... C... DEFINE THE SPATIAL GRID N=20 DX=1./FLOAT(N) X(0)=0. DO 1 I=0,N-1 X(I+1)=X(I)+DX 1 CONTINUE C... C... INITIAL CONDITION (3) PI=4.*ATAN(1.) DO 2 I=0,N UX(I)=SIN(PI*X(I)) 2 CONTINUE C... C... CALCULATE INITIAL DERIVATIVES CALL DERV RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ UX(0:20) 2 /F/UXT(0:20) 3 /C/ U(0:20), X(0:20), DX, N C... C... BOUNDARY CONDITION (2) U(0)=0. C... C... COMPUTE U FROM U BY TRAPEZOIDAL INTEGRATION C... X DO 1 I=1,N U(I)=U(I-1)+(UX(I)+UX(I-1))*DX/2. 1 CONTINUE C... C... EQUATION (1) DO 2 I=1,N UXT(I)=SIN(U(I)) 2 CONTINUE C... C... SINCE U(0,T) = 0 FROM BOUNDARY CONDITION (2), U (0,T) = 0 C... XT UXT(0)=0. RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ UX(0:20) 2 /F/UXT(0:20) 3 /C/ U(0:20), X(0:20), DX, N C... C... DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION DIMENSION XP(21), UP1(11,21), UXP1(11,21), UXTP1(11,21), 1 UP2( 6,21), UXP2( 6,21), UXTP2( 6,21) C... C... INITIALIZE COUNTERS FOR THE PLOTTED SOLUTIONS DATA IP1,IP2/2*0/ C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,1)T,(U(I),I=0,N),(UX(I),I=0,N),(UXT(I),I=0,N) 1 FORMAT(1H ,///,5H T = ,F5.1,//, 1 4H U,6F10.3,/,14X,5F10.3,/,14X,5F10.3,/,14X,5F10.3,//, 2 4H UX,6F10.3,/,14X,5F10.3,/,14X,5F10.3,/,14X,5F10.3,//, 3 4H UXT,6F10.3,/,14X,5F10.3,/,14X,5F10.3,/,14X,5F10.3) C... C... STORE THE SOLUTION FOR PLOTTING C... C... CONTINUOUS PLOTTER. THE SOLUTION IS STORED DURING EACH CALL TO C... SUBROUTINE PRINT IP1=IP1+1 DO 2 I=0,N XP(I+1)=X(I) UP1(IP1,I+1)=U(I) UXP1(IP1,I+1)=UX(I) UXTP1(IP1,I+1)=UXT(I) 2 CONTINUE C... C... DISCRETE (LINE-PRINTER) PLOTTER. THE SOLUTION IS STORED EVERY C... SECOND CALL TO SUBROUTINE PRINT TO PRODUCE A PLOT THAT HAS SIX C... SOLUTON CURVES RATHER THAN THE 11 SOLUTION CURVES FROM THE CON- C... TINUOUS PLOTTER IF(((IP1-1)/2*2).EQ.(IP1-1))THEN IP2=IP2+1 DO 3 I=0,N XP(I+1)=X(I) UP2(IP2,I+1)=U(I) UXP2(IP2,I+1)=UX(I) UXTP2(IP2,I+1)=UXT(I) 3 CONTINUE END IF C... C... PLOT THE SOLUTION IF(IP1.LT.11)RETURN C... C... ****************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND WILL HAVE TO BE CHANGED FOR OTHER COMPUTERS OR CONVERTED TO C... COMMENTS C... C... CALL QIKSAX(3,3) C... CALL QIKPLT(XP,UP1,N+1,3H*X*,8H*U(X,T)*, C... 1 28H*1 + 1 SINE-GORDON EQUATION*,IP1) C... CALL QIKSAX(3,3) C... CALL QIKPLT(XP,UXP1,N+1,3H*X*,9H*UX(X,T)*, C... 1 28H*1 + 1 SINE-GORDON EQUATION*,IP1) C... CALL QIKSAX(3,3) C... CALL QIKPLT(XP,UXTP1,N+1,3H*X*,10H*UXT(X,T)*, C... 1 28H*1 + 1 SINE-GORDON EQUATION*,IP1) C... C... ****************************************************************** C... CALL TPLOTS(IP2,N+1,XP,UP2) WRITE(NO,11) 11 FORMAT(1H ,//,47H U(X,T) VS X WITH T AS A PARAMETER ) CALL TPLOTS(IP2,N+1,XP,UXP2) WRITE(NO,12) 12 FORMAT(1H ,//,47H UX(X,T) VS X WITH T AS A PARAMETER ) CALL TPLOTS(IP2,N+1,XP,UXTP2) WRITE(NO,13) 13 FORMAT(1H ,//,47H UXT(X,T) VS X WITH T AS A PARAMETER ) C... C... REINITIALIZE THE PLOTTING COUNTERS FOR ANOTHER RUN IP1=0 IP2=0 RETURN END 1 + 1 SINE-GORDON EQUATION 0. 10. 1. 21 1000 1 1 REL 0.001 END OF RUNS *APP PHYS04 SUBROUTINE INITAL C... C... LIOUVILLE*S EQUATION C... C... LIOUVILLE*S EQUATION IS DISCUSSED IN ROGERS, C., AND W. F. SHAD- C... WICK, BACKLUND TRANSFORMATIONS AND THEIR APPLICATIONS, ACADEMIC C... PRESS, NY, 1982. C... C... CONSIDER THE FOLLOWING ONE-DIMENSIONAL FORM OF LIOUVILLE*S C... EQUATION C... C... U = EXP(U) (1) C... XT C... C... WITH THE BOUNDARY AND INITIAL CONDITIONS C... C... U(0,T) = 0, T GE 0, U (X,0) = SIN(PI*X), 0 LE X LE 1 (2)(3) C... X C... C... THE FOLLOWING CODE PRODUCES A NUMERICAL SOLUTION TO EQUATIONS C... (1) TO (3) WITH PLOTTING OF THE SOLUTION. C... C... COMMON AREA COMMON/T/ T, NSTOP, NORUN 1 /Y/ UX(0:20) 2 /F/UXT(0:20) 3 /C/ U(0:20), X(0:20), DX, N C... C... DEFINE THE SPATIAL GRID N=20 DX=1./FLOAT(N) X(0)=0. DO 1 I=0,N-1 X(I+1)=X(I)+DX 1 CONTINUE C... C... INITIAL CONDITION (3) PI=4.*ATAN(1.) DO 2 I=0,N UX(I)=SIN(PI*X(I)) 2 CONTINUE C... C... CALCULATE INITIAL DERIVATIVES CALL DERV RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ UX(0:20) 2 /F/UXT(0:20) 3 /C/ U(0:20), X(0:20), DX, N C... C... BOUNDARY CONDITION (2) U(0)=0. C... C... COMPUTE U FROM U BY TRAPEZOIDAL INTEGRATION C... X DO 1 I=1,N U(I)=U(I-1)+(UX(I)+UX(I-1))*DX/2. 1 CONTINUE C... C... EQUATION (1) DO 2 I=1,N UXT(I)=EXP(U(I)) 2 CONTINUE C... C... SINCE U(0,T) = 0 FROM BOUNDARY CONDITION (2), U (0,T) = 0 C... XT UXT(0)=0. RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ UX(0:20) 2 /F/UXT(0:20) 3 /C/ U(0:20), X(0:20), DX, N C... C... DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION DIMENSION XP(21), UP1(11,21), UXP1(11,21), UXTP1(11,21), 1 UP2( 6,21), UXP2( 6,21), UXTP2( 6,21) C... C... INITIALIZE COUNTERS FOR THE PLOTTED SOLUTIONS DATA IP1,IP2/2*0/ C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,1)T,(U(I),I=0,N),(UX(I),I=0,N),(UXT(I),I=0,N) 1 FORMAT(1H ,///,5H T = ,F5.1,//, 1 4H U,6F10.3,/,14X,5F10.3,/,14X,5F10.3,/,14X,5F10.3,//, 2 4H UX,6F10.3,/,14X,5F10.3,/,14X,5F10.3,/,14X,5F10.3,//, 3 4H UXT,6F10.3,/,14X,5F10.3,/,14X,5F10.3,/,14X,5F10.3) C... C... STORE THE SOLUTION FOR PLOTTING C... C... CONTINUOUS PLOTTER. THE SOLUTION IS STORED DURING EACH CALL TO C... SUBROUTINE PRINT IP1=IP1+1 DO 2 I=0,N XP(I+1)=X(I) UP1(IP1,I+1)=U(I) UXP1(IP1,I+1)=UX(I) UXTP1(IP1,I+1)=UXT(I) 2 CONTINUE C... C... DISCRETE (LINE-PRINTER) PLOTTER. THE SOLUTION IS STORED EVERY C... SECOND CALL TO SUBROUTINE PRINT TO PRODUCE A PLOT THAT HAS SIX C... SOLUTON CURVES RATHER THAN THE 11 SOLUTION CURVES FROM THE CON- C... TINUOUS PLOTTER IF(((IP1-1)/2*2).EQ.(IP1-1))THEN IP2=IP2+1 DO 3 I=0,N XP(I+1)=X(I) UP2(IP2,I+1)=U(I) UXP2(IP2,I+1)=UX(I) UXTP2(IP2,I+1)=UXT(I) 3 CONTINUE END IF C... C... PLOT THE SOLUTION IF(IP1.LT.11)RETURN C... C... ****************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND WILL HAVE TO BE CHANGED FOR OTHER COMPUTERS OR CONVERTED TO C... COMMENTS C... C... CALL QIKSAX(3,3) C... CALL QIKPLT(XP,UP1,N+1,3H*X*,8H*U(X,T)*, C... 1 28H* LIOUVILLE'S EQUATION *,IP1) C... CALL QIKSAX(3,3) C... CALL QIKPLT(XP,UXP1,N+1,3H*X*,9H*UX(X,T)*, C... 1 28H* LIOUVILLE'S EQUATION *,IP1) C... CALL QIKSAX(3,3) C... CALL QIKPLT(XP,UXTP1,N+1,3H*X*,10H*UXT(X,T)*, C... 1 28H* LIOUVILLE'S EQUATION *,IP1) C... C... ****************************************************************** C... CALL TPLOTS(IP2,N,XP,UP2) WRITE(NO,11) 11 FORMAT(1H ,//,47H U(X,T) VS X WITH T AS A PARAMETER ) CALL TPLOTS(IP2,N,XP,UXP2) WRITE(NO,12) 12 FORMAT(1H ,//,47H UX(X,T) VS X WITH T AS A PARAMETER ) CALL TPLOTS(IP2,N,XP,UXTP2) WRITE(NO,13) 13 FORMAT(1H ,//,47H UXT(X,T) VS X WITH T AS A PARAMETER ) C... C... REINITIALIZE THE PLOTTING COUNTERS FOR ANOTHER RUN IP1=0 IP2=0 RETURN END LIOUVILLE*S EQUATION 0. 1. 0.1 21 1000 1 1 REL 0.001 END OF RUNS *APP PHYS05 SUBROUTINE INITAL C... C... DYNAMIC RESPONSE OF A MANOMETER C... C... THE DIFFERENTIAL HEIGHT OF LIQUID IN A U-TUBE MANOMETER, H, CAN C... BE CALCULATED AS A FUNCTION OF TIME, T, BY INTEGRATION OF THE C... SECOND-ORDER ORDINARY DIFFERENTIAL EQUATION (ODE) C... C... (TAU**2)D2H/DT2 + 2*TAU*Z*DH/DT + H = F(T) (1) C... C... WHERE C... C... TAU CHACTERISTIC TIME (RECIPROCAL OF THE NATURAL C... FREQUENCY) OF THE MANOMETER (SEC) C... C... Z DAMPING COEFFICIENT (DIMENSIONLESS) C... C... F(T) FORCING FUNCTION CORRESPONDING TO THE PRESSURE APPLIED C... TO ONE LEG OF THE MANOMETER (M) C... C... FOR A MANOMETER, TAU AND Z ARE DEFINED BY C... C... TAU = (L/(2*G))**0.5, Z = (8*VIS/((D**2)*RHO))*(2*L/G)**0.5 C... C... WHERE FOR THE PARTICULAR MANOMETER CONSIDERED IN THE FOLLOWING C... ANALYSIS C... C... L TOTAL LENGTH OF THE LIQUID IN THE MAMOMETER = 0.9915 M C... C... G ACCELERATION OF GRAVITY = 9.8 M/SEC**2 C... C... VIS LIQUID VISCOSITY = 2.5E-03 PA-SEC (1 PA = 1 PASCAL = C... NEWTON/M**2) C... C... D TUBE DIAMETER (M) C... C... RHO LIQUID DENSITY = 13600 KG/M**3 C... C... IN THE FOLLOWING PROGRAMMING, THE TUBE DIAMETER, D, IS SET TO C... THREE VALUES IN THREE RUNS CORRESPONDING TO Z GT 1 (OVERDAMPED), C... Z = 1 (CRITICALLY DAMPED) AND Z LT 1 (UNDERDAMPED) C... C... INITIALLY, THE LIQUID IN THE MAMOMETER IS AT REST. THEREFORE, THE C... INITIAL CONDITIONS FOR EQUATION (1) ARE C... C... H(0) = 0, DH(O)/DT = 0 (2)(3) C... C... AT T = 0, A PRESSURE, P1, OF 2.0E+04 PA IS APPLIED TO ONE OF THE C... LEGS OF THE MAMOMETER. THE OTHER LEG HAS A VACUUM ABOVE IT. THUS C... F(T) IN EQUATION (1) IS A CONSTANT WITH THE VALUE C... C... F(T) = GC*P1/(RHO*G) = 1*2.0E+04/(13600*9.8) = 0.1501 M C... C... WHERE C... C... GC CONVERSION FACTOR = 1 KG-M/N-SEC**2 C... C... EQUATIONS (1) TO (3) CAN BE SOLVED ANALYTICALLY, E.G., BY THE C... LAPLACE TRANSFORM. THE ANALYTICAL SOLUTIONS FOR THE THREE CASES, C... Z GT 1, Z = 1, Z LT 1, ARE PROGRAMMED IN SUBROUTINE PRINT AND C... PRINTED ALONG WITH THE NUMERICAL SOLUTION C... C... THE DERIVATION OF EQUATION (1), AND ITS LAPLACE TRANSFORM SOLUTION C... ARE OUTLINED BELOW C... C... A CONTINUITY BALANCE ON LEG 1 OF THE MANOMETER GIVES C... C... D(H1*S*RHO)/DT = -V*S*RHO C... C... WHERE C... C... H1 HEIGHT OF LIQUID IN LEG 1 (M) C... C... S CROSS SECTIONAL AREA OF THE TUBE (M**2) C... C... V FLUID VELOCITY (M/SEC) C... C... SIMILARLY, FOR LEG 2 C... C... D(H2*S*RHO)/DT = V*S*RHO C... C... A MOMENTUM BALANCE ON THE FLUID IN THE TUBE GIVES C... C... (1/GC)*D(RHO*L*S*V)/DT = (G/GC)*RHO*H1*S - (G/GC)*RHO*H2*S C... C... + P1*S - 0*S - (F*V*ABS(V)/(2*D*GC))*S C... C... F = 64*VIS/(D*ABS(V)*RHO) C... C... THESE EQUATIONS CAN BE SIMPLIFIED (ASSUMING CONSTANT S AND RHO) TO C... C... DH1/DT = -V C... C... DH2/DT = V C... C... DV/DT = (G/L)*(H1 - H2) + P1*GC/(L*RHO) - 32*VIS*V/((D**2)*RHO) C... C... THESE EQUATIONS CAN EASILY BE REARRANGED TO GIVE AN EQUATION IN C... H = H2 - H1 C... C... V = (1/2)*D(H2 - H1)/DT = (1/2)*DH/DT C... C... (1/2)D2H/DT2 = -(G/L)*H + (P1*GC/(L*RHO)) C... C... - (16*VIS/((D**2)*RHO))*DH/DT C... C... WHICH CAN BE REARRANGED TO C... C... (L/(2*G))*D2H/DT2 + (16*VIS*L/((D**2)*RHO*G))*DH/DT + H = C... C... (GC/(RHO*G))*P1 C... C... EQUATION (1) THEN FOLLOWS FROM THIS RESULT WITH TAU AND Z DEFINED C... PREVIOUSLY C... C... EQUATION (1) CAN BE LAPLACE TRANSFORMED TO C... C... ((TAU**2)*(S**2) + 2*TAU*Z*S + 1)*H(S) = (GC*P1/(RHO*G))*(1/S) C... C... WHICH CAN THEN BE SOLVED FOR THE TRANSFORMED HEIGHT, H(S). THE C... CHARACTERISTIC POLYNOMIAL CAN BE WRITTEN IN FACTORED FORM AS C... C... (TAU**2)*(S**2) + 2*TAU*Z*S + 1 = (TAU**2)*(S - R1)*(S - R2) C... C... WHERE R1 AND R2 ARE GIVEN BY THE QUADRATIC FORMULA C... C... R1, R2 = (-2*TAU*Z + (OR -)((2*TAU*Z)**2 - 4*(TAU**2))**0.5)/ C... C... (2*(TAU**2)) C... OR C... C... R1, R2 = (-Z + (OR -)(Z**2 - 1)**0.5)/TAU C... C... R1 AND R2 WILL BE REAL AND DISTINCT, REAL AND REPEATED OR COMPLEX C... FOR Z GT 1, Z = 1 AND Z LT 1. THE ANALYSIS WILL BE CONTINUED FOR C... VALUES OF TAU AND Z CORRESPONDING TO THREE TUBE DIAMETERS C... C... CASE 1 - D = 4.50E-04 M C... C... Z = 3.9691E-07/(4.50E-04)**2 = 1.9600 (OVERDAMPED) C... C... CASE 2 - D = 6.30E-04 M C... C... Z = 3.9691E-07/(6.30E-04)**2 = 1.0000 (CRITICALLY C... DAMPED) C... CASE 3 - D = 1.80E-03 M C... C... Z = 3.9691E-07/(1.80E-03)**2 = 0.1225 (UNDERDAMPED) C... C... SUBSTITUTION OF THESE CONSTANTS IN THE LAPLACE TRANSFORM OF H C... GIVES C... C... H(S) = (0.1501/S)/(0.05058*S**2 + 2*Z*(0.02249)*S + 1) C... C... (0.1501/S)/(0.05085*(S - R1)*(S - R2)) C... C... WHERE FROM THE QUADRATIC FORMULA C... C... CASE 1 - R1 = -1.2196, R2 = -16.210 C... C... CASE 2 - R1 = R2 = R = -4.4464 C... C... CASE 3 - R1 = -0.54468 + 4.41291I, R2 = -0.54468 - 4.41291I C... C... (I = SQRT(-1) SO R1 AND R2 ARE COMPLEX CONJUGATES) C... C... INVERSION OF H(S) CAN BE ACCOMPLISHED BY PARTIAL FRACTIONS C... C... CASE 1 C... C... H(S) = A/S + B/(S - R1) + C/(S - R2) C... C... WHERE A, B AND C ARE CONSTANTS TO BE DETERMINED (THE MULTIPLYING C... CONSTANT 0.1501/0.05058 WILL BE INCLUDED AFTER THE INVERSION). IN C... THE USUAL WAY FOR PARTIAL FRACTION EXPANSIONS C... C... A = 1/(R1*R2) = 1/((-1.12196)*(-16.120)) = 0.05058 C... C... B = 1/(R1*(R1 - R2)) = -0.05470 C... C... C = 1/(R2*(R2 - R1)) = 0.004115 C... C... SUBSTITUTION IN H(S) FOLLOWED BY INVERSION GIVES C... C... H(T) = 0.1501*(1 - 1.0815*EXP(-1.2196*T) C... C... + 0.0815*EXP(-16.210*T) C... C... NOTE BY INSPECTION THAT H(0) = 0 WHICH IS A GOOD CHECK OF THE C... CONSTANTS A, B AND C C... C... CASE 2 C... C... H(S) = A/S + B/(S - R) + C/(S - R)**2 C... C... WHICH IS THE USUAL PARTIAL FRACTION EXPANSION FOR REPEATED ROOTS C... C... A = 1/R**2 = 1/(-4.4464)**2 = 0.05058 C... C... B = D(1/S)/DS WITH S = R = -1/R**2 = -0.05058 C... C... C = 1/R = -0.2249 C... C... SUBSTITUTION IN H(S) FOLLOWED BY INVERSION GIVES C... C... H(T) = 0.1501*(1 - EXP(-4.4464*T) - 4.4464*T*EXP(-4.4464*T)) C... C... AGAIN, BY INSPECTION H(0) = 0 C... C... CASE 3 C... C... H(S) = A/S + B/(S - R1) + C/(S - R2) C... C... WHERE A, B AND C CAN BE EVALUATED AS BEFORE IN CASE 1. THUS C... C... A = 1/(R1*R2) = 1/((-0.54468 + 4.41291I)*(-0.55468 - 4.41291I)) C... C... = 0.05058 C... C... B = 1/(R1*(R1 - R2)) = 1/((-0.55468 + 4.41291I)*(2)*(4.41291I)) C... C... = 1/(-38.9475 - 4.80725I) = 0.02548*EXP(-I*PHI) C... C... WITH PHI = ARCTAN(-4.80725/-38.9475) = 187.04 DEGREES C... C... C = 1/(R2*(R2 - R1)) = 0.02548*EXP(I*PHI) C... C... SUBSTITUTION IN H(S) FOLLOWED BY INVERSION GIVES C... C... H(T) = 0.1501*(1 + 0.02548*EXP(-I*PHI)*EXP(R1*T) C... C... + 0.02548*EXP( I*PHI)*EXP(R2*T)) C... C... WHICH CAN BE WRITTEN IN A FORM MORE CONVENIENT FOR CALCULATION (OF C... COURSE, H(T) MUST BE REAL) C... C... H(T) = 0.1501*(1 C... C... + 0.50376*EXP(-I*PHI)*EXP(RE(R1)*T)*EXP(I*IM(R1)*T)) C... C... + 0.50376*EXP( I*PHI)*EXP(RE(R2)*T)*EXP(I*IM(R2)*T)) C... C... = 0.1501*(1 C... C... + 0.50376*EXP(-0.54468*T)*(EXP( I*(IM(R1)*T - PHI)) + C... C... * EXP(-I*(IM(R1)*T - PHI))) C... C... = 0.1501*(1 C... C... + (2)*(0.50376)*EXP(-0.54468*T)*COS(4.41291*T - PHI)) C... C... WHICH FOLLOWS FROM THE IDENTITY C... C... COS(X) = (EXP(I*X) + EXP(-I*X))/2 C... C... AGAIN, H(0) = 0 C... C... THE THREE ANALYTICAL SOLUTIONS, H(T) FOR THE THREE CASES, ARE PRO- C... GRAMMED IN SUBROUTINE PRINT C... COMMON /T/ T, NOSTOP, NORUN 1 /Y/ Y1, Y2 2 /F/ DY1DT, DY2DT 3 /C/ TAU, Z, F, IP C... C... SET THE MODEL PARAMETERS C... C... TIME CONSTANT C... D(H2*S*RHO)/DT = V*S*RHO XL=0.9915 G=9.8 TAU=SQRT(XL/(2.*G)) C... C... DAMPING COEFFICIENT VIS=1.5E-03 RHO=13600. IF(NORUN.EQ.1)DIA=4.50E-04 IF(NORUN.EQ.2)DIA=6.30E-04 IF(NORUN.EQ.3)DIA=1.80E-03 Z=(8.*VIS/((DIA**2)*RHO))*SQRT(2.*XL/G) C... C... FORCING FUNCTION GC=1. PA=2.0E+04 F=GC*PA/(RHO*G) C... C... INITIAL CONDITIONS Y1=0. Y2=0. C... C... INITIALIZE THE DERIVATIVE CALCULATIONS CALL DERV C... C... SET A COUNTER FOR THE PLOTTED SOLUTION IP=0 RETURN END SUBROUTINE DERV COMMON /T/ T, NOSTOP, NORUN 1 /Y/ Y1, Y2 2 /F/ DY1DT, DY2DT 3 /C/ TAU, Z, F, IP C... C... MODEL DIFFERENTIAL EQUATIONS DY1DT=Y2 DY2DT=(1./TAU**2)*(F-2.*TAU*Z*Y2-Y1) RETURN END SUBROUTINE PRINT(NI,NO) COMMON /T/ T, NOSTOP, NORUN 1 /Y/ Y1, Y2 2 /F/ DY1DT, DY2DT 3 /C/ TAU, Z, F, IP C... C... DIMENSION THE ARRAYS TO STORE THE NUMERICAL SOLUTION FOR PLOTTING DIMENSION TP(51),YP(3,51),YAP(3,51) C... C... PRINT A HEADING FOR THE NUMERICAL SOLUTION IF(IP.EQ.0)WRITE(NO,1)TAU,Z,F 1 FORMAT( 1 23H CHARACTERISTIC TIME = ,F7.3,/, 2 23H DAMPING COEFFICIENT = ,F7.3,/, 3 20H FORCING FUNCTION = , F10.3,//, 4 9X,1HT,3X,7HH (NUM),2X,8HH (ANAL)) C... C... CALCULATE THE ANALYTICAL SOLUTION IF(NORUN.EQ.1)Y1A=F*(1.-1.0815*EXP(-1.2196*T) 1 +0.0815*EXP(-16.210*T)) IF(NORUN.EQ.2)Y1A=F*(1.-EXP(-4.4464*T)-4.4464*T*EXP(-4.4464*T)) IF(NORUN.EQ.3)Y1A=F*(1.+1.00751*EXP(-0.54468*T) 1 *COS(4.41291*T-187.04*3.1415927/180.)) C... C... PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS WRITE(NO,2)T,Y1,Y1A 2 FORMAT(F10.1,2F10.4) C... C... STORE THE NUMERICAL SOLUTION FOR PLOTTING IP=IP+1 TP(IP)=T YP(NORUN,IP)=Y1 YAP(NORUN,IP)=Y1A C... C... PLOT THE NUMERICAL SOLUTION AT THE END OF THE THIRD RUN IF(IP.LT. 51)RETURN IF(NORUN.LT.3)RETURN C... C... PLOT THE NUMERICAL SOLUTIONS FOR THE THREE CASES CALL TPLOTS(NORUN,IP,TP,YP) C... C... LABEL THE PLOT WRITE(NO,3) 3 FORMAT(1H ,//, 1 17H 1 - D = 4.50E-04,/, 2 17H 2 - D = 6.30E-04,/, 3 17H 3 - D = 1.80E-03) C... C... PLOT THE ANALYTICAL SOLUTIONS FOR THE THREE CASES CALL TPLOTS(NORUN,IP,TP,YAP) WRITE(NO,3) RETURN END DYNAMICS OF A MANOMETER, D = 4.50E-04 0. 10. 0.2 2 1000 1 1 REL 0.001 DYNAMICS OF A MANOMETER, D = 6.30E-04 0. 10. 0.2 2 1000 1 1 REL 0.001 DYNAMICS OF A MANOMETER, D = 1.80E-03 0. 10. 0.2 2 1000 1 1 REL 0.001 END OF RUNS *APP PHYS06 SUBROUTINE INITAL C... C... SIMULATION OF A SPACECRAFT MODELED AS A TRIANGULAR TRUSS C... C... THE FOLLOWING EQUATIONS HAVE BEEN PROPOSED BY KANE, LIKINS AND C... LEVINSON* FOR THE DYNAMICS OF A TRIANGULAR TRUSS THAT MOVES IN C... ITS OWN PLANE C... C... (F*L(1) + K(1)*(D(1) - L(1)))*(X(3) - X(1)) C... M(1)*DV(1)/DT - ------------------------------------------- C... D(1) C... C... K(3)*(D(3) - L(3))*(X(1) - X(5)) C... + -------------------------------- = 0 C... D(3) C... C... (F*L(1) + K(1)*(D(1) - L(1)))*(X(4) - X(2)) C... M(1)*DV(2)/DT - ------------------------------------------- C... D(1) C... C... K(3)*(D(3) - L(3))*(X(2) - X(6)) C... + -------------------------------- = 0 C... D(3) C... C... (F*L(2) + K(2)*(D(2) - L(2)))*(X(5) - X(3)) C... M(2)*DV(3)/DT - ------------------------------------------- C... D(2) C... C... K(1)*(D(1) - L(1))*(X(3) - X(1)) C... + -------------------------------- = 0 C... D(1) C... C... (F*L(2) + K(2)*(D(2) - L(2)))*(X(6) - X(4)) C... M(2)*DV(4)/DT - ------------------------------------------- C... D(2) C... C... K(1)*(D(1) - L(1))*(X(4) - X(2)) C... + -------------------------------- = 0 C... D(1) C... C... (F*L(3) + K(3)*(D(3) - L(3)))*(X(1) - X(5)) C... M(3)*DV(5)/DT - ------------------------------------------- C... D(3) C... C... K(2)*(D(2) - L(2))*(X(5) - X(3)) C... + -------------------------------- = 0 C... D(2) C... C... (F*L(3) + K(3)*(D(3) - L(3)))*(X(2) - X(6)) C... M(3)*DV(6)/DT - ------------------------------------------- C... D(3) C... C... K(2)*(D(2) - L(2))*(X(6) - X(4)) C... + -------------------------------- = 0 C... D(2) C... C... WHERE THE DISTANCES ARE GIVEN BY C... C... D(1) = ((X(3) - X(1))**2 + (X(4) - X(2))**2)**0.5 C... C... D(2) = ((X(5) - X(3))**2 + (X(6) - X(4))**2)**0.5 C... C... D(3) = ((X(1) - X(5))**2 + (X(2) - X(6))**2)**0.5 C... C... *KANE, THOMAS R., PETER W. LIKINS AND DAVID A. LEVINSON, C... SPACECRAFT DYNAMICS, MCGRAW-HILL BOOK COMPANY, NEW YORK, C... 1983 C... C... A SOLUTION TO THE PRECEDING EQUATIONS IS TO BE COMPUTED FOR THE C... TWO CASES DISCUSSED BY KANE ET AL C... C... (1) SECTION (B), PAGE 409 C... C... AT T = 0 C... C... X(1) = 0.01 M C... C... X(2) = 0.02 M C... C... X(3) = 4.03 M C... C... X(4) = -0.01 M C... C... X(5) = -0.02 M C... C... X(6) = 3.04 M C... C... V(1) = 0.01 M/SEC C... C... V(2) = -0.01 M/SEC C... C... V(3) = 0.02 M/SEC C... C... V(4) = -0.02 M/SEC C... C... V(5) = 0.03 M/SEC C... C... V(6) = -0.03 M/SEC C... C... (2) SECTION (E), PAGE 411 C... C... AT T = 0 C... C... X(1) = 0.01 M C... C... X(2) = 0.02 M C... C... X(3) = 4.03 M C... C... X(4) = -0.01 M C... C... X(5) = -0.02 M C... C... X(6) = 3.04 M C... C... V(I) = 0 M/SEC, I = 1, 2, 3, 4, 5, 6 C... C... FOR BOTH CASES, THE CONSTANTS IN THE PRECEDING EQUATIONS ARE C... C... F = 1000 N/M, L(1) = 4 M, L(2) = 5 M, L(3) = 3 M, M(1) = C... C... 3.5 KG, M(2) = 4.5 KG, M(3) = 4.0 KG, K(1) = 2.5E+04 N/M, C... C... K(2) = 2.0E+04 N/M, K(3) = 3.33333E+04 N/M. C... C... THE PLOT OF THE ANGLE PHI VS T IN FIGURE P4.21(B) IS OF PARTICULAR C... INTEREST IN PRESENTING THE NUMERICAL SOLUTION. C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(6), V(6) 2 /F/ XT(6), VT(6) 3 /C/ M(3), L(3), K(3), D(3), F, N, IP REAL M, L, K C... C... SET THE MODEL PARAMETERS C... C... NUMBER OF DISPLACEMENT AND VELOCITY ORDINARY DIFFERENTIAL EQUA- C... TIONS (ODES) N=6 C... C... POINT MASSES M(1)=3.5 M(2)=4.5 M(3)=4.0 C... C... SPRING LENGTHS L(1)=4.0 L(2)=5.0 L(3)=3.0 C... C... SPRING CONSTANTS K(1)=2.5E+04 K(2)=2.0E+04 K(3)=(1./3.)*(1.0E+05) C... C... APPLIED FORCE F=1000. C... C... SET THE INITIAL CONDITIONS C... C... CASE 1 (PART (B), PAGE 409) IF(NORUN.EQ.1)THEN C... C... DISPLACEMENTS X(1)= 0.01 X(2)= 0.02 X(3)= 4.03 X(4)=-0.01 X(5)=-0.02 X(6)= 3.04 C... C... VELOCITIES V(1)= 0.01 V(2)=-0.01 V(3)= 0.02 V(4)=-0.02 V(5)= 0.03 V(6)=-0.03 END IF C... C... CASE 2 (PART (E), PAGE 411) IF(NORUN.EQ.2)THEN C... C... DISPLACEMENTS X(1)=0. X(2)=0. X(3)=4. X(4)=0. X(5)=0. X(6)=3. C... C... VELOCITIES DO 1 I=1,N V(I)=0. 1 CONTINUE END IF C... C... CALCULATE THE INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(6), V(6) 2 /F/ XT(6), VT(6) 3 /C/ M(3), L(3), K(3), D(3), F, N, IP REAL M, L, K C... C... MODEL ALGEBRA C... C... INTERPARTICLE DISTANCES CALL DIS C... C... MODEL ODES C... C... POSITION EQUATIONS CALL POS C... C... VELOCITY EQUATIONS CALL VEL RETURN END SUBROUTINE DIS COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(6), V(6) 2 /F/ XT(6), VT(6) 3 /C/ M(3), L(3), K(3), D(3), F, N, IP REAL M, L, K D(1)=(X(3)-X(1))**2+(X(4)-X(2))**2 D(1)=SQRT(D(1)) D(2)=(X(5)-X(3))**2+(X(6)-X(4))**2 D(2)=SQRT(D(2)) D(3)=(X(1)-X(5))**2+(X(2)-X(6))**2 D(3)=SQRT(D(3)) RETURN END SUBROUTINE POS COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(6), V(6) 2 /F/ XT(6), VT(6) 3 /C/ M(3), L(3), K(3), D(3), F, N, IP REAL M, L, K DO 1 I=1,N XT(I)=V(I) 1 CONTINUE RETURN END SUBROUTINE VEL COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(6), V(6) 2 /F/ XT(6), VT(6) 3 /C/ M(3), L(3), K(3), D(3), F, N, IP REAL M, L, K VT(1)=(1./M(1))*((F*L(1)+K(1)*(D(1)-L(1)))*(X(3)-X(1))/D(1) + -K(3)*(D(3)-L(3))*(X(1)-X(5))/D(3)) VT(2)=(1./M(1))*((F*L(1)+K(1)*(D(1)-L(1)))*(X(4)-X(2))/D(1) + -K(3)*(D(3)-L(3))*(X(2)-X(6))/D(3)) VT(3)=(1./M(2))*((F*L(2)+K(2)*(D(2)-L(2)))*(X(5)-X(3))/D(2) + -K(1)*(D(1)-L(1))*(X(3)-X(1))/D(1)) VT(4)=(1./M(2))*((F*L(2)+K(2)*(D(2)-L(2)))*(X(6)-X(4))/D(2) + -K(1)*(D(1)-L(1))*(X(4)-X(2))/D(1)) VT(5)=(1./M(3))*((F*L(3)+K(3)*(D(3)-L(3)))*(X(1)-X(5))/D(3) + -K(2)*(D(2)-L(2))*(X(5)-X(3))/D(2)) VT(6)=(1./M(3))*((F*L(3)+K(3)*(D(3)-L(3)))*(X(2)-X(6))/D(3) + -K(2)*(D(2)-L(2))*(X(6)-X(4))/D(2)) RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(6), V(6) 2 /F/ XT(6), VT(6) 3 /C/ M(3), L(3), K(3), D(3), F, N, IP REAL M, L, K C... C... DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION DIMENSION TP(101), X1P(101), X2P(101), 1 X3P(101), X4P(101), 2 X5P(101), X6P(101), 3 D1P(101), D2P(101), D3P(101), 4 COSP(101), ANGLEP(101), ANGLE1(101) C... C... DIMENSION AND DEFINE AN ARRAY FOR PRINTING THE ODE NUMBERS DIMENSION IS(6) DATA IS/1,2,3,4,5,6/ C... C... INCREMENT THE COUNTER FOR THE PLOTTED POINTS IP=IP+1 C... C... PRINT THE SOLUTION EVERY 10TH CALL TO SUBROUTINE PRINT IF(((IP-1)/10*10).EQ.(IP-1))THEN WRITE(NO,1)T,(IS(I),I=1,N), 1 ( X(I),I=1,N), 2 ( V(I),I=1,N), 3 (VT(I),I=1,N) 1 FORMAT(1H ,//, 5H T = ,F7.3,//, 1 7H I,6I12,/, 2 7H X(I),6E12.3,/, 3 7H V(I),6E12.3,/, 4 7H VT(I),6E12.3) END IF C... C... STORE THE SOLUTION FOR SUBSEQUENT PLOTTING C... C... TIME TP(IP)=T C... C... POSITION X1P(IP)=X(1) X2P(IP)=X(2) X3P(IP)=X(3) X4P(IP)=X(4) X5P(IP)=X(5) X6P(IP)=X(6) C... C... DISTANCES CALL DIS D1P(IP)=D(1) D2P(IP)=D(2) D3P(IP)=D(3) C... C... ANGLE A=D(1) B=D(3) C=D(2) COSINE=(A**2+B**2-C**2)/(2.*A*B) COSP(IP)=COSINE ANGLE=(180./3.1415927)*ACOS(COSINE)-90. ANGLEP(IP)=ANGLE IF(NORUN.EQ.1)ANGLE1(IP)=ANGLE C... C... PLOT THE SOLUTION C... C... POSITION (FOR EACH RUN) IF(IP.NE.101)THEN RETURN ELSE C... C... CALLS TO A POINT PLOTTER CALL TPLOTS(1,IP,X1P,X2P) WRITE(NO,11) 11 FORMAT(1H ,//,40H X(2) VS X(1) ) CALL TPLOTS(1,IP,X3P,X4P) WRITE(NO,12) 12 FORMAT(1H ,//,40H X(4) VS X(3) ) CALL TPLOTS(1,IP,X5P,X6P) WRITE(NO,13) 13 FORMAT(1H ,//,40H X(6) VS X(5) ) CALL TPLOTS(1,IP,TP,D1P) WRITE(NO,14) 14 FORMAT(1H ,//,40H D(1) VS T ) CALL TPLOTS(1,IP,TP,D2P) WRITE(NO,15) 15 FORMAT(1H ,//,40H D(2) VS T ) CALL TPLOTS(1,IP,TP,D3P) WRITE(NO,16) 16 FORMAT(1H ,//,40H D(3) VS T ) CALL TPLOTS(1,IP,TP,COSP) WRITE(NO,17) 17 FORMAT(1H ,//,40H COS(PSI) VS T ) CALL TPLOTS(1,IP,TP,ANGLEP) WRITE(NO,18) 18 FORMAT(1H ,//,40H PHI VS T ) C... C... CALLS TO A CONTINUOUS PLOTTER C... C... ***************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND MUST BE CONVERTED OR CHANGED TO COMMENTS C... C... CALL QIKSAX(3,3) C... CALL QIKSET(6.,-7.5,2.5,6.,-7.5,2.5) C... CALL QIKPLT(X1P,X2P,IP,7H*X1(T)*,7H*X2(T)*,2H**,1) C... CALL QIKSAX(3,3) C... CALL QIKSET(6.,-7.5,2.5,6.,-7.5,2.5) C... CALL QIKPLT(X3P,X4P,IP,7H*X3(T)*,7H*X4(T)*,2H**,1) C... CALL QIKSAX(3,3) C... CALL QIKSET(6.,-7.5,2.5,6.,-7.5,2.5) C... CALL QIKPLT(X5P,X6P,IP,7H*X5(T)*,7H*X6(T)*,2H**,1) C... CALL QIKSAX(3,3) C... CALL QIKSET(6.0,0.0,0.1,8.0,0.0,1.0) C... CALL QIKPLT(TP,D1P,IP,3H*T*,6H*D(1)*,2H**,1) C... CALL QIKSAX(3,3) C... CALL QIKSET(6.0,0.0,0.1,8.0,0.0,1.0) C... CALL QIKPLT(TP,D2P,IP,3H*T*,6H*D(2)*,2H**,1) C... CALL QIKSAX(3,3) C... CALL QIKSET(6.0,0.0,0.1,8.0,0.0,1.0) C... CALL QIKPLT(TP,D3P,IP,3H*T*,6H*D(3)*,2H**,1) C... CALL QIKSAX(3,3) C... CALL QIKSET(6.0,0.0,0.1,8.0,-1.0,0.25) C... CALL QIKPLT(TP,COSP,IP,3H*T*,10H*COS(PHI)*,2H**,1) C... CALL QIKSAX(3,3) C... CALL QIKSET(6.0,0.0,0.1,8.0,-10.0,5.0) C... CALL QIKPLT(TP,ANGLEP,IP,3H*T*,5H*PHI*,2H**,1) END IF C... C... BOTH ANGLES AT THE END OF THE SECOND RUN FOR FIGURE P4.21(B) IF(NORUN.NE.2)THEN RETURN ELSE C... CALL QIKSAX(3,3) C... CALL QIKSET(6.0,0.0,0.1,8.0,-10.0,5.0) C... CALL QIKPLT(TP,ANGLEP1,IP,3H*T*,5H*PHI*,2H**,-101) C... CALL PLOT(-7.0,1.0,-3) C... CALL QLINE(TP,ANGLEP,IP,-11) C... C... ***************************************************************** C... END IF RETURN END KANE, LIKINS AND LEVINSON, SPACECRAFT DYNAMICS, PP. 408-412, CASE (B) 0. 0.6 0.006 1210000 1 1 REL 0.001 KANE, LIKINS AND LEVINSON, SPACECRAFT DYNAMICS, PP. 408-412, CASE (E) 0. 0.6 0.006 1210000 1 1 REL 0.001 END OF RUNS *APP PHYS07 SUBROUTINE INITAL C... C... DYNAMICS OF A DISCRETE-DISTRIBUTED MASS-SPRING SYSTEM C... C... CONSIDER THE FOLLOWING DYNAMICAL SYSTEM* C... C . . G C . . C . A . Y C ............... . C . . + C . . B ---------- C .............. ... ...............................P... + C + . . + . . C . . . . . . C . ............... . . C . . + . . . C Z . S . . W C . . S . . . C . . S . . . C . . S . . . C + . + .+------------- X ------------+. + C ----- . ----- . . ----- C ///// ///// .+--------------- L ---------------+ C C... THE SYSTEM CONSISTS OF A BLOCK, A, A UNIFORM BEAM, B, AND A LIGHT C... SPRING, S. A IS CONSTRAINED TO MOVE IN GUIDE G, ONE END OF B IS C... CLAMPED TO A, AND S CONNECTS A TO A FIXED SUPPORT. C... C... IF W(X,T) DENOTES THE DISPLACEMENT OF A POINT P LOCATED A DISTANCE C... X FROM THE CLAMPED END OF B (W = 0 WHEN BOTH B AND S ARE UNDEFORM- C... ED), EI IS THE FLEXURAL RIGIDITY OF B, KA IS THE MODULUS OF S, C... L IS THE LENGTH OF B, RHO IS THE MASS PER UNIT OF LENGTH B, AND C... MA IS THE MASS OF A, THEN W IS GOVERNED BY THE PARTIAL DIFFEREN- C... TIAL EQUATION (PDE) C... C... EI*W + RHO*W = 0 C... XXXX TT C... C... TOGETHER WITH THE (TIME-DEPENDENT) BOUNDARY CONDITIONS C... C... W (0,T) = W (L,T) = W (L,T) = 0 C... X XX XXX C... C... MA*W (0,T) + EI*W (0,T) + KA*W(0,T) = 0 C... TT XXX C... C... W(X,0) = G(X), W (X,0) = H(X) C... T C... C... TAKING MA = 10 KG, L = 10 M, RHO = 1000 KG/M, EI = 100 N-M**2, C... KA = 394.784 N/M (WHEN KA HAS THIS VALUE, THE SYSTEM FORMED C... BY A AND S ALONE HAS A NATURAL VIBRATION PERIOD OF 1 SEC), AND C... ASSUMING THAT THE SYSTEM STARTS FROM REST WITH B UNDEFORMED AND C... S EXTENDED BY 0.1 M, MAKE TIME-PLOTS, FOR 0 LE T LE 10 SEC, OF C... THE DISPLACEMENT W(L,T) OF THE FREE END OF B. C... C... *KANE, T. R., P. W. LIKINS AND D. A. LEVINSON, SPACECRAFT DY- C... NAMICS, MCGRAW-HILL BOOK COMPANY, NEW YORK, PP. 414-417 C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ W1(21), W2(21) 2 /F/ W1T(21), W2T(21) 3 /S/ W1X(21), W1XX(21), W1XXX(21),W1XXXX(21) 4 /C/ EI, RHO, MA, KA, XL, 5 N, IP REAL MA, KA C... C... SET THE MODEL PARAMETERS EI=100. RHO=1000. MA=10. KA=394.784 XL=10. C... C... INITIAL CONDITIONS N=21 DO 1 I=1,N W1(I)=0.1 W2(I)=0.0 1 CONTINUE C... C... DERIVATIVE W2T IN THE FOURTH BOUNDARY CONDITION IS ALSO INITIALIZ- C... ED TO ZERO W2T(1)=0. C... C... COMPUTE THE INITIAL DERVIATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ W1(21), W2(21) 2 /F/ W1T(21), W2T(21) 3 /S/ W1X(21), W1XX(21), W1XXX(21),W1XXXX(21) 4 /C/ EI, RHO, MA, KA, XL, 5 N, IP REAL MA, KA C... C... COMPUTE THE DERIVATIVE WX CALL DSS004(0.,XL,N,W1,W1X) C... C... BOUNDARY CONDITIONS FOR THE FIRST DERIVATIVE W1X(1)=0. C... C... COMPUTE THE DERIVATIVE UXX CALL DSS004(0.,XL,N,W1X,W1XX) C... C... BOUNDARY CONDITIONS FOR THE SECOND DERIVATIVE W1XX(N)=0. C... C... COMPUTE THE DERIVATIVE WXXX CALL DSS004(0.,XL,N,W1XX,W1XXX) C... C... BOUNDARY CONDITIONS FOR THE THIRD DERIVATIVE W1XXX(N)=0. W1XXX(1)=(1./EI)*(-KA*W1(1)-MA*W2T(1)) C... C... COMPUTE THE DERIVATIVE WXXXX CALL DSS004(0.,XL,N,W1XXX,W1XXXX) C... C... PDE DO 1 I=1,N W1T(I)=W2(I) W2T(I)=-(EI/RHO)*W1XXXX(I) 1 CONTINUE RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ W1(21), W2(21) 2 /F/ W1T(21), W2T(21) 3 /S/ W1X(21), W1XX(21), W1XXX(21),W1XXXX(21) 4 /C/ EI, RHO, MA, KA, XL, 5 N, IP REAL MA, KA C... C... DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION DIMENSION TP(201),W1P(3,201) C... C... INCREMENT THE COUNTER FOR THE PRINTED AND PLOTTED SOLUTION IP=IP+1 C... C... PRINT THE SOLUTION EVERY TWENTIETH CALL TO SUBROUTINE PRINT IF(((IP-1)/20*20).EQ.(IP-1))THEN WRITE(NO,1)T,( W1(I),I=1,N,4),( W1X(I),I=1,N,4), 1 ( W1XX(I),I=1,N,4),( W1XXX(I),I=1,N,4), 2 (W1XXXX(I),I=1,N,4),( W1T(I),I=1,N,4), 3 ( W2T(I),I=1,N,4) 1 FORMAT(1H ,//,' T = ',F6.2,/, 1 ' W1',6E10.3,/, 2 ' W1X',6E10.3,/, 3 ' W1XX',6E10.3,/, 4 ' W1XXX',6E10.3,/, 5 ' W1XXXX',6E10.3,/, 6 ' W1T',6E10.3,/, 7 ' W2T',6E10.3) END IF C... C... STORE THE SOLUTION FOR SUBSEQUENT PLOTTING TP(IP)=T N2=N/2+1 W1P(1,IP)=W1(1) W1P(2,IP)=W1(N2) W1P(3,IP)=W1(N) C... C... PLOT THE SOLUTION AT THE END OF THE RUN IF(IP.LT.201)RETURN CALL TPLOTS(3,IP,TP,W1P) WRITE(NO,2) 2 FORMAT(1H ,/,' 1 - W1(0,T), 2 - W1(L/2,T), 3 - W1(L,T)') C... C... ***************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND IN GENERAL MUST BE CONVERTED FOR THE LOCAL OR COMPUTED OR C... TO COMMENTS C... C... CALL QIKSAX(3,3) C... CALL QIKPLT(TP,W1P,IP,9H*T (SEC)*,24H*W1(X,T), X = 0, L/2, L*, C... 1 2H**,3) C... C... ****************************************************************** C... RETURN END DYNAMICS OF A DISCRETE-DISTRIBUTED MASS-SPRING SYSTEM 0. 10. 0.05 4210000 1 1 REL 0.001 END OF RUNS *APP PHYS08 SUBROUTINE INITAL C... C... NUMERICAL INTEGRATION OF A FOKKER-PLANCK EQUATION C... (VAN DER MEER ANTIPROTON ACCUMULATOR EQUATION) C... C... THE FOLLOWING PARTIAL DIFFERENTIAL EQUATION (PDE) IS DISCUSSED C... BY S. VAN DER MEER (1) WITH APPLICATION TO THE DESIGN OF THE C... ANTIPROTRON ACCUMULATOR (AA) AT CERN C... C... PSI = -(F*PSI) + (D*PSI ) (1) C... T F0 F0 F0 C... WHERE C... C... PSI ANTIPROTON DENSITY C... C... T TIME C... C... F0 REVOLUTION FREQUENCY C... C... F,D SLOWLY VARYING CONSTANTS THAT DEPEND ON VARIOUS C... SYSTEM PARAMETERS AS WELL AS ON THE PARTICLE DIS- C... TRIBUTION C... C... THE SUBSCRIPTS IN T AND F0 DENOTE PARTIAL DERIVATIVES. C... C... (1) VAN DER MEER, S., STOCHASTIC COOLING AND THE ACCUMULATION C... OF ANTIPROTONS, SCIENCE, VOL. 230, PP 900-906, 15 NOVEMBER C... 1985. C... C... THIS PAPER CONTAINED THE ADDITIONAL CITATION C... C... S. VAN DER MEER IS A SENIOR ENGINEER AT CERN (THE EUROPEAN C... LABORATORY FOR PARTICLE PHYSICS), F-01631 CERN CEDEX, CH- C... 1211, GENEVA 23, SWITZERLAND. THIS ARTICLE IS THE LECTURE C... HE DELIVERED IN STOCKHOLM ON 8 DECEMBER 1984 WHEN HE RE- C... CEIVED THE NOBEL PRIZE IN PHYSICS, WHICH HE SHARED WITH C... CARLO RUBBIA. C... C... EQUATION (1) IS AN EXAMPLE OF A FOKKER-PLANCK EQUATION. THE C... INITIAL AND BOUNDARY CONDITIONS WILL BE TAKEN AS C... C... PSI(0,T) = 1, PSI(1,T) = 0, PSI(F0,0) = G(F0) (2)(3)(4) C... F0 C... WHERE C... C... G(F0) INITIAL DISTRIBUTION F0R PSI IN EQUATION (1). C... C... F0R THE PRESENT PROBLEM, TAKE G(F0) = 1 - 2*F0 + F0**2, WHICH C... IS CONSISTENT WITH BOTH BOUNDARY CONDITIONS (2) AND (3). C... C... THE FOLLOWING PROGRAMMING IS FOR A NUMERICAL SOLUTION TO EQUA- C... TIONS (1) TO (4) WITH F = D = 1. C... C AFTER FORMULATION OF THE PRECEDING PROBLEM, THE FOLLOWING LETTER C WAS RECEIVED FROM DR. SIMON VAN DER MEER IN REPLY TO A LETTER C OF INQUIRY ABOUT THE INITIAL AND BOUNDARY CONDITIONS FOR THE C FOKKER-PLANCK EQUATION: C C EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH C EUROPEAN LABORATORY FOR PARTICLE PHYSICS C C CERN DR. W. E. SCHIESSER C CH - 1211 GENEVEE 23 WHITAKER NO. 5 C CERN SITE DE PREVESSIN LEHIGH UNIVERSITY C F - 01631 CERN CEDEX BETHLEHEM, PA 18015 USA C C GENEVA, 9 DECEMBER 1985 C C DEAR DR. SCHIESSER: C C IN REPLY TO YOUR LETTER OF DECEMBER 2, THE FOLLOWING DETAILS. C C EQUATION (6) WAS SOLVED BY DIVIDING THE RANGE OF REVOLUTION C FREQUENCIES F0 COVERED BY THE STACK INTO 50 BINS AND SETTING UP C A DIFFERENCE EQUATION FOR EACH BIN THAT DESCRIBES THE DENSITY C PSI VS TIME. THE DIFFERENTIATION WITH RESPECT TO F0 IS REPLACED C BY CALCULATING DIFFERENCES BETWEEN ADJOINING BINS. THE BEAM- C FEEDBACK EFFECT, WHICH DEPENDS ON THE DENSITY PROFILE, AND THE C PARAMETERS F AND D, ARE NOT CALCULATED FOR EACH TIME STEP, BUT C AT LONGER INTERVALS, SINCE THEY CHANGE SLOWLY. C C THE 50 SIMULTANEOUS FIRST-ORDER DIFFERENCE EQUATIONS WERE SOLVED C BY A RATHER PRIMITIVE ALGORITHM, TO SAVE COMPUTER TIME. THE C INITIAL CONDITION WAS ZERO DENSITY EVERYWHERE. ONE BOUNDARY C CONDITION I ASSUMED WAS THAT THE DENSITY AT A POINT BEYOND THE C LAST BIN (AT THE RIGHT HAND SIDE IN FIG. 11) WAS ZERO, SO THAT C ANY PARTICLES ARRIVING THERE WOULD DISAPPEAR. (IN FACT, NO C PARTICLES EVER GOT THERE). THE OTHER BOUNDARY, AT THE LOW- C FREQUENCY STACK TAIL, WAS GIVEN BY A CONSTANT FLUX ENTERING C THE STACK. THE 'PULSED' NATURE OF THIS FLUX WAS CONSIDERED IN C ANOTHER, INDEPENDENT, CALCULATION. C C THE ACTUAL COMPUTER PROGRAM IS EXTREMELY COMPLEX BECAUSE OF ALL C THE DETAILS CONCERNING FILTERS, PICK-UP RESPONSE, ETC., BUT ALSO C BECAUSE THERE IS A PARTICULAR PROBLEM WITH SPEED. THIS IS BE- C CAUSE IN THE LOW-DENSITY STACK TAIL THE COOLING IS VERY FAST C WHILE IN THE CORE IT IS SEVERAL ORDERS OF MAGNITUDE SLOWER. C THEREFORE, THE STEPS IN TIME FOR INTEGRATING THE TAIL BINS MUST C BE QUITE SHORT, EVEN IF THE DENSITY PROFILE THERE DOES NOT CHANGE C MUCH DURING MOST OF THE TIME; OTHERWISE THE SOLUTION GETS UN- C STABLE. ON THE OTHER HAND, FOR THE CORE MUCH LARGER STEPS ARE C USED, ADJUSTED AS A FUNCTION OF THE DENSITY, TO OBTAIN A REASON- C ABLE SPEED OF EXECUTION. C C I AM AFRAID THERE IS NO PUBLICATION ABOUT THE ACTUAL COMPUTER C PROGRAM. I COULD SEND YOU A LISTING, BUT I AM SURE IT WOULD C TAKE A LOT OF EXPLICATION TO UNDERSTAND HOW IT WORKS. IT IS C ONE OF THOSE PROGRAMS THAT HAVE EVOLVED SLOWLY, ACCUMULATING C ALL KINDS OF DETAILS (E.G., PLOTTING). IN FACT, I USED TWO C DIFFERENT PROGRAMS, ONE THAT WORKED IN AN INTERACTIVE MODE AND C ALLOWED ME TO QUICKLY CHANGE MANY PARAMETERS, BUT THAT ONLY C CALCULATED THE FLUX VS FO FOR A GIVEN DENSITY PROFILE, AND THE C OTHER ONE AS DESCRIBED ABOVE, VERIFYING THE WHOLE SYSTEM BY C CALCULATING DENSITY PROFILES VS TIME. C C THE LAST PROGRAM WAS 'EXPORTED' TO FERMILAB WHERE IT WAS USED C (AND FURTHER IMPROVED), I BELIEVE BY JOHN MARRINER, FOR THE C ANTIPROTON SOURCE THEY HAVE BUILT THERE. C C I HOPE THIS ANSWERS SOME OF YOUR QUESTIONS. THANK YOUR FOR YOUR C CONGRATULATIONS. I AM INDEED AN ENGINEER BY ORIGIN (FROM DELFT, C THE NETHERLANDS). C C YOURS SINCERELY, C C SIMON VAN DER MEER C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ PSI(21) 2 /F/ PSIT(21) 3 /S/ PSIF0(21),PSIFF0(21), F0(21) 4 /C/ F, D 5 /I/ N, IP C... C... SET THE MODEL PARAMETERS F=1.0 D=1.0 C... C... INITIAL CONDITION (4) N=21 DO 1 I=1,N F0(I)=FLOAT(I-1)/FLOAT(N-1) PSI(I)=1.0-2.0*F0(I)+F0(I)**2 1 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ PSI(21) 2 /F/ PSIT(21) 3 /S/ PSIF0(21),PSIFF0(21), F0(21) 4 /C/ F, D 5 /I/ N, IP C... C... BOUNDARY CONDITION (2) PSI(1)=1. C... C... COMPUTE THE DERIVATIVE PHI C... F0 CALL DSS004(0.,1.,N,PSI,PSIF0) C... C... BOUNDARY CONDITION (3) PSIF0(N)=0. C... C... COMPUTE THE DERIVATIVE PHI C... F0F0 CALL DSS004(0.,1.,N,PSIF0,PSIFF0) C... C... EQUATION (1) DO 2 I=1,N PSIT(I)=-F*PSIF0(I)+D*PSIFF0(I) 2 CONTINUE C... C... FROM BOUNDARY CONDITION (2), PHI(0,T) IS INVARIANT WITH T PSIT(1)=0. RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ PSI(21) 2 /F/ PSIT(21) 3 /S/ PSIF0(21),PSIFF0(21), F0(21) 4 /C/ F, D 5 /I/ N, IP C... C... DIMENSION THE ARRAYS FOR PLOTTING DIMENSION PSIP(6,21) C... C... PRINT THE SOLUTION NS=4 WRITE(NO,1)T,(F0(I),I=1,N,NS),(PSI(I),I=1,N,4) 1 FORMAT(1H ,//,' T = ',F5.2,/, 1 ' F0',6F11.3,/, 2 ' PSI',6F11.4) C... C... STORE THE NUMERICAL SOLUTION FOR PLOTTING IP=IP+1 DO 2 I=1,N PSIP(IP,I)=PSI(I) 2 CONTINUE C... C... PLOT THE SOLUTION IF(IP.LT.6)RETURN CALL TPLOTS(IP,N,F0,PSIP) WRITE(NO,10) 10 FORMAT(1H ,//, 1 ' PHI(F0,T) VS F0 WITH T AS A PARAMETER') RETURN END VAN DER MEER ANTIPROTON ACCUMULATOR EQUATION 0. 2.0 0.4 21 1000 1 1 REL 0.001 END OF RUNS *APP PHYS09 SUBROUTINE INITAL C... C... EINSTEIN EQUATION FOR THE SCALE FACTOR OF THE UNIVERSE C... C... THE FOLLOWING ORDINARY DIFFERENTIAL EQUATION (ODE) WAS PROPOSED C... BY EINSTEIN FOR THE SCALE FACTOR, A(T), OF THE EXPANSION OF THE C... UNIVERSE (1) C... C... 2 2 C... D A/DT = 8*PI*G*RHO*A/3 (1) C... C... WHERE C... C... A(T) SCALE FACTOR C... C... T TIME C... C... G NEWTON*S GRAVITATIONAL CONSTANT C... C... RHO VACUUM DENSITY C... C... (1) LAKE, GEORGE, WINDOWS ON A NEW COSMOLOGY, SCIENCE, V. 224, C... NO. 4650, PP 675-682, 18 MAY 85 C... C... EQUATION (1) IS INTEGRATED NUMERICALLY SUBJECT TO THE FOLLOWING C... INITIAL CONDITIONS FOR THE CURRENT EPOCH C... C... A(0) = DA(0)/DT = 1 (2)(3) C... C... AN ANALYTICAL SOLUTION TO EQUATIONS (1) TO (3) CAN EASILY BE C... DERIVED. IF A SOLUTION OF THE FORM C... C... A(T) = D*EXP(C*T) C... C... IS ASSUMED, WHERE C = (8*PI*G*RHO/3)**(1/2), SUBSTUTION IN C... EQUATION (1) CONFIRMS THE SOLUTION. THUS, THE GENERAL SOLUTION C... TO EQUATION (1) IS C... C... A(T) = D1*EXP(C*T) + D2*EXP(-C*T) C... C... WHERE D1 AND D2 ARE CONSTANTS TO BE DETERMINED FROM EQUATIONS (2) C... AND (3). APPLYING THESE INITIAL CONDITIONS, WE HAVE C... C... A(0) = D1 + D2 = 1 C... C... DA(0)/DT = C*D1 - C*D2 = 1 C... C... THUS, D2 = 1 - D1 AND C... C... C*D1 - C*(1 - D1) = 1 C... C... 2*C*D1 = 1 + C C... C... D1 = (1 + C)/(2*C) C... C... AND C... C... D2 = 1 - (1 + C)/(2*C) = (C - 1)/(2*C) C... C... THE SOLUTION IS THEREFORE C... C... A(T) = (1 + C)/(2*C)*EXP(C*T) + (C - 1)/(2*C)*EXP(-C*T) (4) C... C... THE SOLUTION IS VERIFIED AS C... C... 2 2 C... D A/DT = (C/2)*(1 + C)*EXP(C*T) - (C/2)*(1 - C)*EXP(-C*T) C... C... -(C**2)*A = -(C/2)*(1 + C)*EXP(C*T) + (C/2)*(1 - C)*EXP(-C*T) C... C... A(0) = (1 + C)/(2*C) - (1 - C)/(2*C) = 1 C... C... DA(0)/DT = (1/2)*(1 + C) + (1/2)*(1 - C) = 1 C... C... EQUATION (4) CAN THEN BE USED TO CHECK THE NUMERICAL SOLUTION. C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ A, AT 2 /F/ DADT, ATT 3 /C/ C, IP C... C... SET THE MODEL PARAMETERS C=1. C... C... INITIAL CONDITIONS A=1.0 AT=1.0 C... C... COMPUTE THE INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ A, AT 2 /F/ DADT, ATT 3 /C/ C, IP C... C... EQUATION (1) IS PROGRAMMED AS TWO FIRST-ORDER ODES DADT=AT ATT=SQRT(C)*A RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ A, AT 2 /F/ DADT, ATT 3 /C/ C, IP C... C... PRINT A HEADING FOR THE OUTPUT IP=IP+1 IF(IP.EQ.1)THEN WRITE(NO,2) 2 FORMAT(1H ,//,9X,'T',11X,'A(T)',6X,'ANAL A(T)') END IF C... C... COMPUTE THE ANALYTICAL SOLUTION AE=(1.0+C)/(2.0*C)*EXP(C*T)-(1.0-C)/(2.0*C)*EXP(-C*T) C... C... PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS WRITE(NO,1)T,A,AE 1 FORMAT(F10.3,2F15.4) RETURN END EINSTEIN EQUATION FOR THE SCALE FACTOR OF THE UNIVERSE 0. 1.0 0.1 2 1000 1 1 REL 0.0001 END OF RUNS *APP PHYS10 SUBROUTINE INITAL C... C... TIME-DEPENDENT SCHRODINGER EQUATION IN TWO DIMENSIONS C... C... THE TIME-DEPENDENT SCHRODINGER EQUATION IN TWO DIMENSIONS IS C... C... I*U = U + U - K*P(X,Y,T)*U C... T XX YY C... C... WHERE C... C... I = (-1)**(1/2) C... C... SINCE U IS COMPLEX, WE SUBSTITUTE U = V + I*W IN EQUATION (1) TO C... OBTAIN C... C... I*(V + I*W ) = V + I*W + V + I*W - K*P*(V + I*W) C... T T XX XX YY YY C... C... WHICH CAN BE SEPARATED INTO REAL AND IMAGINARY PARTS C... C... -1*W = V + V -K*P*V C... T XX YY C... C... V = W + W -K*P*W C... T XX YY C... C... THE BOUNDARY AND INITIAL CONDITIONS FOR EQUATION (1) WILL BE TAKEN C... AS C... C... U (0,Y,T) = U (1,Y,T) = 0 C... X X C... C... U (X,0,T) = U (X,2,T) = 0 C... Y Y C... C... U(X,Y,0) = (X*(1 - X))**2 + (Y*(2 - Y))**2 C... C... WHICH CAN ALSO BE SEPARATED INTO REAL AND IMAGINARY PARTS C... C... V (0,Y,T) = V (1,Y,T) = W (0,Y,T) = W (1,Y,T) = 0 C... X X X X C... C... V (X,0,T) = V (X,2,T) = W (X,0,T) = W (X,2,T) = 0 C... Y Y Y Y C... C... V(X,Y,0) = (X*(1 - X))**2 + (Y*(2 - Y))**2, W(X,Y,0) = 0 C... C... ALSO, WE WILL TAKE K*P = 1 C... C... A NUMERICAL METHOD OF LINES SOLUTION TO THIS PROBLEM IS PROGRAMMED C... USING AN 11 X 11 POINT X-Y GRID. C... C... THE NUMERICAL PROCEDURE USED IN SUBROUTINE DERV IS TO REPLACE C... THE SECOND-ORDER DERIVATIVES IN X AND Y WITH SECOND-ORDER, C... THREE-POINT CENTERED APPROXIMATIONS. THE DERIVATIVES IN THE C... BOUNDARY CONDITIONS ARE REPLACED WITH SECOND-ORDER, TWO-POINT C... CENTERED APPROXIMATIONS WHICH ARE THEN USED AT THE BOUNDARIES TO C... ELIMINATE FICTITIOUS VALUES OF V(X,Y,T) AND W(X,Y,T) OUTSIDE THE C... BOUNDARIES. C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ V(11,11), W(11,11) 2 /F/ VT(11,11), WT(11,11) 3 /C/ PK, 4 XL, YL, 5 DX, DY, 6 DXS, DYS, 7 NX, NY, 8 IP C... C... SET THE PARAMETERS C... C... X AND Y DIMENSIONS XL=1. YL=2. C... C... X AND Y GRID POINTS NX=11 NY=11 C... C... K*P(X,Y,T) PK=1. C... C... SPATIAL INCREMENTS DX=XL/FLOAT(NX-1) DY=YL/FLOAT(NY-1) DXS=DX**2 DYS=DY**2 C... C... SET THE INITIAL CONDITIONS DO 1 I=1,NX X=FLOAT(I-1)/FLOAT(NX-1)*XL DO 1 J=1,NY Y=FLOAT(J-1)/FLOAT(NY-1)*YL V(I,J)=(X*(1.-X))**2+(Y*(2.-Y))**2 W(I,J)=0. 1 CONTINUE C... C... COMPUTE THE INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ V(11,11), W(11,11) 2 /F/ VT(11,11), WT(11,11) 3 /C/ PK, 4 XL, YL, 5 DX, DY, 6 DXS, DYS, 7 NX, NY, 8 IP COMMON/IO/ NI, NO C... C... INTERIOR POINTS DO 1 I=2,NX-1 DO 1 J=2,NY-1 VT(I,J)= ( W(I+1,J )-2.*W(I ,J )+ W(I-1,J ))/DXS 1 +( W(I ,J+1)-2.*W(I ,J )+ W(I ,J-1))/DYS 2 -PK*W(I ,J ) WT(I,J)=-( V(I+1,J )-2.*V(I ,J )+ V(I-1,J ))/DXS 1 -( V(I ,J+1)-2.*V(I ,J )+ V(I ,J-1))/DYS 2 +PK*V(I ,J ) 1 CONTINUE C... C... BOUNDARY POINTS C... C... Y = 0, X NE 0 OR XL DO 2 I=2,NX-1 VT(I,1)= ( W(I+1, 1)-2.*W(I , 1)+ W(I-1, 1))/DXS 1 +(2.*W(I , 2)-2.*W(I , 1) )/DYS 2 -PK*W(I , 1) WT(I,1)=-( V(I+1, 1)-2.*V(I , 1)+ V(I-1, 1))/DXS 1 -(2.*V(I , 2)-2.*V(I , 1) )/DYS 2 +PK*V(I , 1) 2 CONTINUE C... C... Y = YL, X NE 0 OR XL N=NY DO 3 I=2,NX-1 VT(I,N)= ( W(I+1,N )-2.*W(I ,N )+ W(I-1,N ))/DXS 1 +(2.*W(I ,N-1)-2.*W(I ,N ) )/DYS 2 -PK*W(I ,N ) WT(I,N)=-( V(I+1,N )-2.*V(I ,N )+ V(I-1,N ))/DXS 1 -(2.*V(I ,N-1)-2.*V(I ,N ) )/DYS 2 +PK*V(I ,N ) 3 CONTINUE C... C... X = 0, Y NE 0 OR YL DO 4 J=2,NY-1 VT(1,J)= (2.*W( 2,J )-2.*W( 1,J ) )/DXS 1 +( W( 1,J+1)-2.*W( 1,J )+ W( 1,J-1))/DYS 2 -PK*W( 1,J ) WT(1,J)=-(2.*V( 2,J )-2.*V( 1,J ) )/DXS 1 -( V( 1,J+1)-2.*V( 1,J )+ V( 1,J-1))/DYS 2 +PK*V( 1,J ) 4 CONTINUE C... C... X = XL, Y NE 0 OR YL N=NX DO 5 J=2,NY-1 VT(N,J)= (2.*W(N-1,J )-2.*W(N ,J ) )/DXS 1 +( W(N ,J+1)-2.*W(N ,J ) +W(N ,J-1))/DYS 2 -PK*W(N ,J ) WT(N,J)=-(2.*V(N-1,J )-2.*V(N ,J ) )/DXS 1 -( V(N ,J+1)-2.*V(N ,J ) +V(N ,J-1))/DYS 2 +PK*V(N ,J ) 5 CONTINUE C... C... CORNER POINTS C... C... X = 0 Y = 0 VT(1,1)= (2.*W( 2, 1)-2.*W( 1, 1) )/DXS 1 +(2.*W( 1, 2)-2.*W( 1, 1) )/DYS 2 -PK*W( 1, 1) WT(1,1)=-(2.*V( 2, 1)-2.*V( 1, 1) )/DXS 1 -(2.*V( 1, 2)-2.*V( 1, 1) )/DYS 2 +PK*V( 1, 1) C... C... X = XL, Y = 0 N=NX VT(N,1)= (2.*W(N-1, 1)-2.*W(N , 1) )/DXS 1 +(2.*W(N , 2)-2.*W(N , 1) )/DYS 2 +PK*W(N , 1) WT(N,1)=-(2.*V(N-1, 1)-2.*V(N , 1) )/DXS 1 -(2.*V(N , 2)-2.*V(N , 1) )/DYS 2 -PK*V(N , 1) C... C... X = 0, Y = YL N=NY VT(1,N)= (2.*W( 2,N )-2.*W( 1,N ) )/DXS 1 +(2.*W( 1,N-1)-2.*W( 1,N ) )/DYS 2 -PK*W( 1,N ) WT(1,N)=-(2.*V( 2,N )-2.*V( 1,N ) )/DXS 1 -(2.*V( 1,N-1)-2.*V( 1,N ) )/DYS 2 +PK*V( 1,N ) C... C... X = XL, Y = YL N=NX N=NY VT(N,N)= (2.*W(N-1,N )-2.*W(N ,N ) )/DXS 1 +(2.*W(N ,N-1)-2.*W(N ,N ) )/DYS 2 -PK*W(N ,N ) WT(N,N)=-(2.*V(N-1,N )-2.*V(N ,N ) )/DXS 1 -(2.*V(N ,N-1)-2.*V(N ,N ) )/DYS 2 +PK*V(N ,N ) C... C... THE FOLLOWING WRITE STATEMENT WAS USED TO CHECK FOR ERRORS IN C... SUBROUTINE DERV. IT WAS RETAINED AS COMMENTS TO ILLUSTRATE HOW C... TO LOOK FOR ERRORS IN THE PROGRAMMING OF A RELATIVELY LARGE SET OF C... ODES. NOTE THAT THE OUTPUT UNIT NUMBER, NO, IS ACCESSED THROUGH C... COMMON/IO/ C WRITE(NO,6)((I,J,VT(I,J),WT(I,J),I=1,NX),J=1,NY) C 6 FORMAT(2I10,2E12.3) RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ V(11,11), W(11,11) 2 /F/ VT(11,11), WT(11,11) 3 /C/ PK, 4 XL, YL, 5 DX, DY, 6 DXS, DYS, 7 NX, NY, 8 IP C... C... DIMENSION THE ARRAYS FOR PLOTTING DIMENSION TP(101), YP(2,101) C... C... PRINT THE SOLUTION EVERY 20TH CALL TO SUBROUTINE PRINT IP=IP+1 IF(((IP-1)/20*20).EQ.(IP-1))THEN C... C... PRINT V(X,Y,T) WRITE(NO,1)T 1 FORMAT(1H ,//,' T = ',F5.3,/) WRITE(NO,2) 2 FORMAT(' V(X,Y,T)',/) WRITE(NO,3)((V(I,J),I=1,NX,2),J=1,NY,2) 3 FORMAT(' Y=0.0',6F10.4,/, 1 ' Y=0.4',6F10.4,/, 2 ' Y=0.8',6F10.4,/, 3 ' Y=1.2',6F10.4,/, 4 ' Y=1.6',6F10.4,/, 5 ' Y=2.0',6F10.4) WRITE(NO,4) 4 FORMAT(1H ,9X, 1 ' X=0 X=0.2 X=0.4 X=0.6 X=0.8 X=1.0', 2 //) C... C... PRINT W(X,Y,T) WRITE(NO,5) 5 FORMAT(' W(X,Y,T)',/) WRITE(NO,3)((W(I,J),I=1,NX,2),J=1,NY,2) WRITE(NO,4) C... C... PRINT VT(X,Y,T) WRITE(NO,6) 6 FORMAT(' VT(X,Y,T)',/) WRITE(NO,3)((VT(I,J),I=1,NX,2),J=1,NY,2) WRITE(NO,4) C... C... PRINT WT(X,Y,T) WRITE(NO,7) 7 FORMAT(' WT(X,Y,T)',/) WRITE(NO,3)((WT(I,J),I=1,NX,2),J=1,NY,2) WRITE(NO,4) END IF C... C... STORE THE SOLUTION FOR PLOTTING TP(IP)=T YP(1,IP)=V(NX/2+1,NY/2+1) YP(2,IP)=W(NX/2+1,NY/2+1) C... C... PLOT THE SOLUTION IF(IP.LT.101)RETURN CALL TPLOTS(2,IP,TP,YP) WRITE(NO,8) 8 FORMAT(1H ,/, 1 ' 1 - V(XL/2,YL/2,T) 2 - W(XL/2,YL/2,T) VS T') C... C... PRINT THE PLOTTED SOLUTION WRITE(NO,9) 9 FORMAT(1H1,7X,'T',8X,'V(0.5,1,T)',5X,'W(0.5,1,T)') WRITE(NO,10)(TP(I),YP(1,I),YP(2,I),I=1,IP,4) 10 FORMAT(F10.3,2F15.5) RETURN END TIME-DEPENDENT SCHRODINGER EQUATION IN TWO DIMENSIONS 0. 0.1 0.001 242 1000 3 1 REL 0.001 END OF RUNS *APP PHYS11 SUBROUTINE INITAL C... C... DYNAMICS OF A ROCKET C... C... THE FOLLOWING EQUATIONS HAVE BEEN PROPOSED FOR THE VELOCITY OF C... A ROCKET, VR (1) C... C... DVR/DT + RHO*CD*AR/(2*(M0 - MDOT*T))*VR**2 C... (1) C... - TH*GC/(M0 - MDOT*T) + G = 0 C... C... T = ((PE - P0) + RHO*VR**2)*AE (2) C... C... WHERE C... C... VR VELOCITY OF THE ROCKET C... C... T TIME C... C... RHO AIR DENSITY C... C... CD DRAG COEFFICIENT C... C... AR ROCKET CROSS SECTIONAL AREA C... C... M0 INITIAL ROCKET MASS C... C... MDOT RATE OF CONSUMPTION OF FUEL C... C... TH ROCKET THRUST C... C... GC CONVERSION FACTOR BETWEEN FORCE AND MASS C... C... (1) KREIDER, JAN F., PRINCIPLES OF FLUID MECHANICS, ALLYN AND C... BACON, INC., NEW YORK, PP 50 - 58 C... C... THE ROCKET POSITION, XR, CAN ALSO BE CALCULATED FROM THE EQUATION C... C... DXR/DT = VR (3) C... C... THE FOLLOWING INITIAL CONDITIONS WILL BE USED FOR EQUATIONS (1) C... AND (3) C... C... XR(0) = VR(0) = 0 (4)(5) C... C... THE SOLUTION TO BE COMPUTED IS XR(T) AND VR(T). C... C... COMMON AREA COMMON/T/ T, NSTOP, NORUN 1 /Y/ VR, XR 2 /F/ DVRDT, DXRDT 3 /C/ RHO, CD, AR, M0, MDOT, 4 TH, G, IP REAL M0, MDOT C... C... SET THE MODEL PARAMETERS C... C... AIR DENSITY RHO=0.075 C... C... DRAG COEFFICIENT CD=0.4 C... C... ROCKET FRONTAL AREA AR=20. C... C... INITIAL ROCKET MASS M0=100000. C... C... FUEL BURN RATE MDOT=18500./60. C... C... THRUST TH=800000. C... C... GRAVITATIONAL ACCELERATION G=32.2 C... C... INITIAL CONDITIONS VR=0. XR=0. C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ VR, XR 2 /F/ DVRDT, DXRDT 3 /C/ RHO, CD, AR, M0, MDOT, 4 TH, G, IP REAL M0, MDOT C... C... MODEL ALGEBRA C... C... CHECK IF ALL OF THE SOLID FUEL HAS BEEN CONSUMED. IF SO, THE C... RUN IS TERMINATED IF((M0-MDOT*T).LE.0.)NSTOP=1 C... C... MODEL DIFFERENTIAL EQUATIONS C... C... VELOCITY EQUATION DVRDT=-(RHO*CD*AR)/(2.0*(M0-MDOT*T))*(VR**2) 1 +TH*32.2/(M0-MDOT*T)-G C... C... POSITION EQUATION DXRDT=VR RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ VR, XR 2 /F/ DVRDT, DXRDT 3 /C/ RHO, CD, AR, M0, MDOT, 4 TH, G, IP REAL M0, MDOT C... C... DIMENSION THE ARRAYS FOR PLOTTING DIMENSION TP(51), VRP(51), XRP(51), VRTP(51), XRTP(51) C... C... PRINT A HEADING FOR THE SOLUTION IF(IP.EQ.0)WRITE(NO,1) 1 FORMAT(1H ,//,9X,'T',10X,'VR',10X,'XR',6X,'DVR/DT',6X,'DXR/DT') C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,2)T,VR,XR,DVRDT,DXRDT 2 FORMAT(F10.1,4F12.1) C... C... STORE THE SOLUTION FOR PLOTTING IP=IP+1 TP(IP)=T VRP(IP)=VR XRP(IP)=XR VRTP(IP)=DVRDT XRTP(IP)=DXRDT C... C... PLOT THE NUMERICAL SOLUTION IF(IP.LT.51)RETURN CALL TPLOTS(1,IP,TP,VRP) WRITE(NO,10) 10 FORMAT(1H ,//,' VR(T) (FT/SEC) VS T (SEC)') CALL TPLOTS(1,IP,TP,XRP) WRITE(NO,11) 11 FORMAT(1H ,//,' XR(T) (FT) VS T (SEC)') CALL TPLOTS(1,IP,TP,VRTP) WRITE(NO,12) 12 FORMAT(1H ,//,' DVR/DT (FT/SEC**2) VS T (SEC)') CALL TPLOTS(1,IP,TP,XRTP) WRITE(NO,13) 13 FORMAT(1H ,//,' DXR/DT (FT/SEC) VS T (SEC)') IP=0 RETURN END DYNAMICS OF A ROCKET 0. 200. 4.0 2 1000 1 1 REL 0.001 END OF RUNS *APP PHYS12 PROGRAM EULER C... C... EULER INTEGRATION OF THE DIFFERENTIAL EQUATIONS FOR THE DYNAMICS C... OF A ROCKET C... C... THE FOLLOWING PROGRAM ILLUSTRATES THE EULER INTEGRATION OF A C... SYSTEM OF TWO ORDINARY DIFFERENTIAL EQUATIONS (ODES) WHICH MODEL C... THE DYNAMICS OF A ROCKET. THE DETAILS OF THE CODING ARE GIVEN AS C... COMMENTS IN SUBROUTINE INITAL. C... C... THE MODEL INITIAL CONDITIONS ARE SET IN SUBROUTINE INITAL, AND C... THE MODEL DERIVATIVES ARE PROGRAMMED IN SUBROUTINE DERV. THE C... NUMERICAL SOLUTION IS PRINTED IN SUBROUTINE PRINT. C... C... THE FOLLOWING COMMON AREA LINKS THIS MAIN PROGRAM TO SUBROUTINES C... INITAL, DERV AND PRINT. COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y(15) 2 /F/ F(15) C... C... THE FOLLOWING COMMON AREA PROVIDES THE INPUT/OUTPUT UNIT NUMBERS C... FOR USE IN OTHER SUBROUTINES COMMON/IO/ NI, NO C... C... DEFINE THE INPUT/OUTPUT UNIT NUMBERS NI=5 NO=6 C... C... OPEN INPUT AND OUTPUT FILES OPEN(NI,FILE='INPUT') OPEN(NO,FILE='OUTPUT') C... C... SET THE NUMBER OF RUNS NORUNS=5 C... C... STEP THROUGH THE NORUNS RUNS DO 4 NORUN=1,NORUNS C... C... INITIALIZE THE RUN TERMINATION VARIABLE NSTOP=0 C... C... SET THE INITIAL, FINAL, PRINT INTERVAL VALUES OF TIME T0=0. IF(NORUN.EQ.1)TF=160. IF(NORUN.EQ.2)TF=80. IF(NORUN.EQ.3)TF=40. IF(NORUN.EQ.4)TF=20. IF(NORUN.EQ.5)TF=10. TP=8. C... C... SET THE NUMBER OF DIFFERENTIAL EQUATIONS NEQN=2 C... C... INITIALIZE TIME T=T0 C... C... SET THE INITIAL CONDITIONS CALL INITAL C... C... PRINT THE INITIAL CONDITIONS CALL PRINT(NI,NO) C... C... SET THE INTEGRATION INTERVAL DT=TF/40. C... C... BEGIN THE EULER INTEGRATION 3 TOUT=T+TP C... C... CALCULATE THE DERIVATIVES AT THE BASE POINT 1 CALL DERV C... C... CHECK FOR A RUN TERMINATION IF(NSTOP.NE.0)GO TO 4 C... C... TAKE AN EULER STEP DO 2 J=1,NEQN Y(J)=Y(J)+F(J)*DT 2 CONTINUE T=T+DT C... C... ALL NEQN DEPENDENT VARIABLES HAVE BEEN ADVANCED, SO TAKE THE C... NEXT STEP IN TIME IF(T.LT.(TOUT-0.5*DT))GO TO 1 C... C... ONE PRINT INTERVAL IS COMPLETE, SO PRINT THE SOLUTION CALL PRINT(NI,NO) C... C... CHECK FOR A RUN TERMINATION IF(NSTOP.NE.0)GO TO 4 C... C... CHECK FOR THE END OF THE RUN IF(T.LT.(TF-0.5*TP))GO TO 3 C... C... THE CURRENT RUN IS COMPLETE, SO GO ON TO THE NEXT RUN 4 CONTINUE C... C... ALL OF THE RUNS ARE COMPLETE, SO TERMINATE EXECUTION STOP END SUBROUTINE INITAL C... C... DYNAMICS OF A ROCKET C... C... THE FOLLOWING EQUATIONS HAVE BEEN PROPOSED FOR THE VELOCITY OF C... A ROCKET, VR (1) C... C... DVR/DT + RHO*CD*AR/(2*(M0 - MDOT*T))*VR**2 C... (1) C... - TH*GC/(M0 - MDOT*T) + G = 0 C... C... T = ((PE - P0) + RHO*VR**2)*AE (2) C... C... WHERE C... C... VR VELOCITY OF THE ROCKET C... C... T TIME C... C... RHO AIR DENSITY C... C... CD DRAG COEFFICIENT C... C... AR ROCKET CROSS SECTIONAL AREA C... C... M0 INITIAL ROCKET MASS C... C... MDOT RATE OF CONSUMPTION OF FUEL C... C... TH ROCKET THRUST C... C... GC CONVERSION FACTOR BETWEEN FORCE AND MASS C... C... (1) KREIDER, JAN F., PRINCIPLES OF FLUID MECHANICS, ALLYN AND C... BACON, INC., NEW YORK, PP 50 - 58, PROB C2.5 C... C... THE ROCKET POSITION, XR, CAN ALSO BE CALCULATED FROM THE EQUATION C... C... DXR/DT = VR (3) C... C... THE FOLLOWING INITIAL CONDITIONS WILL BE USED FOR EQUATIONS (1) C... AND (3) C... C... XR(0) = VR(0) = 0 (4)(5) C... C... THE SOLUTION TO BE COMPUTED IS XR(T) AND VR(T). C... C... COMMON AREA COMMON/T/ T, NSTOP, NORUN 1 /Y/ VR, XR 2 /F/ DVRDT, DXRDT 3 /C/ RHO, CD, AR, M0, MDOT, 4 TH, G, IP REAL M0, MDOT C... C... SET THE MODEL PARAMETERS C... C... AIR DENSITY RHO=0.075 C... C... DRAG COEFFICIENT CD=0.4 C... C... ROCKET FRONTAL AREA AR=20. C... C... INITIAL ROCKET MASS M0=100000. C... C... FUEL BURN RATE MDOT=18500./60. C... C... THRUST TH=800000. C... C... GRAVITATIONAL ACCELERATION G=32.2 C... C... INITIAL CONDITIONS VR=0. XR=0. C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ VR, XR 2 /F/ DVRDT, DXRDT 3 /C/ RHO, CD, AR, M0, MDOT, 4 TH, G, IP REAL M0, MDOT C... C... MODEL ALGEBRA C... C... CHECK IF ALL OF THE SOLID FUEL HAS BEEN CONSUMED. IF SO, THE C... RUN IS TERMINATED IF((M0-MDOT*T).LE.0.)NSTOP=1 C... C... MODEL DIFFERENTIAL EQUATIONS C... C... VELOCITY EQUATION DVRDT=-(RHO*CD*AR)/(2.0*(M0-MDOT*T))*(VR**2) 1 +TH*32.2/(M0-MDOT*T)-G C... C... POSITION EQUATION DXRDT=VR RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ VR, XR 2 /F/ DVRDT, DXRDT 3 /C/ RHO, CD, AR, M0, MDOT, 4 TH, G, IP REAL M0, MDOT C... C... PRINT A HEADING FOR THE SOLUTION IF(IP.EQ.0)WRITE(NO,1) 1 FORMAT(1H ,//,9X,'T',10X,'VR',10X,'XR',6X,'DVR/DT',6X,'DXR/DT') C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,2)T,VR,XR,DVRDT,DXRDT 2 FORMAT(F10.1,4F12.1) IP=IP+1 RETURN END *APP PHYS13 SUBROUTINE INITAL C... C... DYNAMIC ANALYSIS BY DIFFERENCE AND DIFFERENTIAL EQUATIONS C... C... THE NUMERICAL SOLUTIONS TO THE FIRST-ORDER DIFFERENCE EQUATION C... C... Y(N+1) = C*Y(N), Y(0) = Y0 (1)(2) C... C... AND THE FIRST-ORDER DIFFERENTIAL EQUATION C... C... DY/DT = R*Y, Y(0) = Y0 (3)(4) C... C... ARE COMPUTED AND COMPARED IN THIS PROGRAM. C... C... EQUATIONS (1) AND (3) HAVE RELATED MATHEMATICAL PROPERTIES AND C... SOLUTIONS. EQUATION (1) IS JUST AN EXPRESSION OF A COMPOUND C... INTEREST CALCULATION IF C... C... C = 1 + R*DT (5) C... C... WHERE C... C... R ANNUAL INTEREST RATE (FRACTION/YEAR) C... C... DT INTEREST PERIOD (YEARS) C... C... Y0 INITIAL PRINCIPAL C... C... TO ILLUSTRATE THIS POINT, CONSIDER THE SOLUTION OF EQUATION (1) C... FOR R = 0.07, DT = 1, Y0 = 100, I.E. ANNUAL COMPOUNDING AT 7 PER C... CENT WITH AN INITIAL PRINCIPAL OF 100. AT THE END OF THE FIRST C... YEAR, EQUATION (1) GIVES C... C... Y(1) = (1 + R*DT)*Y(0) = (1 + 0.07*1)*100 = 100 + 0.07*1*100 C... C... = 100 + 7 = 107 C... C... NOTE THAT THE FIRST YEAR*S INTEREST OF 7 IS ADDED TO THE PRINCIPAL C... FOR A SUM OF 107 AT THE END OF THE FIRST YEAR. EQUATION (1) CAN C... THEN BE APPLIED FOR THE SECOND YEAR AS C... C... Y(2) = (1 + R*DT)*Y(1) = (1 + 0.07*1)*107 = (1.07**2)*100 C... C... BY INDUCTION, THE GENERAL SOLUTION TO EQUATION (1) FOR ANY YEAR C... (ANY N) IS C... C... Y(N+1) = (C**(N+1))*Y(0) = ((1 + R*DT)**(N+1))*Y(0) (6) C... C... (THIS SOLUTION CAN BE VERIFIED BY DIRECT SUBSTITUTION IN EQUATIONS C... (1) AND (2)). C... C... HOWEVER, N CAN BE GENERALIZED SO THAT IT IS NOT NECESSARILY A C... COUNTER FOR YEARS, BUT RATHER FOR INTEREST PERIODS. THUS IF C... INTEREST IS PAID QUARTERLY, FOR THE FIRST YEAR N WOULD HAVE THE C... SUCCESSIVE VALUES N = 1, 2, 3, 4 AND DT WOULD NOW BE DT = 0.25. C... EQUATION (6) CAN STILL BE USED IN THIS MORE GENERAL CASE TO CALCU- C... LATE THE PRINCIPAL PLUS ACCUMULATED INTEREST AT THE END OF ANY C... PERIOD N. C... C... THE USUAL NEXT STEP IN THIS DEVELOPMENT IS TO CONSIDER EQUATION C... (6) IN THE LIMIT AS N APPROACHES INFINITY AND DT APPROACHES ZERO, C... I.E., CONTINUOUS COMPOUNDING. THE FINAL RESULT IS THAT Y INCREASES C... EXPONENTIALLY WITH TIME. HOWEVER, THIS POINT WILL BE DEMONSTRATED C... IN A SLIGHTLY DIFFERENT FASHION. C... C... EQUATION (1) WITH C = 1 + R*DT CAN BE REARRANGED AS C... C... (Y(N+1) - Y(N))/DT = R*Y(N) (7) C... C... IF DISCRETE INDEPENDENT VARIABLE N IS CONSIDERED AS AN INDEX FOR C... A CONTINUOUS VARIABLE T C... C... T = N*DT, (AND T + DT = (N + 1)*DT) (8) C... C... EQUATION (7) CAN BE WRITTEN AS C... C... (Y(T+DT) - Y(T))/DT = R*Y(T) (9) C... C... IF N APPROACHES INFINITY AND DT APPROACHES ZERO, WE RECOGNIZE THE C... LEFT HAND SIDE OF EQUATION (9) AS THE FIRST-ORDER DERIVATIVE OF C... Y WITH RESPECT TO T, DY/DT. THUS THE LIMITING CASE OF CONTINUOUS C... COMPOUNDING CAN BE DESCRIBED MATHEMATICALLY IN TERMS OF THE C... DIFFERENTIAL EQUATION C... C... DY/DT = R*Y (3) C... C... AGAIN, THE ESSENTIAL DIFFERENCE BETWEEN EQUATIONS (1) AND (3) IS C... THE NATURE OF THE RESPECTIVE INDEPENDENT VARIABLES, N AND T, I.E., C... N IS A DISCRETE VARIABLE WHICH CAN ASSUME ONLY INTEGER VALUES C... WHILE T IS A CONTINUOUS VARIABLE. THUS MATHEMATICAL MODELS BASED C... ON EQUATIONS LIKE (1) AND (3) ARE FREQUENTLY TERMED **DISCRETE** C... AND **CONTINUOUS** RESPECTIVELY. C... C... EQUATION (3), WITH THE INITIAL CONDITION (4), HAS THE ANALYTICAL C... (EXACT) SOLUTION C... C... Y(T) = Y0*EXP(R*T) (10) C... C... WHICH CAN EASILY BE VERIFIED BY SUBSTITUTION IN EQUATIONS (3) AND C... (4), AND WHICH CLEARLY INDICATES THAT CONTINUOUS COMPOUNDING C... CORRESPONDS TO EXPONENTIAL GROWTH. C... C... WE SHALL CONSIDER A FEW ADDITIONAL TERMS AND MATHEMATICAL CONCEPTS C... COMMONLY ENCOUNTERED IN THE USE OF DIFFERENCE AND DIFFERENTIAL C... EQUATIONS. EQUATION (1) IS A **FIRST-ORDER LINEAR DIFFERENCE C... EQUATION WITH CONSTANT COEFFICIENTS**. IT IS FIRST-ORDER BECAUSE C... THE MAXIMUM DIFFERENCE IN VALUES OF THE INDEPENDENT VARIABLE N IS C... ONE. THUS THE FOLLOWING EQUATION IS SECOND-ORDER C... C... A2*Y(N+2) + A1*Y(N+1) + A0*Y(N) = 0 (11) C... C... AS IS C... C... B1*Y(N+1) + B0*Y(N) + Y(N-1) = 0 (12) C... C... EQUATION (1) IS LINEAR BECAUSE THE DEPENDENT VARIABLE Y APPEARS C... ONLY TO THE FIRST POWER, I.E., EQUATION (1) IS FIRST-DEGREE. THUS C... THE FOLLOWING EQUATION IS NONLINEAR C... C... Y(N+1) = C*Y(N)**2 (13) C... C... EQUATION (13) IS FIRST-ORDER, BUT SECOND-DEGREE. C... C... EQUATION (1) HAS CONSTANT COEFFICIENTS SINCE THE TERMS MULTIPLYING C... THE DEPENDENT VARIABLE ARE CONSTANT, I.E., ONE FOR Y(N+1) AND C... FOR Y(N). THE FOLLOWING EQUATION HAS VARIABLE COEFFICIENTS (WHICH C... ARE FUNCTIONS OF THE INDEPENDENT VARIABLE N, BUT NOT THE DEPENDENT C... VARIABLE Y) C... C... (N**2)*Y(N+1) = SIN(N)*Y(N) (14) C... C... AS THE ORDER, COMPLEXITY OF THE COEFFICIENTS AND NONLINEARITY OF C... A DIFFERENCE EQUATION ARE INCREASED, IT GENERALLY BECOMES MORE C... DIFFICULT TO SOLVE ANALYTICALLY. IN FACT, BECAUSE EQUATION (1) IS C... THE SIMPLEST POSSIBLE DIFFERENCE EQUATION, WE HAD NO DIFFICULTY IN C... WRITING DOWN ITS ANALYTICAL SOLUTION, EQUATION (6). IN GENERAL, C... HOWEVER, THE DIFFERENCE EQUATIONS IN MATHEMATICAL MODELS ARE TOO C... COMPLICATED, AND TOO MANY IN NUMBER, TO SOLVE ANALYTICALLY. THEY C... CAN, HOWEVER, IN MOST CASES BE SOLVED NUMERICALLY WITH THE AID OF C... A COMPUTER, AS IN THE CASE OF EQUATION (1) WHICH WAS SOLVED NUMER- C... ICALLY BY THE COMPOUND INTEREST CALCULATION (THIS REQUIRED NOTHING C... MORE THAN ARITHMETIC MULTIPLICATION AND ADDITION). C... C... FINALLY EACH DIFFERENCE EQUATION IN A MATHEMATICAL MODEL MUST HAVE C... ONE OR MORE INITIAL CONDITIONS. THE NUMBER OF INITIAL CONDITIONS C... FOR EACH EQUATION EQUALS THE ORDER OF THE EQUATION. THUS FOR C... EQUATION (1), WHICH IS FIRST-ORDER, ONLY ONE INITIAL CONDITION, C... EQUATION (2), IS REQUIRED. EQUATIONS (11) AND (12) WOULD EACH C... REQUIRE TWO INITIAL CONDITIONS. C... C... ALL OF THE PRECEDING STATEMENTS CONCERNING EQUATIONS (1) AND (2) C... CAN ALSO BE APPLIED TO EQUATIONS (3) AND (4). THUS, EQUATION (3) C... IS A **FIRST-ORDER LINEAR DIFFERENTIAL EQUATION WITH CONSTANT C... COEFFICIENTS**. SINCE IT IS FIRST-ORDER, IT REQUIRES ONE INITIAL C... CONDITION, EQUATION (4). C... C... BECAUSE EQUATION (3) IS THE SIMPLEST POSSIBLE DIFFERENTIAL EQUA- C... TION, WE CAN IMMEDIATELY WRITE DOWN ITS ANALYTICAL SOLUTION, C... EQUATION (10). HOWEVER, AS THE ORDER, COMPLEXITY, NONLINEARITY C... AND NUMBER OF DIFFERENTIAL EQUATIONS INCREASE, WE GENERALLY ARE C... UNABLE TO OBTAIN ANALYTICAL SOLUTIONS. THEIR PRACTICAL USE IN A C... CONTINUOUS MATHEMATICAL MODEL THEREFORE USUALLY REQUIRES SOME C... NUMERICAL SOLUTION PROCEDURE. IN FACT, A NUMERICAL INTEGRATION C... ALGORITHM CAN BE USED WHICH IS A DIRECT ANALOG OF THE COMPOUND C... INTEREST CALCULATION APPLIED TO EQUATION (1). C... C... TO EXPLORE THIS POINT FURTHER, IF EQUATION (9) IS WRITTEN AS C... C... Y(T+DT) = Y(T) + R*Y(T)*DT (15) C... C... EQUATION (15) CAN BE USED TO TAKE A STEP ALONG THE SOLUTION, I.E., C... COMPUTE Y(T+DT) USING THE BASE POINT VALUE Y(T). OF COURSE, IF DT C... IS FINITE, AS IT MUST BE IN ANY PRACTICAL COMPUTER CALCULATION, WE C... WOULD NOT EXPECT TO COMPUTE SUCCESSIVE VALUES OF THE SOLUTION, C... Y(DT), Y(2*DT),..., WHICH ARE IN AGREEMENT WITH THE EXACT SOLU- C... TION, EQUATION (10). WE WOULD EXPECT, HOWEVER, THAT AS THE C... CALCULATION IS REPEATED FOR SUCCESSIVELY SMALLER VALUES OF DT, THE C... CORRESPONDING NUMERICAL SOLUTIONS WOULD BE IN BETTER AGREEMENT C... WITH EQUATION (10), I.E., THE NUMERICAL SOLUTION WILL CONVERGE TO C... THE EXACT SOLUTION. THUS WE MIGHT REPEAT THE NUMERICAL CALCULA- C... TION FOR TWO DIFFERENT VALUES OF DT, AND COMPARE THE NUMERICAL C... SOLUTIONS. IF THEY AGREE WITHIN SOME PRESELECTED TOLERANCE, WE C... COULD REASONABLY CONCLUDE THAT THE NUMERICAL SOLUTION IS ACCURATE C... TO THIS TOLERANCE. IF THE TWO SOLUTIONS DO NOT AGREE, WE COULD C... REPEAT THE CALCULATION A THIRD TIME WITH A STILL SMALLER DT, AND C... COMPARE THE NUMERICAL SOLUTIONS FOR THE TWO SMALLER DTS. IF THE C... AGREEMENT IS NOT GOOD ENOUGH, WE COULD AGAIN SELECT A SMALLER C... DT AND COMPARE THE NUMERICAL SOLUTIONS FOR THE TWO SMALLEST DTS. C... HOPEFULLY WE WOULD EVENTUALLY FIND A DT SMALL ENOUGH TO SATISFY C... THE ERROR TOLERANCE. ALTHOUGH THIS TRIAL-AND-ERROR PROCEDURE IS C... TEDIOUS, IT IS WELL SUITED FOR A DIGITAL COMPUTER, AND WE MIGHT C... PROGRAM THE ENTIRE CALCULATION SO THAT THE COMPUTER SELECTS DT TO C... PRODUCE A NUMERICAL SOLUTION SATISFYING THE USER-SPECIFIED ERROR C... CRITERION. THIS CAN INDEED BE DONE AND IN FACT ANY REASONABLY C... PROGRAMMED INTEGRATION ALGORITHM SHOULD HAVE THIS AUTOMATIC ERROR C... MONITORING AND STEP-SIZE ADJUSTMENT FEATURE SO THAT THE USER IS C... RELIEVED OF HAVING TO SELECT A STEP SIZE. THE PROCEDURE JUST C... DESCRIBED IS WORKABLE, BUT CAN BE IMPROVED FOR AT LEAST TWO C... REASONS C... C... (1) THE INDIRECT ERROR ESTIMATION, BASED ON A COMPARISON OF C... SOLUTIONS FOR SUCCESSIVELY SMALLER INTEGRATION STEP SIZES, C... MAY LEAD TO A FALSE CONCLUSION IF THE CONVERGENCE TOWARD C... THE EXACT SOLUTION IS SLOW. FOR EXAMPLE, SUCCESSIVE NUM- C... ERICAL SOLUTIONS MAY AGREE TO WITHIN THE ERROR TOLERANCE, C... BUT HAVE AN ERROR FAR GREATER THAN THE ERROR TOLERANCE C... SIMPLY BECAUSE THEY ARE CONVERGING TO THE EXACT SOLUTION C... VERY SLOWLY. C... C... (2) THE INDIRECT OR IMPLICIT ERROR MONITORING IS RATHER IN- C... EFFICIENT IN THAT AT LEAST TWO COMPLETE SOLUTIONS AT A C... GIVEN T MUST BE COMPUTED FOR COMPARISON. THEN, DEPENDING C... ON THE RESULT OF THE COMPARISON, A THIRD SOLUTION FOR A C... SMALLER DT MAY HAVE TO BE CALCULATED, ETC. THUS THE C... METHOD IS RELATIVELY INEFFICIENT IN TERMS OF THE CALCU- C... LATIONAL EFFORT REQUIRED TO ESTIMATE THE ERROR IN THE C... NUMERICAL SOLUTION. C... C... INTEGRATION ALGORITHMS HAVE BEEN DEVELOPED WHICH LARGELY CIRCUM- C... VENT THESE PROBLEMS BY USING RELIABLE, EXPLICIT ERROR ESTIMATES C... TO ADJUST DT IN ACCORDANCE WITH A USER-PRESCRIBED ERROR CRITERION. C... THE DETAILS OF THESE ALGORITHMS ARE BEYOND THE SCOPE OF THE PRE- C... SENT DISCUSSION, BUT THEY ARE READILY AVAILABLE FOR USE IN THE C... SOLUTION OF CONTINUOUS MATHEMATICAL MODEL EQUATIONS. IN FACT, THE C... USE OF SUCH AN ALGORITHM IS DEMONSTRATED WITH THIS PROGRAM. C... C... EQUATION (15) IS AN EXAMPLE OF THE SIMPLEST OF INTEGRATION ALGO- C... RITHMS, TERMED EULER*S METHOD. IN GENERAL, FOR A DIFFERENTIAL C... EQUATION OF THE FORM C... C... DY/DT = F(Y,T), Y(T0) = Y0 (16)(17) C... C... WHERE F(Y,T) IS A GIVEN FUNCTION, EULER*S METHOD APPLIED TO EQUA- C... TION (16) IS C... C... Y(T+DT) = Y(T) + F(Y(T),T)*DT (18) C... C... EQUATION (18) CAN BE APPLIED TO A BROAD SPECTRUM OF PROBLEMS. THE C... ONLY COMPUTATIONAL REQUIREMENTS ARE (1) THE EVALUATION OF C... F(Y(T),T) AT EACH POINT ALONG THE SOLUTION, (2) THE MULTIPLICATION C... OF F(Y(T),T) BY DT AND (3) THE ADDITION OF F(Y(T),T)*DT TO Y(T). C... THUS, THE DERIVATIVE FUNCTION, F(Y,T), CAN BE OF ANY FORM AND THE C... ONLY REQUIREMENT IS THAT IT BE PROGRAMMABLE FOR COMPUTER EVALUA- C... TION. EQUATION (18) CAN ALSO BE APPLIED DIRECTLY TO SYSTEMS C... OF DIFFERENTIAL EQUATIONS OF THE FORM C... C... DY1/DT = F1(Y1,Y2,...,YM,T) C... C... DY2/DT = F2(Y1,Y2,...,YM,T) C... . . (19) C... . . C... DYM/DT = FM(Y1,Y2,...,YM,T) C... C... WITH THE INITIAL CONDITION VECTOR C... C... Y1(T0) = Y10, Y2(T0) = Y20,..., YM(T0) = YM0 (20) C... C... THE SYSTEM OF M SIMULTANEOUS DIFFERENTIAL EQUATIONS (19) IS COM- C... PLETELY GENERAL, SO THAT NUMERICAL INTEGRATION BASED ON EULER*S C... METHOD, EQUATION (18), OR MORE SOPHISTICATED INTEGRATION ALGO- C... RITHMS IS A VERY POWERFUL TECHNIQUE FOR THE SOLUTION OF CONTINUOUS C... MATHEMATICAL MODELS. MODELS WITH M GT 100 ARE COMMONPLACE. C... C... IN THE REMAINDER OF THIS PROGRAM, THE NUMERICAL SOLUTIONS OF C... EQUATIONS (1) TO (4) ARE COMPUTED FOR R = 0.07, DT = 1, 0.5, 0.25 C... (THREE CASES FOR DIFFERENCE EQUATION (1)), Y0 = 100. IN EACH CASE C... THE SOLUTIONS TO EQUATIONS (1) AND (3) ARE COMPARED WITH EACH C... OTHER AND WITH EQUATION (10). C... C... FINALLY, THREE INTERESTING PROPERTIES OF THE NUMERICAL SOLUTIONS C... SHOULD BE NOTED C... C... (1) IF THE TIME FOR A QUANTITY TO DOUBLE WHICH IS GROWING C... EXPONENTIALLY IS DENOTED AS TD, IT CAN BE EASILY BE C... CALCULATED FROM EQUATION (10) AS C... C... Y(T)/Y0 = 2 = EXP(R*TD) C... C... OR C... C... TD = LN(2)/R = 0.692/R = 70/(100*R) (21) C... C... EQUATION (21) INDICATES THAT THE DOUBLING TIME, TD, IS C... APPROXIMATELY EQUAL TO 70 DIVIDED BY THE RATE OF INCREASE C... R, EXPRESSED AS A PERCENTAGE. THUS IF R = 0.07 AS IN THE C... NUMERICAL EXAMPLE, TD = 70/(100*0.07) = 10 YEARS. THIS C... RESULT CAN BE CHECKED WITH THE NUMERICAL SOLUTIONS. C... C... (2) THE DIFFERENCE BETWEEN THE NUMERICAL SOLUTION TO EQUATION C... (1) FOR DT = 1.0, 0.5 AND 0.25 AND THE SOLUTION OF EQUA- C... TION (3) (I.E., EQUATION (10)) DECREASES LINEARLY WITH C... DT. THUS THE COMPOUND INTEREST CALCULATION, OR MORE C... GENERALLY, EULER*S METHOD, EQUATION (18), IS A FIRST-ORDER C... ALGORITHM. IN OTHER WORDS, THE DISCRETIZATION ERROR DUE C... TO FINITE DT DECREASES LINEARLY WITH DT. A SECOND-ORDER C... ALGORITHM WOULD HAVE A DISCRETIZATION ERROR WHICH DE- C... CREASES AS DT**2, ETC. ALGORITHMS WHICH ARE FIFTH-ORDER C... CORRECT ARE READILY AVAILABLE FOR THE NUMERICAL INTEGRA- C... TION OF ORDINARY DIFFERENTIAL EQUATIONS. C... C... (3) ALL OF THE PRECEDING DISCUSSION IS APPLICABLE TO THE CASE C... R LESS THAN (LT) ZERO. Y(T) GIVEN BY EQUATION (10) NOW C... DECREASES EXPONENTIALLY WITH T, AND APPROACHES ZERO AS T C... APPROACHES INFINITY. THE SOLUTION TO EQUATION (1), C... HOWEVER, (EQUATION (6) OR MORE GENERALLY, EULER*S METHOD, C... EQUATION (18)) GIVES SOME UNEXPE3TED RESULTS, DEPENDING ON C... THE VALUE OF DT. FOR EXAMPLE, IF DT = -1/R, FROM EQUATION C... (6), Y(N) = 0 FOR ALL N. IF - 2/R LT DT LT -1/R, Y(N) C... HAS DECREASING VALUES WITH INCREASING N, BUT THESE VALUES C... WILL ALTERNATE IN SIGN (THE READER MIGHT CONFIRM THIS BY C... AN EXAMPLE COMPOUND INTEREST CALCULATION APPLIED TO EQUA- C... TION (1) FOR DT = -1.5/R). FINALLY, IF DT LT -2/R, Y(N) C... HAS VALUES WHICH INCREASE IN ABSOLUTE MAGNITUDE WITH N, C... AND ALTERNATE IN SIGN (THE READER MIGHT CONFIRM THIS BY C... AN EXAMPLE COMPOUND INTEREST CALCULATION APPLIED TO EQUA- C... TION (1) FOR DT = -3/R). IN THIS CASE, THE SOLUTION TO C... EQUATION (1) IS UNSTABLE (I.E., INCREASES WITHOUT BOUND), C... EVEN THOUGH FOR R LT 0, WE WOULD EXPECT THE SOLUTION TO C... BE STABLE IN MUCH THE SAME WAY AS Y(T) GIVEN BY EQUATION C... (10). THIS DISCUSSION ILLUSTRATES THAT NUMERICAL CALCU- C... LATIONS, AND THEIR UNDERLYING ALGORITHMS, SHOULD BE BOTH C... STABLE AND ACCURATE, SPECIFICALLY, EULER*S METHOD, EQUA- C... TION (18), HAS A STABILITY LIMIT C... C... ABS(R*DT) LT 2 (22) C... C... (THIS FOLLOWS FROM THE PRECEDING DISCUSSION) AND IS FIRST- C... ORDER CORRECT. C... C... COMMON AREAS WHICH LINK THE SUBROUTINES OF THIS PROGRAM COMMON/T/T,NFIN,NORUN/Y/Y/F/DYDT/PARM/R,DT,YN C... C... SET THE PARAMETERS R, DT R=0.07 IF(NORUN.EQ.1)DT=1. IF(NORUN.EQ.2)DT=0.5 IF(NORUN.EQ.3)DT=0.25 C... C... INITIALIZE THE NUMERICAL SOLUTIONS FOR EQUATIONS (1) AND (3), C... I.E., IMPLEMENT INITIAL CONDITIONS (2) AND (4) YN=100. Y =100. RETURN END SUBROUTINE DERV C... C... SUBROUTINE DERV IMPLEMENTS THE DERIVATIVE CALCULATION OF EQUATION C... (3) COMMON/T/T,NFIN,NORUN/Y/Y/F/DYDT/PARM/R,DT,YN C... C... EQUATION (3) DYDT=R*Y RETURN END SUBROUTINE PRINT(NI,NO) C... C... SUBROUTINE PRINT COMPUTES THE EXACT SOLUTION, EQUATION (10), AND C... PRINTS IT ALONG WITH THE TWO NUMERICAL SOLUTIONS FOR COMPARISON C... COMMON/T/T,NFIN,NORUN/Y/Y/F/DYDT/PARM/R,DT,YN C... C... PRINT A HEADING DATA N/0/ IF(N.EQ.0)WRITE(NO,1) 1 FORMAT(20X,5X,10HDIFFERENCE,3X,12HDIFFERENTIAL,7X,8HEQUATION,5X, 110HDIFFERENCE,3X,12HDIFFERENTIAL,/,4X,6HSTEP N,9X,1HT,7X,8HEQUATIO 2N,7X,8HEQUATION,11X,4H(10),5X,10H- EQN (10),5X,10H- EQN (10)) C... C... COMPUTE THE NUMERICAL SOLUTION TO EQUATION (1) YNP1=(1.+R*DT)*YN IF(N.EQ.0)YNP1=YN C... C... COMPUTE THE EXACT SOLUTION FOR EQUATION (3) (EQUATION (10)) YE=100.*EXP(R*T) C... C... COMPUTE THE DIFFERENCES BETWEEN SOLUTIONS YND=YNP1-YE YD=Y-YE C... C... PRINT THE RESULTS WRITE(NO,2)N,T,YNP1,Y,YE,YND,YD 2 FORMAT(I10,F10.3,5F15.4) C... C... UPDATE THE NUMERICAL SOLUTION OF EQUATION (1) YN=YNP1 N=N+1 C... C... RESET THE COUNTER N IF((NORUN.EQ.1).AND.(N.EQ.21))N=0 IF((NORUN.EQ.2).AND.(N.EQ.41))N=0 IF((NORUN.EQ.3).AND.(N.EQ.81))N=0 RETURN END COMPARISON OF DIFFERENCE AND DIFFERENTIAL SOLUTIONS, DT = 1.0 0. 20.0 1.0 1 2000 1 1 ABS 0.001 COMPARISON OF DIFFERENCE AND DIFFERENTIAL SOLUTIONS, DT = 0.5 0. 20.0 0.5 1 2000 1 1 ABS 0.001 COMPARISON OF DIFFERENCE AND DIFFERENTIAL SOLUTIONS, DT = 0.25 0. 20.0 0.25 1 2000 1 1 ABS 0.001 END OF RUNS *APP PHYS14 SUBROUTINE INITAL C... C... INTEGRATION OF THE GIBBS-DUHEM EQUATION C... C... THE ISOTHERMAL GIBBS-DUHEM EQUATION CAN BE USED TO COMPUTE THE C... EQUILIBRIUM VAPOR COMPOSITION OF A BINARY SYSTEM FROM EXPERIMENTAL C... TOTAL PRESSURE VS LIQUID COMPOSITION DATA. FOR THIS STUDY, THE C... FOLLOWING VERSIONS OF THE GIBBS-DUHEM EQUATION ARE INTEGRATED C... C... DY1/DX1 = (1/P)*(DP/DX1)*(1 - Y1)/(1 - X1/Y1) (1) C... C... DY2/DX2 = (1/P)*(DP/DX2)*(1 - Y2)/(1 - X2/Y2) (2) C... C... EQUATIONS (1) AND (2) ARE BASED ON THE ASSUMPTIONS C... C... (1) IDEAL GAS VAPOR PHASE C... C... (2) THE TERM (EXCESS VOL*DP)/RT IN THE EXACT, ISOTHERMAL GIBBS C... DUHEM EQUATION MAY BE NEGLECTED. THIS IS AN EXCELLENT C... ASSUMPTION HERE AND IN MOST CASES AWAY FROM THE CRITICAL C... REGION AND IN THE ABSENCE OF VERY HIGH PRESSURES. C... C... EQUATIONS (1) AND (2) ARE INTEGRATED FOR THE ETHYL ALCOHOL (1) C... CHLOROFORM (2) SYSTEM FOR WHICH C... C... X1 MOLE FRACTION OF ETHYL ALCOHOL IN THE LIQUID PHASE C... C... X2 MOLE FRACTION OF CHLOROFORM IN THE LIQUID PHASE C... C... Y1 MOLE FRACTION OF ETHYL ALCOHOL IN THE EQUILIBRIUM VAPOR C... PHASE C... C... Y2 MOLE FRACTION OF CHLOROFORM IN THE EQUILIBRIUM VAPOR C... PHASE C... C... P TOTAL PRESSURE MEASURED EXPERIMENTALLY AS A FUNCTION OF C... X1 AND X2 C... C... THE FOLLOWING DATA ARE GIVEN IN PRAUSNITZ, J. M., MOLECULAR C... THERMODYNAMICS OF FLUID-PHASE EQUILBRIA, P. 260, PROB. 7, C... PRENTICE-HALL, INC. ENGLEWOOD CLIFFS, N. J. 3 C... C... TOTAL PRESSURE DATA FOR THE SYSTEM C... ETHYL ALCOHOL(1)/CHLOROFORM(2) C... AT 45 C C... C... X2 P(MM HG) X1 P(MM HG) C... 0 172. 0 433.5 C... 0.05 200 0.01 438.5 C... 0.10 233.8 0.02 442.0 C... 0.15 266 0.03 444.9 C... 0.20 298 0.04 447.4 C... 0.25 326 0.05 449.3 C... 0.30 351 0.06 450.9 C... 0.35 372 0.07 452.3 C... 0.45 403.5 0.09 454.5 C... 0.50 416 0.10 455.2 C... 0.55 426.5 C... C... THESE DATA WERE FITTED TO REGRESSION POLYNOMIALS WITH THE FOLLOW- C... ING RESULTS C... C... C... P = 4.335524475527E+02 + C... C... 5.434965034889E+02*X1 - 6.918997668803E+03*X1**2 + (3) C... C... 5.687645687223E+04*X1**3 - 2.309627039370E+05*X1**4 C... C... C... P = 1.723532051343E+02 + C... C... 5.194413688298E+02*X2 + 1.266380147487E+03*X2**2 - (4) C... C... 4.399313649081E+03*X2**3 + 3.469696969444E+03*X2**4 C... C... C... THUS EQUATIONS (3) AND (4) CAN BE SUBSTITUTED IN EQUATIONS (1) AND C... (2) TO ELIMINATE P. THE RESULTING EQUATIONS CAN THEN BE INTEGRAT- C... ED OVER THE INTERVALS 0 LE X1 LE 0.1, 0 LE X2 LE 0.55 TO OBTAIN Y1 C... AND Y2. THE INITIAL CONDITIONS ARE C... C... Y1 = 0, X1 = 0 (5) C... C... Y2 = 0, X2 = 0 (6) C... C... A REFERENCE CONTAINING EXPERIMENTAL VAPOR PHASE COMPOSITIONS FOR C... COMPARISON WITH THE CALCULATED VAPOR PHASE COMPOSITIONS IS GIVEN C... BY PRAUSNITZ. C... COMMON/T/X,NFIN,NRUN/Y/Y/F/DYDX C... C... INITIAL CONDITION (A SMALL NONZERO VALUE OF THE INITIAL CONDITION C... IS USED TO AVOID A DIVISION BY ZERO IN THE CALCULATION OF DY/DX C... IN SUBROUTINE DERV) Y=1.0E-05 RETURN END SUBROUTINE DERV COMMON/T/X,NFIN,NRUN/Y/Y/F/DYDX C... C... DIFFERENTIAL EQUATION FOR DY/DX DYDX=1./P(NRUN)*(DPDX(NRUN))*(1.-Y)/(1.-X/Y) RETURN END FUNCTION P(NRUN) C... C... FUNCTION P COMPUTES THE TOTAL PRESSURE FROM THE REGRESSION C... POLYNOMIALS, EQUATIONS (3) AND (4) C... COMMON/T/X GO TO(1,2),NRUN 1 P = 4.335524475527E+02 1 +5.434965034889E+02*X -6.918997668803E+03*X**2 2 +5.687645687223E+04*X**3 -2.039627039370E+05*X**4 RETURN 2 P = 1.723532051343E+02 1 +5.194413688298E+02*X +1.266380147487E+03*X**2 2 -4.399313649081E+03*X**3 +3.469696969444E+03*X**4 RETURN END FUNCTION DPDX(NRUN) C... C... FUNCTION DPDX COMPUTES THE DERIVATIVE OF THE PRESSURE, DP/DX, FROM C... THE REGRESSION POLYNOMIALS, EQUATIONS (3) AND (4) C... COMMON/T/X GO TO(1,2),NRUN 1 DPDX= 1 +5.434965034889E+02 -6.918997668803E+03*X*2. 2 +5.687645687223E+04*X**2*3.-2.039627039370E+05*X**3*4. RETURN 2 DPDX= 1 +5.194413688298E+02 +1.266380147487E+03*X*2. 2 -4.399313649081E+03*X**2*3.+3.469696969444E+03*X**3*4. RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/X,NFIN,NRUN/Y/Y/F/DYDX C... C... DIMENSION THE ARRAYS FOR PLOTTING DIMENSION XPLOT(45),YPLOT(2,45) DATA IP/0/ C... C... PRINT A HEADING IF(IP.EQ.0)WRITE(NO,2) 2 FORMAT(13X,2HX1,13X,2HY1,13X,2HX2,13X,2HY2,14X,1HP) C... C... PRINT THE NUMERICAL SOLUTION 1 PRESS=P(NRUN) OMX=1.-X OMY=1.-Y GO TO(3,4),NRUN 3 WRITE(NO,5)X,Y,OMX,OMY,PRESS 5 FORMAT(4F15.4,F15.1) GO TO 6 4 WRITE(NO,5)OMX,OMY,X,Y,PRESS C... C... STORE THE NUMERICAL SOLUTION FOR PLOTTING 6 IP=IP+1 GO TO(7,8),NRUN 7 XPLOT(IP)=X YPLOT(1,IP)=Y YPLOT(2,IP)=X GO TO 9 8 XPLOT(IP)=OMX YPLOT(1,IP)=OMY YPLOT(2,IP)=OMX C... C... PLOT THE EQUILIBRIUM DIAGRAM Y1 VS X1 9 IF((NRUN.NE.2).OR.(X.LT.0.549))RETURN CALL TPLOTS(2,IP,XPLOT,YPLOT) WRITE(NO,11) 11 FORMAT(//,16X,87HY1 (MOLE FRACTION ETHYL ALCOHOL IN VAPOR) VS X1 ( 1MOLE FRACTION ETHYL ALCOHOL IN LIQUID)) RETURN END ETHYL ALCOHOL/CHLOROFORM SYSTEM, 0 LE X1 LE 0.10 0. 0.1 0.025 1 1000 3 1 REL 0.0001 ETHYL ALCOHOL/CHLOROFORM SYSTEM, 0 LE X2 LE 0.55 0. 0.55 0.025 1 1000 3 1 REL 0.0001 END OF RUNS *APP PHYS15 SUBROUTINE INITAL C... C... INTEGRATION OF BESSEL*S EQUATION OF ORDER ZERO AND ONE C... C... BESSEL*S EQUATION OF ORDER ZERO WITH INITIAL CONDITIONS IS C... C... (X**2)*D2Y/DT2 + X*DY/DX + (X**2)*Y = 0 (1) C... C... Y(0) = 1 (2) C... C... DY(0)/DX = 0 (3) C... C... THE SOLUTION IS THE WELL-KNOWN BESSEL FUNCTION OF THE FIRST KIND, C... GENERALLY DENOTED AS J0(X) C... C... BESSEL*S EQUATION OF ORDER ONE WITH INITIAL CONDITIONS IS C... C... (X**2)*D2Y/DX2 + X*DY/DX + (X**2 - 1)*Y = 0 (4) C... C... Y(0) = 0 (5) C... C... DY(0)/DX = 1/2 (6) C... C... THE SOLUTION IS GENERALLY DENOTED AS J1(X)) C... C... EQUATIONS (1) AND (4) ARE NUMERICALLY INTEGRATED VIA A SYSTEM OF C... TWO EQUIVALENT FIRST-ORDER EQUATIONS. THE INITIAL VALUE OF THE C... INDEPENDENT VARIABLE (I.E. X = 0) IS REPLACED WITH THE APPROXI- C... MATE VALUE OF X = 1.0E-6 TO AVOID DIFFICULTIES IN THE INITIAL COM- C... PUTATION OF THE DERIVATIVES IN SUBROUTINE DERV. C... C... TABULATED VALUES OF J0(X) AND J1(X) ARE LISTED BELOW FOR COMPARI- C... SON WITH THE NUMERICAL SOLUTIONS C... C... X J0(X) J1(X) C... 1.0 0.76519 0.44005 C... 2.0 0.22389 0.57672 C... 3.0 -0.26005 0.33905 C... 4.0 -0.39714 -0.06604 C... 5.0 -0.17759 -0.32757 C... 6.0 0.15046 -0.27668 C... 7.0 0.30007 -0.00468 C... 8.0 0.17165 0.23463 C... 9.0 -0.09033 0.24531 C... 10.0 -0.24593 0.04347 C... C... REF - ABRAMOWITZ, M., AND I. A. STEGUN, HANDBOOK OF MATHEMATICAL C... FUNCTIONS, NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C... SERIES, NO. 55, SEVENTH PRINTING, PP. 390-392, MAY, 1968. A COM- C... PARISON OF J1(X) WHICH IS TABULATED ABOVE AND DJ0(X)/DX WHICH IS C... COMPUTED NUMERICALLY INDICATES DJ0(X)/DX = -J1(X) C... C... THE PROBLEMS DEFINED BY EQUATIONS (1) TO (6) PROVIDE AN INTER- C... ESTING AND SOMEWHAT UNEXPECTED RESULT. FOR EQUATION (1), NMAX = C... 100 (I.E. THE RATIO OF THE PRINT INTERVAL TO THE MINIMUM ALLOWABLE C... INTEGRATION INTERVAL SET ON THE THIRD DATA LINE) GIVES A NUMERI- C... CAL SOLUTION WITH NO REPORTED INTEGRATION ERRORS AND WHICH AGREES C... WITH THE TABULATED VALUES OF J0(X) TO FOUR FIGURES OR BETTER. C... HOWEVER, FOR EQUATION (4), EVEN WITH NMAX = 99999, THE MAXIMUM C... VALUE WHICH DSS/2 CAN READ UNDER AN I5 FORMAT, THE NUMERICAL C... SOLUTION IS OBVIOUSLY IN ERROR WHEN COMPARED WITH THE TABULATED C... VALUES OF J1(X). FURTHERMORE, AN INTEGRATION ERROR FOR THE C... SECOND DEPENDENT VARIABLE Y2 (= DY/DX) IS REPORTED. C... C... THESE EXAMPLES CLEARLY INDICATE THAT THE PERFORMANCE AND C... EFFICIENCY OF A NUMERICAL INTEGRATION ALGORITHM IS STRONGLY AND C... GENERALLY UNPREDICTABLY RELATED TO THE CHARACTERISTICS OF THE C... PARTICULAR PROBLEM. IN THE CASE OF EQUATION (4), THE ADDITIONAL C... TERM -((N/X)**2)*Y1 WITH N = 1 (V. SUBROUTINE DERV) LEADS TO A C... SIGNIFICANT COMPUTATIONAL ERROR AT X = 0 WHICH CAUSES THE SUB- C... SEQUENT SOLUTION TO VIOLATE THE ERROR CRITERION. FOR EQUATION C... (1), N = 0 AND THIS TERM DROPS OUT OF THE PROBLEM. C... COMMON/T/X,NFIN,NRUN/Y/Y1,Y2/F/DY1DX,DY2DX GO TO(1,2),NRUN C... C... INITIAL CONDITIONS FOR BESSEL*S EQUATION OF ORDER ZERO C... Y(0) = 1, DY(0)/DT = 0 1 Y1=1. Y2=0. RETURN C... C... INITIAL CONDITIONS FOR BESSEL*S EQUATION OF ORDER ONE T C... Y(0) = 0, DY(0)/DT = 1/2 2 Y1=0. Y2=0.5 RETURN END SUBROUTINE DERV COMMON/T/X,NFIN,NRUN/Y/Y1,Y2/F/DY1DX,DY2DX REAL N C... C... BESSEL*S EQUATION OF ORDER N N=FLOAT(NRUN)-1. DY1DX=Y2 DY2DX=-((1./X)*Y2+(1.-(N/X)**2)*Y1) RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/X,NFIN,NRUN/Y/Y1,Y2/F/DY1DX,DY2DX DIMENSION XPLOT(101),YPLOT(2,101) DATA IP/0/ C... C... PRINT A HEADING AT THE TOP OF THE PAGE IF(IP.NE.0)GO TO 3 GO TO(1,2),NRUN 1 WRITE(NO,4) 4 FORMAT(15H X,15H J0(X),15H DJ0(X)/DX) GO TO 3 2 WRITE(NO,5) 5 FORMAT(15H X,15H J1(X),15H DJ1(X)/DX) C... C... PRINT THE SOLUTION EVERY TENTH CALL TO SUBROUTINE PRINT 3 IP=IP+1 IF(((IP-1)/10*10).EQ.(IP-1))WRITE(NO,6)X,Y1,Y2 6 FORMAT(F15.3,2E15.4) C... C... STORE THE SOLUTION FOR SUBSEQUENT PLOTTING XPLOT(IP)=X IF(NRUN.EQ.1)YPLOT(1,IP)=Y1 IF(NRUN.EQ.1)YPLOT(2,IP)=Y2 C... C... AT THE END OF THE SECOND RUN, PLOT THE BESSEL FUNCTION VS X IF(IP.LT.101)RETURN IP=0 IF(NRUN.LT.2)RETURN CALL TPLOTS(2,101,XPLOT,YPLOT) C... C... PRINT A LABEL FOR THE PLOT WRITE(NO,8) 8 FORMAT(10H 1 - J0(X),10H 2 - J1(X)) RETURN END WYLIE, C.R., ADV ENG MATH, PP. 351-357, MCGRAW-HILL, 1966 1.0E-6 10.0 0.1 2 100 3 1 ABS 0.001 WYLIE, C.R., ADV ENG MATH, PP. 351-357, MCGRAW-HILL, 1966 1.0E-6 10.0 0.1 299999 3 1 ABS 0.001 END OF RUNS *APP PHYS16 SUBROUTINE INITAL C... C... BESSEL*S EQUATION OF ORDER ZERO (SECOND APPROACH) C... C... BESSEL*S EQUATION OF ORDER ZERO C... C... T*D2U/DT2 + DU/DT + T*U = 0 (1) C... C... IS INTEGRATED SUBJECT TO THE INITIAL CONDITIONS C... C... U(0) = 1, DU(0)/DT = 0 (2)(3) C... C... THE USUAL PROCEDURE OF STATING A SECOND-ORDER EQUATION AS TWO C... FIRST-ORDER EQUATIONS VIA THE STATE VARIABLES C... C... U1 = U, U2 = DU/DT (4)(5) C... C... GIVES C... C... DU1/DT = U2, DU2/DT = -U2/T - U1 (6)(7) C... C... WITH THE INITIAL CONDITIONS C... C... U1(0) = 1, U2(0) = 0 (8)(9) C... C... A COMPUTATIONAL PROBLEM OCCURS AT T = 0 DURING THE INTEGRATION C... OF EQUATION (7) SINCE THE FIRST TERM ON THE RIGHT HAND SIDE (RHS) C... IS INDETERMINATE, I.E., -U2/T = 0/0 AS A CONSEQUENCE OF INITIAL C... CONDITION (9). APPLICATION OF L*HOSPITAL*S RULE TO THIS TERM C... GIVES C... C... LIM (U2/T) = DU2/DT C... T---+0 C... C... (OBTAINED MERELY BY DIFFERENTIATING THE NUMERATOR AND DENOMINATOR C... OF U2/T WITH RESPECT TO T). THUS EQUATION (7) FOR T = 0 BECOMES C... C... DU2/DT = -U1/2 (10) C... C... IN THE PROGRAMMING IN SUBROUTINE DERV, A TEST FOR T = 0 DETERMINES C... IF EQUATION (7) OR EQUATION (10) IS USED. C... C... THE FOLLOWING TABULATED VALUES OF THE ZERO-ORDER BESSEL FUNCTION C... OF THE FIRST KIND, J0(T), CAN BE COMPARED WITH THE PROGRAM OUTPUT. C... ALSO, FROM THE IDENTITY DJ0(T)/DT = -J1(T), THE FOLLOWING TABU- C... LATED VALUES OF J1(T) CAN BE COMPARED WITH THE NUMERICAL DERIVA- C... TIVE OF THE SOLUTION OF EQUATION (1) PRODUCED BY THE PROGRAM C... (REFERENCE - HANDBOOK OF MATHEMATICAL FUNCTIONS, M. ABRAMOWITZ AND C... I. A. STEGUN (EDS.), NATIONAL BUREAU OF STANDARDS APPLIED MATHE- C... MATICS SERIES NO. 55, PAGE 390, U.S. GOVERNMENT PRINTING OFFICE, C... WASHINGTON, D.C., JUNE, 1964) C... C... T J0(T) J1(T) C... 0.0 1.00000 0.00000 C... 0.1 0.99750 0.04993 C... 0.2 0.99002 0.09950 C... 0.3 0.97762 0.14831 C... 0.4 0.96039 0.19602 C... 0.5 0.93846 0.24226 C... 0.6 0.91200 0.28670 C... 0.7 0.88120 0.32899 C... 0.8 0.84628 0.36884 C... 0.9 0.80752 0.40594 C... 1.0 0.76519 0.44005 C... COMMON/T/T/Y/U1,U2/F/F1,F2 C... C... INITIAL CONDITIONS (8) AND (9) U1=1. U2=0. RETURN END SUBROUTINE DERV COMMON/T/T/Y/U1,U2/F/F1,F2 C... C... EQUATION (6) F1=U2 C... C... EQUATION (7) IF(T.GE.1.0E-15)F2=-U2/T-U1 C... C... EQUATION (10) IF(T.LT.1.0E-15)F2=-U1/2. RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/T/Y/U1,U2/F/F1,F2 C... C... PRINT A HEADING IF(T.LT.0.01)WRITE(NO,2) 2 FORMAT(9X,1HT,10X,5HJ0(T),9X,6H-J1(T)) C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,1)T,U1,U2 1 FORMAT(F10.2,2F15.5) RETURN END T*D2U/DT2 + DU/DT + T*U = 0, U(0) = 1, DU(0)/DT = 0 0. 1.0 0.1 2 100 1 1 ABS 0.0001 END OF RUNS *APP PHYS17 SUBROUTINE INITAL C... C... NUMERICAL INTEGRATION OF THE NORMAL PROBABILITY DISTRIBUTION C... FUNCTION C... C... THE INTEGRAL TO BE EVALUATED IS C... C... 4 C... I = 2/SQRT(PI) INT EXP(-T**2)DT (1) C... 0 C... C... THE INTEGRAND IS THE NORMAL PROBABILITY DISTRIBUTION FUNCTION. C... THE UPPER LIMIT OF THE INTEGRAL (I.E. 4) HAS BEEN SELECTED LARGE C... ENOUGH THAT THE INTEGRAND IS SMALL (I.E. EXP(-16)). THEREFORE, C... THE INTEGRAL I WHICH IS NORMALIZED HAS A VALUE OF UNITY TO SIX C... FIGURES. THE VALUE OF I REPORTED BY SCHEID, FRANCIS, THEORY AND C... PROBLEMS OF NUMERICAL ANALYSIS, SCHAUM*S OUTLINE SERIES, MCGRAW- C... HILL BOOK COMPANY, NEW YORK, P. 137, PROB. 15.31 ARE C... C... N I C... 4 0.986 C... 6 1.000258 C... 8 1.000004 C... 10 1.000000 C... C... WHERE N IS THE NUMBER OF POINTS IN A GAUSSIAN QUADRATURE. C... C... THREE RUNS ARE PROGRAMMED FOR ALGORITHMS 3, 4 AND 7 TO GIVE AN IN- C... DICATION OF THE ACCURACY OF THE ANSWER (I.E. THE VALUE OF I CAN BE C... COMPARED FOR THE THREE RUNS). C... COMMON/T/T/Y/I/F/DIDT/CONST/A REAL I C... C... INITIAL CONDITION AND MULTIPLYING CONSTANT I=0. A=2./SQRT(4.*ATAN(1.)) RETURN END SUBROUTINE DERV COMMON/T/T/Y/I/F/DIDT REAL I C... C... EVALUATE THE INTEGRAND DIDT=EXP(-T**2) RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/T/Y/I/F/DIDT/CONST/A REAL I I=A*I WRITE(NO,1)T,I 1 FORMAT(5H T = ,F3.1,3X,4HI = ,F8.6) RETURN END INTEGRATION OF EXP(-T**2), RUNGE KUTTA MERSON ALGORITHM 0. 4.0 4.0 1 4000 3 1 ABS 0.000001 INTEGRATION OF EXP(-T**2), RUNGE KUTTA TANAKA-4 ALGORITHM 0. 4.0 4.0 1 4000 4 1 ABS 0.000001 INTEGRATION OF EXP(-T**2), RUNGE KUTTA ENGLAND ALGORITHM 0. 4.0 4.0 1 4000 7 1 ABS 0.000001 END OF RUNS *APP PHYS18 SUBROUTINE INITAL C... C... INTEGRATION OF THE SOLAR-RADIATION SPECTRUM C... C... THE INTEGRAL I DEFINED AS C... C... 7.0E-07 C... I = INT R(L)DL (1) C... 3.5E-07 C... C... IS EVALUATED NUMERICALLY FOR THE INTEGRAND R(L) DEFINING THE C... SOLAR-RADIATION SPECTRUM C... C... 3.74E-16 C... R(L) = .................. C... 5 2.52E-6/L C... L (E - 1) C... C... WHERE C... C... L SOLAR-RADIATION WAVELENGTH (METERS) C... C... R(L) SOLAR-RADIATION POWER SPECTRUM (WATTS/METER) C... C... THE INTEGRATION IS TO BE PERFORMED TO A MINIMUM ACCURACY OF THREE C... SIGNIFICANT FUGURES. A COMPLETE STATEMENT OF THE PROBLEM IS GIVEN C... IN THE REFERENCE CITED ON THE FIRST DATA CARD. C... COMMON/T/ L, NFIN, NORUN 1 /Y/ I 2 /F/ DIDL REAL L, I C... C... SET THE INITIAL CONDITION TO A SMALL NONZERO VALUE SO THAT A C... RELATIVE ERROR CAN BE USED I=1.0E-06 RETURN END SUBROUTINE DERV COMMON/T/ L, NFIN, NORUN 1 /Y/ I 2 /F/ DIDL REAL L, I C... C... EVALUATE THE INTEGRAND DIDL=R(L) RETURN END FUNCTION R(L) C... C... THE SOLAR-RADIATION SPECTRUM (WATTS/METER) IS COMPUTED BY FUNCTION C... R(L) C... C... ARGUMENT LIST C... C... L WAVELENGTH OF THE SOLAR RADIATION (METERS) C... REAL L R=3.74E-16/(L**5*(EXP(2.52E-6/L)-1.)) RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ L, NFIN, NORUN 1 /Y/ I 2 /F/ DIDL REAL L, I C... C... PRINT THE FINAL VALUE OF THE INTEGRAL IF(L.LT.6.0E-07)RETURN WRITE(NO,1)I 1 FORMAT(59H SOLAR POWER RADIATED IN THE RANGE OF VISIBLE WAVELENGTH 1 = ,E11.4,6H WATTS) RETURN END EISBERG, R. M., APPLIED MATHEMATICAL PHYSICS, PP. 27-33, ALGORITHM 1 3.5E-07 7.0E-07 3.5E-07 1 1000 1 1 REL 0.0005 EISBERG, R. M., APPLIED MATHEMATICAL PHYSICS, PP. 27-33, ALGORITHM 8 3.5E-07 7.0E-07 3.5E-07 1 1000 8 1 REL 0.0005 END OF RUNS *APP PHYS19 SUBROUTINE INITAL C... C... COUPLED LINEAR HARMONIC OSCILLATORS C... C... THE TWO SIMULTANEOUS LINEAR SECOND-ORDER DIFFERENTIAL EQUATIONS TO C... BE INTEGRATED ARE C... C... D2X1/DT2 = -(A*X1 - B*X2) (1) C... C... D2X2/DT2 = -(A*X2 - B*X1) (2) C... C... WITH THE INITIAL CONDITIONS AND PARAMETERS C... C... X1(0) = 1, DX1(0)/DT = X2 = DX2(0)/DT = 0 (3)(4)(5)(6) C... C... A = 1.25, B = 0.25 C... C... X1(T) AND X2(T) ARE PLOTTED VS T TO DEMONSTRATE THE TRANSFER OF C... ENERGY FROM ONE OSCILLATOR TO THE OTHER C... C... EQUATIONS (1) AND (2) ARE INTEGRATED NUMERICALLY BY DEFINING FOUR C... STATE VARIABLES WHICH IN TURN DEFINE AN EQUIVALENT SYSTEM OF FOUR C... FIRST-ORDER EQUATIONS C... C... X11 = X1, X12 = DX1/DT, X21 = X2, X22 = DX2/DT C... COMMON/T/T,NFIN,NORUN/Y/ X11, X12, X21, X22 1 /F/DX11DT,DX12DT,DX21DT,DX22DT C... C... THE PARAMETERS A AND B ARE SET IN A BLOCK DATA ROUTINE COMMON/PARM/ A, B C... C... INITIAL CONDITIONS (3), (4), (5), AND (6) (IN TERMS OF X11 TO X22) X11=1.0 X12=0. X21=0. X22=0. RETURN END SUBROUTINE DERV COMMON/T/T,NFIN,NORUN/Y/ X11, X12, X21, X22 1 /F/DX11DT,DX12DT,DX21DT,DX22DT COMMON/PARM/ A, B C... C... EQUATION (1) (PROGRAMMED AS TWO FIRST-ORDER EQUATIONS) DX12DT=-(A*X11-B*X21) DX11DT=X12 C... C... EQUATION (2) (PROGRAMMED AS TWO FIRST-ORDER EQUATIONS) DX22DT=-(A*X21-B*X11) DX21DT=X22 RETURN END BLOCK DATA COMMON/PARM/ A, B C... C... THE EQUATION PARAMETERS, A AND B, ARE SET IN THIS ROUTINE DATA A, B 1 / 1.25, 0.25/ END SUBROUTINE PRINT(NI,NO) COMMON/T/T,NFIN,NORUN/Y/ X11, X12, X21, X22 1 /F/DX11DT,DX12DT,DX21DT,DX22DT COMMON/PARM/ A, B C... C... DIMENSION THE ARRAYS REQUIRED FOR PLOTTING THE NUMERICAL SOLUTION DIMENSION X1P(81),X2P(81),TP(81) C... C... INITIALIZE AN INDEX FOR PLOTTING DATA IP/0/ C... C... PRINT A HEADING IF(IP.EQ.0)WRITE(NO,1) 1 FORMAT(9X,1HT,10X,5HX1(T),6X,9HDX1(T)/DT, 2 10X,5HX2(T),6X,9HDX2(T)/DT) C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,2)T,X11,X12,X21,X22 2 FORMAT(F10.1,4E15.3) C... C... STORE THE SOLUTION FOR SUBSEQUENT PLOTTING IP=IP+1 X1P(IP)=X11 X2P(IP)=X21 TP(IP)=T IF(IP.LT.81)RETURN C... C... PLOT THE NUMERICAL SOLUTION AS X1(T) AND X2(T) VS T CALL TPLOTS(1,IP,TP,X1P) WRITE(NO,3) 3 FORMAT(11H X1(T) VS T) CALL TPLOTS(1,IP,TP,X2P) WRITE(NO,4) 4 FORMAT(11H X2(T) VS T) IP=0 RETURN END EISBERG, R. M., APPLIED MATHEMATICAL PHYSICS, PP. 94-97 0. 40.0 0.5 4 150 7 1 REL 0.001 END OF RUNS *APP PHYS20 SUBROUTINE INITAL C... C... SHUTTLE LAUNCH SIMULATION C... C... THE FOLLOWING MODEL AND SIMULATION WERE DEVELOPED BY C... C... PROFESSOR S. H. JOHNSON C... DEPARTMENT OF MECHANICAL ENGINEERING C... AND MECHANICS C... LEHIGH UNIVERSITY C... PACKARD LABORATORY NO. 19 C... BETHLEHEM, PA 18015 C... C... (215) 861-4111 C... C... THE MODEL AND SIMULATION ARE FOR THE VERTICAL, PITCHOVER AND C... GRAVITY TURN PHASES OF A SPACE SHUTTLE FLIGHT. THE MODEL IS C... BASED ON THE CLASSICAL NEWTONIAN EQUATIONS OF MOTION WITH DRAG C... FORCES INCLUDED EMPIRICALLY (IN TABULAR FORM). THE EQUATIONS OF C... MOTION ARE C... C... 2 2 C... D X/DT = (GC/W)*(TLIQ*THROTL+TSOL*REDUC-DRAG)*COS(PSI) (1) C... C... 2 2 C... D Y/DT = (GC/W)*(TLIQ*THROTL*TSOL*REDUC-DRAG)*SIN(PSI) - G (2) C... C... DW/DT = NU (3) C... C... DPSI/DT = DPSIDT (4) C... C... WHERE C... C... X,Y X AND Y COORDINATES (DOWN RANGE AND ALTITUDE) OF THE C... SHUTTLE (FT) C... C... T TIME (SEC) C... C... G ACCELERATION OF GRAVITY (FT/SEC**2) C... C... GC CONVERSION FACTOR = 32.2 LBM-FT/LBF-SEC**2 C... C... TLIQ NOMINAL THRUST OF THE MAIN ENGINES (LBF) C... C... THROTL THROTTLE SETTING FOR THE MAIN ENGINES (DIMENSIONLESS) C... C... TSOL MAXIMUM THRUST OF THE SOLID BOOSTER (LBF) C... C... REDUC SOLID BOOSTER THRUST REDUCTION (DIMENSIONLESS) C... C... DRAG DRAG FORCE (LBF) C... C... W SHUTTLE WEIGHT (LBM) C... C... NU RATE OF SHUTTLE WEIGHT REDUCTION (LBM/SEC) C... C... PSI SHUTTLE TRAJECTORY ANGLE WITH RESPECT TO THE GROUND C... (RADIANS) C... C... DPSIDT TIME DERIVATIVE OF PSI (RADIANS/SEC) C... C... ADDITIONAL DIFFERENTIAL EQUATIONS FOLLOW IMMEDIATELY FROM EQUA- C... TIONS (1) AND (2) C... C... 2 2 C... D X/DT = DVX/DT, DX/DT = VX (5)(6) C... C... 2 2 C... D Y/DT = DVY/DT, DY/DT = VY (7)(8) C... C... THE ASSOCIATED ALGEBRAIC EQUATIONS ARE C... C... VARIATION OF GRAVITY WITH ALTITUDE C... C... G0 = 32.2 C... C... G = G0*(3960.0*5280.0)**2/(3960.0*5280.0+Y)**2 (9) C... C... VARIATION OF AIR DENSITY WITH ALTITUDE C... C... RHO = 0.0234/EXP((Y - 35332.0)/20229.31), Y LE 35332 (10) C... C... RHO = ((145428.0 - Y)/265816.0)**(1.0/0.2347), Y GT 35332 C... C... VARIATION OF TEMPERATURE WITH ALTITUDE C... C... TEMP = 459.67 + 59.0 - 0.003566*Y, Y LE 35332 C... (11) C... TEMP = 459.67 - 67.0, Y GT 35332 C... C... SPEED OF SOUND C... C... C = SQRT(1.4*G*53.3*TEMP) (12) C... C... SCALAR VELOCITY C... C... V = SQRT(DX/DT**2 + DY/DT**2) (13) C... C... MACH NUMBER C... C... MACH = V/C (14) C... C... DRAG COEFFICIENT C... C... CD = CD(MACH) (15) C... C... (THE FUNCTIONAL DEPENDENCY OF CD ON MACH IS GIVEN IN C... TABULAR FORM IN THE BLOCK DATA ROUTINE AS ARRAYS XM C... YCD) C... C... THROTTLE SETTING FOR THE MAIN ENGINES C... C... THROTL = 1.0, 0 LE T LE 6.5 C... C... THROTL = 1.04, 6.5 LT T LE 20 C... C... THROTL = 0.94, 20 LT T LE 36 C... C... THROTL = 0.94 36 LT T LE 45 C... C... - 0.29*(T - 36.0)/9.0 (16) C... C... THROTL = 0.65, 45 LT T LE 52 C... C... THROTL = 0.65 52 LT T LE 67 C... C... + 0.39*(T - 52.0)/15.0 C... C... THROTL = 1.04, 67 LT T C... C... SOLID BOOSTER THRUST REDUCTION C... C... REDUC = 1.0, 0 LE T LE 20 C... C... REDUC = 1.0, 20 LT T LE 45 C... C... - 0.01091*(T - 20.0) C... C... REDUC = 0.7272, 45 LT T LE 50 (17) C... C... REDUC = 0.7272 50 LT T LE 70 C... C... + 0.004545*(T - 50.0) C... C... REDUC = 0.8181 70 LT T C... C... RATE OF WEIGHT REDUCTION OF THE SHUTTLE C... C... NU = NUSOL*REDUC + NULIQ*THROTL (18) C... C... DYNAMIC PRESSURE C... C... Q = 0.5*(RHO/G)*(V**2) (19) C... C... DRAG C... C... DRAG = CD*S*Q (19) C... C... RATE OF CHANGE OF TRAJECTORY ANGLE C... C... VERTICAL ASCENT PORTION C... C... DPSIDT = 0, 0 LE T LE 6.5 C... C... CONSTANT PITCHOVER PORTION C... (20) C... DPSIDT = -0.93*PI/180.0, 6.5 LT T LE 28.01 C... C... GRAVITY TURN PORTION C... C... DPSIDT = -(G/V)*SIN(PI/2.0 - PSI), 28.01 LT T C... C... THE ALGEBRAIC AND DIFFERENTIAL EQUATIONS ARE IMPLEMENTED IN C... SUBROUTINE DERV. C... C... FINALLY, THE INITIAL CONDITIONS FOR THE MODEL, IMPLEMENTED IN C... SUBROUTINE INITAL, ARE C... C... X = Y = VX = VY = 0 (21)(22)(23)(24) C... C.. PSI = PI/2, W = W0 (25)(26) C... C... VARIABLES, PARAMETERS AND CONSTANTS NOT DEFINED PREVIOUSLY ARE C... DEFINED IN THE SUBSEQUENT CODE WHEN THEY ARE FIRST USED. C... C... DSS/2 COMMON AREA COMMON/T/ T, NSTOP, NORUN 1 /Y/ X, Y, VX, VY, W, PSI 2 /F/ DXDT, DYDT, DVXDT, DVYDT, DWDT, DPSIDT 3 /C/ G, G0, RHO, TEMP, C, V, MACH, 4 CD, THROTL, REDUC, Q, DRAG, SUMF, PI, 5 NUSOL, NULIQ, W0, TSOL, TLIQ, S, XM(37), 6 YCD(37), IP C... C... TYPE SELECTED VARIABLES REAL NUSOL, NULIQ, NU, MACH C... C... C... SET THE MODEL CONSTANTS C... C... MAXIMUM RATE OF WEIGHT LOSS OF THE SOLID BOOSTERS (LBM/SEC) NUSOL=-0.0173E+06 C... C... MAXIMUM RATE OF WEIGHT LOSS OF THE MAIN ENGINES (LBM/SEC) NULIQ=-0.0038E+06 C... C... GROUND LEVEL ACCELERATION OF GRAVITY (FT/SEC**2) G0=32.2 C... C... INITIAL SHUTTLE WEIGHT (LBM) W0=4.52E+06 C... C... MAXIMUM THRUST OF THE SOLID BOOSTERS (LBF) TSOL=6.6E+06 C... C... NOMINAL THRUST OF THE MAIN ENGINES (LBF) TLIQ=1.125E+06 C... C... SHUTTLE FRONTAL AREA (FT**2) S=600.0 C... C... PI PI=3.1415927 C... C... MODEL INITIAL CONDITIONS C... C... VELOCITY (FT/SEC) VX=0.0 VY=0.0 C... C... POSITION (FT) X=0.0 Y=0.0 C... C... WEIGHT (LBM) AND ANGLE (RADIANS) W=W0 PSI=PI/2. C... C... INITIAL DERIVATIVES CALL DERV C... C... COUNTER FOR THE PLOTTED OUTPUT IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ X, Y, VX, VY, W, PSI 2 /F/ DXDT, DYDT, DVXDT, DVYDT, DWDT, DPSIDT 3 /C/ G, G0, RHO, TEMP, C, V, MACH, 4 CD, THROTL, REDUC, Q, DRAG, SUMF, PI, 5 NUSOL, NULIQ, W0, TSOL, TLIQ, S, XM(37), 6 YCD(37), IP REAL NUSOL, NULIQ, NU, MACH C... C... MODEL ALGEBRA C... C... VARIATION OF GRAVITY WITH ALTITUDE G=G0*(3960.0*5280.0)**2/(3960.0*5280.0+Y)**2 C... C... VARIATION OF DENSITY WITH ALTITUDE IF(Y.LE.35332.0)THEN RHO=0.0234/EXP((Y-35332.0)/20929.31) ELSE IF((Y.GT.35332.0).AND.(Y.LT.145428.0))THEN RHO=((145428.0-Y)/265816.0)**(1.0/0.2347) ELSE IF(Y.GE.145428.0)THEN RHO=0. END IF C... C... VARIATION OF TEMPERATURE WITH ALTITUDE IF(Y.LE.35332.0)THEN TEMP=459.67+59.0-0.003566*Y ELSE IF(Y.GT.35332.0)THEN TEMP=459.67-67.0 END IF C... C... SPEED OF SOUND C=SQRT(1.4*G*53.3*TEMP) C... C... VELOCITY V=SQRT(VX**2+VY**2) C... C... MACH NUMBER MACH=V/C C... C... VARIATION OF DRAG COEFFICIENT WITH MACH NUMBER CD=0.16 C... C... LINEAR INTERPOLATION OVER 37 POINTS DO 1 I=1,36 IF((MACH.GT.XM(I)).AND.(MACH.LE.XM(I+1))) 1 CD=YCD(I)+(YCD(I+1)-YCD(I))*(MACH-XM(I+1))/(XM(I+1)-XM(I)) 1 CONTINUE C... C... INCREASE IN THE DRAG COEFFICIENT DUE TO DIRTY SHUTTLE STACK CD=1.0*CD C... C... THROTTLE SETTING FOR MAIN ENGINES AS A FUNCTION OF TIME IF((T.GE.0.0).AND.(T.LE.6.5))THROTL=1.0 IF((T.GT.6.5).AND.(T.LE.20.))THROTL=1.04 IF((T.GT.20.).AND.(T.LE.37.))THROTL=0.94 IF((T.GT.37.).AND.(T.LE.45.))THROTL=0.94-0.29*(T-37.0)/9.0 IF((T.GT.45.).AND.(T.LE.52.))THROTL=0.65 IF((T.GT.52.).AND.(T.LE.67.))THROTL=0.65+0.39*(T-52.0)/15.0 IF( T.GT.67.) THROTL=1.04 C... C... SOLID BOOSTER THRUST REDUCTION AS A FUNCTION OF TIME IF((T.GE.0.0).AND.(T.LE.20.))REDUC=1.0 IF((T.GT.20.).AND.(T.LE.45.))REDUC=1.0-0.01091*(T-20.0) IF((T.GT.45.).AND.(T.LE.50.))REDUC=0.7272 IF((T.GT.50.).AND.(T.LE.70.))REDUC=0.7272+0.004545*(T-50.0) IF( T.GT.70.) REDUC=0.8181 C... C... RATE OF SHUTTLE WEIGHT REDUCTION NU=NUSOL*REDUC+NULIQ*THROTL C... C... DYNAMIC PRESSURE Q=0.5*(RHO/G)*(V**2) C... C... DRAG FORCE DRAG=CD*S*Q C... C... SUM OF THE FORCES SUMF=(G/W)*(TLIQ*THROTL+TSOL*REDUC-DRAG) C... C... MODEL DIFFERENTIAL EQUATIONS C... C... HORIZONTAL AND VERTICAL ACCELERATIONS DVXDT=SUMF*COS(PSI) DVYDT=SUMF*SIN(PSI)-G C... C... HORIZONTAL AND VERTICAL VELOCITIES DXDT=VX DYDT=VY C... C... SHUTTLE WEIGHT DWDT=NU C... C... TRAJECTORY ANGLE C... C... VERTICAL PORTION IF((T.GE.0.0).AND.(T.LE.6.5))DPSIDT=0.0 C... C... CONSTANT PITCHOVER PORTION IF((T.GT.6.5).AND.(T.LE.28.01))DPSIDT=-0.93*PI/180. C... C... GRAVITY TURN POSITION IF(T.GT.28.01)DPSIDT=(-G/V)*SIN(PI/2.-PSI) C... C... THIS COMPLETES THE MODEL RETURN END BLOCK DATA C... C... THIS BLOCK DATA ROUTINE DEFINES THE TABULAR VALUES OF XM AND YCD C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ X, Y, VX, VY, W, PSI 2 /F/ DXDT, DYDT, DVXDT, DVYDT, DWDT, DPSIDT 3 /C/ G, G0, RHO, TEMP, C, V, MACH, 4 CD, THROTL, REDUC, Q, DRAG, SUMF, PI, 5 NUSOL, NULIQ, W0, TSOL, TLIQ, S, XM(37), 6 YCD(37), IP REAL NUSOL, NULIQ, NU, MACH C... C... TABULATED VALUES OF MACH NUMBER DATA XM/ 0.0, 0.2, 0.4, 0.6, 0.8, 0.85, 0.9, 1 0.95, 1.0, 1.05, 1.1, 1.15, 1.2, 1.4, 2 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5 5.8, 6.0/ C... C... TABULATED VALUES OF DRAG COEFFICIENT DATA YCD/ 0.16, 0.17, 0.18, 0.22, 0.32, 0.50, 0.58, 1 0.60, 0.62, 0.62, 0.61, 0.60, 0.59, 0.55, 2 0.51, 0.48, 0.45, 0.42, 0.39, 0.38, 0.35, 3 0.34, 0.32, 0.31, 0.30, 0.29, 0.28, 0.27, 4 0.27, 0.27, 0.26, 0.26, 0.25, 0.25, 0.25, 5 0.24, 0.24/ END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ X, Y, VX, VY, W, PSI 2 /F/ DXDT, DYDT, DVXDT, DVYDT, DWDT, DPSIDT 3 /C/ G, G0, RHO, TEMP, C, V, MACH, 4 CD, THROTL, REDUC, Q, DRAG, SUMF, PI, 5 NUSOL, NULIQ, W0, TSOL, TLIQ, S, XM(37), 6 YCD(37), IP REAL NUSOL, NULIQ, NU, MACH C... C... DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION DIMENSION TP(76), HP(76), VP(76), WP(76), TEMPP(76), PSIP(76) C... C... CONVERT THE ALTITUDE FROM FT TO NAUTICAL MILES H=Y/(5280.0*1.151) C... C... CONVERT SHUTTLE WEIGHT AND WEIGHT TIME DERIVATIVE TO TONS WT=W/2000. DWTDT=DWDT/2000. C... C... CONVERT THE TEMPERATURE FROM (R) TO (F) TEMPF=TEMP-459.67 C... C... CONVERT THE ANGLE PSI FROM RADIANS TO DEGREES PSID=PSI*180.0/PI C... C... PRINT THE SOLUTION EVERY 15TH CALL TO SUBROUTINE PRINT IP=IP+1 IF(((IP-1)/15*15).EQ.(IP-1))THEN WRITE(NO,1)T,X,Y,VX,VY,DVXDT,DVYDT 1 FORMAT(1H ,//,' TIME = ',F6.3,//, 1 ' X = ',F7.1,' Y = ',F7.1,/, 2 ' VX = ',F7.1,' VY = ',F7.1,/, 3 ' DVX/DT = ',F7.1,' DVY/DT = ',F7.1) WRITE(NO,2)Q,H,WT,DWTDT,G,RHO,TEMPF,C,MACH,CD,V,PSID 2 FORMAT( 1 ' Q = ',F7.1,' H = ',F7.1,/, 2 ' W = ',F7.1,' DW/DT = ',F7.1,/, 3 ' G = ',F7.1,' RHO = ',F7.4,/, 4 ' TEMP = ',F7.1,' C = ',F7.1,/, 5 ' MACH = ',F7.1,' CD = ',F7.3,/, 6 ' V = ',F7.1,' PSID = ',F7.1) END IF C... C... STORE THE SOLUTION FOR PLOTTING TP(IP)=T HP(IP)=H VP(IP)=MACH WP(IP)=WT TEMPP(IP)=TEMPF PSIP(IP)=PSID C... C... PLOT THE SOLUTION AT THE END OF THE RUN IF(IP.LT.76)RETURN CALL TPLOTS(1,IP,TP,HP) WRITE(NO,11) 11 FORMAT(1H ,/,' ALTITUDE (NAUTICAL MILES) VS TIME (SEC)',/) CALL TPLOTS(1,IP,TP,VP) WRITE(NO,12) 12 FORMAT(1H ,/,' MACH NUMBER VS TIME (SEC)',/) CALL TPLOTS(1,IP,TP,WP) WRITE(NO,13) 13 FORMAT(1H ,/,' WEIGHT (TONS) VS TIME (SEC)',/) CALL TPLOTS(1,IP,TP,TEMPP) WRITE(NO,14) 14 FORMAT(1H ,/,' TEMPERATURE (F) VS TIME (SEC)',/) CALL TPLOTS(1,IP,TP,PSIP) WRITE(NO,15) 15 FORMAT(1H ,/,' ANGLE (DEGREES) VS TIME (SEC)',/) RETURN END SHUTTLE LAUNCH SIMULATION 0. 75. 1. 6 1000 1 1 REL 0.001 END OF RUNS *APP PHYS21 SUBROUTINE INITAL C... C... SETTLING OF A PARTICLE IN A FLUID C... C... THE FOLLOWING ORDINARY DIFFERENTIAL EQUATION (ODE) DEFINES THE C... VELOCITY V(T) OF A PARTICLE OF RADIUS R AND DENSITY RHO WHICH C... IS SETTLING IN A FLUID OF DENSITY RHOF C... C... M*DV/DT = (4/3)*PI*(R**3)*G*(RHO - RHOF) C... (1) C.. - CD*(1/2)*RHOF*(V**2)*PI*(R**2) C... C... WHERE C... C... M MASS OF THE PARTICLE = (4/3)*PI*(R**3)*RHO C... C... CD DRAG COEFFICIENT C... C... = 24/RE, RE LE 0.5 C... C... = (0.545 + (24/RE)**0.5)**2, 0.5 LT RE LE 1.0E+04 C... C... = 0.44, RE GT 1.0E+04 C... C... RE REYNOLDS NUMBER = D*V*RHF/VIS (D = 2*R) C... C... R PARTICLE RADIUS = 0.1 CM C... C... RHO PARTICLE DENSITY = 2000 KG/M**3 C... C... RHOF FLUID DENSITY = 1000 KG/M**3 C... C... G ACCELERATION OF GRAVITY = 9.8 M/SEC**2 C... C... VIS FLUID VISOCITY = 1 CP C... C... THE FOLLOWING PROGRAM COMPUTES THE PARTICLE ACCELERATION, VELCOITY C... AND POSITION AS A FUNCTION OF TIME. C... C... COMMON AREA LINKING THE SUBROUTINES COMMON/T/ T, NSTOP, NORUN 1 /Y/ X, V 2 /F/ DXDT, DVDT 3 /C/ C1, C2, C3, RE, CD, RHOF, VIS, 4 R, IP C... C... SET THE MODEL CONSTANTS (IN CGS UNITS) R=0.1 RHO =2000.0*1000.0/(100.0**3) RHOF=1000.0*1000.0/(100.0**3) VIS=0.01 G=9.8*100.0 PI=3.1415927 VOL=(4./3.)*PI*(R**3) C1=VOL*RHO C2=VOL*G*(RHO-RHOF) C3=(1./2.)*RHOF*PI*(R**2) C... C... INITIAL CONDITIONS X=0. V=0. C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ X, V 2 /F/ DXDT, DVDT 3 /C/ C1, C2, C3, RE, CD, RHOF, VIS, 4 R, IP C... C... REYNOLDS NUMBER RE=(2.0*R)*RHOF*V/VIS C... C... DRAG COEFFICIENT CD=DRAG(RE) C... C... VELOCITY ODE DVDT=(1.0/C1)*(C2-CD*C3*(V**2)) C... C... POSITION ODE DXDT=V RETURN END FUNCTION DRAG(RE) C... C... FUNCTION DRAG COMPUTES THE DRAG COEFFICIENT AS A FUNCTION OF C... THE REYNOLDS NUMBER C... IF( RE.LE.1.0E-10)DRAG=0. IF((RE.GT.1.0E-10).AND.(RE.LE.0.5))DRAG=24.0/RE IF((RE.GT.0.5).AND.(RE.LE.1.0E+04))DRAG= 1 (0.545+SQRT(24.0/RE))**2 IF( RE.GT.1.0E+04)DRAG=0.44 RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ X, V 2 /F/ DXDT, DVDT 3 /C/ C1, C2, C3, RE, CD, RHOF, VIS, 4 R, IP C... C... PRINT A LABEL FOR THE NUMERICAL SOLUTION IF(IP.EQ.0)WRITE(NO,1) 1 FORMAT(9X,'T',11X,'X',11X,'V',10X,'RE',10X,'CD') C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,2)T,X,V,RE,CD 2 FORMAT(F10.4,2F12.4,F12.2,F12.3) IP=IP+1 RETURN END SETTLING OF A PARTICLE IN A FLUID 0. 0.2 0.02 2 1000 1 1 REL 0.001 END OF RUNS *APP PHYS22 SUBROUTINE INITAL C... C... A GOOD TEST PROBLEM IN PARTIAL DIFFERENTIAL EQUATIONS FOR THE C... NUMERICAL METHOD OF LINES C... C... THE FOLLOWING PROBLEM IS DISCUSSED IN RICE, J. R., NUMERICAL C... METHODS, SOFTWARE AND ANALYSIS, MCGRAW-HILL BOOK COMPANY, NY, C... PP 351-352. C... C... THE FOLLOWING TWO SIMULTANEOUS, NONLINEAR, ONE-DIMENSIONAL, C... INITIAL-VALUE PARTIAL DIFFERENTIAL EQUATIONS PROVIDE A GOOD C... INDICATION OF THE VERSATILITY AND FLEXIBILITY OF THE NUMERICAL C... METHOD OF LINES C... C... U = ((V - 1)*U ) + (16*X*T - 2*T - 16*(V - 1))*(U - 1) C... T X X C... (1) C... -4*X C... + 10*X*E C... C... -4*X C... V = V + U + 4*U - 4 + X**2 - 2*T - 10*T*E (2) C... T XX X C... C... E IS THE BASE OF THE NATURAL LOGARITHMS. C... C... THE BOUNDARY CONDITIONS ARE C... C... U(T,0) = V(T,0) = 1 (3)(4) C... C... 3*U(T,1) + U (T,1) = 3, 5*V (T,1) = (E**4)*(U(T,1) - 1) (5)(6) C... X X C... C... AND THE INITIAL CONDITIONS ARE C... C... U(0,X) = V(0,X) = 1 (7)(8) C... C... THE ANALYTICAL SOLUTION TO EQUATIONS (1) TO (8) AS REPORTED IN THE C... REFERENCE CITED ON THE FIRST DATA CARD IS C... C... -4*X C... U(T,X) = 1 + 10*X*T*E , V(T,X) = 1 + (X**2)*T (9)(10) C... C... THUS, EQUATIONS (9) AND (10) CAN BE USED AS AN ABSOLUTE STANDARD C... FOR THE NUMERICAL SOLUTION. C... COMMON/T/T,NFIN,NORUN 1 /Y/ U(21), V(21) 2 /F/ UT(21), VT(21) 3 /SD/ UX(21),UXX(21), VX(21),VXX(21) 4 /X/ X(21) 5 /P/ N, E4 C... SET THE NUMBER OF SPATIAL INTERVALS AND THEIR LENGTH N=20 DX=1./FLOAT(N) C... C... COMPUTE THE VALUES OF X ALONG THE SPATIAL GRID N=N+1 DO 1 I=1,N X(I)=DX*FLOAT(I-1) 1 CONTINUE C... C... SET THE INITIAL CONDITIONS, EQUATIONS (7) AND (8) DO 2 I=1,N U(I)=1. V(I)=1. 2 CONTINUE C... C... SET THE CONSTANT E**4 USED IN ONE OF THE BOUNDARY CONDITIONS E4=EXP(1.)**4 RETURN END SUBROUTINE DERV COMMON/T/T,NFIN,NORUN 1 /Y/ U(21), V(21) 2 /F/ UT(21), VT(21) 3 /SD/ UX(21),UXX(21), VX(21),VXX(21) 4 /X/ X(21) 5 /P/ N, E4 C... C... APPLY THE LEFT BOUNDARY CONDITIONS AS CONSTRAINTS, EQUATIONS (3) C... AND (4) U(1)=1. UT(1)=0. V(1)=1. VT(1)=0. C... C... COMPUTE THE FIRST-ORDER SPATIAL DERIVATIVES XL=X(1) XU=X(N) C... C... THREE-POINT CENTERED DIFFERENCES IF(NORUN.EQ.1)CALL DSS002(XL,XU,N,U,UX) IF(NORUN.EQ.1)CALL DSS002(XL,XU,N,V,VX) C... C... FIVE-POINT CENTERED DIFFERENCES IF(NORUN.EQ.2)CALL DSS004(XL,XU,N,U,UX) IF(NORUN.EQ.2)CALL DSS004(XL,XU,N,V,VX) C... C... APPLY THE RIGHT BOUNDARY CONDITIONS AS CONSTRAINTS, EQUATIONS (5) C... AND (6) UX(N)=3.-3.*U(N) VX(N)=E4*(U(N)-1.)/5. C... C... COMPUTE THE SECOND-ORDER SPATIAL DERIVATIVES C... C... THREE-POINT CENTERED DIFFERENCES IF(NORUN.EQ.1)CALL DSS002(XL,XU,N,VX,VXX) C... C... FIVE-POINT CENTERED DIFFERENCES IF(NORUN.EQ.2)CALL DSS004(XL,XU,N,VX,VXX) C... C... NOTE THAT ARRAY VX IS USED AS TEMPORARY STORAGE IN THE CALCULATION C... OF THE TERM ((V-1)*U ) WHICH IS FINALLY STORED IN ARRAY UXX C... X X DO 1 I=1,N VX(I)=(V(I)-1.)*UX(I) 1 CONTINUE IF(NORUN.EQ.1)CALL DSS002(XL,XU,N,VX,UXX) IF(NORUN.EQ.2)CALL DSS004(XL,XU,N,VX,UXX) C... C... ASSEMBLE THE PDES, EQUATIONS (1) AND (2) DO 3 I=2,N EX=EXP(-4.*X(I)) UT(I)=UXX(I)+(16.*X(I)*T-2.*T-16.*(V(I)-1.))*(U(I)-1.)+10.*X(I)*EX VT(I)=VXX(I)+UX(I)+4.*U(I)-4.+X(I)**2-2.*T-10.*T*EX 3 CONTINUE RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/T,NFIN,NORUN 1 /Y/ U(21), V(21) 2 /F/ UT(21), VT(21) 3 /SD/ UX(21),UXX(21), VX(21),VXX(21) 4 /X/ X(21) 5 /P/ N, E4 6 /A/ UA(21), DU(21), VA(21), DV(21) C... C... PRINT A HEADING FOR THE NUMERICAL AND ANALYTICAL SOLUTIONS IF(T.LT.0.01)WRITE(NO,1) 1 FORMAT(23H U = NUMERICAL U(X,T),/, 1 23H UA = ANALYTICAL U(X,T),/, 2 23H DIFF U = U - UA ,//, 3 23H V = NUMERICAL V(X,T),/, 4 23H VA = ANALYTICAL V(X,T),/, 5 23H DIFF V = V - VA ,//, 6 9X,1HT,9X,1HX,11X,1HU,10X,2HUA,6X,6HDIFF U, 7 11X,1HV,10X,2HVA,6X,6HDIFF V) C... C... COMPUTE ANALYTICAL SOLUTIONS (9) AND (10), DIFFERENCES BETWEEN THE C... NUMERICAL AND ANALYTICAL SOLUTIONS (NOTE THAT THESE CALCULATIONS C... ARE DONE FOR EVERY FOURTH GRID POINT IN X TO REDUCE THE VOLUME OF C... THE FINAL PRINTOUT) DO 2 I=1,N,4 UA(I)=1.+10.*X(I)*T*EXP(-4.*X(I)) DU(I)=U(I)-UA(I) VA(I)=1.+(X(I)**2)*T DV(I)=V(I)-VA(I) 2 CONTINUE C... C... PRINT THE NUMERICAL SOLUTION, ANALYTICAL SOLUTION, DIFFERENCE C... BETWEEN THE SOLUTIONS WRITE(NO,3)(T,X(I),U(I),UA(I),DU(I),V(I),VA(I),DV(I),I=1,N,4) 3 FORMAT(2F10.3,2F12.4,E12.3,2F12.4,E12.3) WRITE(NO,4) 4 FORMAT(1H ,/) RETURN END MADSEN, ET AL, LAPIDUS, ET AL, NUM METH FOR DIFF SYS, PP 238-239, DSS002 0. 2.0 0.4 42 1000 1 1 REL 0.001 MADSEN, ET AL, LAPIDUS, ET AL, NUM METH FOR DIFF SYS, PP 238-239, DSS004 0. 2.0 0.4 42 1000 1 1 REL 0.001 END OF RUNS *APP PHYS23 SUBROUTINE INITAL C... C... LAPLACE TRANSFORM SOLUTION OF SIMULTANEOUS ODES C... C... SOLVE THE FOLLOWING 2X2 SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS C... USING LAPLACE TRANSFORM. COMPUTE A NUMERICAL SOLUTION FOR COM- C... PARISON WITH THE ANALYTICAL SOLUTION. C... C... DY1/DT = 2*Y1 - 3*Y2, Y1(0) = 8 C... C... DY2/DT = -2*Y1 + Y2, Y2(0) = 3 C... C... LAPLACE TRANSFORMATION OF THE TWO ODES, INCLUDING THE INITIAL C... CONDITIONS, GIVES C... C... S*Y1 - 8 = 2*Y1 - 3*Y2 C... C... S*Y2 - 3 = -2*Y1 + Y2 C... C... OR, AFTER REARRANGEMENT, C... C... (S - 2)*Y1 + 3 *Y2 = 8 C... C... 2 *Y1 + (S - 1)*Y2 = 3 C... C... MULTIPLICATION OF THE FIRST EQUATION BY 2, AND THE SECOND EQUATION C... BY (S - 2), FOLLOWED BY SUBTRACTION GIVES C... C... ((S - 1)*(S - 2) - 6)*Y2 = 3*(S - 2) - 16 C... C... (S**2 - 3*S - 4)*Y2 = 3*S - 22 C... C... Y2 = (3*S - 22)/(S**2 - 3*S - 4) C... C... = (3*S - 22)/((S - R1)*(S - R2)) = A/(S - R1) + B/(S - R2) C... C... WHERE A AND B ARE CONSTANTS IN A PARTIAL FRACTIONS EXPANSION OF C... Y2. R1 AND R2 ARE THE EIGENVALUES OF THE ORIGINAL ODE SYSTEM, C... I.E., C... C... R1, R2 = + 3 + (OR -)*((-3)**2 - 4*(1)*-4)**0.5 C... -------------------------------------- C... 2 C... C... = -1, 4 C... C... A AND B CAN NOW BE EVALUATED, C... C... A = (3*S - 22)/(S - R2) = (-3 - 22)/(-1 - 4) = 5 C... S=R1 C... C... B = (3*S - 22)/(S - R1) = (12 - 22)/( 4 + 1) = -2 C... S=R2 C... C... THE INVERSE TRANSFORM FOR Y2 IS THEREFORE C... C... Y2 = 5*EXP(-T) - 2*EXP(4*T) C... C... Y1 CAN BE OBTAINED FROM SUBSTITUTION OF Y2 IN THE SECOND ODE, C... C... -5*EXP(-T) - 8*EXP(4*T) = -2*Y1 + (5*EXP(-T) - 2*EXP(4*T)) C... C... Y1 = 5*EXP(-T) + 3*EXP(4*T) C... C... THIS SOLUTION CAN BE VERIFIED BY SUBSTITUTION IN THE TWO ODES C... C... Y1(0) = 5 + 3 = 8, Y2(0) = 5 - 2 = 3 C... C... FOR THE FIRST ODE, C... C... -5*EXP(-T) + 12*EXP(4*T) = 2*(5*EXP(-T) + 3*EXP(4*T)) C... C... - 3*(5*EXP(-T) - 2*EXP(4*T)) C... C... (-5 - 10 + 15)*EXP(-T) + (12 - 6 - 6)*EXP(4*T) = 0 C... C... FOR THE SECOND ODE, C... C... -5*EXP(-T) - 8*EXP(4*T) = -2*(5*EXP(-T) + 3*EXP(4*T)) C... C... + (5*EXP(-T) - 2*EXP(4*T)) C... C... (-5 + 10 - 5)*EXP(-T) + ( -8 + 6 + 2)*EXP(4*T) = 0 C... C... A NUMERICAL SOLUTION IS COMPUTED IN THE FOLLOWING CODING FOR COM- C... PARISON WITH THE PRECEDING ANALYTICAL SOLUTION. C... C... DSS/2 COMMON COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y1, Y2 2 /F/ DY1DT, DY2DT C... C... INITIAL CONDITIONS Y1=8. Y2=3. C... C... INITIAL CONDITIONS CALL DERV RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y1, Y2 2 /F/ DY1DT, DY2DT C... C... TWO ODES DY1DT= 2.0*Y1-3.0*Y2 DY2DT=-2.0*Y1+1.0*Y2 RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y1, Y2 2 /F/ DY1DT, DY2DT C... C... PRINT A HEADING FOR THE SOLUTION IF(T.EQ.0.)WRITE(NO,1) 1 FORMAT(9X,'T',4X,'Y1 (NUM)',3X,'Y1 (ANAL)', 1 4X,'Y2 (NUM)',3X,'Y2 (ANAL)', 2 4X,'ERROR Y1',4X,'ERROR Y2') C... C... COMPUTE THE ANALYTICAL SOLUTIONS Y1A=5.0*EXP(-T)+3.0*EXP(4.0*T) Y2A=5.0*EXP(-T)-2.0*EXP(4.0*T) C... C... COMPUTE THE PERCENT DIFFERENCE BETWEEN THE NUMERICAL AND ANALYTI- C... CAL SOLUTIONS E1=(Y1-Y1A)/Y1A*100. E2=(Y2-Y2A)/Y2A*100. C... C... PRINT THE ANALYTICAL AND NUMERICAL SOLUTIONS WRITE(NO,2)T,Y1,Y1A,Y2,Y2A,E1,E2 2 FORMAT(F10.1,6F12.4) RETURN END LAPLACE TRANSFORM SOLUTION OF SIMULTANEOUS ODES 0. 1.0 0.1 2 1000 1 1 REL 0.001 END OF RUNS *APP PHYS24 SUBROUTINE INITAL C... C... DAVIDENKO*S METHOD (A DIFFERENTIAL FORM OF NEWTON*S METHOD) C... C... THE FOLLOWING SYSTEM OF SIX NONLINEAR ALGEBRAIC EQUATIONS IS C... SOLVED BY DAVIDENKO*S METHOD FOR THE SIX UNKNOWNS QC, QD, QE, C... P2, P3, P4 C... C... F1 = QC + QD - QE = 0 (1) C... C... F2 = P2 - A1 + B1*QC**2 = 0 (2) C... C... F3 = P3 - A2 + B2*QD**2 = 0 (3) C... C... F4 = 2.31*(P4 - P2) + 8.69E-4*(QC**2)*LC/DC**5 = 0 (4) C... C... F5 = 2.31*(P4 - P3) + 8.69E-4*(QD**2)*LC/DD**5 = 0 (5) C... C... F6 = 2.31*(0 - P4) + 8.69E-4*(QE**2)*LC/DE**5 C... (6) C... + (Z5 - Z4) = 0 C... C... WHERE C... C... QC, QD, QE THREE FLOWS TO BE CALCULATED C... C... P2, P3, P4 THREE PRESSURES TO BE CALCULATED C... C... A1, B1 CONSTANTS FOR PUMP NO. 1 C... C... A2, B2 CONSTANTS FOR PUMP NO. 2 C... C... LC, DC LENGTH AND DIAMETER OF PIPE SECTION C C... C... LD, DD LENGTH AND DIAMETER OF PIPE SECTION D C... C... LE, DE LENGTH AND DIAMETER OF PIPE SECTION E C... C... Z5 - Z4 VERTICAL RISE OF PIPE SECTION E C... C... THE INTERCONNECTION OF PUMPS, PIPES AND TANKS IS DESCRIBED IN THE C... REFERENCE CITED ON THE FIRST DATA LINE. C... C... THE DAVIDENKO DIFFERENTIAL EQUATION FOR EQUATIONS (1) TO (6) CAN C... BE WRITTEN IN MATRIX FORM AS C... C... - - - - - C... J*DX/DT = -F, X(T0) = X0 (7)(8) C... C... WHERE C... C... - C... J JACOBIAN MATRIX FOR EQUATIONS (1) TO (6) C... C... - C... X SOLUTION VECTOR FOR EQUATIONS (1) TO (6) WHEN C... - - C... DX/DT = 0 C... C... T INDEPENDENT VARIABLE IN THE DAVIDENKO METHOD C... C... T0 INITIAL VALUE OF T C... C... - C... X0 INITIAL CONDITION VECTOR FOR EQUATION (7) C... C... - - C... IN THE CODING IN SUBROUTINE DERV, J IS STORED IN ARRAY A, AND F IS C... STORED IN ARRAY B. THE UNCOUPLING OF THE DERIVATIVES IN EQUATION C... (7) IS DONE IN SUBROUTINE DERV BY A CALL TO SUBROUTINES DECOMP AND C... SOLVE (LISTED AND DISCUSSED IN DETAIL IN FORSYTHE, G. E., M. A. C... MALCOLM AND C. B. MOLER, COMPUTER METHODS FOR MATHEMATICAL COMPU- C... TATAIONS, CHAPTER 3, PRENTICE-HALL, ENGLEWOOD-CLIFFS, N.J., 1977). C... ANY OTHER ROUTINE(S) FOR LINEAR ALGEBRAIC EQUATIONS CAN BE USED IN C... -- - C... PLACE OF DECOMP AND SOLVE (THE GENERAL PROBLEM IS AX = B). C... C... - - C... THE PROGRAMMING OF EQUATION (7) IN TERMS OF VECTORS X, DX/DT AND C... - - C... F, AND MATRIX J, IS INTENDED TO ILLUSTRATE THE GENERAL FEATURES OF C... DAVIDENKO*S METHOD, I.E., THESE VECTORS AND MATRIX CAN BE EXPANDED C... TO ACCOMMODATE A PROBLEM OF ANY SIZE. C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(6) 2 /F/ DXDT(6) 3 /CV/ QC, QD, QE, P2, P3, 4 P4, A1, B1, A2, B2, 5 CC, CD, CE, Z5MZ4, A(6,6), 6 B(6), WORK(6), IPVT(6), N C... C... SET THE MODEL CONSTANTS A1=156.5 B1=0.00752 A2=117.1 B2=0.00427 CC=8.69E-4*125./1.278**5 CD=8.69E-4*125./2.067**5 CE=8.69E-4*145./2.469**5 Z5MZ4=70. C... C... SET THE MODEL INITIAL CONDITIONS (TRIAL VALUES OF THE UNKNOWNS) QC=50. QD=100. QE=160. P2=100. P3=80. P4=60. C... C... TRANSFER THE UNKNOWNS TO THE VECTOR X. NOTE THAT THESE STATEMENTS C... CONSTITUTE A DEFINITION OF THE ELEMENTS OF X. THIS IS ALSO AN C... IMPLEMENTATION OF INITIAL CONDITION (8) X(1)=QC X(2)=QD X(3)=QE X(4)=P2 X(5)=P3 X(6)=P4 C... C... INITIALIZE ALL OF THE MODEL CALCULATIONS N=6 CALL DERV RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN C... C... THE CALCULATION OF THE DERIVATIVE VECTOR, DX/DT, IS SLIGHTLY C... DIFFERENT IN SUBROUTINES DERV1 AND DERV2. SEE THE COMMENTS JUST C... BEFORE THE CALL TO SUBROUTINE DECOMP IN DERV1 AND DERV2 WHICH C... EXPLAIN THE DIFFERENCES IN THE CALCULATION OF DX/DT C... IF(NORUN.EQ.1)CALL DERV1 IF(NORUN.EQ.2)CALL DERV2 RETURN END SUBROUTINE DERV1 COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(6) 2 /F/ DXDT(6) 3 /CV/ QC, QD, QE, P2, P3, 4 P4, A1, B1, A2, B2, 5 CC, CD, CE, Z5MZ4, A(6,6), 6 B(6), WORK(6), IPVT(6), N COMMON/IO/ NI, NO C... C... COMPUTE THE ELEMENTS OF THE JACOBIAN MATRIX FOR EQUATIONS (1) TO C... (6) DO 1 I=1,N DO 1 J=1,N A(I,J)=0. 1 CONTINUE A(1,1)=1. A(1,2)=1. A(1,3)=-1. A(2,1)=2.*B1*X(1) A(2,4)=1. A(3,2)=2.*B2*X(2) A(3,5)=1. A(4,1)=2.*CC*X(1) A(4,4)=-2.31 A(4,6)= 2.31 A(5,2)=2.*CD*X(2) A(5,5)=-2.31 A(5,6)= 2.31 A(6,3)=2.*CE*X(3) A(6,6)=-2.31 C... C... COMPUTE THE RIGHT HAND SIDE VECTOR FOR THE DAVIDENKO DIFFERENTIAL C... EQUATIONS. NOTE THAT THESE STATEMENTS CONSTITUTE A DEFINITION OF C... THE ELEMENTS OF B (OR VECTOR F IN EQUATION (7)). THEY FOLLOW C... DIRECTLY FROM THE STATEMENT OF EQUATIONS (1) TO (6) IN SUBROUTINE C... INITAL B(1)=-(X(1)+X(2)-X(3)) B(2)=-(X(4)-A1+B1*X(1)**2) B(3)=-(X(5)-A2+B2*X(2)**2) B(4)=-(2.31*(X(6)-X(4))+CC*X(1)**2) B(5)=-(2.31*(X(6)-X(5))+CD*X(2)**2) B(6)=-(2.31*(0. -X(6))+CE*X(3)**2+Z5MZ4) C... C... SOLVE FOR THE VECTOR OF DERIVATIVES, DX/DT, IN EQUATION (7) (IN C... EFFECT, EQUATION (7) IS PREMULTIPLIED BY THE INVERSE OF JACOBIAN C... MATRIX J). IN DERV1, ARRAY B IS USED TO STORE THE RIGHT HAND SIDE C... APPROACH. IN SUBROUTINE DERV2, THE ARRAY DXDT IS FIRST USED TO C... STORE THE FUNCTION VECTOR, F. ARRAY DXDT IS THEN UPDATED BY THE C... CALL TO SUBROUTINE SOLVE SO THAT IT CONTAINS THE DERIVATIVE VECTOR C... DX/DT IN EQUATION (7) BEFORE THE RETURN FROM DERV2 CALL DECOMP(N,N,A,COND,IPVT,WORK) C... C... CHECK THE CONDITION OF MATRIX A BEFORE COMPUTING THE DERIVATIVE C... VECTOR, DX/DT. IF A IS NEAR SINGULAR, TERMINATE THE CURRENT RUN IF(COND.LT.1.0E+20)GO TO 3 WRITE(NO,4)COND 4 FORMAT(1H ,//, 1 30H THE CONDITION NUMBER OF A IS ,E10.3,33H SO THE CURRENT RUN IS G 1 TERMINATED) NSTOP=1 C... C... COMPUTE THE DERIVATIVE VECTOR DX/DT (INDIRECTLY VIA ARRAY B) 3 CALL SOLVE(N,N,A,B,IPVT) C... C... TRANSFER ARRAY B TO ARRAY DXDT SO THE DERIVATIVE VECTOR, DX/DT, C... IS AVAILABLE IN COMMON/F/ DO 2 I=1,N DXDT(I)=B(I) 2 CONTINUE RETURN END SUBROUTINE DERV2 COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(6) 2 /F/ DXDT(6) 3 /CV/ QC, QD, QE, P2, P3, 4 P4, A1, B1, A2, B2, 5 CC, CD, CE, Z5MZ4, A(6,6), 6 B(6), WORK(6), IPVT(6), N COMMON/IO/ NI, NO C... C... COMPUTE THE ELEMENTS OF THE JACOBIAN MATRIX FOR EQUATIONS (1) TO C... (6) DO 1 I=1,N DO 1 J=1,N A(I,J)=0. 1 CONTINUE A(1,1)=1. A(1,2)=1. A(1,3)=-1. A(2,1)=2.*B1*X(1) A(2,4)=1. A(3,2)=2.*B2*X(2) A(3,5)=1. A(4,1)=2.*CC*X(1) A(4,4)=-2.31 A(4,6)= 2.31 A(5,2)=2.*CD*X(2) A(5,5)=-2.31 A(5,6)= 2.31 A(6,3)=2.*CE*X(3) A(6,6)=-2.31 C... C... COMPUTE THE RIGHT HAND SIDE VECTOR FOR THE DAVIDENKO DIFFERENTIAL C... EQUATIONS (NOTE THE TEMPORARY STORAGE IN ARRAY DXDT) DXDT(1)=-(X(1)+X(2)-X(3)) DXDT(2)=-(X(4)-A1+B1*X(1)**2) DXDT(3)=-(X(5)-A2+B2*X(2)**2) DXDT(4)=-(2.31*(X(6)-X(4))+CC*X(1)**2) DXDT(5)=-(2.31*(X(6)-X(5))+CD*X(2)**2) DXDT(6)=-(2.31*(0. -X(6))+CE*X(3)**2+Z5MZ4) C... C... SOLVE FOR THE VECTOR OF DERIVATIVES, DX/DT, IN EQUATION (7) (IN C... EFFECT, EQUATION (7) IS PREMULTIPLIED BY THE INVERSE OF JACOBIAN C... MATRIX J). IN DERV1, ARRAY B IS USED TO STORE THE RIGHT HAND SIDE C... VECTOR OF FUNCTIONS, F, IN EQUATION (7), WHICH IS THE MOST OBVIOUS C... APPROACH. IN SUBROUTINE DERV2, THE ARRAY DXDT IS FIRST USED TO C... STORE THE FUNCTION VECTOR, F. ARRAY DXDT IS THEN UPDATED BY THE C... CALL TO SUBROUTINE SOLVE SO THAT IT CONTAINS THE DERIVATIVE VECTOR C... DX/DT IN EQUATION (7) BEFORE THE RETURN FROM DERV2 CALL DECOMP(N,N,A,COND,IPVT,WORK) C... C... CHECK THE CONDITION OF MATRIX A BEFORE COMPUTING THE DERIVATIVE C... VECTOR, DX/DT. IF A IS NEAR SINGULAR, TERMINATE THE CURRENT RUN IF(COND.LT.1.0E+20)GO TO 3 WRITE(NO,4)COND 4 FORMAT(1H ,//, 1 30H THE CONDITION NUMBER OF A IS ,E10.3,33H SO THE CURRENT RUN IS G 1 TERMINATED) NSTOP=1 C... C... COMPUTE THE DERIVATIVE VECTOR DX/DT 3 CALL SOLVE(N,N,A,DXDT,IPVT) C... C... AT THIS POINT, THE DERIVATIVE VECTOR, DX/DT, AS COMPUTED BY SOLVE C... IS IN ARRAY DXDT, SO ALL OF THE DERIVATIVES IN COMMON/F/ HAVE BEEN C... COMPUTED RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ X(6) 2 /F/ DXDT(6) 3 /CV/ QC, QD, QE, P2, P3, 4 P4, A1, B1, A2, B2, 5 CC, CD, CE, Z5MZ4, A(6,6), 6 B(6), WORK(6), IPVT(6), N C... C... DIMENSION THE ARRAYS TO STORE THE RESIDUALS OF EQUATIONS (1) TO C... (6) FOR SUBSEQUENT PLOTTING DIMENSION F1R(11), F2R(11), F3R(11), F4R(11), F5R(11), F6R(11), 1 TP(11) C... C... PRINT A HEADING FOR THE NUMERICAL SOLUTION IF(T.GT.0.05)GO TO 3 WRITE(NO,1) 1 FORMAT(9X,1HT,8X,2HQC,8X,2HQD,8X,2HQE,8X,2HP2,8X,2HP3,8X,2HP4,/, 1 12X,8HDX(1)/DT,2X,8HDX(2)/DT,2X,8HDX(3)/DT,2X,8HDX(4)/DT, 2 2X,8HDX(5)/DT,2X,8HDX(6)/DT,/, 3 18X,2HF1,8X,2HF2,8X,2HF3,8X,2HF4,8X,2HF5,8X,2HF6,/) IP=0 C... C... PRINT THE NUMERICAL SOLUTION. THE ELEMENTS OF B MUST FIRST BE C... RECOMPUTED SINCE B CONTAINS THE VECTOR DX/DT AFTER THE CALL TO C... SOLVE IN SUBROUTINE DERV 3 B(1)=-(X(1)+X(2)-X(3)) B(2)=-(X(4)-A1+B1*X(1)**2) B(3)=-(X(5)-A2+B2*X(2)**2) B(4)=-(2.31*(X(6)-X(4))+CC*X(1)**2) B(5)=-(2.31*(X(6)-X(5))+CD*X(2)**2) B(6)=-(2.31*(0. -X(6))+CE*X(3)**2+Z5MZ4) WRITE(NO,2)T,(X(I),I=1,N),(DXDT(I),I=1,N),(B(I),I=1,N) 2 FORMAT(F10.1,6F10.3,/,10X,6F10.3,/,10X,6F10.3,/) C... C... PLOT THE LOG OF THE RESIDUALS OF EQUATION (1) TO (6) VS T TO C... INVESTIGATE THE EXPONENTIAL CONVERGENCE OF DAVIDENKO*S ALGORITHM. C... PLOTTING IS DONE ONLY ONCE IN RUN 2 (NORUN = 2) SINCE THE NUMERI- C... CAL OUTPUT OF RUNS 1 AND 2 IS THE SAME IF(NORUN.EQ.1)RETURN C... C... STORE THE RESIDUALS FOR PLOTTING IP=IP+1 F1R(IP)=ALOG(ABS(B(1))) F2R(IP)=ALOG(ABS(B(2))) F3R(IP)=ALOG(ABS(B(3))) F4R(IP)=ALOG(ABS(B(4))) F5R(IP)=ALOG(ABS(B(5))) F6R(IP)=ALOG(ABS(B(6))) TP(IP)=T C... C... AT THE END OF RUN 2, PLOT THE RESIDUALS IF(T.LT.9.75)RETURN CALL TPLOTS(1,IP,TP,F1R) WRITE(NO,11) CALL TPLOTS(1,IP,TP,F2R) WRITE(NO,12) CALL TPLOTS(1,IP,TP,F3R) WRITE(NO,13) CALL TPLOTS(1,IP,TP,F4R) WRITE(NO,14) CALL TPLOTS(1,IP,TP,F5R) WRITE(NO,15) CALL TPLOTS(1,IP,TP,F6R) WRITE(NO,16) 11 FORMAT(1H ,//,16X,43HLOG OF RESIDUAL OF EQUATION 1 VS T ) 12 FORMAT(1H ,//,16X,43HLOG OF RESIDUAL OF EQUATION 2 VS T ) 13 FORMAT(1H ,//,16X,43HLOG OF RESIDUAL OF EQUATION 3 VS T ) 14 FORMAT(1H ,//,16X,43HLOG OF RESIDUAL OF EQUATION 4 VS T ) 15 FORMAT(1H ,//,16X,43HLOG OF RESIDUAL OF EQUATION 5 VS T ) 16 FORMAT(1H ,//,16X,43HLOG OF RESIDUAL OF EQUATION 6 VS T ) RETURN END SUBROUTINE DECOMP(NDIM,N,A,COND,IPVT,WORK) C INTEGER NDIM,N REAL A(NDIM,N),COND,WORK(N) INTEGER IPVT(N) C C DECOMPOSES A REAL MATRIX BY GAUSSIAN ELIMINATION C AND ESTIMATES THE CONDITION OF THE MATRIX. C C USE SOLVE TO COMPUTE SOLUTIONS TO LINEAR SYSTEMS. C C INPUT.. C C NDIM = DECLARED ROW DIMENSION OF THE ARRAY CONTAINING A. C C N = ORDER OF THE MATRIX. C C A = MATRIX TO BE TRIANGULARIZED. C C OUTPUT.. C C A CONTAINS AN UPPER TRIANGULAR MATRIX U AND A PERMUTED C VERSION OF A LOWER TRIANGULAR MATRIX I-L SO THAT C (PERMUTATION MATRIX)*A = L*U . C C COND = AN ESTIMATE OF THE CONDITION OF A . C FOR THE LINEAR SYSTEM A*X = B, CHANGES IN A AND B C MAY CAUSE CHANGES COND TIMES AS LARGE IN X . C IF COND+1.0 .EQ. COND , A IS SINGULAR TO WORKING C PRECISION. COND IS SET TO 1.0E+32 IF EXACT C SINGULARITY IS DETECTED. C C IPVT = THE PIVOT VECTOR. C IPVT(K) = THE INDEX OF THE K-TH PIVOT ROW C IPVT(N) = (-1)**(NUMBER OF INTERCHANGES) C C WORK SPACE.. THE VECTOR WORK MUST BE DECLARED AND INCLUDED C IN THE CALL. ITS INPUT CONTENTS ARE IGNORED. C ITS OUTPUT CONTENTS ARE USUALLY UNIMPORTANT. C C THE DETERMINANT OF A CAN BE OBTAINED ON OUTPUT BY C DET(A) = IPVT(N) * A(1,1) * A(2,2) * ... * A(N,N). C REAL EK, T, ANORM, YNORM, ZNORM INTEGER NM1, I, J, K, KP1, KB, KM1, M C IPVT(N) = 1 IF (N .EQ. 1) GO TO 80 NM1 = N - 1 C C COMPUTE 1-NORM OF A C ANORM = 0.0 DO 10 J = 1, N T = 0.0 DO 5 I = 1, N T = T + ABS(A(I,J)) 5 CONTINUE IF (T .GT. ANORM) ANORM = T 10 CONTINUE C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C DO 35 K = 1,NM1 KP1= K+1 C C FIND PIVOT C M = K DO 15 I = KP1,N IF (ABS(A(I,K)) .GT. ABS(A(M,K))) M = I 15 CONTINUE IPVT(K) = M IF (M .NE. K) IPVT(N) = -IPVT(N) T = A(M,K) A(M,K) = A(K,K) A(K,K) = T C C SKIP STEP IF PIVOT IS ZERO C IF (T .EQ. 0.0) GO TO 35 C C COMPUTE MULTIPLIERS C DO 20 I = KP1,N A(I,K) = -A(I,K)/T 20 CONTINUE C C INTERCHANGE AND ELIMINATE BY COLUMNS C DO 30 J = KP1,N T = A(M,J) A(M,J) = A(K,J) A(K,J) = T IF (T .EQ. 0.0) GO TO 30 DO 25 I = KP1,N A(I,J) = A(I,J) + A(I,K)*T 25 CONTINUE 30 CONTINUE 35 CONTINUE C C COND = (1-NORM OF A)*(AN ESTIMATE OF 1-NORM OF A-INVERSE) C ESTIMATE OBTAINED BY ONE STEP OF INVERSE ITERATION FOR THE C SMALL SINGULAR VECTOR. THIS INVOLVES SOLVING TWO SYSTEMS C OF EQUATIONS, (A-TRANSPOSE)*Y = E AND A*Z = Y WHERE E C IS A VECTOR OF +1 OR -1 CHOSEN TO CAUSE GROWTH IN Y. C ESTIMATE = (1-NORM OF Z)/(1-NORM OF Y) C C SOLVE (A-TRANSPOSE)*Y = E C DO 50 K = 1, N T = 0.0 IF (K .EQ. 1) GO TO 45 KM1 = K-1 DO 40 I = 1, KM1 T = T + A(I,K)*WORK(I) 40 CONTINUE 45 EK = 1.0 IF (T .LT. 0.0) EK = -1.0 IF (A(K,K) .EQ. 0.0) GO TO 90 WORK(K) = -(EK + T)/A(K,K) 50 CONTINUE DO 60 KB = 1, NM1 K = N - KB T = 0.0 KP1 = K+1 DO 55 I = KP1, N T = T + A(I,K)*WORK(K) 55 CONTINUE WORK(K) = T M = IPVT(K) IF (M .EQ. K) GO TO 60 T = WORK(M) WORK(M) = WORK(K) WORK(K) = T 60 CONTINUE C YNORM = 0.0 DO 65 I = 1, N YNORM = YNORM + ABS(WORK(I)) 65 CONTINUE C C SOLVE A*Z = Y C CALL SOLVE(NDIM, N, A, WORK, IPVT) C ZNORM = 0.0 DO 70 I = 1, N ZNORM = ZNORM + ABS(WORK(I)) 70 CONTINUE C C ESTIMATE CONDITION C COND = ANORM*ZNORM/YNORM IF (COND .LT. 1.0) COND = 1.0 RETURN C C 1-BY-1 C 80 COND = 1.0 IF (A(1,1) .NE. 0.0) RETURN C C EXACT SINGULARITY C 90 COND = 1.0E+32 RETURN END SUBROUTINE SOLVE(NDIM, N, A, B, IPVT) C INTEGER NDIM, N, IPVT(N) REAL A(NDIM,N),B(N) C C SOLUTION OF LINEAR SYSTEM, A*X = B . C DO NOT USE IF DECOMP HAS DETECTED SINGULARITY. C C INPUT.. C C NDIM = DECLARED ROW DIMENSION OF ARRAY CONTAINING A . C C N = ORDER OF MATRIX. C C A = TRIANGULARIZED MATRIX OBTAINED FROM DECOMP . C C B = RIGHT HAND SIDE VECTOR. C C IPVT = PIVOT VECTOR OBTAINED FROM DECOMP . C C OUTPUT.. C C B = SOLUTION VECTOR, X . C INTEGER KB, KM1, NM1, KP1, I, K, M REAL T C C FORWARD ELIMINATION C IF (N .EQ. 1) GO TO 50 NM1 = N-1 DO 20 K = 1, NM1 KP1 = K+1 M = IPVT(K) T = B(M) B(M) = B(K) B(K) = T DO 10 I = KP1, N B(I) = B(I) + A(I,K)*T 10 CONTINUE 20 CONTINUE C C BACK SUBSTITUTION C DO 40 KB = 1,NM1 KM1 = N-KB K = KM1+1 B(K) = B(K)/A(K,K) T = -B(K) DO 30 I = 1, KM1 B(I) = B(I) + A(I,K)*T 30 CONTINUE 40 CONTINUE 50 B(1) = B(1)/A(1,1) RETURN END CARNAHAN, B., ET AL, DIGITAL COMPUTING AND NUM METHODS, WILEY, PP377-378 0. 10.0 1.0 6 1000 1 1 REL 0.001 CARNAHAN, B., ET AL, DIGITAL COMPUTING AND NUM METHODS, WILEY, PP377-378 0. 10.0 1.0 6 1000 1 1 REL 0.001 END OF RUNS *APP PHYS25 SUBROUTINE INITAL C... C... FLOW THROUGH POROUS MEDIA C... C... THE FOLLOWING MODEL FOR THE MOVEMENT OF WATER THROUGH VERY DRY C... SOIL IS DISCUSSED BY FINLAYSON (1) C... C... (-S )*P = (KR(PC)*P ) , PC = - P (1)(2) C... PC T X X C... C... P(0,T) = BP1, P (1,T) = 0 (3)(4) C... X C... C... P(X,0) = BP0 (5) C... C... (1) FINLAYSON, B. A., FLOW THROUGH POROUS MEDIA, IN STIFF C... COMPUTATION (R. C. AIKEN, ED), OXFORD UNIVERSITY PRESS, C... 1985, PP. 130-134 C... C... WHERE C... C... P WATER PRESSURE C... C... PC CAPILLARY PRESSURE C... C... T TIME C... C... X POSITION C... C... S SATURATION C... C... KR RELATIVE PERMEABILITY C... C... BP0, CONSTANT INITIAL AND BOUNDARY VALUES OF P(X,T), C... BP1 RESPECTIVELY C... C... AN ALPHABETIC SUBSCRIPT IN THE PRECEDING EQUATIONS INDICATES A C... PARTIAL DERIVATIVE, E.G., P IS THE FIRST-ORDER PARTIAL DERIVATIVE C... T C... OF P WITH RESPECT TO T. C... C... KR AND S ARE GIVEN BY THE FOLLOWING EQUATIONS C... C... KR = 1/(1 + PC*L/146)**6.65 (6) C... C... S = 0.32 + 0.68/(1 + PC*L/231)**3.65 (7) C... C... WHERE L = 300 CM. C... C... THE NUMERICAL VALUES OF THE MODELS PARAMETERS ARE SET IN SUBROU- C... TINE INITAL, AND THE MODEL EQUATIONS, (1) TO (7), ARE PROGRAMMED C... IN SUBROUTINE DERV. C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ P(51) 2 /F/ PT(51) 3 /S/ PX(51), KPX(51), KPXX(51) 4 /C/ C1, C2, L, N, IP REAL L, K, KPX, KPXX C... C... SET THE MODEL PARAMETERS C... C... THE FOLLOWING VALUES OF C1 AND C2 WERE SELECTED ARBITRARILY. C... FINLAYSON'S FIGURE 3.4 SUGGESTS C1 = 0, C2 = -1, BUT WHEN THESE C... VALUES WERE USED, STRANGE NUMBERS WERE COMPUTED IN FUNCTION K C... WHICH APPEARS TO HAPPEN WHEN P AND PC ARE CLOSE TO 0. THEREFORE, C... VALUES OF C1 AND C2 WERE SELECTED TO KEEP P AND PC AWAY FROM 0 C1=-1.0 C2=-2.0 L=300. N=51 C... C... MODEL INITIAL CONDITIONS DO 1 I=1,N P(I)=C2 1 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ P(51) 2 /F/ PT(51) 3 /S/ PX(51), KPX(51), KPXX(51) 4 /C/ C1, C2, L, N, IP REAL L, K, KPX, KPXX C... C... CONSTRAIN THE PRESSURE TO PHYSICALLY MEANINGFUL VALUES DO 3 I=1,N IF(P(I).GT.0.)P(I)=0. 3 CONTINUE C... C... BOUNDARY CONDITION AT X = 0 P(1)=C1 PT(1)=0. C... C... SPATIAL DERIVATIVE P C... X CALL DSS004(0.,1.,N,P,PX) C... C... BOUNDARY CONDITION AT X = 1 PX(N)=0. C... C... PRODUCT K(P)*P C... X DO 1 I=1,N KPX(I)=K(I)*PX(I) 1 CONTINUE C... C... DERIVATIVE (K(P)*P ) C... X X CALL DSS004(0.,1.,N,KPX,KPXX) C... C... PDE DO 2 I=2,N PT(I)=-(1./SP(I))*KPXX(I) 2 CONTINUE RETURN END REAL FUNCTION K(I) C... C... FUNCTION K COMPUTES THE RELATIVE PERMEABILITY ACCORDING TO C... EQUATION (6) C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ P(51) 2 /F/ PT(51) 3 /S/ PX(51), KPX(51), KPXX(51) 4 /C/ C1, C2, L, N, IP COMMON/IO/ NI, NO REAL L, DUM, KPX, KPXX PC=-P(I) IF(PC.LT.0.)PC=0. ARG=1.+PC*L/146. K=1./ARG**6.65 RETURN END REAL FUNCTION SP(I) C... C... FUNCTION SP COMPUTES THE DERIVATIVE S ACCORDING TO EQUATION C... (7) P C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ P(51) 2 /F/ PT(51) 3 /S/ PX(51), KPX(51), KPXX(51) 4 /C/ C1, C2, L, N, IP COMMON/IO/ NI, NO REAL L, K, KPX, KPXX PC=-P(I) IF(PC.LT.0.)PC=0. ARG=1.+PC*L/231. SP=(0.68)*(-3.65)*(ARG**(-4.65))*(L/231.) RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ P(51) 2 /F/ PT(51) 3 /S/ PX(51), KPX(51), KPXX(51) 4 /C/ C1, C2, L, N, IP REAL L, K, KPX, KPXX, 1 KPXP, KPXXP C... C... DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION DIMENSION PP(3,51), PXP(3,51), 1 KPXP(3,51), KPXXP(3,51), 2 PTP(3,51), XP(51) C... C... PRINT THE NUMERICAL SOLUTION NS=10 WRITE(NO,1)T,( P(I),I=1,N,NS),( PX(I),I=1,N,NS), 1 ( KPX(I),I=1,N,NS),(KPXX(I),I=1,N,NS), 2 ( PT(I),I=1,N,NS) 1 FORMAT(1H ,//,5H T = ,F10.3,/, 1 5H P,6E11.3,/,5H PX,6E11.3,/, 2 5H KPX,6E11.3,/,5H KPXX,6E11.3,/, 3 5H PT,6E11.3,/) C... C... INITIALIZE A SECOND COUNTER FOR THE PLOTTED SOLUTION IF(IP.EQ.0)IP1=0 C... C... STORE THE SOLUTION FOR PLOTTING IP=IP+1 IF(((IP-1)/10*10).NE.(IP-1))RETURN IP1=IP1+1 DO 2 I=1,N XP(I)=FLOAT(I-1)/FLOAT(N-1) PP(IP1,I)=P(I) PXP(IP1,I)=PX(I) KPXP(IP1,I)=KPX(I) KPXXP(IP1,I)=KPXX(I) PTP(IP1,I)=PT(I) 2 CONTINUE C... C... PLOT THE NUMERICAL SOLUTION IF(IP1.LT.3)RETURN C... C... P VS X WITH T AS A PARAMETER CALL TPLOTS(3,IP,XP, PP) WRITE(NO,11) 11 FORMAT(1H ,//,' P(X,T) VS X, T = 0, 1, 2') C... C... PX VS X WITH T AS A PARAMETER CALL TPLOTS(3,IP,XP, PXP) WRITE(NO,12) 12 FORMAT(1H ,//,' PX VS X, T = 0, 1, 2') C... C... KR*PX VS X WITH T AS A PARAMETER CALL TPLOTS(3,IP,XP, KPXP) WRITE(NO,13) 13 FORMAT(1H ,//,' KR*PX VS X, T = 0, 1, 2') C... C... (KR*PX)X VS X WITH T AS A PARAMETER CALL TPLOTS(3,IP,XP,KPXXP) WRITE(NO,14) 14 FORMAT(1H ,//,' (KR*PX)X VS X, T = 0, 1, 2') C... C... PT VS X WITH T AS A PARAMETER CALL TPLOTS(3,IP,XP, PTP) WRITE(NO,15) 15 FORMAT(1H ,//,' PT VS T, T = 0, 1, 2') IP=0 IP1=0 RETURN END MOVEMENT OF WATER IN DRY SOIL 0. 2.0 0.1 51 1000 1 1 ABS 0.01 END OF RUNS *APP PHYS26 SUBROUTINE INITAL C... C... DYNAMICS OF A DOUBLE PENDULUM C... C... THE FOLLOWING SYSTEM OF LINEAR ORDINARY DIFFERENTIAL EQUATIONS C... HAS BEEN REPORTED BY KANE ET AL* FOR THE DYNAMICS OF A PLANAR C... DOUBLE PENDULUM CONSISTING OF TWO IDENTICAL, UNIFORM, PIN- C... CONNECTED RODS, EACH OF LENGTH L AND MASS M C... C... - - - - - C... M*Q + S*Q = 0 (1) C... TT C... C... WHERE THE SUBSCRIPT TT DENOTES A SECOND-ORDER DERIVATIVE WITH C... RESPECT TO TIME, T. C... C... * KANE, THOMAS R., PETER W. LIKINS AND DAVID A LEVINSON, SPACE- C... CRAFT DYNAMICS, MCGRAW-HILL BOOK COMPANY, NEW YORK, 1983, P. C... 405, PROB. 4.17 C... C... EQUATION (1) IS A MATRIX EQUATION WITH THE FOLLOWING ELEMENTS C... C... .... .. .. .. .. C... - .Q1. - .8 3. - .12 -6. C... Q = . ., M = . ., S = (SIGMA/(M*L**2). . C... .Q2. .3 2. .-6 6. C... .... .. .. .. .. C... C... WHERE Q1 AND Q2 ARE ANGULAR DISPLACEMENTS OF THE TWO PENDULUMS C... DEPICTED BELOW C... C... . Q1 Q2 C... . + + C... .... L . .... L . C... .S1.........................S2................... C... .... M .... M C... . C... . C... C... EQUAL LINEAR TORSIONAL SPRINGS, S1 AND S2, OF MODULUS SIGMA ARE C... ATTACHED AS INDICATED AND BOTH SPRINGS ARE UNDEFORMED WHEN ANGLES C... Q1 AND Q2 ARE EQUAL TO ZERO. C... C... (1) CONVERT EQUATION (1) TO A SYSTEM OF FOUR FIRST-ORDER ORDI- C... NARY DIFFERENTIAL EQUATIONS (ODES) BY DEFINING THE FOUR C... STATE VARIABLES C... C... X1 = Q1, X2 = Q1 , X3 = Q2, X4 = Q2 C... T T C... C... TAKE SIGMA/M*L**2 = 1. C... C... (2) WRITE THE SYSTEM OF FIRST-ORDER ODES IN MATRIX FORM AS C... C... - - -- - C... A*X + BX = 0 (2) C... T C... C... - .. ..T - - C... WHERE X = .X1 X2 X3 X4. AND A AND B ARE 4 X 4 MATRICES TO C... .. .. C... C... BE DETERMINED FROM PART (1). A SUPERSCRIPT T DENOTES A MATRIX C... TRANSPOSE (IN PARTICULAR, A ROW MATRIX TRANSPOSES TO A COLUMN C... MATRIX). C... C... (3) NUMERICALLY INTEGRATE EQUATION (2) FOR THE INITIAL CONDITIONS C... C... Q1(0) = 1, Q1 (0) = Q2(0) = Q2 (0) = 0 C... T T C... C... (1) THE SECOND-ORDER ODES ARE C... C... 8*Q1 + 3*Q2 + 12*Q1 - 6*Q2 = 0 C... TT TT C... C... 3*Q1 + 2*Q2 - 6*Q1 + 6*Q2 = 0 C... TT TT C... C... IF THE SUBSTITUTIONS FOR Q1 AND Q2 ARE MADE, THE PRECEDING C... SECOND-ORDER ODES REDUCE TO FOUR FIRST-ORDER ODES C... C... 8*X2 + 3*X4 + 12*X1 - 6*X3 = 0 C... T T C... C... 3*X2 + 2*X4 - 6*X1 + 6*X3 = 0 C... T T C... C... X1 = X2 C... T C... C... X3 = X4 C... T C... C... (2) THE PRECEDING FOUR FIRST-ORDER ODES CAN BE WRITTEN IN MATRIX C... FORM AS C... C... .. .. .. .. .. .. .. .. .. .. C... . . . . . . . . . . C... . 0 8 0 3 . . X1 . . 12 0 -6 0 . . X1 . . 0 . C... . . . T . . . . . . . C... . 0 3 0 2 . . X2 . . -6 0 6 0 . . X2 . . 0 . C... . .*. T . + . .*. . = . . (3) C... . 1 0 0 0 . . X3 . . 0 -1 0 0 . . X3 . . 0 . C... . . . T . . . . . . . C... . 0 0 1 0 . . X4 . . 0 0 0 -1 . . X4 . . 0 . C... . . . T . . . . . . . C... .. .. .. .. .. .. .. .. .. .. C... C... - - C... THE MATRICES A AND B ARE CLEARLY IDENTIFIED FROM THE PRECED- C... ING MATRIX EQUATION. C... C... (3) MATRIX EQUATION (2) DEFINES A SET OF IMPLICIT ODES. THE C... TERM IMPLICIT DENOTES THE COUPLING OF THE ODES THROUGH C... THE DERIVATIVES IN T, I.E., EACH ODE IN GENERAL CONTAINS C... - C... MORE THAN ONE DERIVATIVE DUE TO THE COUPLING MATRIX, A. C... C... MOST, BUT NOT ALL, NUMERICAL INTEGRATORS REQUIRE EXPLICIT C... ODES, I.E., A SYSTEM OF ODES IN WHICH EACH EQUATION CONTAINS C... ONLY ONE DERIVATIVE. THE EXCEPTIONS ARE INTEGRATORS LIKE C... LSODI AND DASSL. THUS, EQUATION (2) MUST BE PREMULTIPLIED C... --1 - C... BY A (THE INVERSE OF A) IN ORDER TO ARRIVE AT A SYSTEM OF C... EXPLICIT ODES). C... C... ANOTHER WAY TO LOOK AT THIS REQUIREMENT IS TO CONSIDER THE C... - C... COUPLING MATRIX TO BE THE IDENTITY MATRIX, I, IN THE CASE OF C... EXPLICIT ODES. THUS, C... C... --1 - - --1 - - - C... A *A*X + A *B*X = 0 C... C... OR C... C... - - --1 - - - C... I*X + A *B*X = 0 C... T C... C... IN ORDER TO ARRIVE AT THIS FINAL FORM IN THE PRESENT CASE, C... THE FIRST TWO EQUATIONS OF MATRIX EQUATION (3) CAN BE SOLVED C... SIMULTANEOUSLY FOR X2 AND X4 . IF THE FIRST EQUATION IS C... T T C... MULTIPLIED BY 2 AND THE SECOND BY 3, AND THE TWO EQUATIONS C... ARE SUBTRACTED, AN EQUATION EXPLICIT IN X2 IS OBTAINED C... T C... 7*X2 + 42*X1 - 30*X3 = 0 C... T C... C... SIMILARLY, IF THE FIRST EQUATION IS MULTIPLIED BY 3 AND C... THE SECOND BY 8, AND THE TWO EQUATIONS ARE SUBTRACTED, AN C... EQUATION EXPLICIT IN X4 IS OBTAINED C... T C... C... -7*X4 + 84*X1 - 66*X3 = 0 C... T C... C... WE NOW HAVE A SET OF FOUR EXPLICIT ODES WHICH CAN BE INTE- C... GRATED SUBJECT TO THE GIVEN INITIAL CONDITIONS. C... C... AS AN ADDITIONAL POINT, NOTE THAT A SET OF LINEARLY COUPLED C... ODES, E.G., EQUATION (3), CAN BE UNCOUPLED, I.E., CONVERTED C... FROM IMPLICIT TO EXPLICIT FORM USING A CONVENTIONAL LINEAR C... ALGEBRAIC EQUATION SOLVER. THIS IS A USEFUL TECHNIQUE IF C... ODE INTEGRATOR DESIGNED SPECIFICALLY FOR IMPLICIT ODES IS C... NOT READILY AVAILABLE. C... C... (3) THE NUMERICAL INTEGRATION OF THE EXPLICIT ODES IS PROGRAMMED C... IN THE SUBSEQUENT SUBROUTINES INITAL, DERV AND PRINT. C... C... DSS/2 COMMON AREA COMMON/T/ T, NOSTOP, NORUN 1 /Y/ X1 , X2 , X3 , X4 2 /F/ X1T, X2T, X3T, X4T 3 /C/ IP C... C... SET THE INITIAL CONDITIONS X1=1. X2=0. X3=0. X4=0. C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NOSTOP, NORUN 1 /Y/ X1 , X2 , X3 , X4 2 /F/ X1T, X2T, X3T, X4T 3 /C/ IP C... C... FOUR FIRST-ORDER EXPLICIT ODES X1T=X2 X2T=( 1.0/7.0)*(-42.0*X1+30.0*X3) X3T=X4 X4T=(-1.0/7.0)*(-84.0*X1+66.0*X3) RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NOSTOP, NORUN 1 /Y/ X1 , X2 , X3 , X4 2 /F/ X1T, X2T, X3T, X4T 3 /C/ IP C... C... DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION DIMENSION X1X3P(2,1001), X2X4P(2,1001), TP(1001) C... C... PRINT A HEADING FOR THE SOLUTION IP=IP+1 IF(IP.EQ.1)WRITE(NO,2) 2 FORMAT(9X,'T',10X,'Q1',9X,'Q1T',10X,'Q2',9X,'Q2T') C... C... PRINT THE SOLUTION EVERY 50TH CALL TO SUBROUTINE PRINT IF(((IP-1)/50*50).EQ.(IP-1))THEN WRITE(NO,1)T,X1,X2,X3,X4 1 FORMAT(F10.2,4F12.4) END IF C... C... STORE THE SOLUTION FOR PLOTTING X1X3P(1,IP)=X1 X1X3P(2,IP)=X3 X2X4P(1,IP)=X2 X2X4P(2,IP)=X4 TP(IP)=T C... C... PLOT THE SOLUTION IF(IP.LT.1001)RETURN C... C... ***************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND MUST BE MODIFIED FOR OTHER COMPUTERS OR CONVERTED TO COMMENTS C... C... PLOT Q1(T), Q2(T) VS T C... CALL QIKSAX(3,3) C... CALL QIKPLT(TP,X1X3P,IP,3H*T*,14H*Q1(T), Q2(T)*,2H**,2) C... C... PLOT Q1T(T), Q2T(T) VS T C... CALL QIKSAX(3,3) C... CALL QIKPLT(TP,X2X4P,IP,3H*T*,16H*Q1T(T), Q2T(T)*,2H**,2) C... C... ***************************************************************** C... RETURN END DYNAMICS OF TWO COUPLED PENDULUMS 0. 10.0 0.01 4 1000 1 1 REL 0.001 END OF RUNS *APP PHYS27 SUBROUTINE INITAL C... C... ONE-DIMENSIONAL, TIME-DEPENDENT MAXWELL'S EQUATIONS C... C... THE TIME-DEPENDENT MAXWELL'S EQUATIONS ARE (1) C... C... 2 C... DEL E = MU*EPS*E + MU*SIGMA*E (1) C... TT T C... C... 2 C... DEL H = MU*EPS*H + MU*SIGMA*H (2) C... TT T C... C... (1) WYLIE, C. R. AND L. C. BARRETT, ADVANCED ENGINEERING C... MATHEMATICS, MCGRAW-HILL BOOK COMPANY, NY, FIFTH EDITION, C... 1982, PP. 817-820 C... C... WHERE C... C... E ELECTRIC INTENSITY C... C... H MAGNETIC INTENSITY C... C... T TIME C... C... DEL THE LAPLACIAN OPERATOR (A SECOND-ORDER SPATIAL C... DIFFERENTIATION OPERATOR) C... C... MU PERMEABILITY C... C... EPS PERMITTIVITY C... C... SIGMA CONDUCTIVITY C... C... AN ALPHABETIC SUBSCRIPT INDICATES A PARTIAL DERIVATIVE, E.G., C... E AND E ARE THE FIRST AND SECOND-ORDER PARTIAL DERIVATIVES C... T TT C... OF E WITH RESPECT TO T, RESPECTIVELY. C... C... FOR ONE SPATIAL DIMENSION, X, IN CARTESIAN COORDINATES, EQUA- C... TIONS (1) AND (2) BECOME C... C... MU*EPS*E + MU*SIGMA*E = E (3) C... TT T XX C... C... MU*EPS*H + MU*SIGMA*H = H (4) C... TT T XX C... C... CONSIDER EQUATION (3) FOR TWO CASES C... C... (1) EPS = 0 (GOOD CONDUCTOR) C... C... COMPUTE A SOLUTION TO EQUATION (3) FOR THE INITIAL CONDITION C... C... E(X,0) = 0 (5) C... C... AND BOUNDARY CONDITIONS C... C... E(0,T) = 1, E (1,T) = 0 (6)(7) C... X C... C... USE SUBROUTINE DSS042 FOR THIS SOLUTION. PLOT E(X,T) VS X FOR C... A SERIES OF VALUES OF T. CHECK THAT YOUR SOLUTION APPROACHES C... THE ANALYTICAL SOLUTION FOR LARGE T. TAKE MU = SIGMA = 1. C... C... (2) EPS = MU = SIGMA = 1 C... C... COMPUTE A SOLUTION TO EQUATION (3) FOR THE INITIAL CONDITIONS C... C... E(X,0) = SIN(PI*X), E (X,0) = 0 (8)(9) C... T C... C... AND BOUNDARY CONDITIONS C... C... E (0,T) = E (1,T) = 0 (10)(11) C... X X C... C... USE SUBROUTINE DSS042 FOR THIS SOLUTION. PLOT E(X,T) VS X FOR C... A SERIES OF VALUES OF T. C... C... DSS/2 COMMON COMMON/T/ T, NSTOP, NORUN 1 /Y/ U(102) 2 /F/UT(102) C... C... INITIAL CONDITION(S) C... C... CASE 1 IF(NORUN.EQ.1)CALL INTAL1 C... C... CASE 2 IF(NORUN.EQ.2)CALL INTAL2 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ U(102) 2 /F/UT(102) C... C... CASE 1 IF(NORUN.EQ.1)CALL DERV1 C... C... CASE 2 IF(NORUN.EQ.2)CALL DERV2 RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ U(102) 2 /F/UT(102) C... C... PRINT AND PLOT THE SOLUTION C... C... CASE 1 IF(NORUN.EQ.1)CALL PRINT1(NI,NO) C... C... CASE 2 IF(NORUN.EQ.2)CALL PRINT2(NI,NO) RETURN END SUBROUTINE INTAL1 C... C... SUBROUTINE INTAL1 IMPLEMENTS INITIAL CONDITION (5) C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E(51) 2 /F/ ET(51) COMMON/C1/ EX(51), EXX(51), 1 EPS, MU, SIGMA, N, IP REAL MU C... C... SET THE PARAMETERS N=51 EPS=0. MU=1. SIGMA=1. C... C... INITIAL CONDITION (5) DO 1 I=1,N E(I)=0. 1 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV1 IP=0 RETURN END SUBROUTINE DERV1 C... C... SUBROUTINE DERV1 IMPLEMENTS EQUATION (3) WITH EPS = 0 (CASE 1) C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E(51) 2 /F/ ET(51) COMMON/C1/ EX(51), EXX(51), 1 EPS, MU, SIGMA, N, IP REAL MU C... C... BOUNDARY CONDITION (6) E(1)=1. ET(1)=0. C... C... BOUNDARY CONDITION (7) EX(N)=0. C... C... COMPUTE UXX CALL DSS042(0.,1.,N,E,EX,EXX,1,2) C... C... EQUATION (3) DO 1 I=2,N ET(I)=EXX(I) 1 CONTINUE RETURN END SUBROUTINE PRINT1(NI,NO) C... C... SUBROUTINE PRINT1 PRINTS AND PLOTS THE NUMERICAL SOLUTION FOR C... CASE 1 C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E(51) 2 /F/ ET(51) COMMON/C1/ EX(51), EXX(51), 1 EPS, MU, SIGMA, N, IP REAL MU C... C... DIMENSION THE ARRAYS FOR PLOTTING COMMON/P1/ XP(51), EP(6,51) C... C... PRINT A HEADING FOR THE SOLUTION IP=IP+1 IF(IP.EQ.1)WRITE(NO,1) 1 FORMAT(9X,'T',/, 1 3X,' E(0,T)',3X,' E(0.2,T)',3X,' E(0.4,T)', 2 3X,' E(0.6,T)',3X,' E(0.8,T)',3X,' E(1.0,T)') C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,2)T,(E(I),I=1,N,10) 2 FORMAT(F10.2,/,6E12.4) C... C... STORE THE SOLUTION FOR PLOTTING DO 3 I=1,N XP(I)=FLOAT(I-1)/FLOAT(N-1) EP(IP,I)=E(I) 3 CONTINUE C... C... PLOT THE SOLUTION IF(IP.LT.6)RETURN CALL TPLOTS(IP,N,XP,EP) WRITE(NO,4) 4 FORMAT(1H ,//,' E(X,T) VS X, T = 0, 0.2,..., 1.0') C... C... ***************************************************************** C... C... THE FOLLOWING CALL TO A CONTINUOUS PLOTTER IS MACHINE DEPENDENT C... AND WILL HAVE TO BE MODIFIED FOR OTHER COMPUTERS OR CONVERTED C... TO A COMMENT C... CALL PLOT1 C... C... ***************************************************************** C... RETURN END SUBROUTINE PLOT1 C... C... SUBROUTINE PLOT1 PRODUCES CONTINUOUS PLOTS OF THE SOLUTION C... FOR CASE 1 C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E(51) 2 /F/ ET(51) COMMON/C1/ EX(51), EXX(51), 1 EPS, MU, SIGMA, N, IP REAL MU COMMON/P1/ XP(51), EP(6,51) C... C... ****************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND WILL HAVE TO BE CONVERTED FOR OTHER COMPUTERS OR CHANGED TO C... COMMENTS C... C... DRAW THE CONTINUOUS PLOTS CALL QIKSAX(3,3) CALL QIKPLT(XP,EP,-N,3H*X*,8H*E(X,T)*, 1 30H*MAXWELL'S EQUATIONS, EPS = 0*,-6) C... C... ****************************************************************** C... RETURN END SUBROUTINE INTAL2 C... C... SUBROUTINE INTAL2 IMPLEMENTS INITIAL CONDITIONS (8) AND (9). C... THE PROGRAMMING IN SUBROUTINES INTAL2, DERV2 AND PRINT2 IS IN C... TERMS OF STATE VARIABLES E1 AND E2 DEFINED AS C... C... E1 = E, E2 = E C... T C... C... THESE STATE VARIABLES, AND THEIR TIME DERIVATIVES, APPEAR IN C... COMMON/Y/ AND COMMON/F/ RESPECTIVELY. THEIR USE IS REQUIRED C... SO THAT ONLY FIRST-ORDER ODES IN T ARE INTEGRATED. C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E1(51), E2(51) 2 /F/ E1T(51), E2T(51) COMMON/C2/ EX(51), EXX(51), 1 EPS, MU, SIGMA, N, IP REAL MU C... C... SET THE PARAMETERS N=51 EPS=1. MU=1. SIGMA=1. C... C... INITIAL CONDITION (8) PI=3.1415927 DO 1 I=1,N X=FLOAT(I-1)/FLOAT(N-1) E1(I)=SIN(PI*X) 1 CONTINUE C... C... INITIAL CONDITION (9) DO 2 I=1,N E2(I)=0. 2 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV2 IP=0 RETURN END SUBROUTINE DERV2 C... C... SUBROUTINE DERV2 IMPLEMENTS EQUATION (3) WITH EPS = 1 (CASE 2) C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E1(51), E2(51) 2 /F/ E1T(51), E2T(51) COMMON/C2/ EX(51), EXX(51), 1 EPS, MU, SIGMA, N, IP REAL MU C... C... BOUNDARY CONDITION (10) EX(1)=0. C... C... BOUNDARY CONDITION (11) EX(N)=0. C... C... COMPUTE UXX CALL DSS042(0.,1.,N,E1,EX,EXX,2,2) C... C... EQUATION (3) DO 1 I=1,N E1T(I)=E2(I) E2T(I)=EXX(I)-E2(I) 1 CONTINUE RETURN END SUBROUTINE PRINT2(NI,NO) C... C... SUBROUTINE PRINT2 PRINTS AND PLOTS THE NUMERICAL SOLUTION FOR C... CASE 2 C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E1(51), E2(51) 2 /F/ E1T(51), E2T(51) COMMON/C2/ EX(51), EXX(51), 1 EPS, MU, SIGMA, N, IP REAL MU C... C... DIMENSION THE ARRAYS FOR PLOTTING COMMON/P2/ XP(51), EP(6,51) IP=IP+1 C... C... PRINT A HEADING FOR THE SOLUTION IF(IP.EQ.1)WRITE(NO,1) 1 FORMAT(9X,'T',/, 1 3X,' E(0,T)',3X,' E(0.2,T)',3X,' E(0.4,T)', 2 3X,' E(0.6,T)',3X,' E(0.8,T)',3X,' E(1.0,T)') C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,2)T,(E1(I),I=1,N,10) 2 FORMAT(F10.2,/,6E12.4) C... C... STORE THE SOLUTION FOR PLOTTING DO 3 I=1,N XP(I)=FLOAT(I-1)/FLOAT(N-1) EP(IP,I)=E1(I) 3 CONTINUE C... C... PLOT THE SOLUTION IF(IP.LT.6)RETURN CALL TPLOTS(IP,N,XP,EP) WRITE(NO,4) 4 FORMAT(1H ,//,' E(X,T) VS X, T = 0, 0.2,..., 1.0') C... C... ***************************************************************** C... C... THE FOLLOWING CALL TO A CONTINUOUS PLOTTER IS MACHINE DEPENDENT C... AND WILL HAVE TO BE MODIFIED FOR OTHER COMPUTERS OR CONVERTED C... TO A COMMENT C... CALL PLOT2 C... C... ***************************************************************** C... RETURN END SUBROUTINE PLOT2 C... C... SUBROUTINE PLOT2 PRODUCES CONTINUOUS PLOTS OF THE SOLUTION C... FOR CASE 2 C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E1(51), E2(51) 2 /F/ E1T(51), E2T(51) COMMON/C2/ EX(51), EXX(51), 1 EPS, MU, SIGMA, N, IP REAL MU COMMON/P2/ XP(51), EP(6,51) C... C... ****************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND WILL HAVE TO BE CONVERTED FOR OTHER COMPUTERS OR CHANGED TO C... COMMENTS C... C... DRAW THE CONTINUOUS PLOTS CALL QIKSAX(3,3) CALL QIKPLT(XP,EP,-N,3H*X*,8H*E(X,T)*, 1 30H*MAXWELL'S EQUATIONS, EPS = 1*,-6) C... C... ****************************************************************** C... RETURN END MAXWELL'S EQUATION, EPS = 0 0. 1.0 0.2 5110000 1 1 ABS 0.01 MAXWELL'S EQUATION, EPS = 1 0. 1.0 0.2 10210000 1 1 ABS 0.01 END OF RUNS *APP PHYS28 SUBROUTINE INITAL C... C... DAMPING IN A SECOND-ORDER LINEAR SYSTEM C... C... CONSIDER THE FOLLOWING THREE CASES OF A SECOND-ORDER, LINEAR, C... CONSTANT COEFFICIENT ORDINARY DIFFERENTIAL EQUATION (ODE) C... C... 2 2 C... (2.1) D Y/DT + 3*DY/DT + Y = 0, Y(0) = 0, DY(0)/DT = 1 C... C... 2 2 C... (2.2) D Y/DT + DY/DT + Y = 0, Y(0) = 0, DY(0)/DT = 1 C... C... 2 2 C... (2.3) D Y/DT + 2*DY/DT + Y = 0, Y(0) = 0, DY(0)/DT = 1 C... C... IN ALL THREE CASES, THE LAPLACE TRANSFORM OF THE ODE, INCLUDING C... THE INITIAL CONDITIONS, IS C... C... (S**2)*Y - 1 + C*S*Y + Y = 0, C = 1, 2, 3 C... C... OR, IN FACTORED FORM C... C... (S - R1)*(S - R2)*Y = 1 C... C... A PARTIAL FRACTIONS EXPANSION THEN GIVES C... C... Y = 1/((S - R1)*(S - R2)) = A/(S - R1) + B/(S - R2) C... C... FOR WHICH THE INVERSE TRANSFORM IS C... C... Y = A*EXP(R1*T) + B*EXP(R2*T) C... C... WITH C... C... A = 1/(R1 - R2), B = 1/(R2 - R1) C... C... THE REMAINING REQUIREMENT, THEN, IS TO EVALUATE R1, R2, A AND B C... FOE THE PARTICULAR VALUES OF C. C... C... FOR CASE (2.1), C... C... -3 + (OR -) (3**2 - 4*1*1)**0.5 C... R1,R2 = ------------------------------- C... 2*1 C... C... = -0.382, -2.618 C... C... A = 1/(-0.382 - (-2.618)) = 0.4472 C... C... B = 1/(-2.618 - (-0.382)) = -0.4472 C... C... TO CHECK THE INITIAL CONDITIONS, C... C... A + B = 0 C... C... A*R1 + B*R2 = (0.4472)*(-0.382) + (-0.4472)*(-2.618) = 0.9999 C... C... FOR CASE (2.2), C... C... -1 + (OR -) (1**3 - 4*1*1)**0.5 C... R1,R2 = ------------------------------- C... 2*1 C... C... = -0.5 + 0.8660*I, -0.5 - 0.8860*I C... C... A = 1/(-0.5 + 0.8660*I - (-0.5 - 0.8660*I)) = 1/1.732*I C... C... B = 1/(-0.5 - 0.8660*I - (-0.5 + 0.8660*I)) = -1/1.732*I C... C... Y = (1/1.732*I)*EXP((-0.5 + 0.8660*I)*T) C... C... - (1/1.732*I)*EXP((-0.5 - 0.8660*I)*T) C... C... = (2/1.732)*EXP(-0.5*T)*(EXP(0.8660*T) - EXP(-0.8660*T)) C... ------------------------------- C... 2*I C... C... = 1.155*EXP(-0.5*T)*SIN(0.8660*T) C... C... TO CHECK THE INITIAL CONDITIONS, C... C... Y(0) = 0 (FROM SIN(0.8660*T) = 0) C... C... DY(0)/DT = 1.155*(-0.5*0 + 0.8660*1*1) = 1.002 C... C... FOR CASE (2.3), C... C... -2 + (OR -) (2**2 - 4*1*1)**0.5 C... R1,R2 = ------------------------------- C... 2*1 C... C... = -1 (REPEATED) C... C... THE SOLUTION IN THE T DOMAIN IS INDETERMINATE SINCE C... C... Y(T) = (1/(R1 - R2))*(EXP(R1*T) - EXP(R2*T)) C... C... = 0/0 WHEN R1 = R2 = -1 C... C... L*HOSPITAL*S RULE COULD BE APPLIED TO Y(T), OR WE CAN RETURN TO C... Y(S). C... C... THUS, IN THE S DOMAIN, C... C... Y = 1/(S - R)**2 C... C... FROM THE TRANSFORM PAIRS C... C... F(T) F(S) C... C... T 1/S**2 C... C... F(T)*EXP(A*T) F(S - A) C... C... Y(T) = T*EXP(-T) C... C... TO CHECK THE INITIAL CONDITIONS, C... C... Y(0) = 0 C... C... DY(0)/DT = EXP(-0) = 1 C... C... THE CHECKING OF THE SOLUTIONS FOR THE DIFFERENTIAL EQUATIONS IS C... LEFT TO THE READER. C... C... DSS/2 COMMON AREA COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y1, Y2 2 /F/ DY1DT, DY2DT 3 /C/ C, IP C... C... SET THE PARAMETER C IF(NORUN.EQ.1)C=1.0 IF(NORUN.EQ.2)C=2.0 IF(NORUN.EQ.3)C=3.0 C... C... SET THE INITIAL CONDITIONS Y1=0. Y2=1. C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y1, Y2 2 /F/ DY1DT, DY2DT 3 /C/ C, IP C... C... ODES DY1DT=Y2 DY2DT = 0.0-C*Y2-Y1 RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y1, Y2 2 /F/ DY1DT, DY2DT 3 /C/ C, IP C... C... DIMENSION THE ARRAYS FOR PLOTTING DIMENSION TP(101), YP(3,101) C... C... COMPUTE THE ANALYTICAL SOLUTION FOR PRINTING C... C... UNDERDAMPED IF(NORUN.EQ.1)THEN YA=1.155*EXP(-0.5*T)*SIN(0.8660*T) C... C... CRITICALLY DAMPED ELSE IF(NORUN.EQ.2)THEN YA=T*EXP(-T) C... C... OVERDAMPED ELSE IF(NORUN.EQ.3)THEN A= 0.4472 B=-0.4472 R1=-0.382 R2=-2.618 YA=A*EXP(R1*T)+B*EXP(R2*T) END IF C... C... PRINT A HEADING FOR THE SOLUTION IF(IP.EQ.0)THEN WRITE(NO,1) 1 FORMAT(9X,'T',11X,'Y',10X,'YA') END IF C... C... PRINT THE SOLUTION EVERY TENTH CALL TO SUBROUTINE PRINT IP=IP+1 IF(((IP-1)/10*10).EQ.(IP-1))THEN WRITE(NO,2)T,Y1,YA 2 FORMAT(F10.2,2F12.4) END IF C... C... STORE THE SOLUTION FOR PLOTTING TP(IP)=T YP(NORUN,IP)=Y1 C... C... PLOT THE SOLUTIONS FOR THE THREE RUNS IF((IP.LT.101).OR.(NORUN.LT.3))RETURN CALL SPLOTS(NORUN,IP,TP,YP) WRITE(NO,3) 3 FORMAT(1H ,/,' Y(T) VS T', 1 ' 1 - UNDERDAMPED 2 - CRITICALLY DAMPED 3 - OVERDAMPED') C... C... **************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND MUST BE CONVERTED FOR OTHER COMPUTERS OR CHANGED TO COMMENTS C... C... CALL QIKSAX(3,3) C... CALL QIKPLT(TP,YP,IP,3H*T*,6H*Y(T)*,9H*DAMPING*,NORUN) C... C... **************************************************************** C... RETURN END DAMPING IN A SECOND-ORDER LINEAR SYSTEM 0. 10.0 0.1 2 1000 1 1 REL 0.001 DAMPING IN A SECOND-ORDER LINEAR SYSTEM 0. 10.0 0.1 2 1000 1 1 REL 0.001 DAMPING IN A SECOND-ORDER LINEAR SYSTEM 0. 10.0 0.1 2 1000 1 1 REL 0.001 END OF RUNS *APP PHYS29 SUBROUTINE INITAL C... C... EFFECT OF DAMPING ON A MASS-SPRING SYSTEM C... C... A PARTICLE P OF MASS 2 GRAMS MOVES ON THE X AXIS AND IS ATTRATCED C... TOWARD THE ORIGIN O WITH A FORCE NUMERICALLY EQUAL TO 8*X. IF IT C... IS INITIALLY AT REST AT X = 10, FIND THE POSITION AT ANY SUBSE- C... QUENT TIME FOR TWO CASES: C... C... (1) NO OTHER FORCES ACT ON THE PARTICLE C... C... (2) A DAMPING FORCE EQUAL TO 8 TIMES THE INSTANTANEOUS C... VELOCITY ACTS. C... C... +...... X ......+ C... ........O...............P.......... C... C... DERIVE THE ANALYTICAL SOLUTION FOR X(T) FOR EACH CASE USING THE C... LAPLACE TRANSFORM. THEN SOLVE THE ODES FOR THE PARTICLE NUMERI- C... CALLY AND COMPARE WITH THE ANALYTICAL SOLUTION. PLOT THE NUMERI- C... CAL SOLUTION, X(T) VS T. C... C... THE EQUATIONS FOR THE TWO CASES ARE (1): C... C... 2 2 C... (1) 2*D X/DT = -8*X, X(0)= 10, DX(0)/DT = 0 C... C... 2 2 C... (2) 2*D X/XT = -8*X - 8*DX/DT, X(0) = 10, DX(0)/DT = 0 C... C... (1) SPIEGEL, M. R., THEORY AND PROBLEMS OF LAPLACE TRANSFORMS, C... SCHAUM PUBLISHING COMPANY, NY, NY, PP 88-89 C... C... FOR CASE (1), THE LAPLACE TRANSFORM IS C... C... (S**2)*X - 10*S + 4*X = 0 OR X = 10*S/(S**2 + 4) C... C... FOR WHICH THE INVERSE LAPLACE TRANSFORM IS C... C... X = 10*COS(2*T) C... C... FOR CASE (2), THE LAPLACE TRANSFORM IS C... C... (S**2)*X - 10*S + 4*(S*X - 10) + 4*X = 0 C... C... OR C... C... X = (10*S + 40)/(S**2 + 4*S + 4) = (10*S + 40)/(S + 2)**2 C... C... = (10*(S + 2) + 20)/(S + 2)**2 C... C... = 10/(S + 2) + 20/(S + 2)**2 C... C... FOR WHICH THE INVERSE LAPLACE TRANSFORM IS C... C... X = 10*EXP(-2*T) + 20*T*EXP(-2*T) = 10*EXP(-2*T)*(1 + 2*T) C... C... THE FOLLOWING PROGRAMMING PRODUCES A NUMERICAL SOLUTION FOR C... COMPARISON WITH THESE ANALYTICAL SOLUTIONS. IF TWO STATE C... VARIABLES ARE DEFINED AS C... C... X1 = X, X2 = DX/DT C... C... THE ORIGINAL EQUATIONS BECOME C... C... (1) DX2/DT + 4*X1 = 0, X1(0) = 10, X2(0) = 0 C... C... (2) DX2/DT + 4*X2 + 4*X1 = 0, X1(0) = 10, X2(0) = 0 C... C... NOW THESE FIRST-ORDER EQUATIONS CAN BE PROGRAMMED IN THE USUAL C... WAY. C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ X1, X2 2 /F/ DX1DT, DX2DT 3 /C/ IP C... C... INITIAL CONDITIONS X1=10. X2=0. C... C... INITIAL DERIVATIVES IP=0 CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ X1, X2 2 /F/ DX1DT, DX2DT 3 /C/ IP C... C... CASE (1) IF(NORUN.EQ.1)THEN DX2DT=-4.0*X1 DX1DT=X2 C... C... CASE (2) ELSE IF(NORUN.EQ.2)THEN DX2DT=-4.0*X1-4.0*X2 DX1DT=X2 END IF RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ X1, X2 2 /F/ DX1DT, DX2DT 3 /C/ IP C... C... DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTIONS REAL TP(101), XP(2,101) C... C... PRINT A LABEL FOR THE SOLUTIONS IP=IP+1 IF(IP.EQ.1)THEN WRITE(NO,1) 1 FORMAT(9X,'T',10X,'X1',10X,'X2',4X,'X (ANAL)',8X,'DIFF') END IF C... C... FOR EVERY TENTH CALL TO PRINT, COMPUTE THE ANALYTICAL SOLUTION, C... AND PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS AND THEIR C... DIFFERENCES IF(((IP-1)/10*10).EQ.(IP-1))THEN IF(NORUN.EQ.1)THEN XANAL=10.0*COS(2.0*T) ELSE IF(NORUN.EQ.2)THEN XANAL=10.0*EXP(-2.0*T)*(1.0+2.0*T) END IF DIFF=XANAL-X1 WRITE(NO,2)T,X1,X2,XANAL,DIFF 2 FORMAT(F10.2,4F12.5) END IF C... C... STORE THE NUMERICAL SOLUTION (TWO CASES) FOR PLOTTING TP(IP)=T XP(NORUN,IP)=X1 C... C... PLOT THE NUMERICAL SOLUTIONS AT THE END OF THE SECOND RUN IF(IP.LT.101)RETURN IP=0 IF(NORUN.LT.2)RETURN CALL TPLOTS(2,101,TP,XP) WRITE(NO,3) 3 FORMAT(1H ,/,' X(T) VS T, 1 - NO DAMPING, 2 - DAMPING') RETURN END EFFECT OF DAMPING ON A MASS-SPRING SYSTEM 0. 5. 0.05 2 1000 1 1 REL 0.001 EFFECT OF DAMPING ON A MASS-SPRING SYSTEM 0. 5. 0.05 2 1000 1 1 REL 0.001 END OF RUNS *APP PHYS30 SUBROUTINE INITAL C... C... EIGENVALUE AND TRANSIENT ANALYSIS OF A RESISTANCE-INDUCTANCE C... CIRCUIT C... C... CONSIDER THE FOLLOWING CIRCUIT DESCRIBED BY EDMINISTER (1): C... C... . C... . C... . ...... C... ....O S O..... R1 .................... C... . ...... . . C... . . ..... C... . ----+ I1 + I2 . . R2. C... .+ . ..... C... ..... ..... . C... . V . . L1. . I2 C... ..... ..... . C... .- . ..... C... . . . L2. C... . . I1 ..... C... . . . C... . . . C... ........................................ C... C... WHERE C... C... V IMPOSED VOLTAGE C... C... I1,I2 CURRENTS THROUGH THE LEFT AND RIGHT HAND LOOPS, C... RESPECTIVELY C... C... R1,R2 RESISTANCES C... C... L1,L2 INDUCTANCES C... C... S SWITCH C... C... (1) EDMINISTER, JOSEPH A., THEORY AND PROBLEMS OF ELECTRIC C... CIRCUITS, SCHAUM'S OUTLINE SERIES, MCGRAW-HILL BOOK CO., C... NEW YORK, P252 C... C... (A) USING KIRCHOFF'S VOLTAGE LAW, DERIVE THE DIFFERENTIAL EQUA- C... TIONS FOR THE DERIVATIVES DI1/DT AND DI2/DT. C... C... (B) DERIVE THE CHARACTERISTIC EQUATION FOR THIS NETWORK FROM C... THE DIFFERENTIAL EQUATIONS C... C... (C) DERIVE THE EIGENVALUES FROM THE CHARACTERISTIC EQUATION C... C... (D) COMPUTE THE EIGENVALUES NUMERICALLY USING IMSL ROUTINE EVLRG C... FOR R1 = 1, R2 = 2, L1 = 1, L2 = 4 AND COMPARE THE RESULTS C... C... (E) COMPUTE THE TRANSIENT SOLUTION TO THE DIFFERENTIAL EQUATIONS C... FOR I1(0) = I2(0) = 0, V = 1. C... C... (F) COMPARE THE NUMERICAL SOLUTION OF (E) WITH THE ANALYTICAL C... SOLUTION. C... C... KIRCHOFF'S VOLTAGE LAW APPLIED TO THE LEFT LOOP GIVES C... C... R1*(I1 + I2) + L1*DI1/DT = V C... C... AND TO THE OUTER LOOP GIVES C... C... R1*(I1 + I2) + R2*I2 + L2*DI2/DT = V C... C... OR C... C... DI1/DT + (R1/L1)*I1 + (R1/L1)*I2 = V/L1 (1) C... C... DI2/DT + (R1/L2)*I1 + (R1 + R2)/L2*I2 = V/L2 (2) C... C... IF SOLUTIONS OF THE FORM I1 = C1*EXP(E*T) AND I2 = C2*EXP(E*T) C... ARE ASSUMED FOR EQUATIONS (1) AND (2), SUBSTITUTION GIVES (WITH C... THE NONHOMOGENEOUS TERM INVOLVING V DROPPED SINCE IT DOES NOT C... BEAR ON THE EIGENVALUE ANALYSIS): C... C... (R1/L1 + E)*I1 + (R1/L1)*I2 = 0 (3) C... C... (R1/L2)*I1 + ((R1 + R2)/L2 + E)*I2 = 0 (4) C... C... EQUATIONS (3) AND (4) CAN BE PUT IN A STANDARD FORM C... C... .. .. .. .. .. .. C... . -(R1/L1) - E - R1/L1 . . I1 . . 0 . C... . . . . = . . C... . -(R1/L2) -(R1 + R2)/L2 - E . . I2 . . 0 . C... .. .. .. .. .. .. C... C... OR C... C... .. .. .. .. .. .. C... . A11 - E A12 . . I1 . . 0 . C... . . . . = . . (5) C... . A21 A22 - E . . I2 . . 0 . C... .. .. .. .. .. .. C... C... WHERE C... C... A11 = -R1/L1, A12 = -R1/L1, A21 = -R1/L2, A22 = -(R1 + R2)/L2. C... C... EQUATION (5) IS A SPECIAL CASE OF THE GENERAL HOMOGENEOUS LINEAR C... ALGEBRAIC EIGENVALUE PROBLEM: C... C... (A - E*I)X = 0 (6) C... C... EQUATION (6) HAS A NONTRIVIAL SOLUTION IF AND ONLY IF THE DETER- C... MINANT OF THE COEFFICIENT MATRIX IS ZERO, I.E., C... C... DET(A - E*I) = 0 (7) C... C... APPLYING CONDITION (7) TO EQUATION (5), C... C... (A11 - E)*(A22 - E) - A12*A21 = 0 C... C... OR (8) C... C... E**2 + (R1*L1 + R2*L1 + R1*L2)/(L1*L2)*E + (R1*R2)/(L1*L2) = 0 C... C... EQUATION (8) IS THE CHARACTERISTIC EQUATION FOR THE CIRCUIT, WHICH C... CAN BE SOLVED FOR THE EIGENVALUES C... C... E1, E2 = (- B + (OR -) (B**2 - 4*A*C)**0.5/(2*A) C... C... WHERE C... C... A 1 C... C... B (R1*L1 + R2*L1 + R1*L2)/(L1*L2) C... C... C (R1*R2)/(L1*L2) C... C... FOR THE GIVEN VALUE OF THE CIRCUIT PARAMETERS, C... C... B = (1*1 + 2*1 + 1*4)/(1*4) = 7/4, C = (1*2)/(1*4) = 1/2, AND C... C... - 7/4 + ((7/4)**2 - 4*1*(1/2))**0.5 C... E1 = ------------------------------------ C... 2 C... C... = (1/8)*(-7 + 17**0.5) = -0.3596118 C... C... - 7/4 - ((7/4)**2 - 4*1*(1/2))**0.5 C... E2 = ------------------------------------ C... 2 C... C... = (1/8)*(-7 - 17**0.5) = -1.3903882 C... C... THUS, THE SYSTEM IS STABLE AND OVERDAMPED (THE EIGENVALUES ARE C... NEGATIVE AND REAL). C... C... TO CHECK THESE EIGENVALUES, WE CAN SUBSTITUTE IN THE CHARACTER- C... ISTIC EQUATION C... C... (-0.3596118)**2 + (7/4)*(-0.3596118) + (1/2) = 0 C... C... (-1.3903882)**2 + (7/4)*(-1.3903882) + (1/2) = 0 C... C... IF EQUATIONS (1) AND (2) ARE LAPLACE TRANSFORMED: C... C... S*I1 - 0 + A11*I1 + A12*I2 = C1/S C... C... S*I2 - 0 + A21*I1 + A22*I2 = C2/S C... C... WHERE C... C... C1 = V/L1, C2 = V/L2 C... C... THROUGH SOME MINOR ALGEBRA C... C... (A11 + S)*I1 + A12*I2 = C1/S (9) C... C... A21*I1 + (A22 + S)*I2 = C2/S (10) C... C... NOTE THE SIMILARITY OF THE LEFT HAND SIDES OF EQUATIONS (3) AND C... (4), AND (9) AND (10). C... C... IF EQUATIONS (9) AND (10) ARE SOLVED FOR I1 AND I2, C... C... (A22 + S)*C1 - A12*C2 C... I1 = -------------------------------------------- (11) C... S*(S**2 + (A11+ A22)*S + (A11*A22 - A12*A21)) C... C... (A11 + S)*C2 - A21*C1 C... I2 = -------------------------------------------- (12) C... S*(S**2 + (A11+ A22)*S + (A11*A22 - A12*A21)) C... C... NOTE THAT EQUATIONS (11) AND (12) HAVE THE SAME DENOMINATOR, AND C... ONE OF THE FACTORS OF THE DENOMINATOR IS THE CHARACTERISTIC C... POLYNOMIAL FOR THE SYSTEM DISCUSSED PREVIOUSLY. C... C... THE ANALYTICAL SOLUTION (I1(T) AND I2(T)) IS OBTAINED BY INVERTING C... EQUATIONS (11) AND (12). THIS CAN BE DONE BY PARTIAL FRACTIONS: C... C... I1 = A/S + B/(S - E1) + C/(S - E2) C... C... I2 = D/S + E/(S - E1) + F/(S - E2) C... C... WHERE A TO F ARE CONSTANTS TO BE DETERMINED AND E1 AND E2 ARE THE C... ROOTS OF THE CHARACTERISTIC POLYNOMIAL, I.E., EQUATION (8). C... C... A CAN BE DETERMINED BY MULTIPLYING EQUATION (11) BY S, THEN C... LETTING S = 0, I.E., C... C... A = (A22*C1 - A12*C2)/(A11*A22 - A12*A21) C... C... SIMILARLY, B CAN BE DETERMINED BY MULTIPLYING EQUATION (11) BY C... (S - E1), THEN LETTING S = E1, C... C... B = ((A22 + E1)*C1 - A12*C2)/(E1*(E1 - E2)) C... C... PROCEEDING IN THE SAME WAY, WE OBTAIN THE OTHER CONSTANTS AS C... C... C = ((A22 + E2)*C1 - A12*C2)/(E2*(E2 - E1)) C... C... D = (A11*C2 - A21*C1)/(A11*A22 - A12*A21) C... C... E = ((A11 + E1)*C2 - A21*C1)/(E1*(E1 - E2)) C... C... F = ((A11 + E2)*C2 - A21*C1)/(E2*(E2 - E1)) C... C... THE INVERSE LAPLACE TRANSFORMS OF EQUATIONS (11) AND (12) ARE C... C... I1(T) = A + B*EXP(E1*T) + C*EXP(E2*T) (13) C... C... I2(T) = D + E*EXP(E1*T) + F*EXP(E2*T) (14) C... C... FINALLY, THE STEADY STATE SOLUTION CAN BE COMPUTED FROM EQUATIONS C... (1) AND (2) BY SETTING DI1/DT = DI2/DT = O, I.E., C... C... A11*I1 + A12*I2 = C1 (15) C... C... A21*I1 + A22*I2 = C2 (16) C... C... OR, SOLVING FOR I1 AND I2, C... C... I1 = (A22*C1 - A12*C2)/(A11*A22 - A12*A21) (17) C... C... I2 = (A11*C2 - A21*C1)/(A11*A22 - A12*A21) (18) C... C... THIS STEADY STATE SOLUTION CAN BE USED TO CHECK THE NUMERICAL AND C... ANALYTICAL SOLUTIONS. C... C... COMMON AREA FOR INPUT/OUTPUT (I/O) UNIT NUMBERS COMMON/IO/ NI, NO C... DSS/2 COMMON AREA COMMON/T/ T 1 /Y/ I1, I2 2 /F/ DI1DT, DI2DT 3 /C/ L1, L2, R1, R2, V, IP REAL L1, L2, I1, I2 C... C... DIMENSION THE ARRAYS REQUIRED BY SUBROUTINE EVLRG REAL A(2,2) COMPLEX E(2) C... C... CIRCUIT PARAMETERS L1=1.0 L2=4.0 R1=1.0 R2=2.0 V =1.0 C... C... ELEMENTS OF THE COEFFICIENT MATRIX A(1,1)=-R1/L1 A(1,2)=-R1/L1 A(2,1)=-R1/L2 A(2,2)=-(R1+R2)/L2 C... C... PRINT THE COEFFICIENT (JACOBIAN) MATRIX WRITE(NO,3)((A(I,J),J=1,N),I=1,N) 3 FORMAT(1H ,//,' COEFFICIENT (JACOBIAN) MATRIX',//, 1 (2F12.4,/)) C... C... PARAMETERS FOR SUBROUTINE EVLRG N=2 LDA=2 C... C... CALL IMSL SUBROUTINE EVLRG TO COMPUTE THE EIGENVALUES CALL EVLRG(N,A,LDA,E) C... C... PRINT THE EIGENVALUES DO 1 I=1,N EREAL= REAL(E(I)) EIMAG=AIMAG(E(I)) WRITE(NO,2)I,EREAL,EIMAG 2 FORMAT(1H ,//, 1 ' EIGENVALUE ',I3,2X,' REAL = ',F17.14,' IMAG = ',F17.14) 1 CONTINUE C... C... INITIAL CONDITIONS FOR EQUATIONS (1) AND (2) I1=0. I2=0. C... C... INITIAL DERIVATIVES CALL DERV IP=0 END SUBROUTINE DERV COMMON/T/ T 1 /Y/ I1, I2 2 /F/ DI1DT, DI2DT 3 /C/ L1, L2, R1, R2, V, IP REAL L1, L2, I1, I2 C... C... EQUATION (1) DI1DT=V/L1-(R1/L1)*I1-(R1/L1)*I2 C... C... EQUATION (2) DI2DT=V/L2-(R1/L2)*I1-(R1+R2)/L2*I2 RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T 1 /Y/ I1, I2 2 /F/ DI1DT, DI2DT 3 /C/ L1, L2, R1, R2, V, IP REAL L1, L2, I1, I2, 1 I1A, I2A, I1SS, I2SS C... C... DIMENSION THE ARRAYS FOR PLOTTING DIMENSION TP(101), PI(2,101) C... C... PRINT A HEADING FOR THE SOLUTION IP=IP+1 IF(IP.EQ.1)THEN WRITE(NO,1) 1 FORMAT(//,9X,'T',10X,'I1',10X,'I2',6X,'DI1/DT',6X,'DI2/DT',/, 1 13X,'I1 (ANAL)',3X,'I2 (ANAL)',5X,'I1 DIFF',5X,'I2 DIFF') END IF C... C... PRINT THE SOLUTION EVERY FIFTH CALL TO PRINT IF(((IP-1)/5*5).EQ.(IP-1))THEN WRITE(NO,2)T,I1,I2,DI1DT,DI2DT 2 FORMAT(F10.2,4F12.4) C... C... COMPUTE AND PRINT THE ANALYTICAL SOLUTION AND THE DIFFERENCE C... BETWEEN THE NUMERICAL AND ANALYTICAL SOLUTIONS E1=-0.3596118 E2=-1.3903882 A11=R1/L1 A12=R1/L1 A21=R1/L2 A22=(R1+R2)/L2 C1=V/L1 C2=V/L2 A=(A22*C1-A12*C2)/(A11*A22-A12*A21) B=((A22+E1)*C1-A12*C2)/(E1*(E1-E2)) C=((A22+E2)*C1-A12*C2)/(E2*(E2-E1)) D=(A11*C2-A21*C1)/(A11*A22-A12*A21) E=((A11+E1)*C2-A21*C1)/(E1*(E1-E2)) F=((A11+E2)*C2-A21*C1)/(E2*(E2-E1)) EXP1=EXP(E1*T) EXP2=EXP(E2*T) I1A=A+B*EXP1+C*EXP2 I2A=D+E*EXP1+F*EXP2 DI1=I1-I1A DI2=I2-I2A WRITE(NO,3)I1A,I2A,DI1,DI2 3 FORMAT(10X,4F12.4,/) END IF C... C... STORE THE SOLUTIONS FOR PLOTTING TP(IP)=T PI(1,IP)=I1 PI(2,IP)=I2 C... C... COMPUTE THE STEADY STATE SOLUTIONS FROM EQUATIONS (17) AND (18) IF(IP.LT.101)RETURN I1SS=(A22*C1-A12*C2)/(A11*A22-A12*A21) I2SS=(A11*C2-A21*C1)/(A11*A22-A12*A21) WRITE(NO,4)I1SS,I2SS 4 FORMAT(//,4X,' I1SS = ',F12.4,3X,'I2SS = ',F12.4) C... C... PLOT THE NUMERICAL SOLUTIONS CALL TPLOTS(2,IP,TP,PI) WRITE(NO,5) 5 FORMAT(1H ,//,' I1 AND I2 VS T') C... C... ****************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND MUST BE CONVERTED FOR OTHER COMPUTERS OR CHANGED TO COMMENTS C... CALL QIKSAX(3,3) CALL QIKPLT(TP,PI,IP,3H*T*,14H*I1(T), I2(T)*,2H**,2) C... C... ****************************************************************** C... RETURN END EIGENVALUE AND TRANSIENT ANALYSIS OF A R-L CIRCUIT 0. 10.0 0.1 2 1000 1 1 REL 0.001 END OF RUNS *APP PHYS31 SUBROUTINE INITAL C... C... LINEAR STABILITY ANALYSIS OF A SYSTEM OF NONLINEAR ORDINARY C... DIFFERENTIAL EQUATIONS C... C... THE FOLLOWING FIRST ORDER, NONLINEAR DIFFERENCE EQUATIONS, C... REPORTED BY HENON, ARE DESCRIBED BY STRANG (1): C... C... N+1 N N 2 C... U = U - (7/5)*(U ) + 1 (1) C... 1 2 1 C... C... N+1 N C... U = (3/10)*U (2) C... 2 1 C... C... (1) STRANG, G., INTRODUCTION TO APPLIED MATHEMATICS, C... WELLESLEY-CAMBRIDGE PRESS, WELLESLEY, MA, PP. 504- C... 505, 1986. C... C... THE SOLUTION TO EQUATIONS (1) AND (2), PLOTTED AS A PHASE PLANE C... (U1 VS U2) EXHIBITS THE INTERESTING PROPERTIES OF HENON'S C... 'STRANGE ATTRACTOR' (E.G., CONVERGENCE TO A COMMON POINT WITH C... INFINITELY FINE DETAIL). C... C... IF EQUATIONS (1) AND (2) ARE REWRITTEN AS: C... C... N+1 N N N 2 N C... U - U = U - (7/5)*(U ) + 1 - U (3) C... 1 1 2 1 1 C... --------- C... (N+1) - N C... C... C... N+1 N N N C... U - U = (3/10)*U - U (4) C... 2 2 1 2 C... --------- C... (N+1) - N C... C... SINCE N IS A DISCRETE VARIABLE WITH SUCCESSIVE VALUES DIFFERING C... BY ONE, THE LEFT HAND SIDES OF EQUATIONS (3) AND (4) CAN BE CON- C... SIDERED DISCRETE (E.G., FINITE DIFFERENCE) APPROXIMATIONS OF THE C... DERIVATIVES DU /DT AND DU /DT. THIS SUGGESTS TWO DIFFERENTIAL C... 1 2 C... EQUATIONS: C... C... 2 C... DU /DT = U - (7/5)*U + 1 - U (5) C... 1 2 1 1 C... C... DU /DT = (3/10)*U - U (6) C... 2 1 2 C... C... THE FOLLOWING CODING NUMERICALLY INTEGRATES EQUATIONS (5) AND C... (6) FOR A SERIES OF INITIAL CONDITIONS, AS EXPLAINED IN THE C... COMMENTS. THE SOLUTION IS PLOTTED AS A PHASE PLANE (U2 VS U1) C... AND THE DERIVATIVES ARE ALSO PLOTTED IN THE SAME FORMAT (DU /DT C... VS DU /DT). 2 C... 1 C... C... THE SOLUTION HAS THE FOLLOWING PROPERTIES: C... C... (1) ALL OF THE SOLUTION CURVES (FOR 25 INITIAL CONDITIONS) C... CONVERGE TO THE SAME POINT, U1 = 0.63135, U2 = 0.18941, C... WHICH APPEARS TO CORRESPOND TO STRANG'S FIGURE 6.14. C... THIS IS AN EQUILIBRIUM POINT FOR EQUATIONS (5) AND (6), C... I.E., DU /DT = DU /DT = 0. C... 1 2 C... C... (2) THE SOLUTION HAS THE 'NEVER ENDING NESTING OF DETAIL' C... DESCRIBED BY STRANG SINCE EACH SMALLER SECTION OF THE C... PHASE PLANE IN THE NEIGHBORHOOD OF THE EQUILIBRIUM POINT C... WILL HAVE THE SAME NUMBER OF SOLUTION CURVES (25 IN THIS C... CASE) CONVERGING TO THE EQUILIBRIUM POINT. C... C... (3) THE EQUILIBRIUM SOLUTION TO THE DIFFERENTIAL EQUATIONS C... IS ALSO THE EQUILIBRIUM SOLUTION TO THE HENON DIFFERNCE C... EQUATIONS (TO WITHIN THE ACCURACY OF THE NUMERICAL INTE- C... GRATION). C... C... THESE RESULTS SUGGEST THAT SIMILAR RESULTS CAN BE PRODUCED FOR C... ANY DIFFERENTIAL EQUATIONS WITH AN EQUILIBRIUM POINT TO WHICH C... ALL SOLUTION CURVES CONVERGE FROM A GIVEN SET OF INITIAL CONDI- C... TIONS. C... C... TO EXPLORE THIS POINT FURTHER USING THE PRESENT SYSTEM OF EQUA- C... TIONS (5) AND (6), IF THE DERIVATIVES ARE SET TO ZERO, THE C... RESULTING SYSTEM OF TWO NONLINEAR ALGEBRAIC EQUATIONS CAN EASILY C... BE SOLVED FOR U1 AND U2: C... C... 2 C... 0 = U2 - (7/5)*U1 + 1 - U1 (7) C... C... 0 = (3/10)*U1 - U2 (8) C... C... SOLVING EQUATION (8) FOR U2 AND SUBSTITUTING IN EQUATION (7) C... GIVES A QUADRATIC EQUATION IN U1 C... C... 2 C... 0 = (3/10)*U1 - (7/5)*U1 + 1 - U1 C... C... WHICH HAS THE ROOTS U1 = 0.63135446 AND -1.1313545. THE CORRES- C... VALUES OF U2 FROM EQUATION (8) ARE U2 = 0.18940634 AND -0.33940 C... 635. C... C... AS THE OUTPUT FROM THE FOLLOWING PROGRAM INDICATES, FOR THE 25 C... INITIAL CONDITIONS, ALL OF THE SOLUTIONS CONVERGE TO THE POINT C... U1 = 0.63135446, U2 =0.18940634. HOWEVER, WHEN THESE INITIAL C... CONDITIONS ARE EXTENDED CLOSER TO THE POINT U1 = -1.1313545, C... U2 = -0.33940635, THE NUMERICAL SOLUTION IMMEDIATELY BECOMES C... UNSTABLE. (NUMERICAL OUTPUT COULD NOT BE COMPUTED FOR THE FIRST C... TIME STEP BEYOND THE INITIAL CONDITION U1(0) = -2, U2(0) = -0.4. C... THIS CASE CAN EASILY BE PROGRAMMED BELOW BY LETTING -2 LE U1(0) C... LE 2 IN INCREMENTS OF 1 AND -0.4 LE U2(0) LE 0.4 IN INCREMENTS C... OF 0.2). C... C... THIS INSTABILITY CAN BE INVESTIGATED IN THE FOLLOWING WAY. THE C... JACOBIAN MATRIX OF EQUATIONS (5) AND (6) IS: C... C... .. .. C... . A11 A12 . C... J = . . C... . A21 A22 . C... .. .. C... C... WHERE C... C... A11 = - 1 - 2*(7/5)*U1, A12 = 1, A21 = 3/10, A22 = -1 C... C... THE EIGENVALUE PROBLEM FOR J IS: C... C... .. .. C... . A11 - E A12 . C... DET. . = 0 C... . A21 A22 - E . C... .. .. C... C... OR C... C... E**2 - (A11 + A12)*E + (A11*A22 - A12*A21) = 0 C... C... SUBSTITUTING FOR A11 TO A22, C... C... E**2 - (- 1 - 2*(7/5)*U1 + 1)*E + C... C... ((- 1 - 2*(7/5)*U1)*(-1) - (1)*(-1)) = 0 C... C... OR C... C... E**2 + A*U1*E + A*U1 = 0 (9) C... C... WHERE A = 2*(7/5). C... C... IF EQUATION (9) IS FACTORED FOR THE EIGENVALUES, C... C... - A*U1 + (OR -) ((A*U1**2) - 4*1*A*U1)**0.5 C... E1, E2 = ------------------------------------------- (10) C... 2*1 C... C... FROM EQUATION (10), WE SEE THAT IF U1 GT 0, E.G., THE PREVIOUS C... EQUILIBRIUM POINT U1 = 0.63135446, E1 AND E2 ARE REAL AND NEGA- C... TIVE. THIS EXPLAINS WHY THE SOLUTIONS IN THE NEIGHBORHOOD OF C... OF U1 = 0.63135446, U2 = 0.18940634 ARE STABLE. C... C... ON THE OTHER HAND, IF U1 LT 0, E.G., THE PREVIOUS EQUILIBRIUM C... POINT U1 = -1.1313545, E1 AND E2 ARE REAL, BUT E1 IS POSITIVE C... WHILE E2 IS NEGATIVE. THEREFORE, E1 IS AN UNSTABLE EIGENVALUE C... WHICH EXPLAINS WHY THE SOLUTIONS IN THE NEIGHBORHOOD OF U1 = C... -1.1313545, U2 = -0.33940635 ARE UNSTABLE. C... C... TO SUMMARIZE, C... C... U1 = 0.63135446, U2 = 0.18940634 C... C... IS A STABLE EQUILIBRIUM POINT, WHILE C... C... U1 = -1.1313545, U2 = -0.33940635 C... C... IS AN UNSTABLE EQUILIBRIUM POINT. C... C... THE BASIS FOR THE PRECEDING STABILTY ANALYSIS IS BRIEFLY DESCRIBED C... BELOW. A GENERAL SYSTEM OF N ORDINARY DIFFERENTIAL EQUATIONS C... (ODES) IN N DEPENDENT VARIABLES IS: C... C... - - - - - C... DY/DT = F(Y,T), Y(TO) = Y0 (11)(12) C... C... WHERE C... C... - T C... Y = (Y1, Y2,...,YN) (13) C... C... - - C... SIMILARLY, F IS A VECTOR OF DERIVATIVE FUNCTIONS, Y0 IS A VECTOR C... OF INITIAL CONDITIONS AND T IS THE INDEPENDENT VARIABLE (EXCEPT C... IN EQUATION (13) WHERE IT DENOTES A TRANSPOSE). C... C... IF A TAYLOR SERIES EXPANSION IS WRITTEN FOR THE RIGHT HAND SIDE C... OF EQUATION (11), FOLLOWED BY TRUNCATION AFTER THE LINEAR TERMS: C... C... DY/DT = F(YE,TE) + J(YE,TE)*(Y - YE) (14) C... C... WHERE E DENOTES AN EQUILIBRIUM POINT ABOUT WHICH THE EXPANSION C... (LINEARIZATION) IS MADE, J IS THE JACOBIAN MATRIX OF F (AND C... OVERBARS DENOTING VECTORS AND MATRICES HAVE BEEN DROPPED). C... C... IF EQUATION (14) IS WRITTEN AS: C... C... D(Y - YE)/DT = J(YE,TE)*(Y - YE) (15) C... C... WHERE WE HAVE USED THE DEFINITION OF THE EQUILIBRIOUM POINT, C... I.E., DYE/DT = F(YE,TE) = 0. EQUATION (15) IS LINEAR IN THE C... PERTURBATION VARIABLE Y - YE, AND THEREFORE HAS THE EXPONENTIAL C... SOLUTION C... C... Y - YE = C*EXP(L*T) (16) C... C... SUBSTITUTION OF EQUATION (16) IN EQUATION (15) GIVES A SYSTEM OF C... HOMOGENEOUS ALGEBRAIC EQUATIONS C... C... (J - L*I)*C = 0 (17) C... C... WHERE I IS THE IDENTITY MATRIX. C... C... EQUATION (17) IS THE WELL KNOWN LINEAR ALGEBRAIC EIGENVALUE C... PROBLEM, WITH NONTRIVIAL SOLUTIONS (C NE 0) IF AND ONLY IF THE C... DETERMINANT OF THE COEFFICIENT MATRIX IS ZERO, C... C... DET(J - L*I) = 0 (18) C... C... SOLUTION OF EQUATION (18) GIVES THE EIGENVALUES OF THE LINEAR C... ODE SYSTEM (EQUATION (14)), L1, L2,..., LN. IF ALL N EIGENVALUES C... HAVE NEGATIVE REAL PARTS, EQUATION (14) HAS A STABLE SOLUTION. C... IF ONE OR MORE EIGENVALUES HAS A POSITIVE REAL PART, THE SOLU- C... TION TO EQUATION (14) IS UNSTABLE. BOTH SITUATIONS OCCURRED C... IN THE PRECEDING EXAMPLE (DEPENDING ON WHETHER THE EQUILIBRIUM C... VALUE OF U1 WAS POSITIVE OR NEGATIVE). C... C... ALSO, THIS IS CLEARLY A LINEAR ANALYSIS, AND THE CONCLUSIONS C... CONCERNING STABILTY MAY NOT APPLY TO THE ORIGINAL NONLINEAR C... PROBLEM, EQUATION (11). NOTE ALSO THAT THE EIGENVALUES COMPUTED C... BY THIS LINEAR ANALYSIS ARE NOT CONSTANT, BUT IN GENERAL, FOR A C... A NONLINEAR PROBLEM, WILL DEPEND ON THE EQUILIBRIUM POINT ABOUT C... WHICH THE LINEARIZATION IS MADE (E.G., THE EQUILIBRIUM VALUE OF C... U1 IN THE PRECEDING EXAMPLE). C... C... THE PROGRAMMING FOR THE STABLE CASE (U1 EQUILIBRIUM GT 0) OF THE C... PRECEDING EXAMPLE FOLLOWS. C... C... WES C... 25OCT87 C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ U1, U2 2 /F/ DU1DT, DU2DT 3 /C/ IP C... C... SET THE INITIAL CONDITIONS SO THAT -1 LE U1(0) LE 1 IN INCREMENTS C... OF 0.5, -0.35 LE U2(0) LE 0.25 IN INCREMENTS OF 0.15, WHICH, IN C... ALL COMBINATIONS, GIVES A TOTAL OF 25 INITIAL CONDITIONS. THESE C... INITIAL CONDITIONS ARE PRINTED AT THE BEGINNING OF EACH RUN (T = C... .00 OUTPUT) DATA U10,U20/-1.5,-0.35/ C... C... RESET U1(0) TO -1.5 AT THE BEGINNING OF EVERY FIFTH RUN, AND INCRE- C... MENT U2(0) BY 0.15 IF(((NORUN-6)*(NORUN-11)*(NORUN-16)*(NORUN-21)).EQ.0)THEN U10=-1.5 U20=U20+0.15 END IF C... C... INCREMENT U1(0) BY 0.5 U10=U10+0.5 C... C... SET THE INITIAL CONDITIONS U1=U10 U2=U20 C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ U1, U2 2 /F/ DU1DT, DU2DT 3 /C/ IP C... C... DIFFERENTIAL EQUATIONS DU1DT=U2-(7.0/5.0)*U1**2+1.0-U1 DU2DT=(3.0/10.0)*U1-U2 RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ U1, U2 2 /F/ DU1DT, DU2DT 3 /C/ IP C... C... DIMENSION THE ARRAYS FOR PLOTTING REAL U1P(25,201), U2P(25,201), 1 DU1DTP(25,201), DU2DTP(25,201), 2 X(201), Y(201) C... C... INITIALIZE THE VARIABLES FOR SAVING THE MINIMUM AND MAXIMUM C... VALUES OF U1, U2, DU1/DT AND DU2/DT DATA U1MIN,U1MAX,DU1MIN,DU1MAX, 1 U2MIN,U2MAX,DU2MIN,DU2MAX 2 /8*0.0/ C... C... PRINT A HEADING FOR THE SOLUTION IP=IP+1 IF(IP.EQ.1)THEN WRITE(NO,1) 1 FORMAT(9X,'T',10X,'U1',10X,'U2',6X,'DU1/DT',6X,'DU2/DT') END IF C... C... PRINT THE SOLUTION EVERY 20TH CALL TO SUBROUTINE PRINT IF(((IP-1)/20*20).EQ.(IP-1))THEN WRITE(NO,2)T,U1,U2,DU1DT,DU2DT 2 FORMAT(F10.2,4F12.5) END IF C... C... STORE THE SOLUTION FOR PLOTTING AS A PHASE PLANE U1P(NORUN,IP)=U1 U2P(NORUN,IP)=U2 DU1DTP(NORUN,IP)=DU1DT DU2DTP(NORUN,IP)=DU2DT C... C... SAVE THE MINUMUM AND MAXIMUM VALUES OF U1, U2, DU1/DT AND DU2/DT C... TO ASSIST IN THE SCALING OF THE PLOTS IF(U1.LT.U1MIN)U1MIN=U1 IF(U1.GT.U1MAX)U1MAX=U1 IF(U2.LT.U2MIN)U2MIN=U2 IF(U2.GT.U2MAX)U2MAX=U2 IF(DU1DT.LT.DU1MIN)DU1MIN=DU1DT IF(DU1DT.GT.DU1MAX)DU1MAX=DU1DT IF(DU2DT.LT.DU2MIN)DU2MIN=DU2DT IF(DU2DT.GT.DU2MAX)DU2MAX=DU2DT C... C... TEST FOR THE LAST RUN WHEN THE SOLUTIONS ARE PLOTTED IF((IP.LT.201).OR.(NORUN.LT.25))RETURN C... C... PRINT THE MIMIMUM AND MAXIMUM VALUES FOR SCALING THE PLOTS WRITE(NO,3) U1MIN, U1MAX, U2MIN, U2MAX, 1 DU1MIN,DU1MAX,DU2MIN,DU2MAX 3 FORMAT(1H1,////, 1 ' MAXIMA AND MINIMA OF U1, U2, DU1/DT AND DU2/DT',//, 2 5X,' U1MIN = ',F7.3,' U1MAX = ',F7.3,/, 3 5X,' U2MIN = ',F7.3,' U2MAX = ',F7.3,/, 4 5X,' DU1MIN = ',F7.3,' DU1MAX = ',F7.3,/, 5 5X,' DU2MIN = ',F7.3,' DU2MAX = ',F7.3) C... C... PRINT AND COMPUTE THE RESIDUAL OF THE HENON DIFFERENCE EQUATIONS C... AT THE EQUILIBRIUM POINT FOR THE DIFFERENTIAL EQUATIONS U1=0.63135446 U2=0.18940634 R1=U1-U2+(7.0/5.0)*U1**2-1.0 R2=U2-(3.0/10.0)*U1 WRITE(NO,10)U1,U2,R1,R2 10 FORMAT(1H ,////, 1 ' RESIDUALS OF THE HENON DIFFERENCE EQUATIONS',//, 2 5X,' U1 = ',F7.5,' U2 = ',F7.5,/, 3 5X,' RES 1 = ',F7.5,' RES 2 = ',F7.5) C... C... CONTINUOUS PLOT OF U2(T) VS U1(T) C... C... **************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND MUST BE CONVERTED FOR OTHER COMPUTERS OR CHANGED TO COMMENTS C... C... SCALE THE AXES C... CALL QIKSET(6.0,-1.5,0.5,8.0,-0.4,0.1) C... CALL QIKSAX(3,3) C... C... FIRST CURVE DO 4 J=1,IP X(J)=U1P(1,J) Y(J)=U2P(1,J) 4 CONTINUE C... CALL QIKPLT(X,Y,IP,7H*U1(T)*,7H*U2(T)*, C... 1 33H*DIFFERENTIAL FORM OF HENON EQNS*,1) C... CALL PLOT(-7.0,1.0,-3) C... C... SUBSEQUENT CURVES DO 6 I=2,NORUN DO 5 J=1,IP X(J)=U1P(I,J) Y(J)=U2P(I,J) 5 CONTINUE C... CALL QLINE(X,Y,IP,M) 6 CONTINUE C... C... CONTINUOUS PLOT OF DU2(T)/DT VS DU1(T)/DT C... C... SCALE THE AXES C... CALL QIKSET(6.0,-2.0,0.5,8.0,-0.8,0.2) C... CALL QIKSAX(3,3) C... C... FIRST CURVE DO 7 J=1,IP X(J)=DU1DTP(1,J) Y(J)=DU2DTP(1,J) 7 CONTINUE C... CALL QIKPLT(X,Y,IP,8H*DU1/DT*,8H*DU2/DT*, C... 1 33H*DIFFERENTIAL FORM OF HENON EQNS*,1) C... CALL PLOT(-7.0,1.0,-3) C... C... SUBSEQUENT CURVES DO 9 I=2,NORUN DO 8 J=1,IP X(J)=DU1DTP(I,J) Y(J)=DU2DTP(I,J) 8 CONTINUE C... CALL QLINE(X,Y,IP,M) 9 CONTINUE C... C... **************************************************************** C... RETURN END DIFFERENTIAL FORM OF THE HENON EQUATIONS 0. 20.0 0.1 2 1000 1 1 REL 0.001 REPEATS 24 END OF RUNS *APP PHYS32 SUBROUTINE INITAL C... C... ONE-DIMENSIONAL, TIME-DEPENDENT MAXWELL'S EQUATIONS - EXPLICIT C... PROGRAMMING OF THE SPATIAL DERIVATIVE C... C... THE TIME-DEPENDENT MAXWELL'S EQUATIONS ARE (1) C... C... 2 C... DEL E = MU*EPS*E + MU*SIGMA*E (1) C... TT T C... C... 2 C... DEL H = MU*EPS*H + MU*SIGMA*H (2) C... TT T C... C... (1) WYLIE, C. R. AND L. C. BARRETT, ADVANCED ENGINEERING C... MATHEMATICS, MCGRAW-HILL BOOK COMPANY, NY, FIFTH EDITION, C... 1982, PP. 817-820 C... C... WHERE C... C... E ELECTRIC INTENSITY C... C... H MAGNETIC INTENSITY C... C... T TIME C... C... DEL THE LAPLACIAN OPERATOR (A SECOND-ORDER SPATIAL C... DIFFERENTIATION OPERATOR) C... C... MU PERMEABILITY C... C... EPS PERMITTIVITY C... C... SIGMA CONDUCTIVITY C... C... AN ALPHABETIC SUBSCRIPT INDICATES A PARTIAL DERIVATIVE, E.G., C... E AND E ARE THE FIRST AND SECOND-ORDER PARTIAL DERIVATIVES C... T TT C... OF E WITH RESPECT TO T, RESPECTIVELY. C... C... FOR ONE SPATIAL DIMENSION, X, IN CARTESIAN COORDINATES, EQUA- C... TIONS (1) AND (2) BECOME C... C... MU*EPS*E + MU*SIGMA*E = E (3) C... TT T XX C... C... MU*EPS*H + MU*SIGMA*H = H (4) C... TT T XX C... C... CONSIDER EQUATION (3) FOR TWO CASES C... C... (1) EPS = 0 (GOOD CONDUCTOR) C... C... COMPUTE A SOLUTION TO EQUATION (3) FOR THE INITIAL CONDITION C... C... E(X,0) = 0 (5) C... C... AND BOUNDARY CONDITIONS C... C... E(0,T) = 1, E (1,T) = 0 (6)(7) C... X C... C... USE SUBROUTINE DSS042 FOR THIS SOLUTION. PLOT E(X,T) VS X FOR C... A SERIES OF VALUES OF T. CHECK THAT YOUR SOLUTION APPROACHES C... THE ANALYTICAL SOLUTION FOR LARGE T. TAKE MU = SIGMA = 1. C... C... (2) EPS = MU = SIGMA = 1 C... C... COMPUTE A SOLUTION TO EQUATION (3) FOR THE INITIAL CONDITIONS C... C... E(X,0) = SIN(PI*X), E (X,0) = 0 (8)(9) C... T C... C... AND BOUNDARY CONDITIONS C... C... E (0,T) = E (1,T) = 0 (10)(11) C... X X C... C... PLOT E(X,T) VS X FOR A SERIES OF VALUES OF T. C... C... DSS/2 COMMON COMMON/T/ T, NSTOP, NORUN 1 /Y/ U(102) 2 /F/UT(102) C... C... INITIAL CONDITION(S) C... C... CASE 1 IF(NORUN.EQ.1)CALL INTAL1 C... C... CASE 2 IF(NORUN.EQ.2)CALL INTAL2 RETURN END SUBROUTINE DERV COMMON/T/ T, NSTOP, NORUN 1 /Y/ U(102) 2 /F/UT(102) C... C... CASE 1 IF(NORUN.EQ.1)CALL DERV1 C... C... CASE 2 IF(NORUN.EQ.2)CALL DERV2 RETURN END SUBROUTINE PRINT(NI,NO) COMMON/T/ T, NSTOP, NORUN 1 /Y/ U(102) 2 /F/UT(102) C... C... PRINT AND PLOT THE SOLUTION C... C... CASE 1 IF(NORUN.EQ.1)CALL PRINT1(NI,NO) C... C... CASE 2 IF(NORUN.EQ.2)CALL PRINT2(NI,NO) RETURN END SUBROUTINE INTAL1 C... C... SUBROUTINE INTAL1 IMPLEMENTS INITIAL CONDITION (5) C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E(51) 2 /F/ ET(51) COMMON/C/ DXS, EPS, MU, SIGMA, N, 1 IP REAL MU C... C... SET THE PARAMETERS N=51 EPS=0. MU=1. SIGMA=1. C... C... GRID SPACING SQUARED DXS=(1.0/FLOAT(N-1))**2 C... C... INITIAL CONDITION (5) DO 1 I=1,N E(I)=0. 1 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV1 IP=0 RETURN END SUBROUTINE DERV1 C... C... SUBROUTINE DERV1 IMPLEMENTS EQUATION (3) WITH EPS = 0 (CASE 1) C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E(51) 2 /F/ ET(51) COMMON/C/ DXS, EPS, MU, SIGMA, N, 1 IP REAL MU C... C... BOUNDARY CONDITION (6) E(1)=1. ET(1)=0. C... C... EQUATION (3) DO 1 I=2,N C... C... X = 1, INCLUDING BOUNDARY CONDITION (7) IF(I.EQ.N)THEN ET(N)=2.0*(E(N-1)-E(N))/DXS C... C... INTERIOR POINTS (X NE 0, X NE 1) ELSE ET(I)=(E(I+1)-2.0*E(I)+E(I-1))/DXS END IF 1 CONTINUE RETURN END SUBROUTINE PRINT1(NI,NO) C... C... SUBROUTINE PRINT1 PRINTS AND PLOTS THE NUMERICAL SOLUTION FOR C... CASE 1 C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E(51) 2 /F/ ET(51) COMMON/C/ DXS, EPS, MU, SIGMA, N, 1 IP REAL MU C... C... DIMENSION THE ARRAYS FOR PLOTTING COMMON/P1/ XP(51), EP(6,51) C... C... PRINT A HEADING FOR THE SOLUTION IP=IP+1 IF(IP.EQ.1)WRITE(NO,1) 1 FORMAT(9X,'T',/, 1 3X,' E(0,T)',3X,' E(0.2,T)',3X,' E(0.4,T)', 2 3X,' E(0.6,T)',3X,' E(0.8,T)',3X,' E(1.0,T)') C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,2)T,(E(I),I=1,N,10) 2 FORMAT(F10.2,/,6E12.4) C... C... STORE THE SOLUTION FOR PLOTTING DO 3 I=1,N XP(I)=FLOAT(I-1)/FLOAT(N-1) EP(IP,I)=E(I) 3 CONTINUE C... C... PLOT THE SOLUTION IF(IP.LT.6)RETURN CALL TPLOTS(IP,N,XP,EP) WRITE(NO,4) 4 FORMAT(1H ,//,' E(X,T) VS X, T = 0, 0.2,..., 1.0') C... C... ***************************************************************** C... C... THE FOLLOWING CALL TO A CONTINUOUS PLOTTER IS MACHINE DEPENDENT C... AND WILL HAVE TO BE MODIFIED FOR OTHER COMPUTERS OR CONVERTED C... TO A COMMENT C... CALL PLOT1 C... C... ***************************************************************** C... RETURN END SUBROUTINE PLOT1 C... C... SUBROUTINE PLOT1 PRODUCES CONTINUOUS PLOTS OF THE SOLUTION C... FOR CASE 1 C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E(51) 2 /F/ ET(51) COMMON/C/ DXS, EPS, MU, SIGMA, N, 1 IP REAL MU COMMON/P1/ XP(51), EP(6,51) C... C... ****************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND WILL HAVE TO BE CONVERTED FOR OTHER COMPUTERS OR CHANGED TO C... COMMENTS C... C... DRAW THE CONTINUOUS PLOTS C... CALL QIKSAX(3,3) C... CALL QIKPLT(XP,EP,-N,3H*X*,8H*E(X,T)*, C... 1 30H*MAXWELL'S EQUATIONS, EPS = 0*,-6) C... C... ****************************************************************** C... RETURN END SUBROUTINE INTAL2 C... C... SUBROUTINE INTAL2 IMPLEMENTS INITIAL CONDITIONS (8) AND (9). C... THE PROGRAMMING IN SUBROUTINES INTAL2, DERV2 AND PRINT2 IS IN C... TERMS OF STATE VARIABLES E1 AND E2 DEFINED AS C... C... E1 = E, E2 = E C... T C... C... THESE STATE VARIABLES, AND THEIR TIME DERIVATIVES, APPEAR IN C... COMMON/Y/ AND COMMON/F/ RESPECTIVELY. THEIR USE IS REQUIRED C... SO THAT ONLY FIRST-ORDER ODES IN T ARE INTEGRATED. C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E1(51), E2(51) 2 /F/ E1T(51), E2T(51) COMMON/C/ DXS, EPS, MU, SIGMA, N, 1 IP REAL MU C... C... SET THE PARAMETERS N=51 EPS=1. MU=1. SIGMA=1. C... C... GRID SPACING SQUARED DXS=(1.0/FLOAT(N-1))**2 C... C... INITIAL CONDITION (8) PI=3.1415927 DO 1 I=1,N X=FLOAT(I-1)/FLOAT(N-1) E1(I)=SIN(PI*X) 1 CONTINUE C... C... INITIAL CONDITION (9) DO 2 I=1,N E2(I)=0. 2 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV2 IP=0 RETURN END SUBROUTINE DERV2 C... C... SUBROUTINE DERV2 IMPLEMENTS EQUATION (3) WITH EPS = 1 (CASE 2) C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E1(51), E2(51) 2 /F/ E1T(51), E2T(51) COMMON/C/ DXS, EPS, MU, SIGMA, N, 1 IP REAL MU C... C... EQUATION (3) DO 1 I=1,N C... C... X = 0, INCLUDING BOUNDARY CONDITION (10) IF(I.EQ.1)THEN E2T(1)=2.0*(E1(2)-E1(1))/DXS-E2(1) C... C... X = 1, INCLUDING BOUNDARY CONDITION (11) ELSE IF(I.EQ.N)THEN E2T(N)=2.0*(E1(N-1)-E1(N))/DXS-E2(N) C... C... INTERIOR POINTS (X NE 0, X NE 1) ELSE E2T(I)=(E1(I+1)-2.0*E1(I)+E1(I-1))/DXS-E2(I) END IF E1T(I)=E2(I) 1 CONTINUE RETURN END SUBROUTINE PRINT2(NI,NO) C... C... SUBROUTINE PRINT2 PRINTS AND PLOTS THE NUMERICAL SOLUTION FOR C... CASE 2 C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E1(51), E2(51) 2 /F/ E1T(51), E2T(51) COMMON/C/ DXS, EPS, MU, SIGMA, N, 1 IP REAL MU C... C... DIMENSION THE ARRAYS FOR PLOTTING COMMON/P2/ XP(51), EP(6,51) IP=IP+1 C... C... PRINT A HEADING FOR THE SOLUTION IF(IP.EQ.1)WRITE(NO,1) 1 FORMAT(9X,'T',/, 1 3X,' E(0,T)',3X,' E(0.2,T)',3X,' E(0.4,T)', 2 3X,' E(0.6,T)',3X,' E(0.8,T)',3X,' E(1.0,T)') C... C... PRINT THE NUMERICAL SOLUTION WRITE(NO,2)T,(E1(I),I=1,N,10) 2 FORMAT(F10.2,/,6E12.4) C... C... STORE THE SOLUTION FOR PLOTTING DO 3 I=1,N XP(I)=FLOAT(I-1)/FLOAT(N-1) EP(IP,I)=E1(I) 3 CONTINUE C... C... PLOT THE SOLUTION IF(IP.LT.6)RETURN CALL TPLOTS(IP,N,XP,EP) WRITE(NO,4) 4 FORMAT(1H ,//,' E(X,T) VS X, T = 0, 0.2,..., 1.0') C... C... ***************************************************************** C... C... THE FOLLOWING CALL TO A CONTINUOUS PLOTTER IS MACHINE DEPENDENT C... AND WILL HAVE TO BE MODIFIED FOR OTHER COMPUTERS OR CONVERTED C... TO A COMMENT C... CALL PLOT2 C... C... ***************************************************************** C... RETURN END SUBROUTINE PLOT2 C... C... SUBROUTINE PLOT2 PRODUCES CONTINUOUS PLOTS OF THE SOLUTION C... FOR CASE 2 C... COMMON/T/ T, NSTOP, NORUN 1 /Y/ E1(51), E2(51) 2 /F/ E1T(51), E2T(51) COMMON/C/ DXS, EPS, MU, SIGMA, N, 1 IP REAL MU COMMON/P2/ XP(51), EP(6,51) C... C... ****************************************************************** C... C... THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT C... AND WILL HAVE TO BE CONVERTED FOR OTHER COMPUTERS OR CHANGED TO C... COMMENTS C... C... DRAW THE CONTINUOUS PLOTS C... CALL QIKSAX(3,3) C... CALL QIKPLT(XP,EP,-N,3H*X*,8H*E(X,T)*, C... 1 30H*MAXWELL'S EQUATIONS, EPS = 1*,-6) C... C... ****************************************************************** C... RETURN END MAXWELL'S EQUATION, EPS = 0 0. 1.0 0.2 5110000 1 1 ABS 0.01 MAXWELL'S EQUATION, EPS = 1 0. 1.0 0.2 10210000 1 1 ABS 0.01 END OF RUNS