C ALGORITHM 627 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.11, NO. 1, C MAR., 1985, P. 58. C PROGRAM TSTVE1 MAN 10 C COMPLETE VE1 PACKAGE LISTING FOR DOCUMENT: MAN 20 C MAN 30 C A FORTRAN SUBROUTINE FOR SOLVING VOLTERRA INTEGRAL EQUATIONS MAN 40 C MAN 50 C BY MAN 60 C MAN 70 C L.APPELBAUM AND J.BOWNDS MAN 80 C MAN 90 C REVISED : OCTOBER,1984. MAN 100 C MAN 110 C MAN 120 C MAN 130 C MAN 140 C ************ START OF DRIVER FOR DEMONSTRATION ****************** MAN 150 C MAN 160 C DRIVER FOR USING SUBROUTINE VE1 MAN 170 C MAN 180 C******************************************************************* MAN 190 C *THIS UPDATED PACKAGE OF SUBROUTINES, REFERRED TO GENERICALLY AS * MAN 200 C * "VE1", SUPERCEDES ALL OTHER VERSIONS AS OF 1 SEPT,1984 * MAN 210 C ****************************************************************** MAN 220 C MAN 230 C THIS DRIVER ILLUSTRATES THE USE OF THE SUBPROGRAM PACKAGE MAN 240 C DESCRIBED IN THE ACCOMPANYING DOCUMENT ENTITLED MAN 250 C " A FORTRAN SUBROUTINE FOR SOLVING VOLTERRA INTEGRAL EQUATIONS ", MAN 260 C BY L.APPELBAUM AND J.BOWNDS. MAN 270 C MAN 280 C THE THEORY, TOGETHER WITH A CERTAIN AMOUNT OF PERFORMANCE DATA MAN 290 C ARE FOUND IN THE ARTICLE, " THEORY AND PERFORMANCE OF A SUBROUTINEMAN 300 C FOR SOLVING VOLTERRA INTEGRAL EQUATIONS ", BY J.M. BOWNDS, MAN 310 C IN THE JOURNAL COMPUTING (VOL. 28(1982),PP. 317-332). MAN 320 C MAN 330 C ** (NOTE) ** MAN 340 C MAN 350 C *** IT IS IMPORTANT TO NOTE THAT THE ARTICLE IN "COMPUTING" HAS CER- MAN 360 C TAIN ERRORS WHICH HAVE BEEN CORRECTED AND INCLUDED IN THE ART- MAN 370 C ICLE BY APPELBAUM AND BOWNDS WHICH ACCOMPANIES THIS PACKAGE. MAN 380 C THESE ERRORS ARE SURELY AN ANNOYANCE, HOWEVER, WHEN CORRECTED, MAN 390 C DO NOT CHANGE THE BASIC STRUCTURE OF ANY OF THE UNDERLYING MAN 400 C ALGORITHMS USED IN PREVIOUS VERSIONS OF THIS PACKAGE. *** MAN 410 C MAN 420 C ** THE AUTHORS ARE INDEBTED TO THE VERY COMPETENT REFEREE ON THE MAN 430 C ACM/TOMS ALGORITHMS STAFF FOR THE DETAILED ATTENTION AND MAN 440 C ASSISTANCE GIVEN THIS PACKAGE AND THE OTHER DOCUMENTS RELATED MAN 450 C TO THE BASIC CONCEPT BEHIND THE ACTUAL PROGRAMS. ** MAN 460 C MAN 470 C MAN 480 C FOR MUCH MORE DETAIL REGARDING EFFORT VS. ACCURACY , THE AUTHOR MAN 490 C JOHN M. BOWNDS MAY BE CONTACTED AT THE ADDRESS: MAN 500 C MAN 510 C DEPARTMENT OF MATHEMATICS MAN 520 C UNIVERSITY OF ARIZONA MAN 530 C TUCSON, ARIZONA 85721 MAN 540 C MAN 550 C MAN 560 C THIS DRIVER INCLUDES SIX PASSES WHICH ARE DESIGNED TO DEMONSTRATE MAN 570 C THE MAIN FEATURES OF SUBROUTINE VE1. THE FOLLOWING TABLE SHOWS MAN 580 C MAN 590 C THE PARAMETER VALUES ASSUMED IN THE PASSES: MAN 600 C MAN 610 C PASS KA NEPS M1 N MAN 620 C ---- -- ---- -- - MAN 630 C MAN 640 C 1 0 0 2 2 MAN 650 C 2 1 0 11 2 MAN 660 C 3 1 1 11 2 MAN 670 C 4 1 1 11 21 MAN 680 C 5 1 0 13 2 MAN 690 C 6 1 1 13 2 MAN 700 C MAN 710 C MAN 720 C MAN 730 C ****************************************************************** MAN 740 C MAN 750 C IMPORTANT. THE FOLLOWING COMMENTS RELATE TO THE SPECIFIC POINTS MAN 760 C WHICH MUST BE CONSIDERED BECAUSE OF DIFFERENT COMPILER MAN 770 C REQUIREMENTS. MAN 780 C MAN 790 C THESE COMMENTS WILL BE REFERRED TO BY NUMBER AT MAN 800 C APPROPRIATE POINTS IN SUBPROGRAMS BELOW. MAN 810 C MAN 820 C 1. THE FORTRAN CODING HERE IS INTENDED TO BE AS PORTABLE AS MAN 830 C POSSIBLE. TOWARD THIS END, THIS CODE HAS BEEN CHECKED WITH MAN 840 C PORTABILITY SOFTWARE, AND, EXCEPT WHERE EXPLICITLY POINTED MAN 850 C OUT BELOW, IT IS THOUGHT THAT THIS SUITE OF CODES SHOULD MAN 860 C BE STANDARD. STATEMENTS WHICH NEED MODIFICATION FOR MAN 870 C OTHER INSTALLATIONS ARE SET OFF WITH COMMENTS. MAN 880 C MAN 890 C 2. COMPILER STATIC STATEMENT. MAN 900 C THIS STATEMENT IS REQUIRED WHERE THIS PARTICULAR CODING MAN 910 C WAS PRODUCED. MANY COMPILERS WILL NEITHER REQUIRE NOR MAN 920 C ACCEPT THIS STATEMENT; THE USER SHOULD BE SO ADVISED. MAN 930 C THE COMPILER STATIC STATEMENT OCCURS IN THE FOLLOWING MAN 940 C SUBPROGRAMS - MAN 950 C MAN 960 C SUBROUTINE COEFFS MAN 970 C SUBROUTINE VE1 MAN 980 C FUNCTION F MAN 990 C SUBROUTINE DIFFEQ MAN 1000 C SUBROUTINE ODE MAN 1010 C SUBROUTINE DE MAN 1020 C SUBROUTINE STEP MAN 1030 C SUBROUTINE INTRP MAN 1040 C MAN 1050 C IT IS IMPORTANT THAT THE INCLUSION OF THE MAN 1060 C MAN 1070 C COMPILER STATIC MAN 1080 C MAN 1090 C STATEMENT IS CRUCIAL ON CERTAIN DATA GENERAL ECLIPSE MAN 1100 C GENERATION MACHINES AND MAY ALSO BE IMPORTANT ON LATER MAN 1110 C VERSIONS IN THE MV CLASS OF 32 BIT SUPERMINIS. MAN 1120 C THE APPROACH TAKEN IN THIS PACKAGE IS TO INCLUDE THE MAN 1130 C COMMENTED FORM FOR EASE OF MODIFICATION WHEN USING MAN 1140 C DG OR NON-DG EQUIPMENT. MAN 1150 C MAN 1160 C MAN 1170 C 3. SIZE OF ARRAYS. MAN 1180 C MAN 1190 C THE SIZE OF MANY OF THE ARRAYS USED IN THE DRIVER AND SUB- MAN 1200 C PROGRAMS DEPEND ON THE NUMBER OF TERMS OCCURRING IN EITHER MAN 1210 C THE KERNEL ITSELF OR ITS APPROXIMATION. MAN 1220 C THE SETTING HERE IS SUCH THAT THE CORRESPONDING SYSTEM MAN 1230 C OF DIFFERENTIAL EQUATIONS SHOULD NOT EXCEED 51 EQUATIONS. MAN 1240 C SOME APPLICATIONS, ON SOME MACHINES, COULD EASILY EXCEED MAN 1250 C THIS OF COURSE, AND THE USER MAY WANT TO MAKE APPROPRIATE MAN 1260 C CHANGES. MAN 1270 C MAN 1280 C IT IS IMPORTANT TO NOTE THAT THE SUBROUTINES MAN 1290 C ODE,DE,STEP,AND INTRP MAN 1300 C USE DIMENSION STATEMENTS WHICH DEPEND ON ALLOCATIONS SET MAN 1310 C IN SUBROUTINE DIFFEQ, THE MANAGER OF THESE SUBROUTINES. AS MAN 1320 C SUCH, THESE PROGRAMS USE DIMENSIONS WHICH INCLUDE TERMS LIKE MAN 1330 C Y(NEQN), WHERE NEQN IS PASSED INTO THE SUBROUTINE. TO THE MAN 1340 C AUTHORS' KNOWLEDGE, THIS IS ACCEPTABLE WITH MANY COMPILERS. MAN 1350 C MAN 1360 C ALSO, THE ABOVE PROGAMS USE AN ARRAY CALLED WORK( ), THE MAN 1370 C DIMENSION OF WHICH MUST BE SET IN SUBROUTINE DIFFEQ AS A MAN 1380 C CONSTANT OF SIZE AT LEAST 21*NEQN + 100. MAN 1390 C MAN 1400 C AS IS THE CUSTOM, THIS ARRAY, WORK( ), IS SIMPLY GIVEN THE MAN 1410 C DIMENSION WORK(1) IN SUBROUTINE ODE. THIS AVOIDS UNNECESSARY MAN 1420 C RECOMPILATIONS OF THE ABOVE SUBROUTINES WHEN THE SIZE OF MAN 1430 C WORK( ) NEEDS TO BE CHANGED. MAN 1440 C MAN 1450 C ON CERTAIN COMPILERS, NAMELY SOME OF THOSE WHICH LOOK FOR MAN 1460 C OUT OF BOUND SUBSCRIPTS, THIS HAS CAUSED ERROR MESSAGES AT MAN 1470 C RUNTIME. IN THIS EVENT, THE SIMPLEST SOLUTION SEEMS TO BE MAN 1480 C THAT WORK( ) SHOULD BE DIMENSIONED EXACTLY AS IT IS IN MAN 1490 C THE CALLING PROGRAM DIFFEQ. REFERENCE IS MADE TO THE MAN 1500 C COMMENTS IN ODE,DE,STEP, AND INTRP (L. SHAMPINE,ET. AL.) MAN 1510 C MAN 1520 C ****************************************************************** MAN 1530 C MAN 1540 C USER SUPPLIED INFORMATION: MAN 1550 C MAN 1560 C THE USER MUST SUPPLY FORTRAN CODING FOR - MAN 1570 C MAN 1580 C 1. A CALLING PROGRAM FOR SUBROUTINE VE1( ), MAN 1590 C 2. THE SUBPROGAMS MAN 1600 C FUNCTION K(X) MAN 1610 C SUBROUTINE CTERMS(X) MAN 1620 C SUBROUTINE GTERMS(X) MAN 1630 C FUNCTION FX(X) MAN 1640 C FUNCTION G(X,U) MAN 1650 C MAN 1660 C DEPENDING ON INDIVIDUAL NEEDS, THE USER MAY WANT TO MODIFY - MAN 1670 C MAN 1680 C SUBROUTINE DIFFEQ( ) MAN 1690 C MAN 1700 C AND THE SUITE OF CODES WHICH ARE USED FOR MAN 1710 C SOLVING A SYSTEM OF O.D.E.S. THE CODES USED HERE, MAN 1720 C NAMELY ODE( ), DE( ), STEP( ), AND INTRP( ), MAN 1730 C ARE DUE TO SHAMPINE,ET.AL.,AND ARE INCLUDED HERE MAN 1740 C DUE TO THEIR EASE OF APPLICATION AND MAN 1750 C THEIR PERFORMANCE CHARACTERISTICS. MAN 1760 C (THESE O.D.E. SUBROUTINES ARE USED HERE WITH THE KIND MAN 1770 C PERMISSION OF DR.SHAMPINE.) MAN 1780 C ********************************************************************* MAN 1790 C MAN 1800 EXTERNAL K, G, FX MAN 1810 DOUBLE PRECISION RELERR, ABSERR, U, DIFF, K, G, FX MAN 1820 DOUBLE PRECISION SOLN, H, XR, XL, XOUT, X, E MAN 1830 DOUBLE PRECISION FLN, FLI, FLM MAN 1840 C MAN 1850 C ************************************************************ MAN 1860 C THE SPECIFIC DIMENSION SETTING OF "51" IN THE NEXT STATEMENT IS, MAN 1870 C OF COURSE, ONLY FOR ILLUSTRATIVE PURPOSES HERE; MAN 1880 C CAUTION: OTHER WORKING ARRAYS USED IN THE CURRENT O.D.E. MAN 1890 C SOLVER WILL LIKELY NEED MODIFICATION IF A CHANGE IS MAN 1900 C MADE HERE. MAN 1910 C MAN 1920 C CONSULT OTHER COMMENTS IN THE SUITE OF CODES USED TO SOLVE THE MAN 1930 C DIFFERENTIAL EQUATIONS HERE. MAN 1940 C MAN 1950 C *************************************************************** MAN 1960 DOUBLE PRECISION C(51), GJ(51) MAN 1970 C MAN 1980 C ************************************************************* MAN 1990 INTEGER PASS MAN 2000 C MAN 2010 C MAN 2020 COMMON /VE/ C, GJ, COUNTK, COUNTF MAN 2030 C MAN 2040 C MAN 2050 C THE FOLLOWING STATEMENT IS REQUIRED ON THE COMPUTER WHERE THE CODING MAN 2060 C WAS ORIGINATED (DATA GENERAL ECLIPSE S 230). ITS PURPOSE IS TO MAN 2070 C DEFINE AN OUTPUT FILE FOR OUTPUT TO DEVICE NUMBER "5". MAN 2080 C OF COURSE OTHER INSTALLATIONS MAY HAVE DIFFERENT MAN 2090 C I/O REQUIREMENTS NECESSITATING EITHER A MODIFICATION OR MAN 2100 C ELIMINATION OF THIS PARTICULAR STATEMENT. THE INCLUSION MAN 2110 C HERE IS IN THE FORM OF A COMMENTED STATEMENT. MAN 2120 C MAN 2130 C MAN 2140 C ******************************************************************* MAN 2150 C OPEN 5, "@LIST", ATT="P" MAN 2160 C******************************************************************* MAN 2170 C MAN 2180 C MAN 2190 C THIS DRIVER WILL COUNT FUNCTION EVALUATIONS AND COMPUTE THE ACCURACYMAN 2200 C AT THE RIGHT HAND ENDPOINT FOR COMPARISON PURPOSES. MAN 2210 C MAN 2220 C MAN 2230 C MAN 2240 C PLEASE SEE THE ABOVE REFERENCED DOCUMENT BY APPELBAUM AND BOWNDS FOR MAN 2250 C EXPLANATION OF VALUES FOR THE VARIOUS PARAMETERS. MAN 2260 C MAN 2270 C IN ALL THE EXAMPLES RUN HERE, THE LEFT-HAND ENDPOINT IS 0.0 AND THE MAN 2280 C RIGHT-HAND ENDPOINT IS 2.0; ALSO, THE INPUT LOCAL ERROR MAN 2290 C TOLERANCES ARE SET AS FOLLOWS. MAN 2300 C MAN 2310 XL = 0.0 MAN 2320 XR = 2.0 MAN 2330 ABSERR = 1.0D-08 MAN 2340 RELERR = ABSERR MAN 2350 C MAN 2360 C MAN 2370 DO 130 II=1,6 MAN 2380 PASS = II MAN 2390 WRITE (5,99999) PASS MAN 2400 WRITE (5,99998) MAN 2410 GO TO (10, 20, 30, 40, 50, 60), PASS MAN 2420 C MAN 2430 C MAN 2440 C THE FOLLOWING TABLE GIVES WHICH ROUTINES ARE/ARE NOT CALLED : MAN 2450 C MAN 2460 C KA CTERMS(X) GTERMS(X,U) K(Z) G(X,U)MAN 2470 C -- --------- ------------ ---- ------MAN 2480 C 0 YES YES NO NO MAN 2490 C 1 NO NO YES YES MAN 2500 C MAN 2510 C MAN 2520 C ************************************************************** MAN 2530 C MAN 2540 C THE PARAMETERS KA,NEPS,M1,AND N ARE DEFINED AS FOLLOWS: MAN 2550 C MAN 2560 C MAN 2570 C KA = 0 => USER WILL SUPPLY KERNEL DECOMPOSITION/APPROX. PER MAN 2580 C EQN. (1B) IN ARTICLE IN "COMPUTING" CITED ABOVE; MAN 2590 C MAN 2600 C KA = 1 => USER WILL EXPECT CODE TO GENERATE CHEBYSHEV APPROX. MAN 2610 C FOR CONVOLUTION KERNEL PER EQN. (1A) IN "COMPUTING" ARTICLE; MAN 2620 C MAN 2630 C NEPS = 0 => USER IS NOT REQUESTING A POINTWISE,SIMULTANEOUS ERROR EST-MAN 2640 C IMATE. THE NUMBER OF DIFFERENTIAL EQUATIONS WHICH MUST BE MAN 2650 C SOLVED IS SIMPLY EQUAL TO THE NUMBER OF TERMS USED IN THE MAN 2660 C EXPANSION OF THE KERNEL( EITHER EQN. (1A) OR (1B) IN ABOVE MAN 2670 C REFERENCE). MAN 2680 C MAN 2690 C MAN 2700 C NEPS = 1 => KERNEL MUST BE OF CONVOLUTION TYPE AND AT LEAST TWO MAN 2710 C TERMS MUST BE USED IN THE CHEBYSHEV EXPANSION. MAN 2720 C MAN 2730 C M1 = TOTAL NUMBER OF TERMS TAKEN IN KERNEL EXPANSION; IN THE CASE MAN 2740 C WHERE A CHEBYSHEV EXPANSION IS BEING USED, M1 IS ONE MORE THAN MAN 2750 C THE HIGHEST DEGREE(M) CHEBYSHEV POLYNOMIAL USED IN THE EXPANSION.MAN 2760 C MAN 2770 C N = THE NUMBER OF EQUALLY SPACED POINTS ON THE BASIC INTERVAL WHERE MAN 2780 C A TABULATED SOLUTION IS DESIRED. NOTE: N = 2 => 1 INTERVAL,ETC. MAN 2790 C MAN 2800 C ********************************************************************* MAN 2810 10 KA = 0 MAN 2820 NEPS = 0 MAN 2830 M1 = 2 MAN 2840 N = 2 MAN 2850 GO TO 70 MAN 2860 20 KA = 1 MAN 2870 NEPS = 0 MAN 2880 M1 = 11 MAN 2890 N = 2 MAN 2900 GO TO 70 MAN 2910 30 KA = 1 MAN 2920 NEPS = 1 MAN 2930 M1 = 11 MAN 2940 N = 2 MAN 2950 GO TO 70 MAN 2960 40 KA = 1 MAN 2970 NEPS = 1 MAN 2980 M1 = 11 MAN 2990 N = 21 MAN 3000 GO TO 70 MAN 3010 50 KA = 1 MAN 3020 NEPS = 0 MAN 3030 M1 = 13 MAN 3040 N = 2 MAN 3050 GO TO 70 MAN 3060 60 KA = 1 MAN 3070 NEPS = 1 MAN 3080 M1 = 13 MAN 3090 N = 2 MAN 3100 70 CONTINUE MAN 3110 COUNTK = 0. MAN 3120 COUNTF = 0. MAN 3130 C MAN 3140 M = M1 - 1 MAN 3150 FLM = M MAN 3160 FLN = N MAN 3170 C MAN 3180 H = XR/(FLN-1.0D00) MAN 3190 X = XL MAN 3200 C MAN 3210 C THE FIRST CALL TO VE1 WITH XOUT = X = XL SERVES TO INITIALIZE MAN 3220 C THE O.D.E. SOLVER AND TABULATE THE INITIAL CONDITIONS; IT MAN 3230 C IS RECOMMENDED THAT THE USER DUPLICATE THIS POLICY, SINCE MAN 3240 C IT FOLLOWS THAT FOUND IN MANY CURRENT O.D.E. CODES. MAN 3250 C MAN 3260 XOUT = XL MAN 3270 DO 100 I=1,N MAN 3280 FLI = I MAN 3290 CALL VE1(K, G, FX, M1, KA, NEPS, RELERR, ABSERR, XL, X, XOUT, MAN 3300 * XR, U, E) MAN 3310 DIFF = SOLN(XOUT) - U MAN 3320 IF (NEPS.EQ.0) GO TO 80 MAN 3330 WRITE (5,99997) XOUT, U, DIFF, E MAN 3340 GO TO 90 MAN 3350 80 WRITE (5,99997) XOUT, U, DIFF MAN 3360 90 XOUT = XL + FLI*H MAN 3370 100 CONTINUE MAN 3380 C MAN 3390 C MAN 3400 COUNTF = (FLM+1D0)*COUNTF MAN 3410 TOT = COUNTF + COUNTK MAN 3420 EFFORT = ALOG10(TOT) MAN 3430 IF (DABS(DIFF).LT.1.0D-15) GO TO 110 MAN 3440 ACC = -DLOG10(DABS(DIFF)) MAN 3450 GO TO 120 MAN 3460 110 ACC = 15. MAN 3470 120 RATIO = EFFORT/ACC MAN 3480 WRITE (5,99996) COUNTK MAN 3490 WRITE (5,99995) COUNTF MAN 3500 WRITE (5,99994) ACC MAN 3510 WRITE (5,99993) EFFORT MAN 3520 WRITE (5,99992) RELERR MAN 3530 WRITE (5,99991) ABSERR MAN 3540 WRITE (5,99990) RATIO MAN 3550 130 CONTINUE MAN 3560 C MAN 3570 C THE FOLLOWING STATEMENT IS NOT NECESSARILY NEEDED SINCE IT RELATES MAN 3580 C TO THE "OPEN" STATEMENT EARLIER. AS BEFORE, THE INCLUSION HERE MAN 3590 C IS AS A COMMENTED STATEMENT. MAN 3600 C ******************************************************************* MAN 3610 C CLOSE 5 MAN 3620 C ****************************************************************** MAN 3630 STOP MAN 3640 C MAN 3650 99999 FORMAT (1H1/3X, 13HPASS NUMBER , I3) MAN 3660 99998 FORMAT (1H0///5X, 2H X, 9X, 12HAPPROX. SOLN, 10X, 12HACTUAL ERROR,MAN 3670 * 5X, 14HCOMPUTED ERROR//) MAN 3680 99997 FORMAT (1H /3X, D9.2, 5X, D11.4, 10X, D11.4, 5X, D11.4) MAN 3690 99996 FORMAT (1H , 20HTOTAL KERNEL COUNT , F12.0) MAN 3700 99995 FORMAT (1H , 36HTOTAL O.D.E. RIGHT HAND SIDE COUNT , F12.0) MAN 3710 99994 FORMAT (1H , 10HACCURACY , E13.5) MAN 3720 99993 FORMAT (1H , 14HTOTAL EFFORT , E13.5) MAN 3730 99992 FORMAT (1H , 32HRELATIVE ERROR INPUT TO O.D.E. , D12.4) MAN 3740 99991 FORMAT (1H , 32HABSOLUTE ERROR INPUT TO O.D.E. , D12.4) MAN 3750 99990 FORMAT (1H , 26HEFFORT TO ACCURACY RATIO , F9.3) MAN 3760 C MAN 3770 C MAN 3780 C ************* END OF DRIVER ***************************** MAN 3790 END MAN 3800 C ********** START OF FUNCTION FX *********************************** FX 10 C FX 20 DOUBLE PRECISION FUNCTION FX(X) FX 30 DOUBLE PRECISION X C C PURPOSE : TO SUPPLY THE "DRIVING" TERM ,F(X), IN THE DOCUMENTED C INTEGRAL EQUATION. C C FOR EXAMPLES IN ACM/TOMS: FX = X + 1D0 - DCOS(X) C RETURN C C ************ END OF FUNCTION SUBPROGRAM FX( ) *************** END C ************* START OF FUNCTION K(Z) *************************** K 10 C K 20 DOUBLE PRECISION FUNCTION K(Z) K 30 DOUBLE PRECISION Z C C RE SIZE OF ARRAYS. SEE COMMENT NO. 3 IN DRIVER. C DOUBLE PRECISION C(51), GJ(51) C COMMON /VE/ C, GJ, COUNTK, COUNTF C C PURPOSE : TO EVALUATE THE KERNEL FUNCTION, K(Z) , AT ANY REQUIRED C POINT. THE KERNEL MUST BE OF "CONVOLUTION" TYPE IN THE C ORIGINAL INTEGRAL EQUATION ( TYPE 1(A) IN DOCUMENTATION ). C THIS FUNCTION IS EVALUATED IF AND ONLY THE PARAMETER KA C IS GREATER THAN ZERO. IF KA IS ZERO THEN THE ONLY CODING C REQUIRED IS RETURN. C C FOR EXAMPLES IN ACM/TOMS : C K = -DCOS(Z) COUNTK = COUNTK + 1.0 C RETURN C *********** END OF FUNCTION SUBPROGRAM K( ) ********* END C ************ START OF FUNCTION GFUN ***************************** G 10 C G 20 DOUBLE PRECISION FUNCTION G(X, U) G 30 DOUBLE PRECISION X, U C C PURPOSE : TO EVALUATE THE TERM G(X,U) IN THE INTEGRAL EQUATION C THE TYPE (1A) IN THE DOCUMENTATION. THIS FUNCTION IS C EVALUATED IF AND ONLY IF THE PARAMETER KA IS GREATER THAN C ZERO. C C IF KA IS ZERO THEN THE ONLY CODING REQUIRED IS RETURN. C C FOR EXAMPLES IN TOMS : C G = U C RETURN C ******** END OF FUNCTION SUBPROGRAM G( ) *********** END C ************* START OF SUBROUTINE CTERMS *********************** CTE 10 C CTE 20 SUBROUTINE CTERMS(X) CTE 30 C C RE SIZE OF ARRAYS; SEE COMMENT NO. 3 IN DRIVER. C DOUBLE PRECISION C(51), GJ(51), X C COMMON /VE/ C, GJ, COUNTK, COUNTF C C C C PURPOSE : TO SUPPLY C(J,X) TERMS FOR USER SUPPLIED DECOMPOSITION. C THE USE OF THIS SUBROUTINE ALSO IMPLIES THE EXISTENCE AND C INCLUSION OF SUBROUTINE GTERMS(X,U). THESE TWO SUBROUTINES C ARE ACTUALLY CALLED IF AND ONLY IF THE PARAMETER KA IS ZERO. C C IF KA IS GREATER THAN ZERO, THEN THE ONLY CODING NEEDED BY C THE USER IS RETURN. C C FOR EXAMPLES IN TOMS: C(1) = -DCOS(X) C(2) = -DSIN(X) C RETURN C ******** END OF SUBROUTINE CTERMS ***************** END C *********** START OF SUBROUTINE GTERMS ************************** GTE 10 C GTE 20 SUBROUTINE GTERMS(X, U) GTE 30 COMMON /VE/ C, GJ, COUNTK, COUNTF C C RE SIZE OF ARRAYS. SEE COMMENT NO. 3 IN DRIVER. C DOUBLE PRECISION X, C(51), GJ(51), U C C C C PURPOSE : TO SUPPLY THE G(J,X,U) TERMS IN A USER SUPPLIED C DECOMPOSITION. THE INCLUSION OF THIS SUBROUTINE IMPLIES C THE EXISTENCE AND INCLUSION OF SUBROUTINE CTERMS(X). THIS C SUBROUTINE IS CALLED IF AND ONLY IF THE PARAMETER KA IS ZERO. C C IF KA IS GREATER THAN ZERO, THE ONLY CODING REQUIRED IS C RETURN. C C FOR EXAMPLES IN TOMS : GJ(1) = (DCOS(X))*U GJ(2) = (DSIN(X))*U RETURN C ******** END OF SUBROUTINE GTERMS ******************+ END C *************** START OF SUBROUTINE VE1 **************************** VE1 10 C VE1 20 SUBROUTINE VE1(K, G, FX, M1, KA, NEPS, RELERR, ABSERR, XL, X, VE1 30 * XOUT, XR, U, E) C C PURPOSE : TO STEP OUT THE SOLUTION TO A GIVEN VOLTERRA INTEGRAL C EQUATION OF THE SECOND KIND BY THE METHODS DISCUSSED IN C THE ASSOCIATED DOCUMENTATION (SEE REFERENCES IN THE DRIVER C SUPPLIED WITH THIS AND OTHER RELEVANT SUBPROGRAMS IN THIS C PACKAGE). AN ATTEMPT IS MADE TO SOLVE FROM X TO XOUT. UPON C RETURNING, VE1 SAVES WHAT IS NECESSARY SO THAT THE USER NEED C ONLY MODIFY XOUT TO PROCEED TO THE NEXT OUTPUT POINT. IT C SHOULD BE NOTED THAT THE SUBROUTINE MAY OR MAY NOT ACTUALLY C COMPUTE SOLUTION VALUES AT OTHER POINTS DUE TO THE USE OF C VARIABLE STEP O.D.E. SOLVERS; HOWEVER, VE1 WILL ATTEMPT TO C PROVIDE A SOLUTION VALUE AT THE USER DEFINED XOUT,IF C POSSIBLE. C THE SUCCESS OR FAILURE OF THE ROUTINE TO REACH XOUT DEPENDS C UPON THE SUCCESS OR FAILURE OF THE O.D.E. SOLVER. THIS COND- C ITION IS HANDLED IN THE SUPPLIED SUBROUTINE DIFFEQ AND IS C LEFT TO THE USER TO DECIDE WHAT ACTION IS TO BE TAKEN. C C C RE THE FOLLOWING COMPILER STATEMENT, C SEE COMMENT NO.2 IN DRIVER. C C COMPILER STATIC C C C THESE SUBROUTINES USE DEVICE NUMBER 5 TO PRINT MESSAGES. THIS MAY C NEED TO BE ALTERED FOR THE SPECIFIC INSTALLATION. C C THE PARAMETERS REPRESENT : C C K SUBROUTINE K(X) TO EVALUATE THE KERNEL AT X. C G FUNCTION G(X,U) TO EVALUATE THE TERMS INVOLVING U IN THE C INTEGRAND. C FX FUNCTION F(X) TO EVALUATE THE F(X) TERM OF EQUATION. C M1 NUMBER OF TERMS IN THE KERNEL DECOMPOSITION. IN THIS C VERSION OF THE PROGRAM, THE MAXIMUM FOR M1 IS 51 IF C THERE IS NO ERROR ESTIMATE, AND 25 IF THERE IS AN ERROR C ESTIMATE. IF M1 IS TO BE GREATER THAN THE GIVEN LIMITS, C THE DIMENSIONS OF C AND D IN THE COMMON STATEMENT AND C OF THE VECTORS IN THE ODE SOLVER SHOULD BE INCREASED C ACCORDINGLY. C KA FLAG, KA = 0 IF EXACT DECOMPOSITION IS USED (EQN 1(B); C KA = 1 IF KERNEL APPROXIMATION IS USED (EQN 1(A)). C NEPS FLAG, NEPS = 0 IF NO SIMULTANEOUS TRUNCATION ERROR C ESTIMATE IS DESIRED. C NEPS = 1 IF SIMULTANEOUS TRUNCATION ERROR C ESTIMATE IS DESIRED. C RELERR RELATIVE ERROR TOLERANCE FOR O.D.E. SOLVER. C ABSERR ABSOLUTE ERROR TOLERANCE FOR O.D.E. SOLVER. C X INDEPENDENT VARIABLE. FOR EACH CALL, X BEGINS C AT THE PREVIOUS OUTPUT POINT AND IS INCREMENTED C TO THE NEW OUPUT POINT IN THE ODE SOLVER. C XOUT OUTPUT POINT AT WHICH SOLUTION IS DESIRED. C CONTROL IS RETURNED TO THE CALLING ROUTINE C WHEN X=XOUT. ON THE FIRST CALL, XOUT=XL INITIALIZES C THE SUBROUTINE;THIS PROCEDURE IS COMMONLY USED WITH C POPULAR O.D.E. SOLVERS AND SO IS ALSO USED HERE. C XL LEFT ENDPOINT OF THE INTERVAL ON WHICH THE C SOLUTION IS DESIRED. REMAINS FIXED THROUHOUT THE C EXECUTION OF THE PROGRAM. C XR RIGHT ENDPOINT OF THE INTERVAL ON WHICH THE C SOLUTION IS DESIRED. REMAINS FIXED THROUGHOUT THE C EXECUTION OF THE PROGRAM. C U COMPUTED SOLUTION AT XOUT. C E ERROR ESTIMATE. IF NEPS = 0 THEN THE VALUE OF E C IS UNPREDICTABLE AND SHOULD BE IGNORED. C C THE USER MUST DECLARE K, G, AND FX IN AN EXTERNAL STATEMENT, C SUPPLY THE RESPECTIVE ROUTINES AND INITIALIZE THE FOLLOWING C PARAMETERS UPON INITIAL ENTRY INTO VE1 : C C M1 NUMBER OF TERMS IN KERNEL DECOMPOSITION. C KA APPROXIMATION METHOD (IF ANY) TO BE USED; C SEE ABOVE FOR DETAILS. C NEPS FLAG FOR SIMULTANEOUS TRUNCATION ERROR CALCULATION. C RELERR SEE DOCUMENTATION OF SPECIFIC O.D.E. SOLVER. C ABSERR SEE DOCUMENTATION OF SPECIFIC O.D.E. SOLVER. C X STARTING POINT. MUST BE A VARIABLE. C XOUT OUTPUT POINT. XOUT = X IS NECESSARY ON THE FIRST CALL C BUT SHOULD NEVER OCCUR AGAIN. C XL LEFT ENDPOINT OF INTERVAL OF SOLUTION. C XR RIGHT ENDPOINT OF INTERVAL OF SOLUTION. C C SUBPROGRAMS CALLED : C C C C COEFFS SUBROUTINE SUPPLIED WITH VE1; CALCULATES THE CHEBYSHEV C COEFFICIENTS IF THE EXPANSION OPTION IS CHOSEN (KA = 1). C C DIFFEQ SUBROUTINE SUPPLIED WITH VE1; CALLS THE ODE SOLVER AND C CHECKS THE FLAGS FOR TROUBLE. THIS SUBROUTINE MAY NEED C TO BE MODIFIED IF THE USER SUPPLIES HIS OWN ODE SOLVER. C CTERMS SUBROUTINE WRITTEN BY USER IF EXACT DECOMPOSITION OPTION C CHOSEN (KA = 0). C SUM SUBROUTINE SUPPLIED WITH VE1; SUMS THE POSITIVE TERMS C AND NEGATIVE TERMS OF THE KERNEL DECOMPOSITION. C C EXTERNAL F, K, G, FX DOUBLE PRECISION K, G, FX, RELERR, ABSERR, X, XOUT, XR, XL, BMA, * U, E C C RE SIZE OF ARRAYS. SEE COMMENT NO.3 IN DRIVER. C DOUBLE PRECISION Z(51), POS, NEG, TERM, TERM1, TERM2, C(51), * GJ(51) C REAL COUNTK, COUNTF INTEGER IFLAG, KA, M, M1, M1P1, NEPS, NEQN C C NOTE THAT THE COUNTING VARIABLES COUNTK AND COUNTF ARE NOT C DOUBLE PRECISION. C COMMON /VE/ C, GJ, COUNTK, COUNTF C C C C C INITIALIZATION - THE CODE DOWN TO STATEMENT 10 IS EXECUTED ONLY C ON THE FIRST CALL TO VE1. THIS PARTICULAR PROCEDURE DUPLICATES C THOSE OF SEVERAL COMMONLY USED O.D.E. SOLVERS, AND IT PROVIDES C A CONVENIENT WAY TO LIST INITIAL DATA AND MANAGE SUCH SOLVERS. C IF (XL.NE.XOUT) GO TO 50 C C IF THE USER SUPPLIES THE DECOMPOSITION OF THE KERNEL(KA=0),THEN C THERE CAN BE NO ERROR ESTIMATE OF THE FORM GIVEN IN THE C ASSOCIATED DOCUMENT. IF NEPS HAD BEEN SET TO 1 FOR THIS C CASE, IT IS SET TO 0 AND A MESSAGE IS PRINTED TO THIS EFFECT. C IF (KA.EQ.0 .AND. NEPS.EQ.1) WRITE (5,99999) IF (KA.EQ.0) NEPS = 0 C C THE FOLLOWING COMMENTS AND LOGIC ARE KINDLY DUE TO A CAREFUL REFEREE. C C IF AN ERROR ESTIMATE IS DESIRED, IT IS ASSUMED IN THIS PACKAGE THAT C THERE ARE AT LEAST TWO TERMS IN THE CHEBYSHEV EXPANSION ASSOCIATED C WITH THE CONVOLUTION KERNEL; HENCE,THE QUANTITY M1 MUST BE 2 OR MORE. C IF (NEPS.NE.1 .OR. M1.GE.2) GO TO 10 WRITE (5,99998) M1 = 2 10 CONTINUE C C (THANKS TO REFEREE.) X = XL M1P1 = M1 + 1 BMA = XR - XL M = M1 - 1 IFLAG = 1 C C NEQN IS THE NUMBER OF DIFFERENTIAL EQUATIONS TO BE SOLVED. C WITHOUT AN ERROR ESTIMATE IT IS THE SAME AS THE NUMBER OF C TERMS IN THE EXPANSION OF THE KERNEL; WITH THE ERROR ESTI- C MATE IT IS TWICE THAT PLUS ONE. C NEQN = M1 IF (NEPS.EQ.1) NEQN = 2*NEQN + 1 C C C C C SET INITIAL CONDITIONS FOR SYSTEM OF DIFFERENTIAL EQUATIONS. C C DO 20 J=1,NEQN Z(J) = 0.0D00 20 CONTINUE C C COMPUTE THE COEFFICIENTS TO THE CHEBYSHEV SERIES FOR THE KA=1 MODE. C C (THANKS TO REFEREE FOR POINTING OUT THE NEED FOR THE FOLLOWING.) C C IF KA = 0, THEN THERE IS NO NEED TO COMPUTE CHEBYSHEV COEFFICIENTS; C IF KA = 1, THEN - C C IF NEPS = 0, THERE IS NEED FOR M1 (= M + 1) COEFFICIENTS; C IF NEPS = 1, THERE IS NEED FOR M2 = M + 2 COEFFICIENTS. C C THE ABOVE FACTS WERE OVERLOOKED IN ALL PREVIOUS C VERSIONS OF THIS PACKAGE. IF M + 2 COEFFS. C ARE NOT STORED IN THE ARRAY C( ) WHEN NEPS=1, C THEN THE O.D.E.S WILL BE ILL-SET. C C C IF (KA.EQ.0) GO TO 40 IF (NEPS.EQ.0) GO TO 30 M2 = M1 + 1 CALL COEFFS(K, M2, BMA) GO TO 40 30 CALL COEFFS(K, M1, BMA) 40 U = FX(XL) CALL DIFFEQ(F, G, FX, NEPS, NEQN, KA, BMA, Z, X, XOUT, RELERR, * ABSERR, IFLAG) GO TO 140 C 50 CONTINUE C C COMPUTE THE SOLUTION U BY SOLVING THE SYSTEM OF DIFFERENTIAL C EQUATIONS ASSOCIATED WITH THE INTEGRAL EQUATION AS DISCUSSED C IN THE PAPER. TO AIDE IN SUMMING THE ALTERNATING SERIES, THE C TWO PARTIAL SUMS ARE KEPT IN THE TWO VARIABLES POS AND NEG; POS C FOR THE SUM OF THE POSITIVE TERMS AND NEG FOR THE NEGATIVE TERMS. C THE ACTUAL CHECKING OF THE SIGN AND THE SUMMING ARE DONE IN THE C SUBROUTINE SUM. C POS = 0.0D00 NEG = 0.0D00 CALL DIFFEQ(F, G, FX, NEPS, NEQN, KA, BMA, Z, X, XOUT, RELERR, * ABSERR, IFLAG) C IF (M.NE.0 .OR. KA.EQ.0) GO TO 60 U = C(1)*Z(1)/2.0D00 + FX(X) GO TO 140 C 60 IF (KA.EQ.0) CALL CTERMS(X) TERM = C(1)*Z(1) IF (KA.EQ.1) TERM = TERM + C(2)*Z(2) CALL SUM(TERM, POS, NEG) IF (M.EQ.1 .AND. KA.EQ.1) GO TO 100 IF (M.EQ.0 .AND. KA.EQ.0) GO TO 110 IF (KA.EQ.1) GO TO 80 DO 70 J=2,M1 TERM1 = C(J)*Z(J) CALL SUM(TERM1, POS, NEG) 70 CONTINUE GO TO 110 C C C REFEREE- NOTE CHANGE; THANK YOU. 80 DO 90 J=3,M1 TERM1 = C(J)*Z(J) TERM2 = -C(J)*Z(J-2) CALL SUM(TERM1, POS, NEG) CALL SUM(TERM2, POS, NEG) C 90 CONTINUE 100 POS = POS/2.0D00 NEG = NEG/2.0D00 C 110 U = POS + NEG + FX(X) C IF (NEPS.EQ.0) GO TO 140 C C IF THE OPTION FOR A SIMULTANEOUS ERROR ESTIMATE IS CHOSEN (NEPS = 1): C COMPUTE THE ERROR ESTIMATE BY SOLVING THE INTEGRAL EQUATION FOR C EPSILON IN A SIMILAR WAY AS THE SOLUTION U. IT IS ASSUMED THAT C THERE ARE AT LEAST TWO TERMS IN THE CHEBYCHEV EXPANSION OF THE C KERNEL. THIS PROCEDURE IS EXPLAINED IN THE DOCUMENT REFERRED TO C ABOVE. C C C REFEREE COMMENTS. C Z(M+2) IS THE SOLUTION OF AN EXTRA DIFFERENTIAL EQUATION REQUIRED C POINTWISE, SIMULTANEOUS ERROR ESTIMATE. THE OTHER COMPONENTS OF Z C ARE AS FOLLOWS (OP.CIT. BOWNDS WITH CORRECTIONS.) C C Z(M+3) = ZETA(0) ( IN ABOVE CITED ARTICLE IN "COMPUTING") C Z(M+4) = ZETA(1), C Z(M+5) = ZETA(2), C * = * , C * = * , C Z(2M+3)= ZETA(M). C POS = 0.0D00 NEG = 0.0D00 TERM = C(1)*Z(M1+2) + C(2)*Z(M1+3) CALL SUM(TERM, POS, NEG) IF (M1.EQ.2) GO TO 130 DO 120 J=3,M1 JJ = M1 + 1 + J TERM1 = C(J)*Z(JJ) TERM2 = -C(J)*Z(JJ-2) CALL SUM(TERM1, POS, NEG) CALL SUM(TERM2, POS, NEG) 120 CONTINUE POS = POS/2.0D00 NEG = NEG/2.0D00 C 130 E = POS + NEG + C(M1+1)*(Z(M1+1)-Z(M))/2.0D00 C 140 RETURN C 99999 FORMAT (50H NEPS RESET TO 0. THERE WILL BE NO ERROR ESTIMATE.) 99998 FORMAT (50H M1 SET TO 2. (MIN.# TERMS IN CASE OF ERROR EST. )) C C C ****************** END OF SUBROUTINE VE1 ************************ END C ************* START OF SUBROUTINE F ******************************** F 10 C F 20 SUBROUTINE F(G, FX, NEPS, NEQN, KA, BMA, X, Z, ZP) F 30 C C PURPOSE : TO COMPUTE THE RIGHT HAND SIDE OF THE APPROPRIATE C SYSTEM OF O.D.E.S PURSUANT TO THE BASIC CONVERSION METHOD C OUTLINED LINED IN THE DOCUMENTS REFERENCED IN THE SUPPLIED C DRIVER. C C C RE THE FOLLOWING COMPILER STATEMENT, C SEE COMMENT NO.2 IN DRIVER. C C COMPILER STATIC C C C THIS SUBROUTINE IS CALLED BY THE O.D.E. SOLVER TO COMPUTE C THE RIGHT HAND SIDE OF THE SYSTEM OF O.D.E.'S. C C THE PARAMETERS REPRESENT : C C G FUNCTION G(X,U) WHICH EVALUATES ALL TERMS DEPENDENT ON (X,U) IN C THE INTEGRAND IF KA = 1. C FX FUNCTION F(X) EVALUATING THE DRIVING TERM. C NEPS FLAG, NEPS = 0 FOR NO ERROR ESTIMATE, C NEPS = 1 FOR ERROR ESTIMATE EVALUATION. C NEQN NUMBER OF EQUATIONS FOR ODE INTEGRATOR. C THE VALUE FOR THIS PARAMETER DEPENDS ON THE QUANTITIES C KA AND NEPS IN THE FOLLOWING WAY- C C IF KA = 0, THEN NEQN = M1 = NO. OF TERMS IN KERNEL C DECOMPOSITION/APPROXIMATION; C C IF KA = 1 AND NEPS = 0, THEN NEQN = M1 = NO. OF TERMS C TAKEN IN CHEBYSHEV EXPANSION FOR THE KERNEL; C C IF KA = 1 AND NEPS = 1, THEN NEQN IS CHANGED BY THE C CODE AND IS RECOMPUTED AS C C NEQN = 2*M1 + 1, C C WHERE M1 IS THE NO. OF TERMS TAKEN IN THE C CHEBYSHEV EXPANSION FOR THE KERNEL, WITH M1 C ASSUMED TO BE AT LEAST 2. C C KA FLAG, KA = 0 IF USER SUPPLIES KERNEL DECOMPOSITION, C KA = 1 IF CHEBYSHEV APPROXIMATION IS USED. C C BMA INTERVAL WIDTH. C C X POINT AT WHICH THE RIGHT-HAND SIDE IS BEING EVALUATED. C C Z SOLUTION VECTOR OF ODE AT X. IF NEPS = 1, THEN THE C VALUE OF NEQN IS RECOMPUTED AS IN THE ABOVE COMMENTS, C THE FIRST M1 VALUES OF Z REPRESENT THE SOLUTION TO THE C ODE CORRESPONDING TO THE INTEGRAL EQUATION FOR U(X). C Z(M1+1) IS THE SOLUTION TO AN EXTRA EQUATION WHICH IS C REQUIRED FOR THE ERROR ESTIMATION, AND Z(M1+2),..., C ...,Z(2*M1+1) REPRESENT THE SOLUTION TO THE ODE C CORRESPONDING TO THE INTEGRAL EQUATION FOR EPSILON C (THE ERROR ESTIMATE). NOTE THAT NEQN = 2*M1 + 1. C C ZP DERIVATIVES OF THE SOLUTION VECTOR AT X. IF NEPS = 1, C THEN THE ZP ARRAY IS DIVIDED IN THE SAME PROPORTIONS C AS THE Z ARRAY DESCRIBED ABOVE. C C SUBPROGRAMS CALLED : C C G FUNCTION G(X,U). C C FX FUNCTION FOR DRIVING TERM. C C SUM SUBROUTINE FOR SUMMING THE POSITIVE AND NEGATIVE OF C THE KERNEL DECOMPOSITION. C C CTERMS SUBROUTINE SUPPLIED BY THE USER TO EVALUATE THE C(J) C TERMS OF THE USER-SUPPLIED KERNEL DECOMPOSITION. C C GTERMS SUBROUTINE SUPPLIED BY THE USER TO EVALUATE THE GJ(X,U) C TERMS OF THE EXACT DECOMPOSITION. C C THE CODE THROUGH THE FIRST RETURN STATEMENT EVALUATES THE RIGHT- C HAND SIDE OF THE SYSTEM OF ODE'S IF THE USER HAS SUPPLIED HIS C OWN KERNEL DECOMPOSITION. C C C RE SIZE OF ARRAYS. SEE COMMENT NO.3 IN DRIVER. C DOUBLE PRECISION G, FX, BMA, X, Z(51), ZP(51), BMA4, FM1 DOUBLE PRECISION TERM, TERM1, TERM2, POS, NEG, U, ST, RI, E DOUBLE PRECISION C(51), GJ(51) C COMMON /VE/ C, GJ, COUNTK, COUNTF C C NOTE THAT THE COUNTING VARIABLES COUNTK AND COUNTF ARE NOT C DOUBLE PRECISION. C IF (KA.EQ.1) GO TO 60 C CALL CTERMS(X) TERM = C(1)*Z(1) IF (TERM.LE.0.0D00) GO TO 10 POS = TERM NEG = 0.0D00 GO TO 20 10 NEG = TERM POS = 0.0D00 20 IF (NEQN.EQ.1) GO TO 40 DO 30 I=2,NEQN TERM = C(I)*Z(I) CALL SUM(TERM, POS, NEG) 30 CONTINUE C 40 U = POS + NEG + FX(X) CALL GTERMS(X, U) C C DO 50 I=1,NEQN ZP(I) = GJ(I) 50 CONTINUE C COUNTF = COUNTF + 1. RETURN C C THE CODE FOM HERE TO THE END EVALUATES THE RIGHT-HAND SIDE OF THE C SYSTEM OF O.D.E.'S IF THE CHEBYSHEV APPROXIMATION IS IN EFFECT C (KA = 1). AS MENTIONED SEVERAL PLACES ELSEWHERE, THE NUMBER OF C EQUATIONS IN THIS SYSTEM IS NEQN; HOWEVER, THE ACTUAL VALUE C DEPENDS ON THE ERROR ESTIMATION SWITCH (NEPS = 0 OR 1). C C IF NEPS = 0 THEN NEQN = M1 = NO. OF TERMS IN EXPANSION; C IF NEPS = 1 THEN NEQN = 2*M1 + 1, WHERE M1 IS THE NUMBER C OF TERMS IN THE EXPANSION. C C 60 M1 = NEQN IF (NEPS.EQ.1) M1 = (NEQN-1)/2 M1M1 = M1 - 1 IF (M1.NE.1) GO TO 70 U = (C(1)*Z(1))*5.D-01 + FX(X) GO TO 120 70 TERM = C(1)*Z(1) + C(2)*Z(2) IF (TERM.GT.0.0D00) GO TO 80 NEG = TERM POS = 0.0D00 GO TO 90 80 POS = TERM NEG = 0.0D00 90 IF (M1.EQ.2) GO TO 110 C DO 100 I=3,M1 TERM1 = C(I)*Z(I) TERM2 = -C(I)*Z(I-2) CALL SUM(TERM1, POS, NEG) CALL SUM(TERM2, POS, NEG) 100 CONTINUE C 110 U = (POS+NEG)/2.0D00 + FX(X) 120 ZP(1) = G(X,U) C REFEREE COMMENTS FOR CLARITY- C AS MENTIONED ELSEWHERE, IT IS ASSUMED THAT THERE ARE AT LEAST C 2 TERMS IN THE CHEBYSHEV EXPANSION OF THE KERNEL IF THE ERROR C ESTIMATE IS IN EFFECT (NEPS = 1); THAT IS, M1>=2 IF NEPS=1. C C C REFEREE: THANK YOU FOR CLARIFICATION AND CORRECTED LOGIC HERE. IF (M1.EQ.1) GO TO 180 C BMA4 = 4.0D00/BMA ZP(2) = BMA4*Z(1) - 2.0D00*ZP(1) C IF (M1.EQ.2) GO TO 140 ST = -2D0*ZP(1) DO 130 I=3,M1 ST = -ST RI = I ZP(I) = ZP(I-2) + BMA4*(RI-1.0D00)*Z(I-1) + ST 130 CONTINUE C C IF THERE IS TO BE NO ERROR ESTIMATE, JUST RETURN; OTHERWISE C THE RIGHT-HAND SIDE OF THE SYSTEM OF O.D.E.S IS EXPANDED TO IN- C CLUDE TERMS FOR ZP(M1+1),ZP(M1+2),...,ZP(2M1+1); HENCE, IF C NEPS = 1, THE SYSTEM CONTAINS 2M1+1 DIFFERENTIAL EQUATIONS. C C 140 IF (NEPS.EQ.0) GO TO 180 C C REFEREE CORRECTIONS: C C EXTRA O.D.E. IN Z(M+2) (=Z(M1+1)) REQUIRED FOR ERROR ESTIMATION; FM1 = M1 ZP(M1+1) = ZP(M1M1) + BMA4*FM1*Z(M1) + 2.0D0*(-1)**M1*ZP(1) C C Z(M+2) IS SOLUTION OF EXTRA O.D.E. IN Z REQUIRED FOR ERROR EST. C Z(M+3),...,Z(2M+3) ARE THE SOLUTIONS OF O.D.E.'S IN THE RESPECTIVE C QUANTITIES ZETA(0),ZETA(1),...,ZETA(M) REFERRED TO IN ARTICLE IN C JOURNAL "COMPUTING" CITED ABOVE. CAUTION: SEE CORRECTIONS FOR THIS C ARTICLE IN MATERIAL ACCOMPANYING THIS PACKAGE. ****************** C C POS = 0.0D00 NEG = 0.0D00 TERM1 = C(1)*Z(M1+2) TERM2 = C(2)*Z(M1+3) CALL SUM(TERM1, POS, NEG) CALL SUM(TERM2, POS, NEG) C IF (M1.EQ.2) GO TO 160 C C REFEREE. NOTE CHANGES IN FOLLOWING LOOP. C DO 150 I=3,M1 IM1 = I + M1 + 1 TERM1 = C(I)*Z(IM1) TERM2 = -C(I)*Z(IM1-2) CALL SUM(TERM1, POS, NEG) CALL SUM(TERM2, POS, NEG) 150 CONTINUE C 160 E = (POS+NEG+C(M1+1)*(Z(M1+1)-Z(M1-1)))/2.0D00 ZP(M1+2) = G(X,E+U) - G(X,U) IF (M1.EQ.1) GO TO 180 ZP(M1+3) = BMA4*Z(M1+2) - 2.D00*ZP(M1+2) C IF (M1.EQ.2) GO TO 180 ST = -2.0D00*ZP(M1+2) DO 170 I=3,M1 ST = -ST RI = I M1I = M1 + 1 + I M1I2 = M1I - 2 M1I1 = M1I - 1 ZP(M1I) = ZP(M1I2) + BMA4*(RI-1.0D00)*Z(M1I1) + ST 170 CONTINUE C 180 COUNTF = COUNTF + 1. RETURN C ******* END OF SUBROUTINE F ******************- END C **************** START OF SUBROUTINE DIFFEQ ********************** DIF 10 C DIF 20 SUBROUTINE DIFFEQ(F, G, FX, NEPS, M1, KA, BMA, Z, X, XOUT, RERR, DIF 30 * AERR, IFLAG) C C PURPOSE : THIS SUBROUTINE IS USED AS A HANDLER FOR THE O.D.E. C SOLVER WHICH IS SUPPLIED BY THE USER. IT DIMENSIONS THE C VARIABLES NEEDED BY THE INTEGRATOR AND HANDLES ERROR FLAGS. C C USER REQUIREMENTS : THIS PARTICULAR VERSION IS SET FOR USE WITH C THE O.D.E. SOLVER "ODE" AS LISTED IN THE BOOK, "COMPUTER C SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS"; THE C INITIAL VALUE PROBLEM", BY L.F. SHAMPINE AND M.K. GORDON, C W.H. FREEMAN AND CO., SAN FRANCISCO, 1975. C C THE AUTHORS OF THE CURRENT CODE WISH TO THANK DR. SHAMPINE C FOR HIS KIND PERMISSION TO INCLUDE THE SUITE OF CODES C (ODE,DE,STEP,AND INTRP) IN THIS PACKAGE. C C OTHER CHOICES FOR O.D.E. SOLVERS MAY REQUIRE DIFFERENT C HANDLING. FURTHERMORE, DECISIONS ON WHAT PROCEDURE TO C FOLLOW DEPENDING ON FLAG CONDITIONS MAY VARY DEPENDING C ON USER REQUIREMENTS RESULTING IN CHANGES IN THIS PARTICULAR C SUBROUTINE. C C C RE THE FOLLOWING COMPILER STATEMENT, C SEE COMMENT NO.2 IN THE DRIVER. C C COMPILER STATIC C EXTERNAL F, G, FX C C NOTE THAT IN THE FOLLOWING DIMENSION - C 1171 = 100 + 21*51. C C THIS SPECIFICATION IS DESCRIBED IN SUBROUTINE ODE(SHAMPINE). C C C RE SIZE OF ARRAYS. SEE COMMENT NO.3 IN THE DRIVER. C DOUBLE PRECISION G, FX, BMA, Z(51), X, XOUT, RERR, AERR, * WORK(1171) DOUBLE PRECISION C(51), GJ(51) C INTEGER NEPS, M1, KA, IFLAG, IWORK(5) COMMON /VE/ C, GJ, COUNTK, COUNTF C 10 CALL ODE(F, G, FX, NEPS, M1, KA, BMA, Z, X, XOUT, RERR, AERR, * IFLAG, WORK, IWORK) C C SEE LISTING FOR SUBROUTINE ODE( ) FOR DEFINITION OF IFLAG OUTPUT. C C THIS TRANSFER DETERMINES HOW THIS INTEGRAL EQN. SOLVER SHOULD REACT C TO CONDITIONS ENCOUNTERED BY THE UNDERLYING O.D.E. SOLVER. C HENCE, THE USER MAY WISH TO MODIFY THESE REACTIONS DEPENDING C ON NEEDS. THE FOLLOWING SEQUENCE DICTATES THAT THE NUMBER OF C TIMES THE RIGHT HAND SIDE OF THE O.D.E. IS EVALUATED NEVER C EXCEEDS 5000 TIMES, UNDER ANY CIRCUMSTANCES. C GO TO (20, 30, 40, 50, 60, 70), IFLAG C C IF IFLAG RETURNS WITH ABNORMAL OUTPUT VALUE (IFLAG = 1) - C 20 IF (X.EQ.XOUT) RETURN WRITE (5,99999) GO TO 80 C C IF IFLAG RETURNS IN NORMAL MODE (IFLAG = 2) - 30 RETURN C C IF IFLAG RETURNS INDICATING THAT TOLERANCES ARE TOO SMALL - 40 WRITE (5,99998) GO TO 90 C C IF IFLAG RETURNS INDICATING TOO MANY FUNCTION EVALUATIONS NEEDED - 50 WRITE (5,99997) GO TO 90 C C IF IFLAG RETURNS INDICATING POSSIBLE STIFFNESS - 60 WRITE (5,99996) GO TO 90 C C IF IFLAG RETURNS INDICATING IMPROPER INPUT - 70 WRITE (5,99995) STOP C 80 WRITE (5,99994) X STOP C 90 CHK = ALOG10(COUNTF+1.) IF (CHK.LT.5000.) GO TO 10 WRITE (5,99993) X C 99999 FORMAT (47H ADAMS TERMINATED WITH IFLAG=1, BUT X .NE. XOUT) 99998 FORMAT (36H TOLS. RESET UNLESS FCOUNT.GT.5000 .) 99997 FORMAT (52H INTEGRATION INTERRUPTED DUE TO HIGH FUNCTION COUNT.) 99996 FORMAT (20H POSSIBLE STIFFNESS.) 99995 FORMAT (43H IMPROPER INPUT TO ODE. FATAL ERROR - STOP.) 99994 FORMAT (42H UNABLE TO CONTINUE. CURRENT VALUE OF X = , D20.12) 99993 FORMAT (42H FUNCTION COUNT EXCEEDS 5000. CURRENT X = , D20.12) C STOP END C ********* END OF SUBROUTINE DIFFEQ ************* ODE 10 C ************* START OF SUBROUTINE ODE (L.SHAMPINE, ET. AL.) ******** ODE 20 C ODE 30 SUBROUTINE ODE(F, G, FX, NEPS, NEQN, KA, BMA, Y, T, TOUT, RELERR, ODE 40 * ABSERR, IFLAG, WORK, IWORK) C C NOTE : C THE ARGUMENT LIST ABOVE AND CALL TO DE HAVE BEEN MODIFIED (LAA). C C SUBROUTINE ODE INTEGRATES A SYSTEM OF NEQN FIRST ORDER C ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM C DY(I)/DT = F(T,Y(1),Y(2),...,Y(NEQN)) C Y(I) GIVEN AT T. C THE SUBROUTINE INTEGRATES FROM T TO TOUT . ON RETURN THE C PARAMETERS IN THE CALL LIST ARE SET FOR CONTINUING THE INTEGRATION. C THE USER HAS ONLY TO DEFINE A NEW VALUE TOUT AND CALL ODE AGAIN. C C THE DIFFERENTIAL EQUATIONS ARE ACTUALLY SOLVED BY A SUITE OF CODES C DE , STEP , AND INTRP . ODE ALLOCATES VIRTUAL STORAGE IN THE C ARRAYS WORK AND IWORK AND CALLS DE . DE IS A SUPERVISOR WHICH C DIRECTS THE SOLUTION. IT CALLS ON THE ROUTINES STEP AND INTRP C TO ADVANCE THE INTEGRATION AND TO INTERPOLATE AT OUTPUT POINTS. C STEP USES A MODIFIED DIVIDED DIFFERENCE FORM OF THE ADAMS C FORMULAS AND LOCAL EXTRAPOLATION. IT ADJUSTS THE ORDER AND STEP C SIZE TO CONTROL THE LOCAL EERROR PER UNIT STEP IN A GENERALIZED C SENSE. NORMALLY EACH CALL TO STEP ADVANCES THE SOLUTION ONE STEP C IN THE DERECTION OF TOUT . FOR REASONS OF EFFICIENCY DE C INTEGRATES BEYOND TOUT INTERNALLY, THOUGH NEVER BEYOND C T+10*(TOUT-T), AND CALLS INTRP TO INTERPOLATE THE SOLUTION AT C TOUT . AN OPTION IS PROVIDED TO STOP THE INTEGRATION AT TOUT BUT C IT SHOULD BE USED ONLY IF IT IS IMPOSSIBLE TO CONTINUE THE C INTEGRATION BEYOND TOUT . C C THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT, C COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS: THE INITIAL C VALUE PROBLEM BY L. F. SHAMPINE AND M. K. GORDON. C C THE PARAMETERS REPRESENT: C F -- SUBROUTINE F(T,Y,YP) TO EVALUATE DERIVATIVES YP(I)=DY(I)/DT C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED C Y(*) -- SOLUTION VECTOR AT T C T -- INDEPENDENT VARIABLE C TOUT -- POINT AT WHICH SOLUTION IS DESIRED C RELERR,ABSERR -- RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR LOCAL C ERROR TEST. AT EACH STEP THE CODE REQUIRES C DABS(LOCAL ERROR) .LE. DABS(Y)*RELERR + ABSERR C FOR EACH COMPONENT OF THE LOCAL ERROR AND SOLUTION VECTORS C IFLAG -- INDICATES STATUS OF INTEGRATION C WORK(*),IWORK(*) -- ARRAYS TO HOLD INFORMATION INTERNAL TO CODE C WHICH IS NECESSARY FOR SUBSEQUENT CALLS C C FIRST CALL TO ODE -- C C THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR THE ARRAYS C IN THE CALL LIST C Y(NEQN), WORK(100+21*NEQN), IWORK(5) C DECLARE F IN AN EXTERNAL STATEMENT, SUPPLY THE SUBROUTINE C F(T,Y,YP) TO EVALUATE C DY(I)/DT = YP(I) = F(T,Y(1),Y(2),...,Y(NEQN)) C AND INITIALIZE THE PARAMETERS: C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED C Y(*) -- VECTOR OF INITIAL CONDITIONS C T -- STARTING POINT OF INTEGRATION C TOUT -- POINT AT WHICH SOLUTION IS DESIRED C RELERR,ABSERR -- RELATIVE AND ABSOLUTE LOCAL ERROR TOLERANCES C IFLAG -- +1,-1. INDICATOR TO INITIALIZE THE CODE. NORMAL INPUT C IS +1. THE USER SHOULD SET IFLAG = -1 ONLY IF IT IS C IMPOSSIBLE TOR CONTINUE THE INTEGRATION BEYOND TOUT . C ALL PARAMETERS EXCEPT F . NEQN AND TOUT MAY BE ALTERED BY THE C CODE ON OUTPUT SO MUST BE VARIABLES IN THE CALLING PROGRAM. C C OUTPUT FROM ODE -- C C NEQN -- UNCHANGED C Y(*) -- SOLUTION AT TOUT C T -- LAST POINT REACHED IN INTEGRATION. NORMAL RETURN HAS C T = TOUT . C TOUT -- UNCHANGED C RELERR,ABSERR -- NORMAL RETURN HAS TOLERANCES UNCHANGED. IFLAG=3 C SIGNALS TOLERANCES INCREASED C IFLAG = 2 -- NORMAL RETURN. INTEGRATION REACHED TOUT C = 3 -- INTEGRATION DID NOT REACH TOUT BECAUSE ERROR C TOLERANCES TOO SMALL. RELERR , ABSERR INCREASED C APPROPRIATELY FOR CONTINUING C = 4 -- INTEGRATION DID NOT REACH TOUT BECAUSE EQUATIONS C 500 STEPS NEEDED C = 5 -- INTEGRATION DID NOT REACH TOUT BECAUSE EQUATIONS C APPEAR TO BE STIFF C = 6 -- INVALID INPUT PARAMETERS (FATAL ERROR) C THE VALUE OF IFLAG IS RETURNED NEGATIVE WHEN THE INPUT C VALUE IS NEGATIVE AND THE INTEGRATION DOES NOT REACH TOUT, C I.E. , -3, -4, -5. C WORK(*),IWORK(*) -- INFORMATION GENERALLY OF NO INTEREST TO THE C USER BUT NECESSARY FOR SUBSEQUENT CALLS. C C SUBSEQUENT CALLS TO ODE -- C C SUBROUTINE ODE RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE C THE INTEGRATION. IF THE INTEGRATION REACHED TOUT , THE USER NEED C ONLY DEFINE A NEW TOUT AND CALL AGAIN. IF THE INTEGRATION DID NOT C REACH TOUT AND THE USER WANTS TO CONTINUE, HE JUST CALLS AGAIN. C THE OUTPUT VALUE OF IFLAG IS THE APPROPRIATE INPUT VALUE FOR C SUBSEQUENT CALLS. THE ONLY SITUATION IN WHICH IT SHOULD BE ALTERED C IS TO STOP THE INTEGRATION INTERNALLY AT THE NEW TOUT , I.E., C CHANGE OUTPUT IFLAG=2 TO INPUT IFLAG=-2 . ERROR TOLERANCES MAY C BE CHANGED BY THE USER BEFORE CONTINUING. ALL OTHER PARAMETERS MUST C REMAIN UNCHANGED. C C********************************************************************** C* SUBROUTINES DE AND STEP CONTAIN MACHINE DEPENDENT CONSTANTS. * C* BE SURE THEY ARE SET BEFORE USING ODE . * C********************************************************************** C C C RE THE FOLLOWING COMPILER STATEMENT, C SEE COMMENT NO.2 IN THE DRIVER. C C COMPILER STATIC C DOUBLE PRECISION Y, WORK, F, G, FX, BMA, T DOUBLE PRECISION TOUT, RELERR, ABSERR LOGICAL START, PHASE1, NORND C C RE SIZE OF ARRAYS. SEE COMMENT NO.3 IN THE DRIVER. C ALSO, REGARDING THE SIZE OF WORK( ), SEE COMMENT IN SUB ODE. C DIMENSION Y(NEQN), WORK(1) C C DIMENSION IWORK(5) EXTERNAL F, G, FX DATA IALPHA, IBETA, ISIG, IV, IW, IG, IPHASE, IPSI, IX, IH, * IHOLD, ISTART, ITOLD, IDELSN /1,13,25,38,50,62,75,76,88,89,90,91, * 92,93/ C C CHECK FOR FIRST ENTRY INTO ROUTINE (IFLAG=1); IF T = TOUT, (BOWNDS) C INITIAL CONDITIONS ARE SIMPLY RETURNED. THIS CHANGE (BOWNDS) C MAKES NORMAL MODE (IFLAG=1) USE MORE COMPATIBLE (BOWNDS) C WITH OTHER GENERAL PURPOSE O.D.E. SOLVERS. (BOWNDS) C C************************************************************** (BOWNDS) IF (IFLAG.EQ.1 .AND. T.EQ.TOUT) RETURN C************************************************************** (BOWNDS) IYY = 100 IWT = IYY + NEQN IP = IWT + NEQN IYP = IP + NEQN IYPOUT = IYP + NEQN IPHI = IYPOUT + NEQN IF (IABS(IFLAG).EQ.1) GO TO 10 START = WORK(ISTART).GT.0D0 PHASE1 = WORK(IPHASE).GT.0D0 N1 = -1 NORND = IWORK(2).NE.N1 10 CALL DE(F, G, FX, NEPS, NEQN, KA, BMA, Y, T, TOUT, RELERR, * ABSERR, IFLAG, WORK(IYY), WORK(IWT), WORK(IP), WORK(IYP), * WORK(IYPOUT), WORK(IPHI), WORK(IALPHA), WORK(IBETA), WORK(ISIG), * WORK(IV), WORK(IW), WORK(IG), PHASE1, WORK(IPSI), WORK(IX), * WORK(IH), WORK(IHOLD), START, WORK(ITOLD), WORK(IDELSN), * IWORK(1), NORND, IWORK(3), IWORK(4), IWORK(5)) WORK(ISTART) = -1.0D0 IF (START) WORK(ISTART) = 1.0D0 WORK(IPHASE) = -1.0D0 IF (PHASE1) WORK(IPHASE) = 1.0D0 IWORK(2) = -1 IF (NORND) IWORK(2) = 1 RETURN C ****** END OF SUBROUTINE ODE ************************* END C ********** START OF SUBROUTINE DE (L.SHAMPINE, ET. AL.) *********** DE 10 C DE 20 SUBROUTINE DE(F, GFUN, FX, NEPS, NEQN, KA, BMA, Y, T, TOUT, DE 30 * RELERR, DABERR, IFLAG, YY, WT, P, YP, YPOUT, PHI, ALPHA, BETA, * SIG, V, W, G, PHASE1, PSI, X, H, HOLD, START, TOLD, DELSGN, NS, * NORND, K, KOLD, ISNOLD) C C ODE MERELY ALLOCATES STORAGE FOR DE TO RELIEVE THE USER OF THE C INCONVENIENCE OF A LONG CALL LIST. CONSEQUENTLY DE IS USED AS C DESCRIBED IN THE COMMENTS FOR ODE . C C NOTE : C THE ABOVE ARGUMENT LIST AND CALLS TO STEP HAVE BEEN MODIFIED (LAA). C C THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT, C COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS: THE INITIAL C VALUE PROBLEM BY L. F. SHAMPINE AND M. K. GORDON. C C C RE THE FOLLOWING COMPILER STATEMENT, C SEE COMMENT NO.2 IN THE DRIVER. C C COMPILER STATIC C LOGICAL START, CRASH, STIFF, PHASE1, NORND DOUBLE PRECISION ALPHA(12), BETA(12), SIG(12), V(12), W(12), G(13) C C RE SIZE OF ARRAYS. SEE COMMENT NO.3 IN THE DRIVER. C DOUBLE PRECISION Y(NEQN), PSI(12), YP(NEQN), YPOUT(NEQN) DOUBLE PRECISION YY(NEQN), WT(NEQN), PHI(NEQN,16), P(NEQN) C DOUBLE PRECISION T, TOUT, RELERR, DABERR, EPS, DEL, DABDEL, BMA, * TEND, RELEPS, DABEPS, X, DELSGN, H, TOLD, HOLD, FOURU, GFUN, FX, * A C EXTERNAL F, GFUN, FX C C********************************************************************** C* THE ONLY MACHINE DEPENDENT CONSTANT IS BASED ON THE MACHINE UNIT * C* ROUNDOFF ERROR U WHICH IS THE SMALLEST POSITIVE NUMBER SUCH THAT* C* 1.0+U .GT. 1.0 . U MUST BE CALCULATED AND FOURU=4.0*U INSERTED * C* IN THE FOLLOWING DATA STATEMENT BEFORE USING DE . THE ROUTINE * C* MACHIN CALCULATES U . FOURU AND TWO=2.0*U MUST ALSO BE * C* INSERTED IN SUBROUTINE STEP BEFORE CALLING DE . * C********************************************************************** C C* THE FOLLOWING VALUE FOR FOURU IS SPECIFICALLY FOR DATA GENERAL, C* S-230 ECLIPSE. OTHER COMPUTERS REQUIRE OTHER VALUES : C ************************************************************* DATA FOURU /.8881784197001252D-15/ C ************************************************************ C THE CONSTANT MAXNUM IS THE MAXIMUM NUMBER OF STEPS ALLOWED IN ONE C CALL TO DE . THE USER MAY CHANGE THIS LIMIT BY ALTERING THE C FOLLOWING STATEMENT DATA MAXNUM /500/ C C *** *** *** C TEST FOR IMPROPER PARAMETERS C IF (NEQN.LT.1 .OR. NEQN.GT.51) GO TO 10 IF (T.EQ.TOUT) GO TO 10 IF (RELERR.LT.0D0 .OR. DABERR.LT.0D0) GO TO 10 EPS = DMAX1(RELERR,DABERR) IF (EPS.LE.0D0) GO TO 10 IF (IFLAG.EQ.0) GO TO 10 ISN = ISIGN(1,IFLAG) IFLAG = IABS(IFLAG) IF (IFLAG.EQ.1) GO TO 20 C IF(T .NE. TOLD) GO TO 10 IF (IFLAG.GE.2 .AND. IFLAG.LE.5) GO TO 20 10 IFLAG = 6 RETURN C C ON EACH CALL SET INTERVAL OF INTEGRATION AND COUNTER FOR NUMBER OF C STEPS. ADJUST INPUT ERROR TOLERANCES TO DEFINE WEIGHT VECTOR FOR C SUBROUTINE STEP C 20 DEL = TOUT - T DABDEL = DABS(DEL) TEND = T + 1.0D1*DEL IF (ISN.LT.0) TEND = TOUT NOSTEP = 0 KLE4 = 0 STIFF = .FALSE. RELEPS = RELERR/EPS DABEPS = DABERR/EPS IF (IFLAG.EQ.1) GO TO 30 IF (ISNOLD.LT.0) GO TO 30 IF (DELSGN*DEL.GT.0D0) GO TO 50 C C ON START AND RESTART ALSO SET WORK VARIABLES X AND YY(*), STORE THE C DIRECTION OF INTEGRATION AND INITIALIZE THE STEP SIZE C 30 START = .TRUE. X = T DO 40 L=1,NEQN YY(L) = Y(L) 40 CONTINUE DELSGN = DSIGN(1.0D0,DEL) H = DSIGN(DMAX1(DABS(TOUT-X),FOURU*DABS(X)),TOUT-X) C C IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND RETURN C 50 IF (DABS(X-T).LT.DABDEL) GO TO 60 CALL INTRP(X, YY, TOUT, Y, YPOUT, NEQN, KOLD, PHI, PSI) IFLAG = 2 T = TOUT TOLD = T ISNOLD = ISN RETURN C C IF CANNOT GO PAST OUTPUT POINT AND SUFFICIENTLY CLOSE, C EXTRAPOLATE AND RETURN C 60 IF (ISN.GT.0 .OR. DABS(TOUT-X).GE.FOURU*DABS(X)) GO TO 80 H = TOUT - X CALL F(GFUN, FX, NEPS, NEQN, KA, BMA, X, YY, YP) DO 70 L=1,NEQN Y(L) = YY(L) + H*YP(L) 70 CONTINUE IFLAG = 2 T = TOUT TOLD = T ISNOLD = ISN RETURN C C TEST FOR TOO MUCH WORK C 80 IF (NOSTEP.LT.MAXNUM) GO TO 100 IFLAG = ISN*4 IF (STIFF) IFLAG = ISN*5 DO 90 L=1,NEQN Y(L) = YY(L) 90 CONTINUE T = X TOLD = T ISNOLD = 1 RETURN C C LIMIT STEP SIZE, SET WEIGHT VECTOR AND TAKE A STEP C 100 H = DSIGN(DMIN1(DABS(H),DABS(TEND-X)),H) DO 110 L=1,NEQN WT(L) = RELEPS*DABS(YY(L)) + DABEPS 110 CONTINUE CALL STEP(F, GFUN, FX, NEPS, NEQN, KA, BMA, YY, X, H, EPS, WT, * START, HOLD, K, KOLD, CRASH, PHI, P, YP, PSI, ALPHA, BETA, SIG, * V, W, G, PHASE1, NS, NORND) C C TEST FOR TOLERANCES TOO SMALL C IF (.NOT.CRASH) GO TO 130 IFLAG = ISN*3 RELERR = EPS*RELEPS DABERR = EPS*DABEPS DO 120 L=1,NEQN Y(L) = YY(L) 120 CONTINUE T = X TOLD = T ISNOLD = 1 RETURN C C AUGMENT COUNTER ON WORK AND TEST FOR STIFFNESS C 130 NOSTEP = NOSTEP + 1 KLE4 = KLE4 + 1 IF (KOLD.GT.4) KLE4 = 0 IF (KLE4.GE.50) STIFF = .TRUE. GO TO 50 C ******* END OF SUBROUTINE DE ************************ END C *********** START OF SUBROUTINE STEP (L.SHAMPINE, ET. AL.) ******** STE 10 C STE 20 SUBROUTINE STEP(F, GFUN, FX, NEPS, NEQN, KA, BMA, Y, X, H, EPS, STE 30 * WT, START, HOLD, K, KOLD, CRASH, PHI, P, YP, PSI, ALPHA, BETA, * SIG, V, W, G, PHASE1, NS, NORND) C C NOTE : C THE ABOVE ARGUMENT LIST AND CALLS TO F HAVE BEEN MODIFIED (LAA). C C SUBROUTINE STEP INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY C DIFFERENTIAL EQUATIONS ONE STEP NORMALLY FROM X TO X+H, USING A C MODIFIED DIVIDED DIFFENCE FORM OF THE ADAMS PECE FORMULAS. LOCAL C EXTRAPOLATION IS USED TO IMPROVE ABSOLUTE STABILITY AND ACCURACY. C THE CODE ADJUSTS ITS ORDER AND STEP SIZE TO CONTROL THE LOCAL ERROR C PER UNIT STEP IN A GENERALIZED SENSE. SPECIAL DEVICES ARE INCLUDED C TO CONTROL ROUNDOFF ERROR AND TO DETECT WHEN THE USER IS REQUESTING C TOO MUCH ACCURACY. C C THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT, C COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS: THE INITIAL C VALUE PROBLEM BY L. F. SHAMPINE AND M. K. GORDON. C C C THE PARAMETERS REPRESENT: C X -- INDEPENDENT VARIABLE C Y(*) -- SOLUTION VECTOR AT X C YP(*) -- DERIVATIVE OF SOLUTION VECTOR AT X AFTER SUCCESSFUL C STEP C NEQM -- NUMBER OF EQUATIONS TO BE INTEGRATED C H -- APPROPRIATE STEP SIZE FOR NEXT STEP. NORMALLY DETERMINED BY C CODE C EPS -- LOCAL ERROR TOLERANCE. MUST BE VARIABLE C WT(*) -- VECTOR OF WEIGHTS FOR ERROR CRITERION C START -- LOGICAL VARIABLE SET .TRUE. FOR FIRST STEP, .FALSE. C OTHERWISE C HOLD -- STEP SIZE USED FOR LAST SUCCESSFUL STEP C K -- APPROPRIATE ORDER FOR NEXT STEP (DETERMINED BY CODE) C KOLD -- ORDER USED FOR LAST SUCCESSFUL STEP C CRASH -- LOGICAL VARIABLE SET .TRUE. WHEN NO STEP CAN BE TAKEN, C .FALSE. OTHERWISE C THE ARRAYS PHI, PSI ARE REQUIRED FOR THE INTERPOLATION SUBROUTINE C INTRP . THE ARRAY P IS INTERNAL TO THE CODE. C C INPUT TO STEP C C FIRST CALL -- C C THE USER MUST PROVIDE STORAGE IN HIS DRIVER PROGRAM FOR ALL ARRAYS C IN THE CALL LIST, NAMELY C C DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12) C C THE USER MUST ALSO DECLARE START AND CRASH LOGICAL VARIABLES C AND F AN EXTERNAL SUBROUTINE, SUPPLY THE SUBROUTINE F(X,Y,YP) C TO EVALUATE C DY(I)/DX = YP(I) = F(X,Y(1),Y(2),...,Y(NEQN)) C AND INITIALIZE ONLY THE FOLLOWING PARAMETERS: C X -- INITIAL VALUE OF THE INDEPENDENT VARIABLE C Y(*) -- VECTOR OF INITIAL VALUES OF DEPENDENT VARIABLES C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED C H -- NOMINAL STEP SIZE INDICATING DIRECTION OF INTEGRATION C AND MAXIMUM SIZE OF STEP. MUST BE VARIABLE C EPS -- LOCAL ERROR TOLERANCE PER STEP. MUST BE VARIABLE C WT(*) -- VECTOR OF NON-ZERO WEIGHTS FOR ERROR CRITERION C C -- .TRUE. C C STEP REQUIRES THE L2 NORM OF THE VECTOR WITH COMPONENTS C LOCAL ERROR(L)/WT(L) BE LESS THAN EPS FOR A SUCCESSFUL STEP. THE C ARRAY T ALLOWS THE USER TO SPECIFY AN ERROR TEST APPROPRIATE C FOR HIS PROBLEM. FOR EXAMPLE, C WT(L) = 1.0 SPECIFIES DABSOLUTE ERROR C = DABS(Y(L)) ERROR RELATIVE TO THE MOST RECENT VALUE OF THE C L-TH COMPONENT OF THE SOLUTION, C = DABS(YP(L)) ERROR RELATIVE TO THE MOST RECENT VALUE OF C THE L-TH COMPONENT OF THE DERIVATIVE, C = DMAX1(WT(L),DABS(Y(L)) ERROR RELATIVE TO THE LARGEST C MAGNITUDE OF L-TH COMPONENT OBTINED SO FAR, C = DABS(Y(L))*RELERR/EPS SPECIFIES A MIXED C RELATIVE-DABSOLUTE TEST WHERE RELERR IS RELATIVE C ERROR, DABSERR IS DABSOLUTE ERROR AND EPS = C DMAX1(RELERR,DABSERR) C C SUBSEQUENT CALLS -- C C SUBROUTINE STEP IS DESIGNED SO THAT ALL INFORMATION NEEDED TO C CONTINUE THE INTEGRATION, INCLUDING THE STEP SIZE H AND THE ORDER C K , IS RETURNED WITH EACH STEP. WITH THE EXCEPTION OF THE STEP C SIZE, THE ERROR TOLERANCE, AND THE WEIGHTS, NONE OF THE PARAMETERS C SHOULD BE ALTERED. THE ARRAY WY MUST BE UPDATED AFTER EACH STEP C TO MAINTAIN RELATIVE ERROR TESTS LIKE THOSE ABOVE. NORMALLY THE C INTEGRATION IS CONTINUED JUST BEYOND THE DESIRED ENDPOINT AND THE C SOLUTION INTERPOLATED THERE WITH SUBROUTINE INTRP . IF IT IS C IMPOSSIBLE TO INTEGRATE BEYOND THE ENDPOINT, THE STEP SIZE MAY BE C REDUCED TO HIT THE ENDPOINT SINCE THE CODE WILL NOT TAKE A STEP C LARGER THAN THE H INPUT. CHANGING THE DIRECTION OF INTEGRATION, C I.E., THE DSIGN OF H , REQUIRES THE USER SET START = .TRUE. BEFORE C CALLING STEP AGAIN. THIS IS THE ONLY SITUATION IN WHICH START C SHOULD BE ALTERED. C C OUTPUT FROM STEP C C SUCCESSFUL STEP -- C C THE SUBROUTINE RETURNS AFTER EACH SUCCESSFUL STEP WITH START AND C CRASH SET .FALSE. . X REPRESENTS THE INDEPENDENT VARIABLE C ADVANCED ONE STEP OF LENGTH HOLD FROM ITS VALUE ON INPUT AND Y C THE SOLUTION VECTOR AT THE NEW VALUE OF X . ALL OTHER PARAMETERS C REPRESENT INFOMATION CORRESPONDING TO THE NEW X NEEDED TO C CONTINUE THE INTEGRATION. C C UNSUCCESSFUL STEP -- C C WHEN THE ERROR TOLERANCE IS TOO SMALL FOR THE MACHINE PRECISION, C THE SUBROUTINE RETURNS WITHOUT TAKING A STEP AND CRASH = .TRUE. C AN APPROPRIATE STEP SIZE AND ERROR TOLERANCE FOR CONTINUING ARE C ESTIMATED AND ALL OTHER INFORMATION IS RESTORED A UPON INPUT C BEFORE RETURNING. TO CONTINUE WITH THE LARGER TOLERANCE, THE USER C JUST CALLS THE CODE AGAIN. A RESTART IS NEITHER REQUIRED NOR C DESIRABLE. C C C RE THE FOLLOWING COMPILER STATEMENT, C SEE COMMENT NO.2 IN THE DRIVER. C C COMPILER STATIC C LOGICAL START, CRASH, PHASE1, NORND C C RE ARRAY SIZE. SEE COMMENT NO.3 IN THE DRIVER. C DOUBLE PRECISION Y(NEQN), WT(NEQN), PHI(NEQN,16), P(NEQN), * YP(NEQN), PSI(12) C DOUBLE PRECISION ALPHA(12), BETA(12), SIG(13), W(12), V(12), * G(13), GSTR(13), TWO(13), GFUN, BMA, FX, A DOUBLE PRECISION X, H, EPS, HOLD, XOLD, P5EPS, ROUND, SUM, DABSH, * TEMP1, TEMP2, TEMP3, TEMP4, TEMP5, TEMP6, ERR, ERK, ERKM1, * ERKM2, ERKP1, HNEW, R, TAU, RHO, REALI, REALNS, FOURU, TWOU C EXTERNAL F, GFUN, FX C********************************************************************* C THE ONLY MACHINE DEPENDENT CONSTANTS ARE BASED ON THE MACHINE UNIT * C ROUNDOFF ERROR U WHICH IS THE SMALLEST POSITIVE NUMBER SUCH THAT * C 1.0+U .GT. 1.0 . THE USER MUST CALCULATE U AND INSERT * C TWOU=2.0D0*U AND FOURU=4.0*U IN THE DATA STATEMENT BEFORE A CALL* C TO THE CODE. * C * C ******************************************************************** C THE FOLLOWING CONSTANTS ARE SPECIFICALLY FOR USE WITH THE C DATA GENERAL S-230 ECLIPSE. OTHER COMPUTERS WILL REQUIRE C DIFFERENT VALUES : DATA TWOU, FOURU /.4440892098500626D-15,.8881784197001252D-15/ C C ################################################################# C C THE FOLLOWING CONSTANTS ARE NOT MACHINE DEPENDENT. C C DATA TWO(1), TWO(2), TWO(3), TWO(4), TWO(5) /2.0D0,4.0D0,8.0D0, * 1.6D1,3.2D1/ DATA TWO(6), TWO(7), TWO(8), TWO(9), TWO(10) /6.4D1,1.28D2,2.56D2, * 5.12D2,1.024D3/ DATA TWO(11), TWO(12), TWO(13) /2.048D3,4.096D3,8.192D3/ DATA GSTR(1), GSTR(2), GSTR(3), GSTR(4), GSTR(5) * /0.5D0,0.833D-1,0.417D-1,0.264D-1,0.188D-1/ DATA GSTR(6), GSTR(7), GSTR(8), GSTR(9), GSTR(10) * /0.143D-1,0.114D-1,0.936D-2,0.789D-2,0.679D-2/ DATA GSTR(11), GSTR(12), GSTR(13) /0.592D-2,0.594D-2,0.468D-2/ C C C *** BEGIN BLOCK 0 *** C CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE C PRECISION. IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A C STARTING STEP SIZE. C *** C C IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE C CRASH = .TRUE. IF (DABS(H).GE.FOURU*DABS(X)) GO TO 10 H = DSIGN(FOURU*DABS(X),H) RETURN 10 P5EPS = 0.5D0*EPS C C IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE C ROUND = 0D0 DO 20 L=1,NEQN ROUND = ROUND + (Y(L)/WT(L))**2 20 CONTINUE ROUND = TWOU*DSQRT(ROUND) IF (P5EPS.GE.ROUND) GO TO 30 EPS = 2.0D0*ROUND*(1.0D0+FOURU) RETURN 30 CRASH = .FALSE. G(1) = 1.0D0 G(2) = 0.5D0 SIG(1) = 1.0D0 IF (.NOT.START) GO TO 60 C C INITIALIZE. COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP C CALL F(GFUN, FX, NEPS, NEQN, KA, BMA, X, Y, YP) SUM = 0D0 DO 40 L=1,NEQN PHI(L,1) = YP(L) PHI(L,2) = 0D0 SUM = SUM + (YP(L)/WT(L))**2 40 CONTINUE SUM = DSQRT(SUM) DABSH = DABS(H) IF (EPS.LT.1.6D1*SUM*H*H) DABSH = 0.25D0*DSQRT(EPS/SUM) H = DSIGN(DMAX1(DABSH,FOURU*DABS(X)),H) HOLD = 0D0 K = 1 KOLD = 0 START = .FALSE. PHASE1 = .TRUE. NORND = .TRUE. IF (P5EPS.GT.100D0*ROUND) GO TO 60 NORND = .FALSE. DO 50 L=1,NEQN PHI(L,15) = 0.0 50 CONTINUE 60 IFAIL = 0 C *** END BLOCK 0 *** C C *** BEGIN BLOCK 1 *** C COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP. AVOID COMPUTING C THOSE QUANTITNIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED. C *** C 70 KP1 = K + 1 KP2 = K + 2 KM1 = K - 1 KM2 = K - 2 C C NS IS THE NUMBER OF STEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT C ONE. WHEN K.LT.NS, NO COEFFICIENTS CHANGE. C IF (H.NE.HOLD) NS = 0 NS = MIN0(NS+1,KOLD+1) NSP1 = NS + 1 IF (K.LT.NS) GO TO 180 C C COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH C ARE CHANGED C BETA(NS) = 1.0D0 REALNS = NS ALPHA(NS) = 1.0D0/REALNS TEMP1 = H*REALNS SIG(NSP1) = 1.0D0 IF (K.LT.NSP1) GO TO 90 DO 80 I=NSP1,K IM1 = I - 1 TEMP2 = PSI(IM1) PSI(IM1) = TEMP1 BETA(I) = BETA(IM1)*PSI(IM1)/TEMP2 TEMP1 = TEMP2 + H ALPHA(I) = H/TEMP1 REALI = I SIG(I+1) = REALI*ALPHA(I)*SIG(I) 80 CONTINUE 90 PSI(K) = TEMP1 C C COMPUTE COEFFICIENTS G(*) C C INITIALIZE V(*) AND SET W(*). G(2) IS SET IN DATA STATEMENT C IF (NS.GT.1) GO TO 110 DO 100 IQ=1,K TEMP3 = IQ*(IQ+1) V(IQ) = 1.0D0/TEMP3 W(IQ) = V(IQ) 100 CONTINUE GO TO 150 C C IF ORDER WAS RAISED, UPDATE DIAGONAL PART OD V(*) C 110 IF (K.LE.KOLD) GO TO 130 TEMP4 = K*KP1 V(K) = 1.0D0/TEMP4 NSM2 = NS - 2 IF (NSM2.LT.1) GO TO 130 DO 120 J=1,NSM2 I = K - J V(I) = V(I) - ALPHA(J+1)*V(I+1) 120 CONTINUE C C UPDATE V(*) AND SET W(*) C 130 LIMIT1 = KP1 - NS TEMP5 = ALPHA(NS) DO 140 IQ=1,LIMIT1 V(IQ) = V(IQ) - TEMP5*V(IQ+1) W(IQ) = V(IQ) 140 CONTINUE G(NSP1) = W(1) C C COMPUTE THE G(*) IN THE WORK VECTOR W(*) C 150 NSP2 = NS + 2 IF (KP1.LT.NSP2) GO TO 180 DO 170 I=NSP2,KP1 LIMIT2 = KP2 - I TEMP6 = ALPHA(I-1) DO 160 IQ=1,LIMIT2 W(IQ) = W(IQ) - TEMP6*W(IQ+1) 160 CONTINUE G(I) = W(1) 170 CONTINUE 180 CONTINUE C *** END BLOCK 1 *** C C *** BEGIN BLOCK 2 *** C PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED C SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K, C K-1, K-2, AS IF CONSTANT STEP SIZE WERE USED. C C C CHANGE PHI TO PHI STAR C IF (K.LT.NSP1) GO TO 210 DO 200 I=NSP1,K TEMP1 = BETA(I) DO 190 L=1,NEQN PHI(L,I) = TEMP1*PHI(L,I) 190 CONTINUE 200 CONTINUE C C PREDICT SOLUTION AND DIFFERENCES C 210 DO 220 L=1,NEQN PHI(L,KP2) = PHI(L,KP1) PHI(L,KP1) = 0D0 P(L) = 0D0 220 CONTINUE DO 240 J=1,K I = KP1 - J IP1 = I + 1 TEMP2 = G(I) DO 230 L=1,NEQN P(L) = P(L) + TEMP2*PHI(L,I) PHI(L,I) = PHI(L,I) + PHI(L,IP1) 230 CONTINUE 240 CONTINUE IF (NORND) GO TO 260 DO 250 L=1,NEQN TAU = H*P(L) - PHI(L,15) P(L) = Y(L) + TAU PHI(L,16) = (P(L)-Y(L)) - TAU 250 CONTINUE GO TO 280 260 DO 270 L=1,NEQN P(L) = Y(L) + H*P(L) 270 CONTINUE 280 XOLD = X X = X + H DABSH = DABS(H) CALL F(GFUN, FX, NEPS, NEQN, KA, BMA, X, P, YP) C C ESTIMATE ERRORS AT ORDERS K,K-1,K-2 C ERKM2 = 0D0 ERKM1 = 0D0 ERK = 0D0 DO 320 L=1,NEQN TEMP3 = 1.0/WT(L) TEMP4 = YP(L) - PHI(L,1) IF (KM2) 310, 300, 290 290 ERKM2 = ERKM2 + ((PHI(L,KM1)+TEMP4)*TEMP3)**2 300 ERKM1 = ERKM1 + ((PHI(L,K)+TEMP4)*TEMP3)**2 310 ERK = ERK + (TEMP4*TEMP3)**2 320 CONTINUE IF (KM2) 350, 340, 330 330 ERKM2 = DABSH*SIG(KM1)*GSTR(KM2)*DSQRT(ERKM2) 340 ERKM1 = DABSH*SIG(K)*GSTR(KM1)*DSQRT(ERKM1) 350 TEMP5 = DABSH*DSQRT(ERK) ERR = TEMP5*(G(K)-G(KP1)) ERK = TEMP5*SIG(KP1)*GSTR(K) KNEW = K C C TEST IF ORDER SHOULD BE LOWERED C IF (KM2) 380, 370, 360 360 IF (DMAX1(ERKM1,ERKM2).LE.ERK) KNEW = KM1 GO TO 380 370 IF (ERKM1.LE.0.5D0*ERK) KNEW = KM1 C C TEST IF STEP SUCCESSFUL C 380 IF (ERR.LE.EPS) GO TO 470 C *** END BLOCK 2 *** C C *** BEGIN BLOCK 3 *** C THE STEP IS UNSUCCESSFUL. RESTORE X, PHI(*,*), PSI(*) . C IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE. IF STEP FAILS MORE C THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE. DOUBLE ERROR C TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINE C PRECISION. C *** C C RESTORE X, PHI(*,*) AND PSI(*) C PHASE1 = .FALSE. X = XOLD DO 400 I=1,K TEMP1 = 1.0D0/BETA(I) IP1 = I + 1 DO 390 L=1,NEQN PHI(L,I) = TEMP1*(PHI(L,I)-PHI(L,IP1)) 390 CONTINUE 400 CONTINUE IF (K.LT.2) GO TO 420 DO 410 I=2,K PSI(I-1) = PSI(I) - H 410 CONTINUE C C ON THIRD FAILURE, SET ORDER TO ONE. THEREAFTER, USE OPTIMAL STEP C SIZE C 420 IFAIL = IFAIL + 1 TEMP2 = 0.5D0 IF (IFAIL-3) 450, 440, 430 430 IF (P5EPS.LT.0.25D0*ERK) TEMP2 = DSQRT(P5EPS/ERK) 440 KNEW = 1 450 H = TEMP2*H K = KNEW IF (DABS(H).GE.FOURU*DABS(X)) GO TO 460 CRASH = .TRUE. H = DSIGN(FOURU*DABS(X),H) EPS = EPS + EPS RETURN 460 GO TO 70 C *** END BLOCK 3 *** C C *** BEGIN BLOCK 4 *** C THE STEP IS SUCCESSFUL. CORRECT THE PREDICTED SOLUTION, EVALUATE C THE DERIVATIVES USING THE CORRECTED SOLUTION AND UPDATE THE C DIFFERENCES. DETERMINE BEST ORDER AND STEP SIZE FOR NEXT STEP. C 470 KOLD = K HOLD = H C C CORRECT AND EVALUATE C TEMP1 = H*G(KP1) IF (NORND) GO TO 490 DO 480 L=1,NEQN RHO = TEMP1*(YP(L)-PHI(L,1)) - PHI(L,16) Y(L) = P(L) + RHO PHI(L,15) = (Y(L)-P(L)) - RHO 480 CONTINUE GO TO 510 490 DO 500 L=1,NEQN Y(L) = P(L) + TEMP1*(YP(L)-PHI(L,1)) 500 CONTINUE 510 CALL F(GFUN, FX, NEPS, NEQN, KA, BMA, X, Y, YP) C C UPDATE DIFFERENCES FOR NEXT STEP C DO 520 L=1,NEQN PHI(L,KP1) = YP(L) - PHI(L,1) PHI(L,KP2) = PHI(L,KP1) - PHI(L,KP2) 520 CONTINUE DO 540 I=1,K DO 530 L=1,NEQN PHI(L,I) = PHI(L,I) + PHI(L,KP1) 530 CONTINUE 540 CONTINUE C C ESTIMATE ERROR AT ORDER K+1 UNLESS: C IN FIRST PHASE WHEN ALWAYS RAISE ORDER, C ALREADY DECIDED TO LOWER ORDER, C STEP SIZE NOT CONSTANT SO ESTIMATE UNRELIABLE C ERKP1 = 0.0 IF (KNEW.EQ.KM1 .OR. K.EQ.12) PHASE1 = .FALSE. IF (PHASE1) GO TO 570 IF (KNEW.EQ.KM1) GO TO 580 IF (KP1.GT.NS) GO TO 590 DO 550 L=1,NEQN ERKP1 = ERKP1 + (PHI(L,KP2)/WT(L))**2 550 CONTINUE ERKP1 = DABSH*GSTR(KP1)*DSQRT(ERKP1) C C USING ESTIMATED ERROR AT ORDER K+1, DETERMINE APPROPRIATE ORDER C FOR NEXT STEP C IF (K.GT.1) GO TO 560 IF (ERKP1.GE.0.5D0*ERK) GO TO 590 GO TO 570 560 IF (ERKM1.LE.DMIN1(ERK,ERKP1)) GO TO 580 IF (ERKP1.GE.ERK .OR. K.EQ.12) GO TO 590 C C HERE ERKP1 .LT. ERK .LT. DMAX1(ERK1,ERK2) ELSE ORDER WOULD HAVE C BEEN LOWERED IN BLOCK 2. THUS ORDER IS TO BE RAISED C C RAISE ORDER C 570 K = KP1 ERK = ERKP1 GO TO 590 C C LOWER ORDER C 580 K = KM1 ERK = ERKM1 C C WITH NEW ORDER DETERMINE APPROPRIAAATE STEP SIZE FOR NEXT STEP C 590 HNEW = H + H IF (PHASE1) GO TO 600 IF (P5EPS.GE.ERK*TWO(K+1)) GO TO 600 HNEW = H IF (P5EPS.GE.ERK) GO TO 600 TEMP2 = K + 1 R = (P5EPS/ERK)**(1.0D0/TEMP2) HNEW = DABSH*DMAX1(0.5D0,DMIN1(0.9D0,R)) HNEW = DSIGN(DMAX1(HNEW,FOURU*DABS(X)),H) 600 H = HNEW RETURN C *** END BLOCK 4 *** C ******* END OF SUBROUTINE STEP ************************* END C ************ START OF SUBROUTINE INTRP (L. SHAMPINE, ET. AL.) ****** INT 10 C INT 20 SUBROUTINE INTRP(X, Y, XOUT, YOUT, YPOUT, NEQN, KOLD, PHI, PSI) INT 30 C C THE METHODS IN SUBROUTINE STEP APPROXIMATE THE SOLUTION NEAR X C BY A POLYNOMIAL. SUBROUTINE INTRP APPROXIMATES THE SOLUTION AT C XOUT BY EVALUATING THE POLYNOMIAL THERE. INFORMATION DEFINING THIS C POLYNOMIAL IS PASSED FROM STEP SO INTRP CINNOT BE USED ALONE. C C AS MENTIONED ELSEWHERE HERE, C THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT, C COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS: THE INITIAL C VALUE PROBLEM BY L. F. SHAMPINE AND M. K. GORDON. C C INPUT TO INTRP -- C C THE USER PROVIDES STORAGE IN THE CALLING PROGRAM FOR THE ARRAYS IN C THE CALL LIST C AND DEFINES C XOUT -- POINT AT WHICH SOLUTION IS DESIRED. C THE REMAINING PARAMETERS ARE DEFINED IN STEP AND PASSED TO INTRP C FROM THAT SUBROUTINE C C OUTPUT FROM INTRP -- C C YOUT(*) -- SOLUTION AT XOUT C YPOUT(*) -- DERIVATIVE OF SOLUTION AT XOUT C THE REMAINING PARAMETERS ARE RETURNED UNALTERED FROM THEIR INPUT C VALUES. INTEGRATION WITH STEP MAY BE CONTINUED. C C C RE THE FOLLOWING COMPILER STATEMENT, C SEE COMMENT NO.2 IN THE DRIVER. C C COMPILER STATIC C C C RE ARRAY SIZE. SEE COMMENT NO.3 IN THE DRIVER. C DOUBLE PRECISION Y(NEQN), YOUT(NEQN), YPOUT(NEQN), PHI(NEQN,16), * PSI(12) C DOUBLE PRECISION G(13), W(13), RHO(13) DOUBLE PRECISION X, HI, TEMP1, TEMP2, TEMP3, PSIJM1, GAMMA, ETA, * TERM, XOUT DATA G(1) /1.0D0/, RHO(1) /1.0D0/ C HI = XOUT - X KI = KOLD + 1 KIP1 = KI + 1 C C INITIALIZE W(*) FOR COMPUTING G(*) C DO 10 I=1,KI TEMP1 = I W(I) = 1.0D0/TEMP1 10 CONTINUE TERM = 0D0 C C COMPUTE G(*) C DO 30 J=2,KI JM1 = J - 1 PSIJM1 = PSI(JM1) GAMMA = (HI+TERM)/PSIJM1 ETA = HI/PSIJM1 LIMIT1 = KIP1 - J DO 20 I=1,LIMIT1 W(I) = GAMMA*W(I) - ETA*W(I+1) 20 CONTINUE G(J) = W(1) RHO(J) = GAMMA*RHO(JM1) TERM = PSIJM1 30 CONTINUE C C C INTERPOLATE C DO 40 L=1,NEQN YPOUT(L) = 0D0 YOUT(L) = 0D0 40 CONTINUE DO 60 J=1,KI I = KIP1 - J TEMP2 = G(I) TEMP3 = RHO(I) DO 50 L=1,NEQN YOUT(L) = YOUT(L) + TEMP2*PHI(L,I) YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I) 50 CONTINUE 60 CONTINUE DO 70 L=1,NEQN YOUT(L) = Y(L) + HI*YOUT(L) 70 CONTINUE RETURN C ******* END OF SUBROUTINE INTRP ************** END C *********** START OF SUBROUTINE SUM ******************************* SUM 10 C SUM 20 SUBROUTINE SUM(T, POS, NEG) SUM 30 C C PURPOSE : TO COMPUTE ALTERNATING SUMS BY SUMMING POSITIVE AND C NEGATIVE TERMS SEPARATELY. C DOUBLE PRECISION T, POS, NEG C IF (T.LE.0.0D0) GO TO 10 POS = POS + T GO TO 20 10 NEG = NEG + T 20 RETURN C ******* END OF SUBROUTINE SUM *************** END C ********** START OF SUBROUTINE COEFFS ***************************** COE 10 C COE 20 SUBROUTINE COEFFS(K, M1, BMA) COE 30 C C PURPOSE : TO COMPUTE THE FOURIER/CHEBYSHEV COEFFICIENTS FOR A C KERNEL OF THE FORM K(X-T) ON THE INTERVAL Õ0,BMAå. THE C COEFFICIENTS ARE STORED IN THE ARRAY C(51) AND ARE PASSED C TO OTHER NEEDING SUBPROGRAMS THROUGH THE LABELLED COMMON C BELOW. THIS SUBROUTINE IS EXECUTED IF AND ONLY IF THE C PARAMETER KA IS A POSITIVE INTEGER. C C PARAMETERS: C K FUNCTION SUBPROGRAM DEFINING THE KERNEL. C M1 NO. OF TERMS IN CHEBY. APPROXIMATION. C BMA LENGTH OF INTERVAL OVER WHICH SOLN IS REQUESTED. C C C RE THE FOLLOWING COMPILER STATEMENT. C SEE COMMENT NO.2 IN DRIVER. C C COMPILER STATIC C C DOUBLE PRECISION K, BMA, PI, RM1, BMAD2, TEMP1, TEMP2, RJ C C RE SIZE OF ARRAYS. SEE COMMENT NO.3. IN DRIVER. C DOUBLE PRECISION TERM, NEG, POS, RNU, Z, W, C(51), GJ(51) C COMMON /VE/ C, GJ, COUNTK, COUNTF C C NOTE THAT THE COUNTING VARIABLES COUNTK AND COUNTF ARE NOT DOUBLE C PRECISION. C C ALSO, INSTALLATIONS WITH A LARGER WORD SIZE MAY NEED MORE C PRECISION FOR PI. C DATA PI /3.14159265358979323D0/ C RM1 = M1 BMAD2 = BMA/2.0D00 C TEMP1 = K(BMA) TEMP2 = K(0.0D00) DO 50 J=1,M1 RJ = J TERM = (TEMP1+TEMP2*DCOS(PI*(RJ-1D0)))*5.0D-01 IF (TERM.GT.0.0D00) GO TO 10 NEG = TERM POS = 0.0D00 GO TO 20 10 POS = TERM NEG = 0.0D00 20 IF (M1.EQ.1) GO TO 40 DO 30 NU=2,M1 RNU = NU Z = BMAD2*(1.0D00+DCOS(PI*(RNU-1.0D00)/RM1)) W = DCOS((PI*(RJ-1.0D00)*(RNU-1.0D00))/RM1) TERM = K(Z)*W IF (TERM.GT.0.0D00) POS = POS + TERM IF (TERM.LT.0.0D00) NEG = NEG + TERM 30 CONTINUE C 40 C(J) = (2.0D00/RM1)*(POS+NEG) 50 CONTINUE C RETURN C ******* END OF SUBROUTINE COEFFS ******************* END C *********** START OF FUNCTION SOLN ******************************** SOL 10 C SOL 20 DOUBLE PRECISION FUNCTION SOLN(X) SOL 30 DOUBLE PRECISION X C C C THIS FUNCTION IS SUPPLIED STRICTLY FOR DEMONSTRATION PURPOSES FOR C USE WITH THE SUPPLIED DRIVER. OF COURSE THIS FUNCTION WOULD C NEVER HAVE ANY SERIOUS USE OTHER THAN DEMONSTRATION OR TESTING. C C FOR EXAMPLES IN ACM/TOMS : SOLN = X C RETURN C ******* END OF FUNCTION SOLN *****************+ END