Mercurial > octave-nkf
view libcruft/odessa/odessa_cfode.f @ 5018:1c65a8e44ef9 ss-2-1-59
[project @ 2004-09-22 03:33:29 by jwe]
author | jwe |
---|---|
date | Wed, 22 Sep 2004 03:33:29 +0000 |
parents | 258c1d15ad78 |
children |
line wrap: on
line source
SUBROUTINE ODESSA_CFODE (METH, ELCO, TESCO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION ELCO(13,12), TESCO(3,12) C----------------------------------------------------------------------- C ODESSA_CFODE IS CALLED BY THE INTEGRATOR ROUTINE TO SET COEFFICIENTS C NEEDED THERE. THE COEFFICIENTS FOR THE CURRENT METHOD, AS C GIVEN BY THE VALUE OF METH, ARE SET FOR ALL ORDERS AND SAVED. C THE MAXIMUM ORDER ASSUMED HERE IS 12 IF METH = 1 AND 5 IF METH = 2. C (A SMALLER VALUE OF THE MAXIMUM ORDER IS ALSO ALLOWED.) C ODESSA_CFODE IS CALLED ONCE AT THE BEGINNING OF THE PROBLEM, C AND IS NOT CALLED AGAIN UNLESS AND UNTIL METH IS CHANGED. C C THE ELCO ARRAY CONTAINS THE BASIC METHOD COEFFICIENTS. C THE COEFFICIENTS EL(I), 1 .LE. I .LE. NQ+1, FOR THE METHOD OF C ORDER NQ ARE STORED IN ELCO(I,NQ). THEY ARE GIVEN BY A GENETRATING C POLYNOMIAL, I.E., C L(X) = EL(1) + EL(2)*X + ... + EL(NQ+1)*X**NQ. C FOR THE IMPLICIT ADAMS METHODS, L(X) IS GIVEN BY C DL/DX = (X+1)*(X+2)*...*(X+NQ-1)/FACTORIAL(NQ-1), L(-1) = 0. C FOR THE BDF METHODS, L(X) IS GIVEN BY C L(X) = (X+1)*(X+2)* ... *(X+NQ)/K, C WHERE K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ). C C THE TESCO ARRAY CONTAINS TEST CONSTANTS USED FOR THE C LOCAL ERROR TEST AND THE SELECTION OF STEP SIZE AND/OR ORDER. C AT ORDER NQ, TESCO(K,NQ) IS USED FOR THE SELECTION OF STEP C SIZE AT ORDER NQ - 1 IF K = 1, AT ORDER NQ IF K = 2, AND AT ORDER C NQ + 1 IF K = 3. C----------------------------------------------------------------------- DIMENSION PC(12) PARAMETER (ONE=1.0D0,ZERO=0.0D0) C GO TO (100, 200), METH C 100 ELCO(1,1) = ONE ELCO(2,1) = ONE TESCO(1,1) = ZERO TESCO(2,1) = 2.0D0 TESCO(1,2) = ONE TESCO(3,12) = ZERO PC(1) = ONE RQFAC = ONE DO 140 NQ = 2,12 C----------------------------------------------------------------------- C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL C P(X) = (X+1)*(X+2)*...*(X+NQ-1). C INITIALLY, P(X) = 1. C----------------------------------------------------------------------- RQ1FAC = RQFAC RQFAC = RQFAC/DBLE(NQ) NQM1 = NQ - 1 FNQM1 = DBLE(NQM1) NQP1 = NQ + 1 C FORM COEFFICIENTS OF P(X)*(X+NQ-1). ---------------------------------- PC(NQ) = ZERO DO 110 IB = 1,NQM1 I = NQP1 - IB 110 PC(I) = PC(I-1) + FNQM1*PC(I) PC(1) = FNQM1*PC(1) C COMPUTE INTEGRAL, -1 TO 0, OF P(X) AND X*P(X). ----------------------- PINT = PC(1) XPIN = PC(1)/2.0D0 TSIGN = ONE DO 120 I = 2,NQ TSIGN = -TSIGN PINT = PINT + TSIGN*PC(I)/DBLE(I) 120 XPIN = XPIN + TSIGN*PC(I)/DBLE(I+1) C STORE COEFFICIENTS IN ELCO AND TESCO. -------------------------------- ELCO(1,NQ) = PINT*RQ1FAC ELCO(2,NQ) = ONE DO 130 I = 2,NQ 130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/DBLE(I) AGAMQ = RQFAC*XPIN RAGQ = ONE/AGAMQ TESCO(2,NQ) = RAGQ IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/DBLE(NQP1) TESCO(3,NQM1) = RAGQ 140 CONTINUE RETURN C 200 PC(1) = ONE RQ1FAC = ONE DO 230 NQ = 1,5 C----------------------------------------------------------------------- C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL C P(X) = (X+1)*(X+2)*...*(X+NQ). C INITIALLY, P(X) = 1. C----------------------------------------------------------------------- FNQ = DBLE(NQ) NQP1 = NQ + 1 C FORM COEFFICIENTS OF P(X)*(X+NQ). ------------------------------------ PC(NQP1) = ZERO DO 210 IB = 1,NQ I = NQ + 2 - IB 210 PC(I) = PC(I-1) + FNQ*PC(I) PC(1) = FNQ*PC(1) C STORE COEFFICIENTS IN ELCO AND TESCO. -------------------------------- DO 220 I = 1,NQP1 220 ELCO(I,NQ) = PC(I)/PC(2) ELCO(2,NQ) = ONE TESCO(1,NQ) = RQ1FAC TESCO(2,NQ) = DBLE(NQP1)/ELCO(1,NQ) TESCO(3,NQ) = DBLE(NQ+2)/ELCO(1,NQ) RQ1FAC = RQ1FAC/FNQ 230 CONTINUE RETURN C----------------------- END OF SUBROUTINE ODESSA_CFODE ----------------------- END