00001 SUBROUTINE CFODE (METH, ELCO, TESCO) 00002 CLLL. OPTIMIZE 00003 INTEGER METH 00004 INTEGER I, IB, NQ, NQM1, NQP1 00005 DOUBLE PRECISION ELCO, TESCO 00006 DOUBLE PRECISION AGAMQ, FNQ, FNQM1, PC, PINT, RAGQ, 00007 1 RQFAC, RQ1FAC, TSIGN, XPIN 00008 DIMENSION ELCO(13,12), TESCO(3,12) 00009 C----------------------------------------------------------------------- 00010 C CFODE IS CALLED BY THE INTEGRATOR ROUTINE TO SET COEFFICIENTS 00011 C NEEDED THERE. THE COEFFICIENTS FOR THE CURRENT METHOD, AS 00012 C GIVEN BY THE VALUE OF METH, ARE SET FOR ALL ORDERS AND SAVED. 00013 C THE MAXIMUM ORDER ASSUMED HERE IS 12 IF METH = 1 AND 5 IF METH = 2. 00014 C (A SMALLER VALUE OF THE MAXIMUM ORDER IS ALSO ALLOWED.) 00015 C CFODE IS CALLED ONCE AT THE BEGINNING OF THE PROBLEM, 00016 C AND IS NOT CALLED AGAIN UNLESS AND UNTIL METH IS CHANGED. 00017 C 00018 C THE ELCO ARRAY CONTAINS THE BASIC METHOD COEFFICIENTS. 00019 C THE COEFFICIENTS EL(I), 1 .LE. I .LE. NQ+1, FOR THE METHOD OF 00020 C ORDER NQ ARE STORED IN ELCO(I,NQ). THEY ARE GIVEN BY A GENETRATING 00021 C POLYNOMIAL, I.E., 00022 C L(X) = EL(1) + EL(2)*X + ... + EL(NQ+1)*X**NQ. 00023 C FOR THE IMPLICIT ADAMS METHODS, L(X) IS GIVEN BY 00024 C DL/DX = (X+1)*(X+2)*...*(X+NQ-1)/FACTORIAL(NQ-1), L(-1) = 0. 00025 C FOR THE BDF METHODS, L(X) IS GIVEN BY 00026 C L(X) = (X+1)*(X+2)* ... *(X+NQ)/K, 00027 C WHERE K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ). 00028 C 00029 C THE TESCO ARRAY CONTAINS TEST CONSTANTS USED FOR THE 00030 C LOCAL ERROR TEST AND THE SELECTION OF STEP SIZE AND/OR ORDER. 00031 C AT ORDER NQ, TESCO(K,NQ) IS USED FOR THE SELECTION OF STEP 00032 C SIZE AT ORDER NQ - 1 IF K = 1, AT ORDER NQ IF K = 2, AND AT ORDER 00033 C NQ + 1 IF K = 3. 00034 C----------------------------------------------------------------------- 00035 DIMENSION PC(12) 00036 C 00037 GO TO (100, 200), METH 00038 C 00039 100 ELCO(1,1) = 1.0D0 00040 ELCO(2,1) = 1.0D0 00041 TESCO(1,1) = 0.0D0 00042 TESCO(2,1) = 2.0D0 00043 TESCO(1,2) = 1.0D0 00044 TESCO(3,12) = 0.0D0 00045 PC(1) = 1.0D0 00046 RQFAC = 1.0D0 00047 DO 140 NQ = 2,12 00048 C----------------------------------------------------------------------- 00049 C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL 00050 C P(X) = (X+1)*(X+2)*...*(X+NQ-1). 00051 C INITIALLY, P(X) = 1. 00052 C----------------------------------------------------------------------- 00053 RQ1FAC = RQFAC 00054 RQFAC = RQFAC/DBLE(NQ) 00055 NQM1 = NQ - 1 00056 FNQM1 = DBLE(NQM1) 00057 NQP1 = NQ + 1 00058 C FORM COEFFICIENTS OF P(X)*(X+NQ-1). ---------------------------------- 00059 PC(NQ) = 0.0D0 00060 DO 110 IB = 1,NQM1 00061 I = NQP1 - IB 00062 110 PC(I) = PC(I-1) + FNQM1*PC(I) 00063 PC(1) = FNQM1*PC(1) 00064 C COMPUTE INTEGRAL, -1 TO 0, OF P(X) AND X*P(X). ----------------------- 00065 PINT = PC(1) 00066 XPIN = PC(1)/2.0D0 00067 TSIGN = 1.0D0 00068 DO 120 I = 2,NQ 00069 TSIGN = -TSIGN 00070 PINT = PINT + TSIGN*PC(I)/DBLE(I) 00071 120 XPIN = XPIN + TSIGN*PC(I)/DBLE(I+1) 00072 C STORE COEFFICIENTS IN ELCO AND TESCO. -------------------------------- 00073 ELCO(1,NQ) = PINT*RQ1FAC 00074 ELCO(2,NQ) = 1.0D0 00075 DO 130 I = 2,NQ 00076 130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/DBLE(I) 00077 AGAMQ = RQFAC*XPIN 00078 RAGQ = 1.0D0/AGAMQ 00079 TESCO(2,NQ) = RAGQ 00080 IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/DBLE(NQP1) 00081 TESCO(3,NQM1) = RAGQ 00082 140 CONTINUE 00083 RETURN 00084 C 00085 200 PC(1) = 1.0D0 00086 RQ1FAC = 1.0D0 00087 DO 230 NQ = 1,5 00088 C----------------------------------------------------------------------- 00089 C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL 00090 C P(X) = (X+1)*(X+2)*...*(X+NQ). 00091 C INITIALLY, P(X) = 1. 00092 C----------------------------------------------------------------------- 00093 FNQ = DBLE(NQ) 00094 NQP1 = NQ + 1 00095 C FORM COEFFICIENTS OF P(X)*(X+NQ). ------------------------------------ 00096 PC(NQP1) = 0.0D0 00097 DO 210 IB = 1,NQ 00098 I = NQ + 2 - IB 00099 210 PC(I) = PC(I-1) + FNQ*PC(I) 00100 PC(1) = FNQ*PC(1) 00101 C STORE COEFFICIENTS IN ELCO AND TESCO. -------------------------------- 00102 DO 220 I = 1,NQP1 00103 220 ELCO(I,NQ) = PC(I)/PC(2) 00104 ELCO(2,NQ) = 1.0D0 00105 TESCO(1,NQ) = RQ1FAC 00106 TESCO(2,NQ) = DBLE(NQP1)/ELCO(1,NQ) 00107 TESCO(3,NQ) = DBLE(NQ+2)/ELCO(1,NQ) 00108 RQ1FAC = RQ1FAC/FNQ 00109 230 CONTINUE 00110 RETURN 00111 C----------------------- END OF SUBROUTINE CFODE ----------------------- 00112 END
1.7.1
