C LBVP 10
C ..................................................................LBVP 20
C LBVP 30
C SUBROUTINE LBVP LBVP 40
C LBVP 50
C PURPOSE LBVP 60
C TO SOLVE A LINEAR BOUNDARY VALUE PROBLEM, WHICH CONSISTS OF LBVP 70
C A SYSTEM OF NDIM LINEAR FIRST ORDER DIFFERENTIAL EQUATIONS LBVP 80
C DY/DX=A(X)*Y(X)+F(X) LBVP 90
C AND NDIM LINEAR BOUNDARY CONDITIONS LBVP 100
C B*Y(XL)+C*Y(XU)=R. LBVP 110
C LBVP 120
C USAGE LBVP 130
C CALL LBVP (PRMT,B,C,R,Y,DERY,NDIM,IHLF,AFCT,FCT,DFCT,OUTP, LBVP 140
C AUX,A) LBVP 150
C PARAMETERS AFCT,FCT,DFCT,OUTP REQUIRE AN EXTERNAL STATEMENT.LBVP 160
C LBVP 170
C DESCRIPTION OF PARAMETERS LBVP 180
C PRMT - AN INPUT AND OUTPUT VECTOR WITH DIMENSION GREATER LBVP 190
C OR EQUAL TO 5, WHICH SPECIFIES THE PARAMETERS OF LBVP 200
C THE INTERVAL AND OF ACCURACY AND WHICH SERVES FOR LBVP 210
C COMMUNICATION BETWEEN OUTPUT SUBROUTINE (FURNISHED LBVP 220
C BY THE USER) AND SUBROUTINE LBVP. LBVP 230
C THE COMPONENTS ARE LBVP 240
C PRMT(1)- LOWER BOUND XL OF THE INTERVAL (INPUT), LBVP 250
C PRMT(1)- UPPER BOUND XU OF THE INTERVAL (INPUT), LBVP 260
C PRMT(3)- INITIAL INCREMENT OF THE INDEPENDENT VARIABLE LBVP 270
C (INPUT), LBVP 280
C PRMT(4)- UPPER ERROR BOUND (INPUT). IF RELATIVE ERROR IS LBVP 290
C GREATER THAN PRMT(4), INCREMENT GETS HALVED. LBVP 300
C IF INCREMENT IS LESS THAN PRMT(3) AND RELATIVE LBVP 310
C ERROR LESS THAN PRMT(4)/50, INCREMENT GETS DOUBLED.LBVP 320
C THE USER MAY CHANGE PRMT(4) BY MEANS OF HIS LBVP 330
C OUTPUT SUBROUTINE. LBVP 340
C PRMT(5)- NO INPUT PARAMETER. SUBROUTINE LBVP INITIALIZES LBVP 350
C PRMT(5)=0. IF THE USER WANTS TO TERMINATE LBVP 360
C SUBROUTINE LBVP AT ANY OUTPUT POINT, HE HAS TO LBVP 370
C CHANGE PRMT(5) TO NON-ZERO BY MEANS OF SUBROUTINE LBVP 380
C OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE LBVP 390
C FEASIBLE IF ITS DIMENSION IS DEFINED GREATER LBVP 400
C THAN 5. HOWEVER SUBROUTINE LBVP DOES NOT REQUIRE LBVP 410
C AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL LBVP 420
C FOR HANDING RESULT VALUES TO THE MAIN PROGRAM LBVP 430
C (CALLING LBVP) WHICH ARE OBTAINED BY SPECIAL LBVP 440
C MANIPULATIONS WITH OUTPUT DATA IN SUBROUTINE OUTP. LBVP 450
C B - AN NDIM BY NDIM INPUT MATRIX. (DESTROYED) LBVP 460
C IT IS THE COEFFICIENT MATRIX OF Y(XL) IN LBVP 470
C THE BOUNDARY CONDITIONS. LBVP 480
C C - AN NDIM BY NDIM INPUT MATRIX (POSSIBLY DESTROYED). LBVP 490
C IT IS THE COEFFICIENT MATRIX OF Y(XU) IN LBVP 500
C THE BOUNDARY CONDITIONS. LBVP 510
C R - AN INPUT VECTOR WITH DIMENSION NDIM. (DESTROYED) LBVP 520
C IT SPECIFIES THE RIGHT HAND SIDE OF THE LBVP 530
C BOUNDARY CONDITIONS. LBVP 540
C Y - AN AUXILIARY VECTOR WITH DIMENSION NDIM. LBVP 550
C IT IS USED AS STORAGE LOCATION FOR THE RESULTING LBVP 560
C VALUES OF DEPENDENT VARIABLES COMPUTED AT LBVP 570
C INTERMEDIATE POINTS. LBVP 580
C DERY - INPUT VECTOR OF ERROR WEIGHTS. (DESTROYED) LBVP 590
C ITS MAXIMAL COMPONENT SHOULD BE EQUAL TO 1. LBVP 600
C LATERON DERY IS THE VECTOR OF DERIVATIVES, WHICH LBVP 610
C BELONG TO FUNCTION VALUES Y AT INTERMEDIATE POINTS.LBVP 620
C NDIM - AN INPUT VALUE, WHICH SPECIFIES THE NUMBER OF LBVP 630
C DIFFERENTIAL EQUATIONS IN THE SYSTEM. LBVP 640
C IHLF - AN OUTPUT VALUE, WHICH SPECIFIES THE NUMBER OF LBVP 650
C BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS LBVP 660
C GREATER THAN 10, SUBROUTINE LBVP RETURNS WITH LBVP 670
C ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM. LBVP 680
C ERROR MESSAGE IHLF=12 OR IHLF=13 APPEARS IN CASE LBVP 690
C PRMT(3)=0 OR IN CASE SIGN(PRMT(3)).NE.SIGN(PRMT(2)-LBVP 700
C PRMT(1)) RESPECTIVELY. FINALLY ERROR MESSAGE LBVP 710
C IHLF=14 INDICATES, THAT THERE IS NO SOLUTION OR LBVP 720
C THAT THERE ARE MORE THAN ONE SOLUTION OF THE LBVP 730
C PROBLEM. LBVP 740
C A NEGATIVE VALUE OF IHLF HANDED TO SUBROUTINE OUTP LBVP 750
C TOGETHER WITH INITIAL VALUES OF FINALLY GENERATED LBVP 760
C INITIAL VALUE PROBLEM INDICATES, THAT THERE WAS LBVP 770
C POSSIBLE LOSS OF SIGNIFICANCE IN THE SOLUTION OF LBVP 780
C THE SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS FOR LBVP 790
C THESE INITIAL VALUES. THE ABSOLUTE VALUE OF IHLF LBVP 800
C SHOWS, AFTER WHICH ELIMINATION STEP OF GAUSS LBVP 810
C ALGORITHM POSSIBLE LOSS OF SIGNIFICANCE WAS LBVP 820
C DETECTED. LBVP 830
C AFCT - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT LBVP 840
C COMPUTES THE COEFFICIENT MATRIX A OF VECTOR Y ON LBVP 850
C THE RIGHT HAND SIDE OF THE SYSTEM OF DIFFERENTIAL LBVP 860
C EQUATIONS FOR A GIVEN X-VALUE. ITS PARAMETER LIST LBVP 870
C MUST BE X,A. SUBROUTINE AFCT SHOULD NOT DESTROY X. LBVP 880
C FCT - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT LBVP 890
C COMPUTES VECTOR F (INHOMOGENEOUS PART OF THE LBVP 900
C RIGHT HAND SIDE OF THE SYSTEM OF DIFFERENTIAL LBVP 910
C EQUATIONS) FOR A GIVEN X-VALUE. ITS PARAMETER LIST LBVP 920
C MUST BE X,F. SUBROUTINE FCT SHOULD NOT DESTROY X. LBVP 930
C DFCT - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT LBVP 940
C COMPUTES VECTOR DF (DERIVATIVE OF THE INHOMOGENEOUSLBVP 950
C PART ON THE RIGHT HAND SIDE OF THE SYSTEM OF LBVP 960
C DIFFERENTIAL EQUATIONS) FOR A GIVEN X-VALUE. ITS LBVP 970
C PARAMETER LIST MUST BE X,DF. SUBROUTINE DFCT LBVP 980
C SHOULD NOT DESTROY X. LBVP 990
C OUTP - THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED. LBVP1000
C ITS PARAMETER LIST MUST BE X,Y,DERY,IHLF,NDIM,PRMT.LBVP1010
C NONE OF THESE PARAMETERS (EXCEPT, IF NECESSARY, LBVP1020
C PRMT(4),PRMT(5),...) SHOULD BE CHANGED BY LBVP1030
C SUBROUTINE OUTP. IF PRMT(5) IS CHANGED TO NON-ZERO,LBVP1040
C SUBROUTINE LBVP IS TERMINATED. LBVP1050
C AUX - AN AUXILIARY STORAGE ARRAY WIRH 20 ROWS AND LBVP1060
C NDIM COLUMNS. LBVP1070
C A - AN NDIM BY NDIM MATRIX, WHICH IS USED AS AUXILIARY LBVP1080
C STORAGE ARRAY. LBVP1090
C LBVP1100
C REMARKS LBVP1110
C THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM, IF LBVP1120
C (1) MORE THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE LBVP1130
C NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE LBVP1140
C IHLF=11), LBVP1150
C (2) INITIAL INCREMENT IS EQUAL TO 0 OR IF IT HAS WRONG SIGN LBVP1160
C (ERROR MESSAGES IHLF=12 OR IHLF=13), LBVP1170
C (3) THERE IS NO OR MORE THAN ONE SOLUTION OF THE PROBLEM LBVP1180
C (ERROR MESSAGE IHLF=14), LBVP1190
C (4) THE WHOLE INTEGRATION INTERVAL IS WORKED THROUGH, LBVP1200
C (5) SUBROUTINE OUTP HAS CHANGED PRMT(5) TO NON-ZERO. LBVP1210
C LBVP1220
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED LBVP1230
C SUBROUTINE GELG SYSTEM OF LINEAR EQUATIONS. LBVP1240
C THE EXTERNAL SUBROUTINES AFCT(X,A), FCT(X,F), DFCT(X,DF), LBVP1250
C AND OUTP(X,Y,DERY,IHLF,NDIM,PRMT) MUST BE FURNISHED LBVP1260
C BY THE USER. LBVP1270
C LBVP1280
C METHOD LBVP1290
C EVALUATION IS DONE USING THE METHOD OF ADJOINT EQUATIONS. LBVP1300
C HAMMINGS FOURTH ORDER MODIFIED PREDICTOR-CORRECTOR METHOD LBVP1310
C IS USED TO SOLVE THE ADJOINT INITIAL VALUE PROBLEMS AND FI- LBVP1320
C NALLY TO SOLVE THE GENERATED INITIAL VALUE PROBLEM FOR Y(X).LBVP1330
C THE INITIAL INCREMENT PRMT(3) IS AUTOMATICALLY ADJUSTED. LBVP1340
C FOR COMPUTATION OF INTEGRAL SUM, A FOURTH ORDER HERMITEAN LBVP1350
C INTEGRATION FORMULA IS USED. LBVP1360
C FOR REFERENCE, SEE LBVP1370
C (1) LANCE, NUMERICAL METHODS FOR HIGH SPEED COMPUTERS, LBVP1380
C ILIFFE, LONDON, 1960, PP.64-67. LBVP1390
C (2) RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL LBVP1400
C COMPUTERS, WILEY, NEW YORK/LONDON, 1960, PP.95-109. LBVP1410
C (3) RALSTON, RUNGE-KUTTA METHODS WITH MINIMUM ERROR BOUNDS, LBVP1420
C MTAC, VOL.16, ISS.80 (1962), PP.431-437. LBVP1430
C (4) ZURMUEHL, PRAKTISCHE MATHEMATIK FUER INGENIEURE UND LBVP1440
C PHYSIKER, SPRINGER, BERLIN/GOETTINGEN/HEIDELBERG, 1963, LBVP1450
C PP.227-232. LBVP1460
C LBVP1470
C ..................................................................LBVP1480
C LBVP1490
0SUBROUTINE LBVP(PRMT,B,C,R,Y,DERY,NDIM,IHLF,AFCT,FCT,DFCT,OUTP, LBVP1500
1AUX,A) LBVP1510
C LBVP1520
DIMENSION PRMT(1),B(1),C(1),R(1),Y(1),DERY(1),AUX(20,1),A(1) LBVP1530
C LBVP1540
C ERROR TEST LBVP1550
IF(PRMT(3)*(PRMT(2)-PRMT(1)))2,1,3 LBVP1560
1 IHLF=12 LBVP1570
RETURN LBVP1580
2 IHLF=13 LBVP1590
RETURN LBVP1600
C LBVP1610
C SEARCH FOR ZERO-COLUMNS IN MATRICES B AND C LBVP1620
3 KK=-NDIM LBVP1630
IB=0 LBVP1640
IC=0 LBVP1650
DO 7 K=1,NDIM LBVP1660
AUX(15,K)=DERY(K) LBVP1670
AUX(1,K)=1. LBVP1680
AUX(17,K)=1. LBVP1690
KK=KK+NDIM LBVP1700
DO 4 I=1,NDIM LBVP1710
II=KK+I LBVP1720
IF(B(II))5,4,5 LBVP1730
4 CONTINUE LBVP1740
IB=IB+1 LBVP1750
AUX(1,K)=0. LBVP1760
5 DO 6 I=1,NDIM LBVP1770
II=KK+I LBVP1780
IF(C(II))7,6,7 LBVP1790
6 CONTINUE LBVP1800
IC=IC+1 LBVP1810
AUX(17,K)=0. LBVP1820
7 CONTINUE LBVP1830
C LBVP1840
C DETERMINATION OF LOWER AND UPPER BOUND LBVP1850
IF(IC-IB)8,11,11 LBVP1860
8 H=PRMT(2) LBVP1870
PRMT(2)=PRMT(1) LBVP1880
PRMT(1)=H LBVP1890
PRMT(3)=-PRMT(3) LBVP1900
DO 9 I=1,NDIM LBVP1910
9 AUX(17,I)=AUX(1,I) LBVP1920
II=NDIM*NDIM LBVP1930
DO 10 I=1,II LBVP1940
H=B(I) LBVP1950
B(I)=C(I) LBVP1960
10 C(I)=H LBVP1970
C LBVP1980
C PREPARATIONS FOR CONSTRUCTION OF ADJOINT INITIAL VALUE PROBLEMS LBVP1990
11 X=PRMT(2) LBVP2000
CALL FCT(X,Y) LBVP2010
CALL DFCT(X,DERY) LBVP2020
DO 12 I=1,NDIM LBVP2030
AUX(18,I)=Y(I) LBVP2040
12 AUX(19,I)=DERY(I) LBVP2050
C LBVP2060
C POSSIBLE BREAK-POINT FOR LINKAGE LBVP2070
C LBVP2080
C THE FOLLOWING PART OF SUBROUTINE LBVP UNTIL NEXT BREAK-POINT FOR LBVP2090
C LINKAGE HAS TO REMAIN IN CORE DURING THE WHOLE REST OF THE LBVP2100
C COMPUTATIONS LBVP2110
C LBVP2120
C START LOOP FOR GENERATING ADJOINT INITIAL VALUE PROBLEMS LBVP2130
K=0 LBVP2140
KK=0 LBVP2150
100 K=K+1 LBVP2160
IF(AUX(17,K))108,108,101 LBVP2170
C LBVP2180
C INITIALIZATION OF ADJOINT INITIAL VALUE PROBLEM LBVP2190
101 X=PRMT(2) LBVP2200
CALL AFCT(X,A) LBVP2210
SUM=0. LBVP2220
GL=AUX(18,K) LBVP2230
DGL=AUX(19,K) LBVP2240
II=K LBVP2250
DO 104 I=1,NDIM LBVP2260
H=-A(II) LBVP2270
DERY(I)=H LBVP2280
AUX(20,I)=R(I) LBVP2290
Y(I)=0. LBVP2300
IF(I-K)103,102,103 LBVP2310
102 Y(I)=1. LBVP2320
103 DGL=DGL+H*AUX(18,I) LBVP2330
104 II=II+NDIM LBVP2340
XEND=PRMT(1) LBVP2350
H=.0625*(XEND-X) LBVP2360
ISW=0 LBVP2370
GOTO 400 LBVP2380
C THIS IS BRANCH TO ADJOINT LINEAR INITIAL VALUE PROBLEM LBVP2390
C LBVP2400
C THIS IS RETURN FROM ADJOINT LINEAR INITIAL VALUE PROBLEM LBVP2410
105 IF(IHLF-10)106,106,117 LBVP2420
C LBVP2430
C UPDATING OF COEFFICIENT MATRIX B AND VECTOR R LBVP2440
106 DO 107 I=1,NDIM LBVP2450
KK=KK+1 LBVP2460
H=C(KK) LBVP2470
R(I)=AUX(20,I)+H*SUM LBVP2480
II=I LBVP2490
DO 107 J=1,NDIM LBVP2500
B(II)=B(II)+H*Y(J) LBVP2510
107 II=II+NDIM LBVP2520
GOTO 109 LBVP2530
108 KK=KK+NDIM LBVP2540
109 IF(K-NDIM)100,110,110 LBVP2550
C LBVP2560
C GENERATION OF LAST INITIAL VALUE PROBLEM LBVP2570
110 X=PRMT(4) LBVP2580
CALL GELG(R,B,NDIM,1,X,I) LBVP2590
IF(I)111,112,112 LBVP2600
111 IHLF=14 LBVP2610
RETURN LBVP2620
C LBVP2630
112 PRMT(5)=0. LBVP2640
IHLF=-I LBVP2650
X=PRMT(1) LBVP2660
XEND=PRMT(2) LBVP2670
H=PRMT(3) LBVP2680
DO 113 I=1,NDIM LBVP2690
113 Y(I)=R(I) LBVP2700
ISW=1 LBVP2710
114 ISW2=12 LBVP2720
GOTO 200 LBVP2730
115 ISW3=-1 LBVP2740
GOTO 300 LBVP2750
116 IF(IHLF)400,400,117 LBVP2760
C THIS WAS BRANCH INTO INITIAL VALUE PROBLEM LBVP2770
C LBVP2780
C THIS IS RETURN FROM INITIAL VALUE PROBLEM LBVP2790
117 RETURN LBVP2800
C LBVP2810
C THIS PART OF LINEAR BOUNDARY VALUE PROBLEM COMPUTES THE RIGHT LBVP2820
C HAND SIDE DERY OF THE SYSTEM OF ADJOINT LINEAR DIFFERENTIAL LBVP2830
C EQUATIONS (IN CASE ISW=0) OR OF THE GIVEN SYSTEM (IN CASE ISW=1). LBVP2840
200 CALL AFCT(X,A) LBVP2850
IF(ISW)201,201,205 LBVP2860
C LBVP2870
C ADJOINT SYSTEM LBVP2880
201 LL=0 LBVP2890
DO 203 M=1,NDIM LBVP2900
HS=0. LBVP2910
DO 202 L=1,NDIM LBVP2920
LL=LL+1 LBVP2930
202 HS=HS-A(LL)*Y(L) LBVP2940
203 DERY(M)=HS LBVP2950
204 GOTO(502,504,506,407,415,418,608,617,632,634,421,115),ISW2 LBVP2960
C LBVP2970
C GIVEN SYSTEM LBVP2980
205 CALL FCT(X,DERY) LBVP2990
DO 207 M=1,NDIM LBVP3000
LL=M-NDIM LBVP3010
HS=0. LBVP3020
DO 206 L=1,NDIM LBVP3030
LL=LL+NDIM LBVP3040
206 HS=HS+A(LL)*Y(L) LBVP3050
207 DERY(M)=HS+DERY(M) LBVP3060
GOTO 204 LBVP3070
C LBVP3080
C THIS PART OF LINEAR BOUNDARY VALUE PROBLEM COMPUTES THE VALUE OF LBVP3090
C INTEGRAL SUM, WHICH IS A PART OF THE OUTPUT OF ADJOINT INITIAL LBVP3100
C VALUE PROBLEM (IN CASE ISW=0) OR RECORDS RESULT VALUES OF THE LBVP3110
C FINAL INITIAL VALUE PROBLEM (IN CASE ISW=1). LBVP3120
300 IF(ISW)301,301,305 LBVP3130
C LBVP3140
C ADJOINT PROBLEM LBVP3150
301 CALL FCT(X,R) LBVP3160
GU=0. LBVP3170
DGU=0. LBVP3180
DO 302 L=1,NDIM LBVP3190
GU=GU+Y(L)*R(L) LBVP3200
302 DGU=DGU+DERY(L)*R(L) LBVP3210
CALL DFCT(X,R) LBVP3220
DO 303 L=1,NDIM LBVP3230
303 DGU=DGU+Y(L)*R(L) LBVP3240
SUM=SUM+.5*H*((GL+GU)+.1666667*H*(DGL-DGU)) LBVP3250
GL=GU LBVP3260
DGL=DGU LBVP3270
304 IF(ISW3)116,422,618 LBVP3280
C LBVP3290
C GIVEN PROBLEM LBVP3300
305 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) LBVP3310
IF(PRMT(5))117,304,117 LBVP3320
C LBVP3330
C POSSIBLE BREAK-POINT FOR LINKAGE LBVP3340
C LBVP3350
C THE FOLLOWING PART OF SUBROUTINE LBVP SOLVES IN CASE ISW=0 THE LBVP3360
C ADJOINT INITIAL VALUE PROBLEM. IT COMPUTES INTEGRAL SUM AND LBVP3370
C THE VECTOR Y OF DEPENDENT VARIABLES AT THE LOWER BOUND PRMT(1). LBVP3380
C IN CASE ISW=1 IT SOLVES FINALLY GENERATED INITIAL VALUE PROBLEM. LBVP3390
400 N=1 LBVP3400
XST=X LBVP3410
IHLF=0 LBVP3420
DO 401 I=1,NDIM LBVP3430
AUX(16,I)=0. LBVP3440
AUX(1,I)=Y(I) LBVP3450
401 AUX(8,I)=DERY(I) LBVP3460
ISW1=1 LBVP3470
GOTO 500 LBVP3480
C LBVP3490
402 X=X+H LBVP3500
DO 403 I=1,NDIM LBVP3510
403 AUX(2,I)=Y(I) LBVP3520
C LBVP3530
C INCREMENT H IS TESTED BY MEANS OF BISECTION LBVP3540
404 IHLF=IHLF+1 LBVP3550
X=X-H LBVP3560
DO 405 I=1,NDIM LBVP3570
405 AUX(4,I)=AUX(2,I) LBVP3580
H=.5*H LBVP3590
N=1 LBVP3600
ISW1=2 LBVP3610
GOTO 500 LBVP3620
C LBVP3630
406 X=X+H LBVP3640
ISW2=4 LBVP3650
GOTO 200 LBVP3660
407 N=2 LBVP3670
DO 408 I=1,NDIM LBVP3680
AUX(2,I)=Y(I) LBVP3690
408 AUX(9,I)=DERY(I) LBVP3700
ISW1=3 LBVP3710
GOTO 500 LBVP3720
C LBVP3730
C TEST ON SATISFACTORY ACCURACY LBVP3740
409 DO 414 I=1,NDIM LBVP3750
Z=ABS(Y(I)) LBVP3760
IF(Z-1.)410,411,411 LBVP3770
410 Z=1. LBVP3780
411 DELT=.06666667*ABS(Y(I)-AUX(4,I)) LBVP3790
IF(ISW)413,413,412 LBVP3800
412 DELT=AUX(15,I)*DELT LBVP3810
413 IF(DELT-Z*PRMT(4))414,414,429 LBVP3820
414 CONTINUE LBVP3830
C LBVP3840
C SATISFACTORY ACCURACY AFTER LESS THAN 11 BISECTIONS LBVP3850
X=X+H LBVP3860
ISW2=5 LBVP3870
GOTO 200 LBVP3880
415 DO 416 I=1,NDIM LBVP3890
AUX(3,I)=Y(I) LBVP3900
416 AUX(10,I)=DERY(I) LBVP3910
N=3 LBVP3920
ISW1=4 LBVP3930
GOTO 500 LBVP3940
C LBVP3950
417 N=1 LBVP3960
X=X+H LBVP3970
ISW2=6 LBVP3980
GOTO 200 LBVP3990
418 X=XST LBVP4000
DO 419 I=1,NDIM LBVP4010
AUX(11,I)=DERY(I) LBVP4020
4190Y(I)=AUX(1,I)+H*(.375*AUX(8,I)+.7916667*AUX(9,I) LBVP4030
1-.2083333*AUX(10,I)+.04166667*DERY(I)) LBVP4040
420 X=X+H LBVP4050
N=N+1 LBVP4060
ISW2=11 LBVP4070
GOTO 200 LBVP4080
421 ISW3=0 LBVP4090
GOTO 300 LBVP4100
422 IF(N-4)423,600,600 LBVP4110
423 DO 424 I=1,NDIM LBVP4120
AUX(N,I)=Y(I) LBVP4130
424 AUX(N+7,I)=DERY(I) LBVP4140
IF(N-3)425,427,600 LBVP4150
C LBVP4160
425 DO 426 I=1,NDIM LBVP4170
DELT=AUX(9,I)+AUX(9,I) LBVP4180
DELT=DELT+DELT LBVP4190
426 Y(I)=AUX(1,I)+.3333333*H*(AUX(8,I)+DELT+AUX(10,I)) LBVP4200
GOTO 420 LBVP4210
C LBVP4220
427 DO 428 I=1,NDIM LBVP4230
DELT=AUX(9,I)+AUX(10,I) LBVP4240
DELT=DELT+DELT+DELT LBVP4250
428 Y(I)=AUX(1,I)+.375*H*(AUX(8,I)+DELT+AUX(11,I)) LBVP4260
GOTO 420 LBVP4270
C LBVP4280
C NO SATISFACTORY ACCURACY. H MUST BE HALVED. LBVP4290
429 IF(IHLF-10)404,430,430 LBVP4300
C LBVP4310
C NO SATISFACTORY ACCURACY AFTER 10 BISECTIONS. ERROR MESSAGE. LBVP4320
430 IHLF=11 LBVP4330
X=X+H LBVP4340
IF(ISW)105,105,114 LBVP4350
C LBVP4360
C THIS PART OF LINEAR INITIAL VALUE PROBLEM COMPUTES LBVP4370
C STARTING VALUES BY MEANS OF RUNGE-KUTTA METHOD. LBVP4380
500 Z=X LBVP4390
DO 501 I=1,NDIM LBVP4400
X=H*AUX(N+7,I) LBVP4410
AUX(5,I)=X LBVP4420
501 Y(I)=AUX(N,I)+.4*X LBVP4430
C LBVP4440
X=Z+.4*H LBVP4450
ISW2=1 LBVP4460
GOTO 200 LBVP4470
502 DO 503 I=1,NDIM LBVP4480
X=H*DERY(I) LBVP4490
AUX(6,I)=X LBVP4500
503 Y(I)=AUX(N,I)+.2969776*AUX(5,I)+.1587596*X LBVP4510
C LBVP4520
X=Z+.4557372*H LBVP4530
ISW2=2 LBVP4540
GOTO 200 LBVP4550
504 DO 505 I=1,NDIM LBVP4560
X=H*DERY(I) LBVP4570
AUX(7,I)=X LBVP4580
505 Y(I)=AUX(N,I)+.2181004*AUX(5,I)-3.050965*AUX(6,I)+3.832865*X LBVP4590
C LBVP4600
X=Z+H LBVP4610
ISW2=3 LBVP4620
GOTO 200 LBVP4630
506 DO 507 I=1,NDIM LBVP4640
5070Y(I)=AUX(N,I)+.1747603*AUX(5,I)-.5514807*AUX(6,I) LBVP4650
1+1.205536*AUX(7,I)+.1711848*H*DERY(I) LBVP4660
X=Z LBVP4670
GOTO(402,406,409,417),ISW1 LBVP4680
C LBVP4690
C POSSIBLE BREAK-POINT FOR LINKAGE LBVP4700
C LBVP4710
C STARTING VALUES ARE COMPUTED. LBVP4720
C NOW START HAMMINGS MODIFIED PREDICTOR-CORRECTOR METHOD. LBVP4730
600 ISTEP=3 LBVP4740
601 IF(N-8)604,602,604 LBVP4750
C LBVP4760
C N=8 CAUSES THE ROWS OF AUX TO CHANGE THEIR STORAGE LOCATIONS LBVP4770
602 DO 603 N=2,7 LBVP4780
DO 603 I=1,NDIM LBVP4790
AUX(N-1,I)=AUX(N,I) LBVP4800
603 AUX(N+6,I)=AUX(N+7,I) LBVP4810
N=7 LBVP4820
C LBVP4830
C N LESS THAN 8 CAUSES N+1 TO GET N LBVP4840
604 N=N+1 LBVP4850
C LBVP4860
C COMPUTATION OF NEXT VECTOR Y LBVP4870
DO 605 I=1,NDIM LBVP4880
AUX(N-1,I)=Y(I) LBVP4890
605 AUX(N+6,I)=DERY(I) LBVP4900
X=X+H LBVP4910
606 ISTEP=ISTEP+1 LBVP4920
DO 607 I=1,NDIM LBVP4930
0DELT=AUX(N-4,I)+1.333333*H*(AUX(N+6,I)+AUX(N+6,I)-AUX(N+5,I)+ LBVP4940
1AUX(N+4,I)+AUX(N+4,I)) LBVP4950
Y(I)=DELT-.9256198*AUX(16,I) LBVP4960
607 AUX(16,I)=DELT LBVP4970
C PREDICTOR IS NOW GENERATED IN ROW 16 OF AUX, MODIFIED PREDICTOR LBVP4980
C IS GENERATED IN Y. DELT MEANS AN AUXILIARY STORAGE. LBVP4990
C LBVP5000
ISW2=7 LBVP5010
GOTO 200 LBVP5020
C DERIVATIVE OF MODIFIED PREDICTOR IS GENERATED IN DERY. LBVP5030
C LBVP5040
608 DO 609 I=1,NDIM LBVP5050
0DELT=.125*(9.*AUX(N-1,I)-AUX(N-3,I)+3.*H*(DERY(I)+AUX(N+6,I)+ LBVP5060
1AUX(N+6,I)-AUX(N+5,I))) LBVP5070
AUX(16,I)=AUX(16,I)-DELT LBVP5080
609 Y(I)=DELT+.07438017*AUX(16,I) LBVP5090
C LBVP5100
C TEST WHETHER H MUST BE HALVED OR DOUBLED LBVP5110
DELT=0. LBVP5120
DO 616 I=1,NDIM LBVP5130
Z=ABS(Y(I)) LBVP5140
IF(Z-1.)610,611,611 LBVP5150
610 Z=1. LBVP5160
611 Z=ABS(AUX(16,I))/Z LBVP5170
IF(ISW)613,613,612 LBVP5180
612 Z=AUX(15,I)*Z LBVP5190
613 IF(Z-PRMT(4))614,614,628 LBVP5200
614 IF(DELT-Z)615,616,616 LBVP5210
615 DELT=Z LBVP5220
616 CONTINUE LBVP5230
C LBVP5240
C H MUST NOT BE HALVED. THAT MEANS Y(I) ARE GOOD. LBVP5250
ISW2=8 LBVP5260
GOTO 200 LBVP5270
617 ISW3=1 LBVP5280
GOTO 300 LBVP5290
618 IF(H*(X-XEND))619,621,621 LBVP5300
619 IF(ABS(X-XEND)-.1*ABS(H))621,620,620 LBVP5310
620 IF(DELT-.02*PRMT(4))622,622,601 LBVP5320
621 IF(ISW)105,105,117 LBVP5330
C LBVP5340
C LBVP5350
C H COULD BE DOUBLED IF ALL NECESSARY PRECEEDING VALUES ARE LBVP5360
C AVAILABLE. LBVP5370
622 IF(IHLF)601,601,623 LBVP5380
623 IF(N-7)601,624,624 LBVP5390
624 IF(ISTEP-4)601,625,625 LBVP5400
625 IMOD=ISTEP/2 LBVP5410
IF(ISTEP-IMOD-IMOD)601,626,601 LBVP5420
626 H=H+H LBVP5430
IHLF=IHLF-1 LBVP5440
ISTEP=0 LBVP5450
DO 627 I=1,NDIM LBVP5460
AUX(N-1,I)=AUX(N-2,I) LBVP5470
AUX(N-2,I)=AUX(N-4,I) LBVP5480
AUX(N-3,I)=AUX(N-6,I) LBVP5490
AUX(N+6,I)=AUX(N+5,I) LBVP5500
AUX(N+5,I)=AUX(N+3,I) LBVP5510
AUX(N+4,I)=AUX(N+1,I) LBVP5520
DELT=AUX(N+6,I)+AUX(N+5,I) LBVP5530
DELT=DELT+DELT+DELT LBVP5540
6270AUX(16,I)=8.962963*(Y(I)-AUX(N-3,I))-3.361111*H*(DERY(I)+DELT LBVP5550
1+AUX(N+4,I)) LBVP5560
GOTO 601 LBVP5570
C LBVP5580
C LBVP5590
C H MUST BE HALVED LBVP5600
628 IHLF=IHLF+1 LBVP5610
IF(IHLF-10)630,630,629 LBVP5620
629 IF(ISW)105,105,114 LBVP5630
630 H=.5*H LBVP5640
ISTEP=0 LBVP5650
DO 631 I=1,NDIM LBVP5660
0Y(I)=.00390625*(80.*AUX(N-1,I)+135.*AUX(N-2,I)+40.*AUX(N-3,I)+ LBVP5670
1AUX(N-4,I))-.1171875*(AUX(N+6,I)-6.*AUX(N+5,I)-AUX(N+4,I))*H LBVP5680
0AUX(N-4,I)=.00390625*(12.*AUX(N-1,I)+135.*AUX(N-2,I)+ LBVP5690
1108.*AUX(N-3,I)+AUX(N-4,I))-.0234375*(AUX(N+6,I)+18.*AUX(N+5,I)- LBVP5700
29.*AUX(N+4,I))*H LBVP5710
AUX(N-3,I)=AUX(N-2,I) LBVP5720
631 AUX(N+4,I)=AUX(N+5,I) LBVP5730
DELT=X-H LBVP5740
X=DELT-(H+H) LBVP5750
ISW2=9 LBVP5760
GOTO 200 LBVP5770
632 DO 633 I=1,NDIM LBVP5780
AUX(N-2,I)=Y(I) LBVP5790
AUX(N+5,I)=DERY(I) LBVP5800
633 Y(I)=AUX(N-4,I) LBVP5810
X=X-(H+H) LBVP5820
ISW2=10 LBVP5830
GOTO 200 LBVP5840
634 X=DELT LBVP5850
DO 635 I=1,NDIM LBVP5860
DELT=AUX(N+5,I)+AUX(N+4,I) LBVP5870
DELT=DELT+DELT+DELT LBVP5880
0AUX(16,I)=8.962963*(AUX(N-1,I)-Y(I))-3.361111*H*(AUX(N+6,I)+DELT LBVP5890
1+DERY(I)) LBVP5900
635 AUX(N+3,I)=DERY(I) LBVP5910
GOTO 606 LBVP5920
C LBVP5930
C END OF INITIAL VALUE PROBLEM LBVP5940
END LBVP5950