SUBROUTINE INVR(A,N,NI) CCC THIS SUBROUTINE CALLS DECOMP ONCE AND SOLVE SEVERAL TIMES CCC TO CONSTRUCT THE INVERSE. IT SOLVES AX=I WHERE EACH MATRIX CCC IS N BY N. CCC INPUT CCC A(N,N) - AN N BY N ARRAY, STORED IN MATRIX WITH DIMENSIONS CCC NI BY NI CCC N - THE SIZE OF THE MATRIX TO BE INVERTED, <=20 CCC NI - THE SIZE OF THE DIMENSION OF A CCC OUTPUT CCC A(N,N) - ON OUTPUT THIS IS THE INVERSE OF THE ORIGINAL A CCC THE ORIGINAL A IS DESTROYED DIMENSION A(NI,NI), B(20), C(400) CCC PACK A DENSELY INTO C DO 5 J=1,N IND = N*(J-1) DO 5 I=1,N C(IND+I) = A(I,J) 5 CONTINUE CCC PERFORM AN LU DECOMPOSITION ON C CALL INVERT (N,1,C,B,1) CCC SOLVE AX=I DO 20 J=1,N DO 10 I=1,N 10 B(I) = 0. B(J) = 1. CALL INVSW (N,1,C,B,1) CCC PUT B INTO INVERSE A DO 15 I=1,N A(I,J) = B(I) 15 CONTINUE 20 CONTINUE RETURN END SUBROUTINE INVERT(N,NE,A,B,ITYPE) DIMENSION A(1000,3), B(N), A1(1000), B1(1000), C1(1000) CCC THIS SUBROUTINE CALLS DECOMP, LUDEC, OR INVTRI TO PREFORM LOWER CCC DECOMPOSITION DEPENDING ON WHETHER THE MATRIX A IS CCC ITYPE A CCC 1 DENSE CCC 2 BLOCK TRIDIAGONAL CCC 3 TRIDIAGONAL CCC IF A IS DENSE THE MATRIX IS STORED AS A( N,N ) CCC IF A IS BLOCK DIAGONAL THE MATRIX IS STORED AS A ( N,N,NE ) CCC IF A IS TRIDIAGONAL THE ELEMENTS A,B,C ARE STORED AS A(N,1) = A(N), CCC A(N,2) = B(N), A(N,3) = C(N) CCC N IS LIMITED TO 201 WITHOUT CHANGING DIMENSION STATEMENT FOR A(N), CCC B(N), C(N) CCC THE MATRICES ARE DIMENSIONED AS A(N,N), A(N,N,NE) OR A(N,3) GO TO (5,10,15), ITYPE 5 CALL DECOMP (N,A) RETURN 10 NP = (N-1)/NE+1 CALL LUDEC (NP,NE,A) RETURN 15 DO 20 K=1,N A1(K) = A(K,1) B1(K) = A(K,2) C1(K) = A(K,3) 20 CONTINUE CALL INVTRI (N,A1,B1,C1) DO 25 K=1,N A(K,1) = A1(K) A(K,2) = B1(K) A(K,3) = C1(K) 25 CONTINUE RETURN ENTRY INVSW CCC THIS PORTION OF SUBROUTINE DOES THE FORWARD AND BACKWARD SWEEP TO CCC SOLVE AX = B, WITH THE X STORED IN B AND INVR MUST BE CALLED FIRST GO TO (30,35,40), ITYPE 30 CALL SOLVE (N,A,B) RETURN 35 CALL FAS (NP,NE,N,A,B) RETURN 40 CALL SWEEP (N,A1,B1,C1,B) RETURN END SUBROUTINE SING(I) GO TO (5,15), I 5 WRITE (6,10) 10 FORMAT (////,' MATRIX WITH ZERO ROW IN DECOMPOSE '////) RETURN 15 WRITE (6,20) 20 FORMAT (////,' SINGULAR MATRIX IN DECOMPOSE. ZERO DIVIDE IN SOLVE $ . '////) RETURN END SUBROUTINE DECOMP(N,A) DIMENSION A(N,N) COMMON /DENSE/ IPS(201),SC(201) CCC PAGE 68, FORSYTH AND MOLER CCC INITIALIZE IPS, A AND SCALES IF (N.EQ.1) RETURN DO 25 I=1,N IPS(I) = I ROWNRM = 0.0 DO 10 J=1,N IF (ROWNRM-ABS(A(I,J))) 5,10,10 5 ROWNRM = ABS(A(I,J)) 10 CONTINUE IF (ROWNRM) 15,20,15 15 SC(I) = 1./ROWNRM GO TO 25 C 20 CALL SING (1) SC(I) = 0. 25 CONTINUE CCC GAUSIAN ELIMINATION WITH PARTIAL PIVOTING NM1 = N-1 DO 65 K=1,NM1 BIG = 0. DO 35 I=K,N IP = IPS(I) SIZE = ABS(A(IP,K)*SC(IP)) IF (SIZE-BIG) 35,35,30 30 BIG = SIZE IDXPIV = I 35 CONTINUE IF (BIG) 45,40,45 40 CALL SING (2) GO TO 65 C 45 IF (IDXPIV-K) 50,55,50 50 J = IPS(K) IPS(K) = IPS(IDXPIV) IPS(IDXPIV) = J 55 KP = IPS(K) PIVOT = A(KP,K) KP1 = K+1 DO 60 I=KP1,N IP = IPS(I) EM = -A(IP,K)/PIVOT A(IP,K) = -EM DO 60 J=KP1,N A(IP,J) = A(IP,J)+EM*A(KP,J) CCC INNER LOOP 60 CONTINUE 65 CONTINUE KP = IPS(N) IF (A(KP,N)) 75,70,75 70 CALL SING (2) 75 CONTINUE RETURN END SUBROUTINE SOLVE(N,A,B) DIMENSION A(N,N), B(N) COMMON /DENSE/ IPS(201),SC(201) CCC PAGE 69, FORSYTH AND MOLER CCC FORWARD SWEEP DENSE MATRIX IF (N.GT.1) GO TO 5 B(1) = B(1)/A(1,1) RETURN 5 CONTINUE NP1 = N+1 IP = IPS(1) SC(1) = B(IP) DO 15 I=2,N IP = IPS(I) IM1 = I-1 SUM = 0. DO 10 J=1,IM1 SUM = SUM+A(IP,J)*SC(J) 10 CONTINUE SC(I) = B(IP)-SUM 15 CONTINUE CCC BACK SUBSTITUTION IP = IPS(N) SC(N) = SC(N)/A(IP,N) DO 25 IBACK=2,N I = NP1-IBACK IP = IPS(I) IP1 = I+1 SUM = 0. DO 20 J=IP1,N SUM = SUM+A(IP,J)*SC(J) 20 CONTINUE SC(I) = (SC(I)-SUM)/A(IP,I) 25 CONTINUE DO 30 I=1,N 30 B(I) = SC(I) RETURN END SUBROUTINE FAS(NP,NE,NT,A,B) DIMENSION A(NP,NP,NE), B(NT) CCC FORWARD SWEEP BLOCK DIAGONAL MATRIX NP1 = NP DO 10 L=1,NE DO 10 I=2,NP1 I2 = (L-1)*(NP-1)+I S = 0. I1 = I-1 DO 5 J=1,I1 J2 = I2-I+J S = S+A(I,J,L)*B(J2) 5 CONTINUE B(I2) = B(I2)-S 10 CONTINUE CCC BACK SUBSTITUTION DO 25 L1=1,NE L = NE-L1+1 IF (L.NE.NE) GO TO 15 B(NT) = B(NT)/A(NP,NP,NE) 15 N1 = NP-1 DO 25 K=1,N1 I = N1+1-K I2 = (L-1)*(NP-1)+I M = I+1 N2 = N1+1 S = 0. DO 20 J=M,N2 J2 = (L-1)*(NP-1)+J 20 S = S+A(I,J,L)*B(J2) B(I2) = (B(I2)-S)/A(I,I,L) 25 CONTINUE RETURN END SUBROUTINE LUDEC(NP,NE,A) DIMENSION A(NP,NP,NE) CCC LOWER DECOMPOSITION BLOCK DIAGONAL MATRIX N1 = NP-1 DO 10 L=1,NE DO 5 K=1,N1 K1 = K+1 DO 5 I=K1,NP S = A(I,K,L)/A(K,K,L) A(I,K,L) = S DO 5 J=K1,NP A(I,J,L) = A(I,J,L)-S*A(K,J,L) 5 CONTINUE IF (L.EQ.NE) RETURN A(1,1,L+1) = A(NP,NP,L) 10 CONTINUE END SUBROUTINE INVTRI(N,A,B,C) DIMENSION A(N), B(N), C(N) CCC LOWER DECOMPOSITION TRIDIAGONAL MATRIX CCC SOLVES A(I-1)*T(I-1)+B(I)*T(I)+C(I+1)*T(I+1) = D(I) DO 5 L=2,N S = A(L)/B(L-1) B(L) = B(L)-S*C(L-1) 5 A(L) = S RETURN END SUBROUTINE SWEEP(N,A,B,C,D) DIMENSION A(N), B(N), C(N), D(N) CCC FORWARD SWEEP TRIDIAGONAL MATRIX DO 5 L=2,N 5 D(L) = D(L)-A(L)*D(L-1) CCC BACK SUBSTITUTION D(N) = D(N)/B(N) DO 10 L=2,N K = N-L+1 10 D(K) = (D(K)-C(K)*D(K+1))/B(K) RETURN END SUBROUTINE COLL(A,B,Q,X,W,ND,N,AA) DIMENSION A(ND,ND), B(ND,ND), Q(ND,ND), X(ND), W(ND) DIMENSION QINV(7,7), Z(7), C(7,7), D(7,7) CCC THIS SUBROUTINE COMPUTES THE MATRICES FOR ORTHOGONAL CCC COLLOCATION USING SYMMETRIC POLYNOMIALS, TABLE 4-5,4-6 CCC CCC INPUT VARIABLES CCC N = NUMBER OF INTERIOR COLLOCATION POINTS CCC ND = ARRAY DIMENSION OF MATRICES IN CALLING PROGRAM CCC AA = GEOMETRY FACTOR CCC = 1 PLANAR CCC = 2 CYLINDRICAL CCC = 3 SPHERICAL CCC CCC OUTPUT VARIABLES CCC A = MATRIX FOR FIRST DERIVATIVE, EQ. 4-205 CCC B = MATRIX FOR LAPLACIAN, EQ. 4-205 CCC Q = MATRIX FOR Q INVERSE, EQ. 4-203 CCC X = VECTOR OF COLLOCATION POINTS, TABLE 4-5 CCC W = VECTOR OF WEGHTS, EQ. 4-207 CCC DIMENSION X1(36), X2(36), X3(36) DATA X1/.4472135955,0.0000000000,0.0000000000,0. $0000000000,0.0000000000,0.0000000000,.3399810436,.8611363116,0. $0000000000,0.0000000000,0.0000000000,0.0000000000,.2386191861,. $6612093865,.9324695142,0.0000000000,0.0000000000,0.0000000000,. $1834346425,.5255324099,.7966664774,.9602898565,0.0000000000,0. $0000000000,.1488743390,.4333953941,.6794095683,.8650633667,. $9739065285,0.0000000000,.1252334085,.3678314990,.5873179543,. $7699026742,.9041172564,.9815606342/ DATA X2/.5773502692,0.0000000000,0.0000000000,0. $0000000000,0.0000000000,0.0000000000,.4597008434,.8880738340,0. $0000000000,0.0000000000,0.0000000000,0.0000000000,.3357106870,. $7071067812,.9419651451,0.0000000000,0.0000000000,0.0000000000,. $2634992300,.5744645143,.8185294874,.9646596062,0.0000000000,0. $0000000000,.2165873427,.4803804169,.7071067812,.8770602346,. $9762632447,0.0000000000,.1837532119,.4115766111,.6170011402,. $7869622564,.9113751660,.9829724091/ DATA X3/.6546536707,0.0000000000,0.0000000000,0. $0000000000,0.0000000000,0.0000000000,.5384693101,.9061798459,0. $0000000000,0.0000000000,0.0000000000,0.0000000000,.4058451514,. $7415311856,.9491079123,0.0000000000,0.0000000000,0.0000000000,. $3242534234,.6133714327,.8360311073,.9681602395,0.0000000000,0. $0000000000,.2695431560,.5190961292,.7301520056,.8870625998,. $9782286581,0.0000000000,.2304583160,.4484927510,.6423493394,. $8015780907,.9175983992,.9841830547/ WRITE (6,5) 5 FORMAT (' THE COLLOCATION POINTS ARE'/) N1 = N+1 DO 25 LL=1,N IA = AA+0.001 LLN = 6*(N-1)+LL GO TO (10,15,20), IA 10 X(LL) = X1(LLN) GO TO 25 C 15 X(LL) = X2(LLN) GO TO 25 C 20 X(LL) = X3(LLN) 25 CONTINUE X(N1) = 1.0 DO 30 I=1,N1 AI = I Z(I) = 1./(2.*AI+AA-2.) 30 CONTINUE DO 35 I=1,N1 Q(I,1) = 1. QINV(I,1) = 1. DO 35 J=2,N1 Q(I,J) = X(I)**(2*J-2) 35 QINV(I,J) = Q(I,J) DO 40 J=1,N1 CA = 2.*J-2. DA = (2.*J-2.)*(2.*J+AA-4.) DO 40 I=1,N1 C(I,J) = CA*X(I)**(2*J-3) D(I,J) = DA*X(I)**(2*J-4) 40 CONTINUE CALL INVR (QINV,N1,7) DO 50 I=1,N1 W(I) = 0.0 DO 50 J=1,N1 A(I,J)=0. B(I,J)=0. DO 45 K=1,N1 A(I,J) = A(I,J)+C(I,K)*QINV(K,J) B(I,J) = B(I,J)+D(I,K)*QINV(K,J) 45 CONTINUE Q(I,J) = QINV(I,J) W(I) = W(I)+Z(J)*QINV(J,I) 50 CONTINUE WRITE (6,55) (X(I),I=1,N1) 55 FORMAT (7E15.5//) RETURN END SUBROUTINE PLANAR(A,B,Q,X,W,N,NX) DIMENSION A(N,N), B(N,N), Q(N,N), X(N), W(N), ZU(7) DIMENSION R(30,30), S(30) DIMENSION POINT(15) CCC THIS SUBROUTINE COMPUTES THE MATRICES FOR COLLOCATION WITHOUT CCC SYMMETRY, USING COLLOCATION POINTS IN TABLE 4-3. CCC INPUT VARIABLES CCC NX = NUMBER OF COLLOCATION POINTS, INCLUDING TWO END POINTS CCC N = ARRAY DIMENSIONS OF MATRICES IN CALLING PROGRAM CCC CCC OUTPUT VARIABLES CCC A = MATRIX FOR FIRST DERIVATIVE, EQ. 4-103 CCC B = MATRIX FOR SECOND DERIVATIVE, EQ. 4-103 CCC Q = MATRIX FOR Q INVERSE, EQ. 4-101 CCC X = VECTOR OF COLLOCATIONS POINTS, FROM TABLE 4-3 CCC W = VECTOR OF WEIGHTS, EQ. 4-106 CCC DATA (POINT(I),I=1,15)/0.,0.577350269189626,0.774596669241483,0. $861136311594053,0.906179845938664,3*0.,0.339981043584856,0. $538469310105683,5*0./ NCOL = NX-2 IF (NX.GT.7) GO TO 35 JRT = (NCOL+1)/2 DO 5 J=1,JRT 5 ZU(J) = POINT((J-1)*5+NCOL) J = 1 DO 10 I=1,JRT X(J) = (1.0-ZU(I))/2. X(J+1) = (1.0+ZU(I))/2. 10 J = J+2 C BUBBLE SORT ON COLLOCATION POINTS NC1 = NCOL-1 DO 20 J=1,NC1 I = J 15 IF (X(I+1).GT.X(I)) GO TO 20 STOR = X(I) X(I) = X(I+1) X(I+1) = STOR IF (I.EQ.1) GO TO 20 I = I-1 GO TO 15 C 20 CONTINUE WRITE (6,25) (I,X(I),I=1,NCOL) 25 FORMAT (3(/),10X,'COLLOCATION PTS',/,5X,'POINT',5X,'ORDINATE',/, $(5X,I5,5X,E15.8)) NC1 = NCOL+1 DO 30 I=2,NC1 K = NC1-I+2 30 X(K) = X(K-1) X(1) = 0.0 X(NX) = 1.0 35 DO 50 I=1,NX R(I,I) = 0.0 A(I,I) = 0.0 S(I) = 1.0 B(I,I) = 0.0 DO 40 J=1,NX IF (I.EQ.J) GO TO 40 R(I,J) = 1.0/(X(I)-X(J)) S(I) = S(I)*R(I,J) 40 CONTINUE DO 45 J=1,NX JX = NX-J+1 IF (JX.LT.J) GO TO 50 IF (JX.EQ.J) A(I,I) = A(I,I)+R(I,J) IF (JX.GT.J) A(I,I) = A(I,I)+R(I,J)+R(I,JX) 45 CONTINUE 50 CONTINUE DO 60 I=1,NX DO 55 J=1,NX IF (I.EQ.J) GO TO 55 A(I,J) = S(J)*R(I,J)/S(I) B(I,J) = 2.0*A(I,J)*(A(I,I)-R(I,J)) B(I,I) = B(I,I)+R(I,J)*(A(I,I)-R(I,J)) 55 CONTINUE 60 CONTINUE DO 85 I=1,NX Q(1,I) = S(I) K = 1 W(I) = 0.0 DO 75 J=1,NX IF (J.EQ.I) GO TO 75 L = K K = K+1 Q(K,I) = Q(L,I) 65 IF (L.EQ.1) GO TO 70 M = L-1 Q(L,I) = Q(M,I)-X(J)*Q(L,I) L = M GO TO 65 C 70 Q(1,I) = -X(J)*Q(1,I) 75 CONTINUE DO 80 J=1,NX 80 W(I) = W(I)+Q(J,I)/FLOAT(J) 85 CONTINUE RETURN END PROGRAM OCRXN DIMENSION A(7,7), B(7,7), Q(7,7), XC(7), W(7), AA(49), D(7), TH(7) COMMON /RXN/ PAR(8) CCC THIS PROGRAM USES ORTHOGONAL COLLOCATION TO SOLVE CCC CCC DEL**2 C = PHI*PHI*R(C) CCC CCC -DC/DR = BIM*(C - 1 ) AT R = 1 CCC CCC VARIABLES CCC CCC N = NUMBER OF INTERIOR COLLOCATION POINTS CCC AS = GEOMETRY FACTOR CCC = 1 PLANAR CCC = 2 CYLINDRICAL CCC = 3 SPHERICAL CCC PHI = THIELE MODULUS CCC BIM = BIOT NUMBER FOR MASS TRANSFER CCC CGUESS = INITIAL GUESS FOR C(R), A CONSTANT CCC INTERACTIVE VERSION OPEN(UNIT=1,FILE='INPUT.DAT',STATUS='OLD') DO 5 KOUNT=1,3 C 10 FORMAT (' ENTER N, A, PHI,BIM,CGUESS ') C read *,N,AS,PHI,BIM,CGUESS C WRITE (6,15) C 15 FORMAT (' ENTER FOUR REACTION RATE PARAMETERS ') C read *,PAR(1),PAR(2),PAR(3),PAR(4) read (1,*)N,AS,PHI,BIM,CGUESS read (1,*)PAR(1) N1 = N+1 N2 = N*N DELTA = PHI*PHI CCC SET INITIAL CONDITION DO 20 I=1,N1 20 TH(I) = CGUESS CCC CALCULATE MATRICES CALL COLL (A,B,Q,XC,W,7,N,AS) CCC BEGIN ITERATION DO 50 ITER=1,20 CCC SET THE MATRICES DO 25 I=1,N2 25 AA(I) = 0. DO 35 J=1,N CALL RXN (TH(J),RATE,DR) D(J) = DELTA*(RATE-DR*TH(J)) D(J) = D(J)-BIM*B(J,N1)/(A(N1,N1)+BIM) DO 30 I=1,N KN = N*(I-1)+J AA(KN) = B(J,I)-B(J,N1)*A(N1,I)/(A(N1,N1)+BIM) IF (I.EQ.J) AA(KN) = AA(KN)-DELTA*DR 30 CONTINUE 35 CONTINUE CCC DO THE LU DECOMPOSITION CALL INVERT (N,1,AA,D,1) CCC SOLVE FOR THE RIGHT HAND SIDE CALL INVSW (N,1,AA,D,1) CCC FIND MAXIMUM CHANGE IN SOLUTION ER = 0. DO 40 I=1,N ERR = ABS(TH(I)-D(I)) TH(I) = D(I) IF (ERR.GT.ER) ER = ERR 40 CONTINUE WRITE (6,45) ITER,ER 45 FORMAT (' ITERATION ',I3,' ERROR IS ',E15.4) IF (ER.LT.1.E-6) GO TO 55 50 CONTINUE CCC CALCULATE THE EFFECTIVENESS FACTOR CCC ETA1 USES EQ. 4-228 55 SUM = 0. DO 60 I=1,N 60 SUM = SUM+A(N1,I)*D(I) D(N1) = (BIM-SUM)/(A(N1,N1)+BIM) SUM = 0. SUM1 = 0. DO 65 I=1,N1 CALL RXN (D(I),RATE,DR) SUM1 = SUM1+W(I) 65 SUM = SUM+W(I)*RATE c=1 call rxn(c,r,dr) ccc divide by r(c=1) to calculate eta; ETA1 = SUM/SUM1/r WRITE (6,70) ETA1 70 FORMAT (' EFF. FACTOR ',F20.15) WRITE (6,80) (D(I),I=1,N1) 80 FORMAT (25(4F10.6,/)) 5 CONTINUE CLOSE(UNIT=1) STOP END SUBROUTINE RXN(C,R,DR) COMMON /RXN/ PAR(8) CCC THIS SUBROUTINE COMPUTES R AND DR/DC, GIVEN C denom=1./(par(1)+c) r=c*denom dr=par(1)*denom*denom RETURN END PROGRAM FDRXN DIMENSION A(1000), B(1000), C(1000), TH(1000), D(1000), AA(1000,3) DIMENSION R(1000) EQUIVALENCE (A(1),AA(1,1)), (B(1),AA(1,2)), (C(1),AA(1,3)) COMMON /RXN/ PAR(8) OPEN(UNIT=1,FILE='INPUT.DAT',STATUS='OLD') DO 5 KOUNT=1,3 read (1,*)N,AS,PHI,BIM,CGUESS read (1,*)PAR(1) N1 = N+1 DELX = 1./N DELTA = PHI*PHI BB = DELTA*DELX**2 CCC SET INITIAL CONDITION DO 20 I=1,N1 R(I) = DELX*(I-1) 20 TH(I) = CGUESS CCC BEGIN ITERATION DO 40 ITER=1,20 CCC SET THE MATRICES C(1) = 2.*AS Dadd = -BIM*DELX*(2.+DELX*(AS-1.)) Badd = Dadd DO 25 I=1,N1 CALL RXN (TH(I),RATE,DR) B(I) = -2.-BB*DR D(I) = BB*(RATE-DR*TH(I)) IF (I.EQ.1) B(1) = B(1)+2.-2.*AS IF (I.EQ.1) GO TO 25 C(I) = (AS-1.)*DELX/(2.*R(I)) A(I) = 1.-C(I) C(I) = 1.+C(I) if (I.LT.N1) goto 25 A(N1) = 2. b(i)=b(i)+badd d(i)=d(i)+dadd 25 CONTINUE CCC DO THE LU DECOMPOSITION CALL INVERT (N1,1,AA,D,3) CCC SOLVE FOR THE RIGHT HAND SIDE CALL INVSW (N1,1,AA,D,3) CCC FIND MAXIMUM CHANGE IN SOLUTION ER = 0. DO 30 I=1,N1 ERR = ABS(TH(I)-D(I)) TH(I) = D(I) IF (ERR.GT.ER) ER = ERR 30 CONTINUE WRITE (6,35) ITER,ER 35 FORMAT (' ITERATION ',I3,' ERROR IS ',E15.4) IF (ER.LT.1.E-6) GO TO 45 40 CONTINUE CCC CALCULATE THE EFFECTIVENESS FACTOR CCC ETA1 USES EQ. 4-144,4-49 45 SL1 = (D(N-1)-4.*D(N)+3.*D(N1))/(DELX*2.) cc=1 call rxn(cc,rr,dr) ccc divide by r(conc.=1) to calculate eta; ETA1 = AS*SL1/DELTA/rr WRITE (6,50) ETA1 50 FORMAT (' EFF. FACTOR ',F20.15) WRITE (6,60) (D(I),I=1,N1) 60 FORMAT (25(4F10.6,/)) 5 CONTINUE CLOSE(UNIT=1) STOP end SUBROUTINE RXN(C,R,DR) COMMON /RXN/ PAR(8) CCC THIS SUBROUTINE COMPUTES R AND DR/DC, GIVEN C denom=1./(par(1)+c) r=c*denom dr=par(1)*denom*denom RETURN END PROGRAM IVRXN COMMON /RXN/ PAR(8) DIMENSION XI(4),STORE(30),ISTORE(5) EXTERNAL RHS CCC READ GEOMETRY AA= 1, PLANAR CCC 2, CYLINDRICAL CCC 3, SPHERICAL CCC NUMBER OF CASES, ANUM CCC ACCURACY = ACCUR CCC REACTION PARAMETERS (8F10.0) CCC ANUM TIMES: PHI,S0 = INITIAL GUESS OF C(X=0) WRITE (5,5) 5 FORMAT(' ENTER AA,NUMBER OF CASES, RELERR,ABSERR') read *, AA,ANUM,ACCUR IF (RELERR.EQ.0.) RELERR = 1.E-6 IF (ABSERR.EQ.0.)ABSERR=1.E-6 WRITE(5,6) 6 FORMAT(' ENTER FOUR PARAMETERS FOR RATE EQUATIONS') read *, (PAR(I),I=1,4) NUM = ANUM DO 35 KK=1,NUM WRITE (5,7) 7 FORMAT(' ENTER PHI AND GUESS FOR CENTER CONCENTRATION ') read *, PHI,S0 PAR(8) = PHI*PHI ITER = 0 S = S0 PAR(7) = AA CCC SET RUNGE-KUTTA PARAMETERS CCC SET ON ENTRY CCC RHS = NAME OF SUBROUTINE CALCULATING RIGHT HAND SIDE CCC NV = NUMBER OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS CCC X = INITIAL "TIME" CCC XP = FINAL "TIME" CCC XI(1-4) = VECTOR OF SOLUTION AT X CCC RELERR = RELATIVE ERROR TOLERANCE CCC ABSERR = ABSOLUTE ERROR TOLERANCE CCC IRR=1 TO START CCC RETURNED CCC XI(1-4) = VECTOR OF SOLUTION AT XP CCC IRR = 0 - NO ERROR CCC 1 - ERROR TEST COULD NOT BE SATISFIED CCC 2 - EXCESSIVE COMPUTATION TIME WILL RESULT 10 CONTINUE NV = 4 X = 0. IRR = 1 XP = 1. IF (S.GE.0.) GO TO 15 S0 = S0/10. S = S0 15 CONTINUE XI(1) = S XI(2) = 0. XI(3) = 1. XI(4) = 0. CCC XI(1) = U=C, XI(2) = V = DC/DX CCC XI(3) = DU/DS, XI(4) = DV/DS CCC SEE EQ. 4-365 ITER = ITER+1 CALL RKF45(RHS,NV,XI,X,XP, 1 RELERR,ABSERR,IRR,STORE,ISTORE) PSI = XI(1)-1. DPSI = XI(3) SN = S-PSI/DPSI IF (ABS(SN-S).LT.RELERR) GO TO 25 WRITE (5,20) S,SN,PSI,DPSI 20 FORMAT (4E15.4) S = SN IF (ITER.GE.100) GO TO 25 GO TO 10 C 25 AE=1. CALL RXN(AE,R,DR) ETA=AA*XI(2)/(PAR(8)*R) WRITE (5,30) S,PHI,XI(2),ETA 30 FORMAT (//,28(1H+),/,4E15.4,//,28(1H+),/) 35 CONTINUE STOP END SUBROUTINE RHS(X,XI,XF) COMMON /RXN/ PAR(8) DIMENSION XI(4), XF(4) CCC THIS SUBROUTINE COMPUTES F(XI) IN DXI/DX=F(XI,X) CCC PAR(8) = PHI**2 PAR(7) = AA = GEOMETRY FACTOR CCC SEE EQ. 4-365 CALL RXN (XI(1),R,DR) DIV = 1. IF (X.EQ.0.) DIV = PAR(7) XF(2) = PAR(8)*R/DIV XF(4) = PAR(8)*DR*XI(3)/DIV XF(1) = XI(2) XF(3) = XI(4) IF (X.EQ.0.) GO TO 5 XF(2) = XF(2)-(PAR(7)-1.)*XI(2)/X XF(4) = XF(4)-(PAR(7)-1.)*XI(4)/X 5 RETURN END SUBROUTINE RXN(C,R,DR) COMMON /RXN/ PAR(8) CCC THIS SUBROUTINE COMPUTES R(C) AND DR/DC, GIVEN C. T = 1.+PAR(1)*(1.-C) E = EXP(PAR(2)*(1.-1./T)) R = C*E DR = E*(1.-PAR(1)*C*PAR(2)/T**2) RETURN END PROGRAM PDE COMMON /SUBFD/ DELX,XX(102) COMMON /SUBOC/ NP,A(7,7),B(7,7),F(7,7),W(7),XC(7) COMMON /GEAR9/ HUSED,NQUSED,NSTEP,NFE,NJE COMMON /DE/ PAR(7),METH CCC CCC CCC THIS PROGRAM SOLVES THE PARTIAL DIFFERENTIAL EQUATION CC CCC DC/DT = F(C) CCC CCC WHERE F(C) IS A DIFFERENTIAL OPERATOR. IT USES FINITE CCC DIFFERENCE OR ORTHOGONAL COLLOCATION, AS CHOSEN BY METH. CCC CCC METH - 1 FINITE DIFFERENCE CCC - 2 ORTHOGONAL COLLOCATION, SYMMETRIC, TABLES 4.5,4.6 CCC - 3 ORTHOGONAL COLLOCATION, UNSYMMETRIC, TABLES 4.3,4.4 CCC CCC DATA INPUT CCC CCC 1. N,METH,AG 2I5,F10.0 CCC N = NUMBER OF INTERIOR FINITE DIFFERENCE POINTS CCC OR NUMBER OF INTERIOR COLLOCATION POINTS CCC (<=6 IF METH = 2 AND <=5 IF METH = 3) CCC METH = METHOD INDICATOR, AS LISTED ABOVE CCC AG = GEOMETRY FACTOR CCC 1-PLANAR, 2-CYLINDRICAL, 3-SPHERICAL CCC 2. PAR(1-7) 7E10.0 CCC PARAMETERS THAT CAN BE USED IN THE DIFFERENTIAL EQUATION CCC 3. MF,EPS CCC MF = METHOD PARAMETER FOR GEARB, USUALLY 22 CCC EPS = ERROR CRITERION USED IN GEARB CCC 4. C(1-N) USER SUPPLIES FORMAT IN INITIAL CCC INITIAL CONDITIONS CCC 5. TF E10.0 CCC TIME YOU WISH TO SAMPLE THE SOLUTION CCC REPEAT THIS DATA INPUT FOR AS MANY TIMES (UP TO 100) CCC THAT YOU WANT THE SOLUTION PRINTED OUT. CCC THE LAST DATA INPUT SHOULD BE ZERO OR A NEGATIVE NUMBER CCC TO STOP. CCC CCC CCC THE USER ALSO SUPPLIES SUBROUTINES AS FOLLOWS CCC CCC SUBROUTINE BC(C,N,METH,C1,CN) CCC INPUT - N,METH,C(1-N) CCC OUTPUT - C1,CN CCC C1=C(X=0), USED IF METH = 1 OR 3 CCC CN=C(X=1), USED FOR ALL METH CCC CCC SUBROUTINE DIFFUN(N,T,C,CDOT) CCC INPUT - N,T,C(1-N) CCC T - TIME CCC C = SOLUTION AT THAT TIME CCC OUTPUT - CDOT(1,N) CCC CDOT = RIGHT HAND SIDE OF DIFFERENTIAL EQUATION, CCC EVALUATED AT T FOR C(1-N) CCC CCC SUBROUTINE INITIAL(C,N) CCC INPUT - N CCC OUTPUT - C(1,N) = INITIAL CONDITIONS CCC THIS IS WHERE DATA INPUT 4. IS USED CCC IF NO DATA IS READ IN INITIAL THEN SKIP DATA INPUT 4. CCC CCC CCC DATA OUTPUT CCC CCC METHOD AND NUMBER OF POINTS CCC COLLOCATION MATRICS IF METH = 2 OR 3 CCC PARAMETERS IN DIFFERENTIAL EQUATIONS CCC MF,EPS CCC TF, HUSED, NQUSED, NSTEP, NFE,NJE CCC TF = TIME AT WHICH THE FOLLOWING SOLUTION APPLIES CCC HUSED = LAST STEP SIZE USED CCC NQUSED = LAST ORDER USED CCC NSTEP = NUMBER OF STEPS TAKEN (TOTAL) CCC NFE = NUMBER OF FUNCTION EVALUATIONS (TOTAL), INCLUDIN CCC THOSE TO EVALUATE THE JACOBIAN CCC NJE = NUMBER OF JACOBIANS EVALUATED (TOTAL) CCC I,C(I), I=1,N, THE SOLUTION CCC FD - C(1) = C(X=DELX) CCC C(N) = C(1.-DELX) CCC OC, SYM - C(1) = C(FIRST INTERIOR COLLOCATION POINT) CCC C(N) = C(LAST INTERIOR COLLOCATIN POINT) CCC C(N+1) = C(X=1.) CCC OC, UNSYM - C(1) = C(0.) CCC C(2) = C(FIRST INTERIOR COLLOCATION POINT) CCC C(N+1) = C(LAST INTERIOR COLLOCATION POINT CCC C(N+2) = C(X=1.0) CCC CCC DIMENSION C(100) READ (5,5) N,METH,AG 5 FORMAT (2I5,F10.0) IF (METH.GT.1) GO TO 20 CCC FINITE DIFFERENCE WRITE (6,10) N 10 FORMAT (' FINITE DIFFERENCE METHOD WITH ',I5,' GRID POINTS') DELX = 1./(N+1.) XX(1) = 0. N1 = N+1 c MU = ML = 1 mu = 1 ml = mu DO 15 I=1,N1 15 XX(I+1) = DELX*I GO TO 90 CCC ORTHOGONAL COLLOCATION 20 WRITE (6,25) N 25 FORMAT (' ORTHOGONAL COLLOCATION WITH ',I3,' INTERIOR',' COLLOCATI $ON POINTS ') c MU = ML = MAX0(1,N/2) mu = max0(1,n/2) ml = mu IF (METH.GE.3) GO TO 35 WRITE (6,30) 30 FORMAT (' SYMMETRIC POLYNOMIALS, TABLES 4.5,4.6 ') NP = N+1 CALL COLL (A,B,F,XC,W,7,N,AG) GO TO 45 C 35 WRITE (6,40) 40 FORMAT (' UNSYMMETRIC POLYNOMIALS, TABLES 4.3,4.4 ') NP = N+2 CALL PLANAR (A,B,F,XC,W,7,NP) 45 WRITE (6,50) 50 FORMAT (//,' A-MATRIX',/) DO 55 I=1,NP 55 WRITE (6,85) (A(I,J),J=1,NP) WRITE (6,60) 60 FORMAT (//,' B-MATRIX',/) DO 65 I=1,NP 65 WRITE (6,85) (B(I,J),J=1,NP) WRITE (6,70) 70 FORMAT (//,' Q-MATRIX',/) DO 75 I=1,NP 75 WRITE (6,85) (F(I,J),J=1,NP) WRITE (6,80) 80 FORMAT (//,' W-MATRIX',/) WRITE (6,85) (W(J),J=1,NP) 85 FORMAT (7E17.8) 90 WRITE (6,95) AG 95 FORMAT (' GEOMETRY FACTOR IS ',F5.0) READ (5,100) (PAR(I),I=1,7) 100 FORMAT (7E10.0) WRITE (6,105) (PAR(I),I=1,7) 105 FORMAT (//,' PARAMETERS IN DIFFERENTIAL EQUATIONS ARE ',/,4E15.7,/ $,3E15.7,//) T = 0. READ (5,115) MF,EPS INDEX = 1 H0 = EPS*EPS HUSED = 0. c NQUSED = NSTEP = NFE = NJE = 0 nqused = 0 nstep = nqused nfe = nqused nje = nqused WRITE (6,110) MF,EPS 110 FORMAT (' GEARB IS USED WITH MF = ',I3,5X,' AND EPS = ',E10.3) 115 FORMAT (I5,E10.0) CALL INITIAL (C,N) CALL PRINT (C,N,T) DO 125 KK=1,100 READ (5,120) TF 120 FORMAT (E10.0) IF (TF.LE.0.) GO TO 130 CALL DRIVEB (N,T,H0,C,TF,EPS,MF,INDEX,MU,ML) CALL PRINT (C,N,TF) 125 CONTINUE 130 STOP END SUBROUTINE PRINT(C,N,TF) COMMON /SUBFD/ DELX,XX(102) COMMON /SUBOC/ NP,A(7,7),B(7,7),F(7,7),W(7),XC(7) COMMON /GEAR9/ HUSED,NQUSED,NSTEP,NFE,NJE COMMON /DE/ PAR(7),METH DIMENSION C(100), CX(102) CALL BC (C,N,METH,C1,CN,TF) IF (METH.EQ.2) GO TO 10 CCC FINITE DIFFERENCE AND OC, UNSYMMTRIC NX = N+2 CX(1) = C1 DO 5 I=1,N 5 CX(I+1) = C(I) GO TO 20 CCC OC,SYM 10 NX = N+1 DO 15 I=1,N 15 CX(I) = C(I) CCC ALL METHODS 20 CX(NX) = CN IF (METH.EQ.1) GO TO 30 DO 25 I=1,NX 25 XX(I) = XC(I) 30 WRITE (6,35) TF 35 FORMAT (//,' TIME = ',E12.5) WRITE (6,40) HUSED,NQUSED,NSTEP,NFE,NJE 40 FORMAT (/,' LAST STEP SIZE USED = ',E10.3,' FOR ',I3,' -TH', $' ORDER METHOD',/,' NUMBER OF STEPS = ',I10,/,' NUMBER OF FUNCTION $ EVALUATIONS = ',I10,/,' NUMBER OF JACOBIAN EVALUATIONS = ',I10,/) DO 50 I=1,NX,2 J = I+1 IF (J.GT.NX) GO TO 45 WRITE (6,55) I,XX(I),CX(I),J,XX(J),CX(J) GO TO 50 C 45 WRITE (6,55) I,XX(I),CX(I) 50 CONTINUE 55 FORMAT (I5,2X,F10.6,F12.9,10X,I5,2X,F10.6,F12.9) RETURN END SUBROUTINE INITIAL(C,N) DIMENSION C(N) DO 5 I=1,N 5 C(I) = 0. RETURN END SUBROUTINE BC(C,N,METH,C1,CN,T) DIMENSION C(N) CCC SET C1=C(X=0.) IF METH = 1 OR 3 CCC SET CN=C(X=1.) FOR ALL METH IF (METH.EQ.2) GO TO 5 C1 = 1. 5 CN = 0. RETURN END SUBROUTINE DIFFUN(N,T,C,CDOT) COMMON /SUBFD/ DELX,XX(102) COMMON /SUBOC/ NP,A(7,7),B(7,7),F(7,7),W(7),XC(7) COMMON /GEAR9/ HUSED,NQUSED,NSTEP,NFE,NJE COMMON /DE/ PAR(7),METH DIMENSION C(N), CDOT(N), TH(10) CCC SOLVE FOR CCC CCC DC/DT = D((1.+PAR(1)*C)DC/DX)/DX CCC CCC IF PAR(2) DIFFERS FROM ZERO ADD TO THE EQUATION CCC CCC /PAR(2) + 2*X*DC/DX CCC CCC CALL BC (C,N,METH,C1,CN,T) GO TO (5,15,30), METH CCC FD 5 DO 10 I=1,N DELX = 1./(N+1.) X = I*DELX CI = C(I) C0 = C1 IF (I.NE.1) C0 = C(I-1) CL = CN IF (I.NE.N) CL = C(I+1) AA = (1.+0.5*PAR(1)*(CL+CI))*(CL-CI) AA = AA-(1.+0.5*PAR(1)*(CI+C0))*(CI-C0) AA = AA/DELX**2 IF (PAR(2).EQ.0.) GO TO 10 AA = AA/PAR(2)+X*(CL-C0)/DELX 10 CDOT(I) = AA GO TO 55 CCC ORTHOGONAL COLLOCATION, SYMMETRIC CCC THIS IS NOT AN APPROPRIATE CHOICE FOR DIFFUSION PROBLEM. CCC SET UP IS PROVIDED FOR DC/DT = DEL**2 C + R 15 NX = N+1 N0 = 0 N1 = 1 R = 1. DO 25 J=1,N SUM1 = 0. DO 20 I=1,NX 20 SUM1 = SUM1+B(J,I)*C(I) 25 CDOT(J) = SUM1+R GO TO 55 CCC ORTHOGONAL COLLOCATION, UNSYMMETRIC 30 NX = N+2 N0 = 1 N1 = 2 DO 35 I=1,N 35 TH(I+N0) = C(I) IF (METH.EQ.3) TH(1) = C1 TH(N+1+N0) = CN DO 50 J=1,N SUM2 = 0. SUM3 = 0. JJ = J+N0 DO 45 K=1,NX SUM2 = SUM2+A(JJ,K)*TH(K) SUM1 = 0. DO 40 I=1,NX 40 SUM1 = SUM1+A(K,I)*TH(I) 45 SUM3 = SUM1*(1.+TH(K)*PAR(1))*A(JJ,K)+SUM3 IF (PAR(2).EQ.0.) GO TO 50 SUM2 = SUM2*2.*XC(JJ) SUM3 = SUM3/PAR(2)+SUM2 50 CDOT(J) = SUM3 55 RETURN END