1SECTION : 5.2.1.1.1.2.C (OCTOBER 1974) PAGE 1
AUTHOR: K.DEKKER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 1973/09/01.
BRIEF DESCRIPTION:
LINIGER1VS SOLVES INITIAL VALUE PROBLEMS , GIVEN AS AN AUTONOMOUS
SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS , BY MEANS OF AN
IMPLICIT,FIRST ORDER ACCURATE, EXPONENTIALLY FITTED ONESTEP METHOD.
AUTOMATIC STEPSIZE CONTROL IS PROVIDED.
IN PARTICULAR THIS METHOD IS SUITABLE FOR THE INTEGRATION OF STIFF
DIFFERENTIAL EQUATIONS.
KEYWORDS:
DIFFERENTIAL EQUATIONS,
INITIAL VALUE PROBLEMS,
STIFF EQUATIONS,
EXPONENTIAL FITTING,
IMPLICIT ONESTEP METHODS.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE LINIGER1VS READS:
"PROCEDURE" LINIGER1VS (X,XE,M,Y,SIGMA,DERIVATIVE,J,JACOBIAN,
ITMAX,HMIN,HMAX,AETA,RETA,INFO,OUTPUT);
"VALUE" M;
"INTEGER" M,ITMAX;
"REAL" X,XE,SIGMA,HMIN,HMAX,AETA,RETA;
"ARRAY" Y,J,INFO;
"PROCEDURE" DERIVATIVE,JACOBIAN,OUTPUT;
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <VARIABLE>;
THE INDEPENDENT VARIABLE X;
ENTRY: THE INITIAL VALUE X0;
EXIT : THE FINAL VALUE XE;
XE: <ARITHMETIC EXPRESSION>;
THE FINAL VALUE OF X (XE>=X);
M: <ARITHMETIC EXPRESSION>;
THE NUMBER OF EQUATIONS;
Y: <ARRAY IDENTIFIER>;
"ARRAY" Y[1:M];
THE DEPENDENT VARIABLE;
ENTRY: THE INITIAL VALUES OF THE SYSTEM OF DIFFERENTIAL
EQUATIONS: Y[I] AT X=X0;
EXIT : THE FINAL VALUES OF THE SOLUTION: Y[I] AT X=XE;
1SECTION : 5.2.1.1.1.2.C (OCTOBER 1974) PAGE 2
SIGMA: <ARITHMETIC EXPRESSION>;
THE MODULUS OF THE POINT AT WHICH EXPONENTIAL FITTING IS
DESIRED, FOR EXAMPLE THE LARGEST NEGATIVE EIGENVALUE OF THE
JACOBIAN OF THE SYSTEM OF DIFFERENTIAL EQUATIONS;
DERIVATIVE: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" DERIVATIVE(Y); "ARRAY" Y;
THIS PROCEDURE SHOULD DELIVER THE RIGHT HAND SIDE OF THE
I-TH DIFFERENTIAL EQUATION AT THE POINT (Y) AS Y[I];
J: <ARRAY IDENTIFIER>;
"ARRAY" J[1:M,1:M];
THE JACOBIAN MATRIX OF THE SYSTEM;
THE ARRAY J SHOULD BE UPDATED IN THE PROCEDURE JACOBIAN;
JACOBIAN: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:
"PROCEDURE" JACOBIAN(J,Y); "ARRAY" J,Y;
IN THIS PROCEDURE (AN APPROXIMATION OF) THE JACOBIAN HAS TO
BE ASSIGNED TO THE ARRAY J;
ITMAX: <ARITHMETIC EXPRESSION>;
AN UPPERBOUND FOR THE NUMBER OF ITERATIONS IN NEWTON'S
PROCESS, USED TO SOLVE THE IMPLICIT EQUATIONS;
HMIN: <ARITHMETIC EXPRESSION>;
MINIMAL STEPSIZE BY WHICH THE INTEGRATION IS PERFORMED;
HMAX: <ARITHMETIC EXPRESSION>;
MAXIMAL STEPSIZE BY WHICH THE INTEGRATION IS PERFORMED;
AETA: <ARITHMETIC EXPRESSION>;
REQUIRED ABSOLUTE PRECISION IN THE INTEGRATION PROCESS;
RETA: <ARITHMETIC EXPRESSION>;
REQUIRED RELATIVE PRECISION IN THE INTEGRATION PROCESS;
IF BOTH AETA AND RETA HAVE NEGATIVE VALUES , INTEGRATION
WILL BE PERFORMED WITH A STEPSIZE EQUAL TO HMAX , WHICH MAY
BE VARIATED BY USER ; IN THIS CASE THE ABSOLUTE VALUES OF
AETA AND RETA WILL CONTROL THE NEWTON ITERATION;
INFO: <ARRAY IDENTIFIER>;
"ARRAY" INFO[1:9];
DURING INTEGRATION AND UPON EXIT THIS ARRAY CONTAINS THE
FOLLOWING INFORMATION:
INFO[1]: NUMBER OF INTEGRATION STEPS TAKEN;
INFO[2]: NUMBER OF DERIVATIVE EVALUATIONS USED;
INFO[3]: NUMBER OF JACOBIAN EVALUATIONS USED;
INFO[4]: NUMBER OF INTEGRATION STEPS EQUAL TO HMIN TAKEN ;
INFO[5]: NUMBER OF INTEGRATION STEPS EQUAL TO HMAX TAKEN ;
INFO[6]: MAXIMAL NUMBER OF ITERATIONS TAKEN IN THE NEWTON
PROCESS;
INFO[7]: LOCAL ERROR TOLERANCE;
INFO[8]: ESTIMATED LOCAL ERROR;
INFO[9]: MAXIMUM VALUE OF THE ESTIMATED LOCAL ERROR;
OUTPUT: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS;
"PROCEDURE" OUTPUT;
THIS PROCEDURE IS CALLED AT THE END OF EACH INTEGRATION
STEP ; THE USER CAN ASK FOR OUTPUT OF THE PARAMETERS , FOR
EXAMPLE X, Y, INFO.
1SECTION : 5.2.1.1.1.2.C (OCTOBER 1974) PAGE 3
DATA AND RESULTS: SEE EXAMPLE OF USE, AND REF[2].
PROCEDURES USED:
INIVEC= CP31010,
MULVEC= CP31020,
MULROW= CP31021,
DUPVEC= CP31030,
MATVEC= CP34011,
ELMVEC= CP34020,
VECVEC= CP34010,
DEC = CP34300,
SOL = CP34051.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH : CIRCA 20 + M * (5 + M).
RUNNING TIME: DEPENDS STRONGLY ON THE DIFFERENTIAL EQUATION TO SOLVE.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
LINIGER1VS: INTEGRATES THE SYSTEM OF DIFFERENTIAL EQUATIONS FROM X0
UNTIL XE, BY MEANS OF A FIRST ORDER FORMULA.
THE INTEGRATION METHOD IS BASED ON THE F(1) FORMULA DESCRIBED BY
LINIGER AND WILLOUGHBY (SEE REF[1]). ERROR ESTIMATES AND A STEPSIZE
STRATEGY FOR THIS METHOD ARE DESCRIBED IN [2] , AND A VARIABLE STEP
METHOD IS CONSTRUCTED FOR THE CONVENIENCE OF THE USER. HOWEVER, THE
STEPSIZE STRATEGY REQUIRES MANY EXTRA ARRAY OPERATIONS. THE USER
MAY AVOID THIS EXTRA WORK BY GIVING AETA AND RETA A NEGATIVE VALUE
AND PRESCRIBING A STEPSIZE (HMAX) HIMSELF.
REFERENCES:
[1]. W.LINIGER AND R.A.WILLOUGHBY.
EFFICIENT INTEGRATION METHODS FOR STIFF SYSTEMS OF ORDINARY
DIFFERENTIAL EQUATIONS.
SIAM J. NUM. ANAL. 7 (1970) PAGE 47.
[2]. K.DEKKER.
ERROR ESTIMATES AND STEPSIZE STRATEGIES FOR THE LINIGER-
WILLOUGHBY FORMULAE.
(TO APPEAR IN 1974).
1SECTION : 5.2.1.1.1.2.C (OCTOBER 1974) PAGE 4
EXAMPLE OF USE:
CONSIDER THE SYSTEM OF DIFFERENTIAL EQUATIONS:
DY[1]/DX = -Y[1] + Y[1] * Y[2] + .99 * Y[2]
DY[2]/DX = -1000 * ( -Y[1] + Y[1] * Y[2] + Y[2] )
WITH THE INITIAL CONDITIONS AT X = 0:
Y[1] = 1 AND Y[2] = 0.
THE SOLUTION AT X = 50 IS APPROXIMATELY:
Y[1] = .765 878 320 2487 AND Y[2] = .433 710 353 5768.
THE FOLLOWING PROGRAM SHOWS INTEGRATION OF THIS PROBLEM WITH
VARIABLE AND CONSTANT STEPSIZES:
"BEGIN" "COMMENT" TEST LINIGER1VS;
"PROCEDURE" LINIGER1VS(X,XE,M,Y,SIGMA,F,J,JACOBIAN,
ITMAX,HMIN,HMAX,AETA,RETA,INFO,OUTPUT);
"CODE" 33132;
"INTEGER" ITMAX;
"REAL" X,SIGMA,RETA,TIME;
"REAL" "ARRAY" Y[1:2],J[1:2,1:2],INFO[1:9];
"PROCEDURE" F(A); "ARRAY" A;
"BEGIN" "REAL" A1,A2; A1:=A[1]; A2:=A[2];
A[1]:=(A1+.99)*(A2-1)+.99;
A[2]:=1000*((1+A1)*(1-A2)-1);
"END";
"PROCEDURE" JACOBIAN(J,Y); "ARRAY" J,Y;
"BEGIN" J[1,1]:=Y[2]-1; J[1,2]:=.99+Y[1];
J[2,1]:=1000*(1-Y[2]); J[2,2]:=-1000*(1+Y[1]);
SIGMA:=ABS(J[2,2]+J[1,1]-SQRT((J[2,2]-J[1,1])**2+
4*J[2,1]*J[1,2]))/2;
"END" JACOBIAN;
"PROCEDURE" OUT;
"IF" X=50 "THEN"
OUTPUT(61,"("6(3ZDB),2BD"-ZD,2(2B+.3DB3D),-3ZD.3D,/")",
INFO[1],INFO[2],INFO[3],INFO[4],INFO[5],INFO[6],INFO[9],Y[1],
Y[2],CLOCK-TIME);
"FOR" RETA:="-2,"-4,"-6 "DO"
"BEGIN" X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;
LINIGER1VS(X,50,2,Y,SIGMA,F,J,JACOBIAN,10,.1,50,RETA,
RETA,INFO,OUT);
"END"; OUTPUT(61,"("//")");
"FOR" RETA:=-"-2,-"-4,-"-6 "DO"
"BEGIN" X:=Y[2]:=0; Y[1]:=1; TIME:=CLOCK;
LINIGER1VS(X,50,2,Y,SIGMA,F,J,JACOBIAN,10,.1,1,RETA,
RETA,INFO,OUT);
"END";
"END"
17 21 8 2 0 2 2" -2 +.772 017 +.435 672 0.525
13 25 23 2 0 2 2" -2 +.767 414 +.434 202 0.717
105 210 105 2 0 2 2" -2 +.766 027 +.433 758 4.741
50 52 1 0 50 2 0" 0 +.766 670 +.433 081 0.549
50 104 3 0 50 3 0" 0 +.766 183 +.433 811 1.158
50 152 12 0 50 4 0" 0 +.766 185 +.433 809 1.653
1SECTION : 5.2.1.1.1.2.C (OCTOBER 1974) PAGE 5
SOURCE TEXT(S):
"CODE" 33132;
"PROCEDURE" LINIGER1VS(X,XE,M,Y,SIGMA,DERIVATIVE,J,
JACOBIAN,ITMAX,HMIN,HMAX,AETA,RETA,INFO,OUTPUT);
"VALUE" M;
"INTEGER" M,ITMAX;
"REAL" X,XE,SIGMA,AETA,RETA,HMIN,HMAX;
"ARRAY" Y,J,INFO;
"PROCEDURE" DERIVATIVE,JACOBIAN,OUTPUT;
"BEGIN" "INTEGER" I,ST,LASTJAC;
"REAL" H,HNEW,MU,MU1,BETA,P,E,E1,ETA,ETA1,DISCR;
"BOOLEAN" LAST,FIRST,EVALJAC,EVALCOEF;
"INTEGER" "ARRAY" PI[1:M];
"REAL" "ARRAY" DY,YL,YR,F[1:M],A[1:M,1:M],AUX[1:3];
"PROCEDURE" INIVEC(L,U,A,X); "CODE" 31010;
"PROCEDURE" MULVEC(AL,U,SHIFT,A,B,X); "CODE" 31020;
"PROCEDURE" MULROW(L,U,I,J,A,B,X); "CODE" 31021;
"PROCEDURE" DUPVEC(L,U,SHIFT,A,B); "CODE" 31030;
"REAL" "PROCEDURE" VECVEC(L,U,SHIFT,A,B); "CODE" 34010;
"REAL" "PROCEDURE" MATVEC(A,B,C,D,E); "CODE" 34011;
"PROCEDURE" ELMVEC(L,U,SHIFT,A,B,X); "CODE" 34020;
"PROCEDURE" DEC(A,N,AUX,P); "CODE" 34300;
"PROCEDURE" SOL(A,N,P,B); "CODE" 34051;
"REAL" "PROCEDURE" NORM(A); "ARRAY" A;
NORM:=SQRT(VECVEC(1,M,0,A,A));
"PROCEDURE" COEFFICIENT;
"BEGIN" "REAL" B,E; B:=ABS(H*SIGMA);
"IF" B>40 "THEN"
"BEGIN" MU:=1/B; BETA:=1; P:=2+2/(B-2)
"END" "ELSE"
"IF" B<.04 "THEN"
"BEGIN" E:=B*B/30; P:=3-E;
MU:=.5-B/12*(1-E/2);
BETA:=.5+B/6*(1-E)
"END" "ELSE"
"BEGIN" E:=EXP(B)-1;
MU:=1/B-1/E;
BETA:=(1-B/E)*(1+1/E);
P:=(BETA-MU)/(.5-MU)
"END";
MU1:=H*(1-MU);
LUDECOMP
"END" COEFFICIENT; "COMMENT"
1SECTION : 5.2.1.1.1.2.C (OCTOBER 1974) PAGE 6
;
"PROCEDURE" LUDECOMP;
"BEGIN" "INTEGER" I;
"FOR" I:=1 "STEP" 1 "UNTIL" M "DO"
"BEGIN" MULROW(1,M,I,I,A,J,-MU1);
A[I,I]:=A[I,I]+1
"END";
AUX[2]:=0; DEC(A,M,AUX,PI)
"END" LUDECOMP;
"PROCEDURE" STEPSIZE;
"IF" ETA<0 "THEN"
"BEGIN" "REAL" HL; HL:=H;
H:=HNEW:=HMAX; INFO[5]:=INFO[5]+1;
"IF" 1.1*HNEW>XE-X "THEN"
"BEGIN" LAST:="TRUE"; H:=HNEW:=XE-X
"END";
EVALCOEF:=H^=HL;
"END" "ELSE"
"IF" FIRST "THEN"
"BEGIN" H:=HNEW:=HMIN; FIRST:="FALSE"; INFO[4]:=INFO[4]+1
"END" "ELSE"
"BEGIN" "REAL" B,HL;
B:=DISCR/ETA; HL:=H; "IF" B<.01 "THEN" B:=.01;
HNEW:= "IF" B>0 "THEN" H*B**(-1/P) "ELSE" HMAX;
"IF" HNEW<HMIN "THEN"
"BEGIN" HNEW:=HMIN; INFO[4]:=INFO[4]+1
"END" "ELSE"
"IF" HNEW>HMAX "THEN"
"BEGIN" HNEW:=HMAX; INFO[5]:=INFO[5]+1 "END";
"IF" 1.1*HNEW>=XE-X "THEN"
"BEGIN" LAST:="TRUE"; H:=HNEW:=XE-X
"END" "ELSE"
"IF" ABS(H/HNEW-1)>.1 "THEN" H:=HNEW;
EVALCOEF:=H^=HL
"END" STEPSIZE;
"PROCEDURE" LINEARITY(ERROR); "REAL" ERROR;
"BEGIN" "INTEGER" K;
"FOR" K:=1 "STEP" 1 "UNTIL" M "DO"
DY[K]:=Y[K]-MU1*F[K];
SOL(A,M,PI,DY);
ELMVEC(1,M,0,DY,Y,-1);
ERROR:=NORM(DY)
"END" LINEARITY; "COMMENT"
1SECTION : 5.2.1.1.1.2.C (OCTOBER 1974) PAGE 7
;
"PROCEDURE" ITERATION(I); "INTEGER" I;
"IF" RETA<0 "THEN"
"BEGIN" "INTEGER" K;
"IF" I=1 "THEN"
"BEGIN" MULVEC(1,M,0,DY,F,H);
"FOR" K:=1 "STEP" 1 "UNTIL" M "DO" YL[K]:=Y[K]+MU*DY[K];
SOL(A,M,PI,DY); E:=1;
"END" "ELSE"
"BEGIN" "FOR" K:=1 "STEP" 1 "UNTIL" M "DO"
DY[K]:=YL[K]-Y[K]+MU1*F[K];
"IF" E*NORM(Y)>E1*E1 "THEN"
"BEGIN" EVALJAC:=I>=3;
"IF" I>3 "THEN"
"BEGIN" INFO[3]:=INFO[3]+1; JACOBIAN(J,Y);
LUDECOMP
"END";
"END";
SOL(A,M,PI,DY)
"END";
E1:=E; E:=NORM(DY);
ELMVEC(1,M,0,Y,DY,1);
ETA:=NORM(Y)*RETA+AETA;
DISCR:=0;
DUPVEC(1,M,0,F,Y);
DERIVATIVE(F);
INFO[2]:=INFO[2]+1;
"END" "ELSE"
"BEGIN" "INTEGER" K;
"IF" I=1 "THEN"
"BEGIN" LINEARITY(E);
"IF" E*(ST-LASTJAC)>ETA "THEN"
"BEGIN" JACOBIAN(J,Y); LASTJAC:=ST;
INFO[3]:=INFO[3]+1;
H:=HNEW; COEFFICIENT;
LINEARITY(E)
"END";
EVALJAC:= E*(ST+1-LASTJAC)>ETA;
MULVEC(1,M,0,DY,F,H);
"FOR" K:=1 "STEP" 1 "UNTIL" M "DO" YL[K]:=Y[K]+MU*DY[K];
SOL(A,M,PI,DY);
"FOR" K:=1 "STEP" 1 "UNTIL" M "DO"
YR[K]:=H*BETA*MATVEC(1,M,K,J,DY);
SOL(A,M,PI,YR);
ELMVEC(1,M,0,YR,DY,1); "COMMENT"
1SECTION : 5.2.1.1.1.2.C (OCTOBER 1974) PAGE 8
;
"END" "ELSE"
"BEGIN" "FOR" K:=1 "STEP" 1 "UNTIL" M "DO"
DY[K]:=YL[K]-Y[K]+MU1*F[K];
"IF" E>ETA1 "AND" DISCR>ETA1 "THEN"
"BEGIN" INFO[3]:=INFO[3]+1; JACOBIAN(J,Y);
LUDECOMP
"END";
SOL(A,M,PI,DY);
E:=NORM(DY)
"END";
ELMVEC(1,M,0,Y,DY,1);
ETA:=NORM(Y)*RETA+AETA;
ETA1:=ETA/SQRT(RETA);
DUPVEC(1,M,0,F,Y);
DERIVATIVE(F);
INFO[2]:=INFO[2]+1;
"FOR" K:=1 "STEP" 1 "UNTIL" M "DO" DY[K]:=YR[K]-H*F[K];
DISCR:=NORM(DY)/2
"END" ITERATION;
FIRST:=EVALJAC:="TRUE"; LAST:=EVALCOEF:="FALSE";
INIVEC(1,9,INFO,0);
ETA:=RETA*NORM(Y)+AETA;
ETA1:=ETA/SQRT(ABS(RETA));
DUPVEC(1,M,0,F,Y);
DERIVATIVE(F);
INFO[2]:=1;
"FOR" ST:=1,ST+1 "WHILE" ^LAST "DO"
"BEGIN" STEPSIZE; INFO[1]:=INFO[1]+1;
"IF" EVALJAC "THEN"
"BEGIN" JACOBIAN(J,Y);
INFO[3]:=INFO[3]+1;
H:=HNEW;
COEFFICIENT;
EVALJAC:="FALSE"; LASTJAC:=ST
"END" "ELSE"
"IF" EVALCOEF "THEN" COEFFICIENT;
"FOR" I:=1,I+1 "WHILE" E>ABS(ETA) "AND" DISCR>1.3*ETA
"AND" I<=ITMAX "DO"
"BEGIN" ITERATION(I); "IF" I>INFO[6] "THEN" INFO[6]:=I;
"END"; INFO[7]:=ETA; INFO[8]:=DISCR;
X:=X+H;
"IF" DISCR>INFO[9] "THEN" INFO[9]:=DISCR;
OUTPUT;
"END";
"END" LINIGER1VS;
"EOP"
1SECTION : 5.2.1.1.1.2.F � (OCTOBER 1975) PAGE 1
AUTHOR: B.LINDBERG.
CONTRIBUTOR: K.DEKKER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 741101.
BRIEF DESCRIPTION:�
IMPEX SOLVES AN INITIAL VALUE PROBLEM,GIVEN AS AN AUTONOMOUS SYSTEM
OF FIRST ORDER DIFFERENTIAL EQUATIONS , BY MEANS OF THE IMPLICIT
MIDPOINT RULE WITH SMOOTHING AND EXTRAPOLATION.
AUTOMATIC STEPSIZE CONTROL IS PROVIDED.
IN PARTICULAR THIS METHOD IS SUITABLE FOR THE INTEGRATION OF STIFF
DIFFERENTIAL EQUATIONS.
KEYWORDS:
DIFFERENTIAL EQUATIONS,
INITIAL VALUE PROBLEMS,
STIFF EQUATIONS,
IMPLICIT MIDPOINT RULE,
SMOOTHING,
EXTRAPOLATION.�
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE IMPEX READS:
"PROCEDURE" IMPEX (N, T0, TEND, Y0, DERIV, AVAILABLE, H0, HMAX,
PRESCH, EPS, WEIGHTS, UPDATE, FAIL, CONTROL);
"VALUE" N;
"INTEGER" N;
"REAL" T0,TEND,H0,HMAX,EPS;
"BOOLEAN" PRESCH,FAIL;
"ARRAY" Y0,WEIGHTS;
"BOOLEAN" "PROCEDURE" AVAILABLE;
"PROCEDURE" DERIV,UPDATE,CONTROL;
IMPEX: INTEGRATES THE SYSTEM OF DIFFERENTIAL EQUATIONS , GIVEN AS
DY/DT=F(T,Y), FROM T0 UNTIL TEND.
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE NUMBER OF EQUATIONS;
T0: <ARITHMETIC EXPRESSION>;
THE INITIAL VALUE OF THE INDEPENDENT VARIABLE;�
1SECTION : 5.2.1.1.1.2.F � (OCTOBER 1975) PAGE 2
TEND: <ARITHMETIC EXPRESSION>;
THE FINAL VALUE OF THE INDEPENDENT VARIABLE;
Y0: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" Y0[1:N];
ENTRY: THE INITIAL VALUES OF THE SYSTEM OF DIFFERENTIAL
EQUATIONS: Y0[I] AT T=T0;
DERIV: <PROCEDURE" IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:�
"PROCEDURE" DERIV(T,Y,F,N);
"INTEGER" N; "REAL" T; "ARRAY" Y,F;
THIS PROCEDURE SHOULD DELIVER THE VALUE OF F(T,Y) IN THE
ARRAY F[1:N];
AVAILABLE: <PROCEDURE IDENTIFIER>;�
THE HEADING OF THIS PROCEDURE READS:�
"BOOLEAN" "PROCEDURE" AVAILABLE(T,Y,A,N);
"INTEGER" N; "REAL" T; "ARRAY" Y,A;
IF AN ANALYTIC EXPRESSION OF THE JACOBIAN MATRIX AT THE
POINT (T,Y) IS NOT AVAILABLE, THIS PROCEDURE SHOULD DELIVER
THE VALUE "FALSE";
OTHERWISE THE PROCEDURE SHOULD DELIVER THE VALUE "TRUE",
AND THE JACOBIAN MATRIX SHOULD BE ASSIGNED TO THE ARRAY
A[1:N,1:N];
H0: <ARITHMETIC EXPRESSION>;
THE INITIAL STEPSIZE;
HMAX: <ARITHMETIC EXPRESSION>;
MAXIMAL STEPSIZE BY WHICH THE INTEGRATION IS PERFORMED;
PRESCH: <BOOLEAN EXPRESSION>;
INDICATOR FOR CHOICE OF STEPSIZE;
THE STEPSIZE IS AUTOMATICALLY CONTROLLED IF PRESCH="FALSE";
OTHERWISE THE STEPSIZE HAS TO BE PRESCRIBED , EITHER ONLY
INITIALLY OR ALSO BY THE PROCEDURE CONTROL;
EPS: <ARITHMETIC EXPRESSION>;
BOUND FOR THE ESTIMATE OF THE LOCAL ERROR;
WEIGHTS: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" WEIGHTS[1:N];
WEIGHTS FOR THE COMPUTATION OF THE WEIGHTED EUCLIDEAN NORM
OF THE ERRORS;
ENTRY: INITIAL WEIGHTS;
NOTE THAT THE CHOICE WEIGHTS[I] = 1 IMPLIES AN ESTIMATION
OF THE ABSOLUTE ERROR , WHEREAS WEIGTHS[I] = Y[I] DEFINES A
RELATIVE ERROR;
UPDATE: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:�
"PROCEDURE" UPDATE(WEIGHTS,Y2,N);
"INTEGER" N; "ARRAY" WEIGHTS,Y2;
THIS PROCEDURE MAY CHANGE THE ARRAY WEIGHTS , ACCORDING TO
THE VALUE OF AN APPROXIMATION FOR Y(T) , GIVEN IN THE ARRAY
Y2[1:N];
1SECTION : 5.2.1.1.1.2.F � (OCTOBER 1975) PAGE 3
FAIL: <BOOLEAN EXPRESSION>;
EXIT :
IF THE PROCEDURE FAILS TO SOLVE THE SYSTEM OF EQUATIONS,
FAIL WILL HAVE THE VALUE "TRUE" UPON EXIT;
THIS MAY OCCUR UPON DIVERGENCE OF THE ITERATION PROCESS,
USED IN THE MIDPOINT RULE , WHILE INTEGRATION IS PERFORMED
WITH A USER DEFINED PRESCRIBED STEPSIZE;
CONTROL: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS:�
"PROCEDURE" CONTROL(TPRINT,T,H,HNEW,Y,ERROR,N);
"INTEGER" N; "REAL" TPRINT,T,H,HNEW; "ARRAY" Y,ERROR;
CONTROL IS CALLED ON ENTRY OF IMPEX, AND FURTHER AS SOON AS
THE INEQUALITY TPRINT <= T HOLDS;
DURING A CALL OF CONTROL PRINTING OF RESULTS AND
CHANGE OF STEPSIZE (IF PRESCH = "TRUE") IS THEN POSSIBLE;
THE MEANING OF THE FORMAL PARAMETERS IS:
TPRINT: <VARIABLE>;
ENTRY: THE VALUE OF THE INDEPENDENT VARIABLE AT
WHICH A CALL OF CONTROL WAS DESIRED;
EXIT: A NEW VALUE (TPRINT>T) AT WHICH A CALL OF
CONTROL IS DESIRED;
T: <VARIABLE>;
THE ACTUAL VALUE OF THE INDEPENDENT VARIABLE, UP TO
WHICH INTEGRATION HAS BEEN PERFORMED;
H: <VARIABLE>;
HALVE THE ACTUAL STEPSIZE;
HNEW: <VARIABLE>;
THE NEW STEPSIZE;
IF PRESCH="TRUE", THEN THE USER MAY PRESCRIBE A NEW
STEPSIZE BY CHANGING HNEW;
Y: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" Y[1:5,1:N];
THE VALUE OF THE DEPENDENT VARIABLE AND ITS FIRST
FOUR DIVIDED DIFFERENCES AT THE POINT T ARE GIVEN
IN THIS ARRAY;
ERROR: <ARRAY IDENTIFIER>;
"REAL" "ARRAY" ERROR[1:3];
THE ELEMENTS OF THIS ARRAY CONTAIN THE FOLLOWING
ERRORS:
ERROR[1]: THE LOCAL ERROR;
ERROR[2]: THE GLOBAL ERROR OF SECOND ORDER IN H;�
ERROR[3]: THE GLOBAL ERROR OF FOURTH ORDER IN H;�
N: <VARIABLE>;
THE NUMBER OF EQUATIONS;
EXAMPLE OF USE: SEE EXAMPLE OF USE OF THE PROCEDURE IMPEX;
DATA AND RESULTS:
FOR DATA, SEE REF[1].
THE RESULTS OF THE INTEGRATION ARE ATTAINABLE THROUGH THE PROCEDURE
CONTROL , WHICH IS CALLED AT SPECIFIED, USER DEFINED, VALUES OF THE
INDEPENDENT VARIABLE . IN PARTICULAR , THE VALUES OF THE DEPENDENT
VARIABLE AT THE ENDPOINT OF INTEGRATION ARE OBTAINED BY A CALL OF
CONTROL WITH TPRINT=TEND.
1SECTION : 5.2.1.1.1.2.F � (OCTOBER 1975) PAGE 4
PROCEDURES USED:
INIVEC = CP31010,
INIMAT = CP31011,
MULVEC = CP31020,
MULROW = CP31021,
DUPVEC = CP31030,
DUPROWVEC= CP31032,
DUPMAT = CP31035,
VECVEC = CP34010,
MATVEC = CP34011,
MATMAT = CP34013,
ELMVEC = CP34020,
ELMROW = CP34024,
DEC = CP34300,
SOL = CP34051.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: CIRCA N * ( 23 + 2 * N ) (DECIMAL).
RUNNING TIME: DEPENDS STRONGLY ON THE DIFFERENTIAL EQUATION TO SOLVE.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
THE INTEGRATION METHOD (REF[1]) IS BASED ON THE COMPUTATION OF TWO
INDEPENDENT SOLUTIONS Y(T,H) AND Y(T,H/2) BY THE IMPLICIT MIDPOINT
RULE. PASSIVE SMOOTHING AND PASSIVE EXTRAPOLATION IS PERFORMED TO
OBTAIN STABILITY AND HIGH ACCURACY . THE ALGORITHM USES FOR EACH
STEP AT LEAST THREE FUNCTION EVALUATIONS, AND ON CHANGE OF STEPSIZE
OR AT SLOW CONVERGENCE IN THE ITERATION PROCESS AN APPROXIMATION OF
THE JACOBIAN MATRIX (COMPUTED BY DIVIDED DIFFERENCES OR EXPLICITLY
SPECIFIED BY THE USER) . IF THE COMPUTED LOCAL ERROR EXCEEDS THE
TOLERANCE , THE LAST STEP IS REJECTED. MOREOVER , TWO GLOBAL ERRORS
ARE COMPUTED.
REFERENCES:
[1]. B.LINDBERG.
IMPEX 2 , A PROCEDURE FOR THE SOLUTION OF SYSTEMS OF STIFF
DIFFERENTIAL EQUATIONS.
ROYAL INSTITUTE OF TECHNOLOGY,STOCKHOLM. TRITA-NA-7303 (1973).
[2]. (TO APPEAR)
COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 3 (DUTCH).
M.C. SYLLABUS 15.3 (1974), MATHEMATICAL CENTRE.
1SECTION : 5.2.1.1.1.2.F � (OCTOBER 1975) PAGE 5
EXAMPLE OF USE:
CONSIDER THE AUTONOMOUS SYSTEM OF DIFFERNETIAL EQUATIONS:
DY[1]/DX = .2 * ( Y[2] - Y[1] ),
DY[2]/DX = 10 * Y[1] - ( 60 - Y[3]/8 ) * Y[2] + Y[3]/8,
DY[3]/DX = 1,
WITH INITIAL CONDITIONS AT X=0 : Y[1]=Y[2]=Y[3]=0 (SEE REF[2]).
THE SOLUTION AT SEVERAL POINTS IN THE INTERVAL [0, 400] MAY BE
OBTAINED BY THE FOLLOWING PROGRAM:�
(THE SOLUTION AT X=400 IS: Y[1]=22.24222011, Y[2]=27.11071335)
"BEGIN" "INTEGER" N,NFE,NJE,POINT;
"REAL" T,TEND,EPS,HMAX,L,H2,TIME;
"ARRAY" Y,SW[1:3],PRINT[1:5];
"BOOLEAN" FAIL;
"PROCEDURE" IMPEX(N,T0,TEND,Y0,DERIV,AVAIL,H0,HMAX,PRESCH,EPS,
WEIGHTS,UPDATE,FAIL,CONTROL); "CODE" 33135;
"PROCEDURE" LIPEST(L,Y,EPS,T,F,N);�
"REAL" T,L,EPS; "ARRAY" Y; "INTEGER" N; "PROCEDURE" F;�
"BEGIN" "REAL" N1,N2; "INTEGER" I,IT; "ARRAY" F1,F2,Z,X[1:N];
"PROCEDURE" DUPVEC(L,U,SHIFT,A,B); "CODE" 31030;
"REAL" "PROCEDURE" VECVEC(L,U,SHIFT,A,B); "CODE" 34010;
"PROCEDURE" ELMVEC(L,U,SHIFT,A,B,X); "CODE" 34020;�
"REAL" "PROCEDURE" NORM(Y); "ARRAY" Y;
NORM:=SQRT(VECVEC(1,N,0,Y,Y));�
DUPVEC(1,N,0,Z,Y);
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
X[I]:="IF" Y[I]=0 "THEN" EPS "ELSE" (1+EPS)*Y[I];
N1:=NORM(X)*EPS; F(T,X,F1,N);
"FOR" IT:=1 "STEP" 1 "UNTIL" 5 "DO"
"BEGIN" F(T,Z,F2,N);�
ELMVEC(1,N,0,F2,F1,-1);
N2:=N1/NORM(F2);�
DUPVEC(1,N,0,Z,X); ELMVEC(1,N,0,Z,F2,N2)
"END";
F(T,Z,F2,N);
ELMVEC(1,N,0,F2,F1,-1);
L:=NORM(F2)/N1
"END" LIPEST;
"PROCEDURE" F(T,Y,F1,N); "VALUE" T; "REAL" T; "ARRAY" Y,F1;
"INTEGER" N;
"BEGIN" NFE:=NFE+1;
F1[1]:=0.2*(Y[2]-Y[1]);
F1[2]:=10*Y[1]-(60-.125*Y[3])*Y[2]+.125*Y[3];
F1[3]:=1
"END";
"BOOLEAN" "PROCEDURE" AVAILABLE(T,Y,A,N);
"INTEGER" N; "REAL" T; "ARRAY" Y,A;
"BEGIN" NJE:=NJE+1; AVAILABLE:="TRUE";
A[1,1]:=-.2; A[1,2]:=.2; A[1,3]:=A[3,1]:=A[3,2]:=A[3,3]:=0;
A[2,1]:=10; A[2,2]:=.125*Y[3]-60; A[2,3]:=.125*(1+Y[2])
"END"
1SECTION : 5.2.1.1.1.2.F � (OCTOBER 1975) PAGE 6
;
"PROCEDURE" UPDATE(SW,R1,N); "INTEGER" N; "ARRAY" SW,R1;
"BEGIN" "REAL" S1,S2; "INTEGER" I;�
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
"BEGIN" S1:=1/SW[I]; S2:=ABS(R1[I]);
"IF" S1<S2 "THEN" SW[I]:=1/S2
"END"
"END";
"PROCEDURE" CONTROL(TP,T,H,HNEW,Y,ERR,N);
"REAL" TP,T,H,HNEW; "ARRAY" Y,ERR; "INTEGER" N;
"BEGIN" "INTEGER" I;
"ARRAY" C[3:5],X[1:N];
"REAL" S,S2,S3,S4,C1;
NEXT: S:=(T-TP)/H;
S2:=S*S; S3:=S2*S; S4:=S3*S;
C[3]:=(S2-S)/2;
C[4]:=-S3/6+S2/2-S/3;
C[5]:=S4/24-S3/4+11*S2/24-S/4;�
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
X[I]:=Y[1,I]-S*Y[2,I]+C[3]*Y[3,I]+C[4]*Y[4,I]+C[5]*Y[5,I];
OUTPUT(61,"("3ZD.2D2B,D.D"-D2B,2(+D.8D"-D2B),2(4ZD),3ZD.2D,
/")",TP,ERR[3],X[1],X[2],NFE,NJE,CLOCK-TIME);
"IF" TP<TEND "THEN"
"BEGIN" POINT:=POINT+1; TP:=PRINT[POINT];
"IF" TP<=T "THEN" "GOTO" NEXT
"END"
"END" CONTROL;�
N:=3; NJE:=NFE:=0; T:=0; TEND:=400; EPS:="-5; HMAX:=400;
Y[1]:=Y[2]:=Y[3]:=0; SW[1]:=SW[2]:=SW[3]:=1;�
PRINT[1]:=.1; PRINT[2]:=1; PRINT[3]:=10; PRINT[4]:=100;
PRINT[5]:=400;�
LIPEST(L,Y,"-5,T,F,N);
H2:=(EPS*320)**(1/5)/(4*L);
OUTPUT(61,"(""("EPS=")",D.2D"-D,/,"("INTERVAL OF INTEGRATION=(")",
3ZD,"(",")",3ZD,"(")")",/,"("MAXIMALLY ALLOWED STEPSIZE=")",
D.2D"-D,//")",EPS,T,TEND,HMAX);
OUTPUT(61,"(""("LIPSCHCONST=")",BD.3D"+D,/,"("STARTING STEPSIZE")"
"("=")",BD.2D"+D,/,"("FUNCTIONAL EVAL=")",4ZD,//")",L,H2,NFE);
TIME:=CLOCK;
OUTPUT(61,"(""(" X ERROR Y[1] Y[2]")"
"(" NFE NJE TIME")",/")");
IMPEX(N,T,TEND,Y,F,AVAILABLE,H2,HMAX,"FALSE",EPS,SW,UPDATE,FAIL,�
CONTROL);
OUTPUT(61,"("/"("NO OF FUNCTIONAL EVALUATIONS= ")",3ZD,/,
"("NO OF JACOBEAN EVALUATIONS= ")",3ZD,/")",NFE,NJE)
"END"
1SECTION : 5.2.1.1.1.2.F � (OCTOBER 1975) PAGE 7
THIS PROGRAM DELIVERS:
EPS=1.00"-5
INTERVAL OF INTEGRATION=( 0, 400)
MAXIMALLY ALLOWED STEPSIZE=4.00" 2
LIPSCHCONST= 6.003"+1
STARTING STEPSIZE= 1.32"-3
FUNCTIONAL EVAL= 7
X ERROR Y[1] Y[2] NFE NJE TIME
0.00 0.0" 0 +0.00000000" 0 +0.00000000" 0 7 0 0.01
0.10 6.3"-7 +1.49614151"-6 +1.74013792"-4 46 4 0.72
1.00 1.5"-6 +1.91041887"-4 +2.08361269"-3 85 8 1.48
10.00 8.7"-7 +1.30147663"-2 +2.34487800"-2 119 9 1.99
100.00 1.3"-5 +3.06302487"-1 +3.27552180"-1 225 13 3.47
400.00 1.4"-5 +2.22406546" 1 +2.71090507" 1 556 30 7.51
NO OF FUNCTIONAL EVALUATIONS= 556
NO OF JACOBEAN EVALUATIONS= 30
1SECTION : 5.2.1.1.1.2.F (OCTOBER 1975) PAGE 8
SOURCE TEXT(S):
0"CODE" 33135;
"PROCEDURE" IMPEX (N, T0, TEND, Y0, DERIV, AVAILABLE, H0, HMAX,
PRESCH, EPS, WEIGHTS, UPDATE, FAIL, CONTROL);
"VALUE" N;
"INTEGER" N;
"REAL" T0,TEND,H0,HMAX,EPS;
"BOOLEAN" PRESCH,FAIL;
"ARRAY" Y0,WEIGHTS;
"BOOLEAN" "PROCEDURE" AVAILABLE;
"PROCEDURE" DERIV,UPDATE,CONTROL;
"BEGIN" "INTEGER" I,K,ECI;
"REAL" T,T1,T2,T3,TP,H,H2,HNEW,ALF,LQ;
"ARRAY" Y,Z,S1,S2,S3,U1,U3,W1,W2,W3,EHR[1:N],R,RF[1:5,1:N],
ERR[1:3],A1,A2[1:N,1:N];
"INTEGER" "ARRAY" PS1,PS2[1:N];
"BOOLEAN" START,TWO,HALV;
"PROCEDURE" INIVEC(L,U,A,X); "CODE" 31010;
"PROCEDURE" INIMAT(LR,UR,LC,UC,A,X); "CODE" 31011;
"PROCEDURE" MULVEC(L,U,SHIFT,A,B,X); "CODE" 31020;
"PROCEDURE" MULROW(L,U,I,J,A,B,X); "CODE" 31021;
"PROCEDURE" DUPVEC(L,U,SHIFT,A,B); "CODE" 31030;
"PROCEDURE" DUPROWVEC(L,U,I,A,B); "CODE" 31032;
"PROCEDURE" DUPMAT(L,U,I,J,A,B); "CODE" 31035;
"REAL" "PROCEDURE" VECVEC(L,U,SHIFT,A,B);"CODE" 34010;
"REAL" "PROCEDURE" MATVEC(L,U,I,A,B); "CODE" 34011;
"REAL" "PROCEDURE" MATMAT(L,U,I,J,A,B); "CODE" 34013;
"PROCEDURE" ELMVEC(L,U,SHIFT,A,B,X); "CODE" 34020;
"PROCEDURE" ELMROW(L,U,I,J,A,B,X); "CODE" 34024;
"PROCEDURE" DEC(A,N,AUX,P); "CODE" 34300;
"PROCEDURE" SOL(A,N,P,B); "CODE" 34051;
"PROCEDURE" DFDY(T,Y,A); "REAL" T; "ARRAY" Y,A;
"BEGIN" "INTEGER" I,J; "REAL" SL; "ARRAY" F1,F2[1:N];
DERIV(T,Y,F1,N);
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
"BEGIN"
SL:="-6*Y[I]; "IF" ABS(SL)<"-6 "THEN" SL:="-6;
Y[I]:=Y[I]+SL; DERIV(T,Y,F2,N);
"FOR" J:=1 "STEP" 1 "UNTIL" N "DO"
A[J,I]:=(F2[J]-F1[J])/SL;
Y[I]:=Y[I]-SL;
"END"
"END" DFDY;
"PROCEDURE" STARTV(Y,T); "VALUE" T; "REAL" T; "ARRAY" Y;
"BEGIN" "REAL" A,B,C;
A:=(T-T1)/(T1-T2); B:=(T-T2)/(T1-T3);
C:=(T-T1)/(T2-T3)*B; B:=A*B;
A:=1+A+B; B:=A+C-1;
MULVEC(1,N,0,Y,S1,A); ELMVEC(1,N,0,Y,S2,-B);
ELMVEC(1,N,0,Y,S3,C)
"END" STARTV;
1SECTION : 5.2.1.1.1.2.F (OCTOBER 1975) PAGE 9
;
"PROCEDURE" ITERATE(Z,Y,A,H,T,WEIGHTS,FAIL,PS);
"ARRAY" Z,Y,A,WEIGHTS; "REAL" H,T; "LABEL" FAIL;
"INTEGER" "ARRAY" PS;
"BEGIN" "INTEGER" IT,LIT; "REAL" MAX,MAX1,CONV; "ARRAY" DZ,F1[1:N];
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO" Z[I]:=(Z[I]+Y[I])/2;
IT:=LIT:=1; CONV:=1;
ATER: DERIV(T,Z,F1,N);
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
F1[I]:=DZ[I]:=Z[I]-H*F1[I]/2-Y[I];
SOL(A,N,PS,DZ);
ELMVEC(1,N,0,Z,DZ,-1);
MAX:=0;
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
MAX:=MAX+(WEIGHTS[I]*DZ[I])**2;
MAX:=SQRT(MAX);
"IF" MAX*CONV<EPS/10 "THEN" "GOTO" OUT;
IT:=IT+1; "IF" IT=2 "THEN" "GOTO" ASS;
CONV:=MAX/MAX1;
"IF" CONV>.2 "THEN"
"BEGIN" "IF" LIT=0 "THEN" "GOTO" FAIL;
LIT:=0; CONV:=1; IT:=1;
RECOMP(A,H,T,Z,FAIL,PS);
"END";
ASS: MAX1:=MAX;
"GOTO" ATER;
OUT: "FOR" I:=1 "STEP" 1 "UNTIL" N "DO" Z[I]:=2*Z[I]-Y[I];
"END" ITERATE;
"PROCEDURE" RECOMP(A,H,T,Y,FAIL,PS);
"REAL" H,T; "ARRAY" A,Y; "LABEL" FAIL; "INTEGER" "ARRAY" PS;
"BEGIN" "REAL" SL; "ARRAY" AUX[1:3];
SL:=H/2;
"IF" "NOT" AVAILABLE(T,Y,A,N) "THEN" DFDY(T,Y,A);
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
"BEGIN" MULROW(1,N,I,I,A,A,-SL); A[I,I]:=1+A[I,I]
"END";
AUX[2]:="-14;
DEC(A,N,AUX,PS);
"IF" AUX[3]<N "THEN" "GOTO" FAIL
"END" RECOMP;
"PROCEDURE" INITIALIZATION;
"BEGIN" H2:=HNEW; H:=H2/2;
DUPVEC(1,N,0,S1,Y0); DUPVEC(1,N,0,S2,Y0); DUPVEC(1,N,0,S3,Y0);
DUPVEC(1,N,0,W1,Y0); DUPROWVEC(1,N,1,R,Y0);
INIVEC(1,N,U1,0); INIVEC(1,N,W2,0);
INIMAT(2,5,1,N,R,0); INIMAT(1,5,1,N,RF,0);
T:=T1:=T0; T2:=T0-2*H-"6; T3:=2*T2+1;
RECOMP(A1,H,T,S1,MISS,PS1);RECOMP(A2,H2,T,W1,MISS,PS2);
"END"
1SECTION : 5.2.1.1.1.2.F (OCTOBER 1975) PAGE 10
;
"PROCEDURE" ONE LARGE STEP;
"BEGIN" STARTV(Z,T+H);
ITERATE(Z,S1,A1,H,T+H/2,WEIGHTS,MISS,PS1);
DUPVEC(1,N,0,Y,Z);
STARTV(Z,T+H2);
ITERATE(Z,Y,A1,H,T+3*H/2,WEIGHTS,MISS,PS1);
DUPVEC(1,N,0,U3,U1); DUPVEC(1,N,0,U1,Y);
DUPVEC(1,N,0,S3,S2); DUPVEC(1,N,0,S2,S1);
DUPVEC(1,N,0,S1,Z);
ELMVEC(1,N,0,Z,W1,1); ELMVEC(1,N,0,Z,S2,-1);
ITERATE(Z,W1,A2,H2,T+H,WEIGHTS,MISS,PS2);
T3:=T2; T2:=T1; T1:=T+H2;
DUPVEC(1,N,0,W3,W2); DUPVEC(1,N,0,W2,W1); DUPVEC(1,N,0,W1,Z);
"END";
"PROCEDURE" CHANGE OF INFORMATION;
"BEGIN" "REAL" ALF1,C1,C2,C3; "ARRAY" KOF[2:4,2:4],E,D[1:4];
C1:=HNEW/H2; C2:=C1*C1; C3:=C2*C1;
KOF[2,2]:=C1; KOF[2,3]:=(C1-C2)/2; KOF[2,4]:=C3/6-C2/2+C1/3;
KOF[3,3]:=C2; KOF[3,4]:=C2-C3; KOF[4,4]:=C3;
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
U1[I]:=R[2,I]+R[3,I]/2+R[4,I]/3;
ALF1:=MATVEC(1,N,1,RF,U1)/VECVEC(1,N,0,U1,U1);
ALF:=(ALF+ALF1)*C1;
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
"BEGIN"
E[1]:=RF[1,I]-ALF1*U1[I];
E[2]:=RF[2,I]-ALF1*2*R[3,I];
E[3]:=RF[3,I]-ALF1*4*R[4,I];
E[4]:=RF[4,I];
D[1]:=R[1,I]; RF[1,I]:=E[1]:=E[1]*C2;
"FOR" K:=2 "STEP" 1 "UNTIL" 4 "DO"
"BEGIN" R[K,I]:=D[K]:=MATMAT(K,4,K,I,KOF,R);
RF[K,I]:=E[K]:=C2*MATVEC(K,4,K,KOF,E)
"END" K;
S1[I]:=D[1]+E[1];W1[I]:=D[1]+4*E[1];
S2[I]:=S1[I]-(D[2]+E[2]/2);
S3[I]:=S2[I]-(D[2]+E[2])+(D[3]+E[3]/2);
"END" I;
T3:=T-HNEW; T2:=T-HNEW/2; T1:=T;
H2:=HNEW; H:=H2/2; ERR[1]:=0;
"IF" HALV "THEN"
"BEGIN" DUPVEC(1,N,0,PS2,PS1); DUPMAT(1,N,1,N,A2,A1) "END";
"IF" TWO "THEN"
"BEGIN" DUPVEC(1,N,0,PS1,PS2); DUPMAT(1,N,1,N,A1,A2)
"END" "ELSE" RECOMP(A1,HNEW/2,T,S1,MISS,PS1);
"IF" ^HALV "THEN" RECOMP(A2,HNEW,T,W1,MISS,PS2);
"END";
1SECTION : 5.2.1.1.1.2.F (OCTOBER 1975) PAGE 11
;
"PROCEDURE" BACKWARD DIFFERENCES;
"FOR"I:=1 "STEP" 1 "UNTIL" N "DO"
"BEGIN" "REAL" B0,B1,B2,B3;
B1:=(U1[I]+2*S2[I]+U3[I])/4;
B2:=(W1[I]+2*W2[I]+W3[I])/4;
B3:=(S3[I]+2*U3[I]+S2[I])/4;
B2:=(B2-B1)/3; B0:=B1-B2;
B2:=B2-(S1[I]-2*S2[I]+S3[I])/16;
B1:=2*B3-(B2+RF[1,I])-(B0+R[1,I])/2;
B3:=0;
"FOR" K:=1 "STEP" 1 "UNTIL" 4 "DO"
"BEGIN" B1:=B1-B3; B3:=R[K,I]; R[K,I]:=B0; B0:=B0-B1
"END"; R[5,I]:=B0;
"FOR" K:=1 "STEP" 1 "UNTIL" 4 "DO"
"BEGIN" B3:=RF[K,I]; RF[K,I]:=B2; B2:=B2-B3 "END";
RF[5,I]:=B2;
"END";
"PROCEDURE" ERROR ESTIMATES;
"BEGIN" "REAL" C0,C1,C2,C3,B0,B1,B2,B3,W,SL1,SN,LR;
C0:=C1:=C2:=C3:=0;
"FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
"BEGIN" W:=WEIGHTS[I]**2;
B0:=RF[4,I]/36; C0:=C0+B0*B0*W; LR:=ABS(B0);
B1:=RF[1,I]+ALF*R[2,I]; C1:=C1+B1*B1*W;
B2:=RF[3,I]; C2:=C2+B2*B2*W;
SL1:=ABS(RF[1,I]-RF[2,I]);
SN:="IF" SL1<"-10 "THEN"1"ELSE"ABS(RF[1,I]-R[4,I]/6)/SL1;
"IF" SN>1 "THEN" SN:=1;
"IF" START "THEN" "BEGIN" SN:=SN**4; LR:=LR*4 "END";
EHR[I]:=B3:=SN*EHR[I]+LR; C3:=C3+B3*B3*W;
"END" I;
B0:=ERR[1];
ERR[1]:=B1:=SQRT(C0); ERR[2]:=SQRT(C1);
ERR[3]:=SQRT(C3)+SQRT(C2)/2;
LQ:=EPS/("IF" B0<B1 "THEN" B1"ELSE" B0);
"IF" B0<B1 "AND" LQ>=80 "THEN" LQ:=10;
"END";
"PROCEDURE" REJECT;
"IF" START "THEN"
"BEGIN" HNEW:=LQ**(1/5)*H/2; "GOTO" INIT
"END" "ELSE"
"BEGIN" "FOR" K:=1,2,3,4,1,2,3 "DO" ELMROW(1,N,K,K+1,R,R,-1);
"FOR" K:=1,2,3,4 "DO" ELMROW(1,N,K,K+1,RF,RF,-1);
T:=T-H2; HALV:="TRUE"; HNEW:=H; "GOTO" MSTP
"END";
1SECTION : 5.2.1.1.1.2.F (OCTOBER 1975) PAGE 12
;
"PROCEDURE" STEPSIZE;
"IF" LQ<2 "THEN"
"BEGIN" HALV:="TRUE"; HNEW:=H "END" "ELSE"
"BEGIN" "IF" LQ>80 "THEN"
HNEW:=("IF" LQ>5120 "THEN" (LQ/5)**(1/5) "ELSE" 2)*H2;
"IF" HNEW>HMAX "THEN" HNEW:=HMAX;
"IF" TEND>T "AND" TEND-T<HNEW "THEN" HNEW:=TEND-T;
TWO:=HNEW=2*H2;
"END";
"IF" PRESCH "THEN" H:=H0 "ELSE"
"BEGIN" "IF" H0>HMAX "THEN" H:=HMAX "ELSE" H:=H0;
"IF" H>(TEND-T0)/4 "THEN" H:=(TEND-T0)/4;
"END";
HNEW:=H;
ALF:=0; T:=TP:=T0;
INIVEC(1,3,ERR,0); INIVEC(1,N,EHR,0);
DUPROWVEC(1,N,1,R,Y0);
CONTROL(TP,T,H,HNEW,R,ERR,N);
INIT: INITIALIZATION; START:="TRUE";
"FOR" ECI:=0,1,2,3 "DO"
"BEGIN" ONE LARGE STEP; T:=T+H2;
"IF" ECI>0 "THEN"
"BEGIN" BACKWARD DIFFERENCES; UPDATE(WEIGHTS,S2,N) "END"
"END";
ECI:=4;
MSTP: "IF" HNEW^=H2 "THEN"
"BEGIN" ECI:=1; CHANGE OF INFORMATION;
ONE LARGE STEP; T:=T+H2; ECI:=2;
"END";
ONE LARGE STEP;
BACKWARD DIFFERENCES;
UPDATE(WEIGHTS,S2,N);
ERROR ESTIMATES;
"IF" ECI<4 "AND" LQ>80 "THEN" LQ:=20;
HALV:=TWO:="FALSE";
"IF" PRESCH "THEN" "GOTO" TRYCK;
"IF" LQ<1 "THEN" REJECT "ELSE" STEPSIZE;
TRYCK: "IF" TP<=T "THEN" CONTROL(TP,T,H,HNEW,R,ERR,N);
"IF" START "THEN" START:="FALSE";
"IF" HNEW=H2 "THEN" T:=T+H2; ECI:=ECI+1;
"IF" T<TEND+H2 "THEN" "GOTO" MSTP "ELSE" "GOTO" END;
MISS: FAIL:=PRESCH;
"IF" ^ FAIL "THEN"
"BEGIN" "IF" ECI>1 "THEN" T:=T-H2;
HALV:=TWO:="FALSE"; HNEW:=H2/2;
"IF" START "THEN" "GOTO" INIT "ELSE" "GOTO" TRYCK
"END";
END:
"END" IMPEX;
"EOP"
1SECTION : 5.2.1.2.2.1.2 (DECEMBER 1979) PAGE 1
AUTHORS: T.M.T.COOLEN AND R.PLOEGER.
INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM.
RECEIVED: 740301.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS TWO PROCEDURES :
RICHARDSON SOLVES A SYSTEM OF LINEAR EQUATIONS WITH A COEFFICIENT
MATRIX HAVING POSITIVE REAL EIGENVALUES BY MEANS OF A NON-
STATIONARY SECOND ORDER ITERATIVE METHOD: RICHARDSON'S METHOD.
SINCE RICHARDSON'S METHOD IS PARTICULARLY SUITABLE FOR SOLVING
A SYSTEM OF LINEAR EQUATIONS THAT IS OBTAINED BY DISCRETIZING A
TWO-DIMENSIONAL ELLIPTIC BOUNDARY VALUE PROBLEM, THE PROCEDURE
RICHARDSON IS PROGRAMMED IN SUCH A WAY THAT THE SOLUTION VECTOR
IS GIVEN AS A TWO-DIMENSIONAL ARRAY U[J,L], LJ<=J<=UJ, LL<=L<=UL.
THE COEFFICIENT MATRIX IS NOT STORED, BUT EACH ROW CORRESPONDING
TO A PAIR (J,L) IS GENERATED WHEN NEEDED.
RICHARSON CAN ALSO BE USED TO DETERMINE THE EIGENVALUE OF THE
COEFFICIENT MATRIX CORRESPONDING TO THE DOMINANT EIGENFUNCTION.
ELIMINATION, USED IN CONNECTION WITH THE PROCEDURE RICHARDSON,
(THIS SECTION) SOLVES A SYSTEM OF LINEAR EQUATIONS WITH
A COEFFICIENT MATRIX HAVING POSITIVE REAL
EIGENVALUES BY MEANS OF A NON-STATIONARY SECOND ORDER ITERATIVE
METHOD, WHICH IS AN ACCELERATION OF RICHARDSON'S METHOD. IN
GENERAL, ELIMINATION CANNOT BE USED BY ITSELF IN A SENSIBLE WAY.
SINCE RICHARDSON'S METHOD AND ITS ACCELERATION ARE PARTICULARLY
SUITABLE FOR SOLVING A SYSTEM OF LINEAR EQUATIONS THAT IS OBTAINED
BY DISCRETIZING A TWO-DIMENSIONAL ELLIPTIC BOUNDARY VALUE PROBLEM,
THE PROCEDURES RICHARDSON AND ELIMINATION ARE PROGRAMMED IN SUCH A
WAY THAT THE SOLUTION VECTOR IS GIVEN AS A TWO-DIMENSIONAL ARRAY
U[J,L], LJ<=J<=UJ, LL<=L<=UL. THE COEFFICIENT MATRIX IS NOT STORED,
BUT EACH ROW CORSPONDING TO A PAIR (J,L) IS GENERATED WHEN NEEDED.
KEYWORDS:
DIFFERENTIAL EQUATION,
TWO-DIMENSIONAL BOUNDARY VALUE PROBLEM,
SYSTEM OF LINEAR EQUATIONS,
COEFFICIENT MATRIX HAVING POSITIVE REAL EIGENVALUES,
NON-STATIONARY SECOND ORDER ITERATIVE METHOD,
RICHARDSON'S METHOD.
ACCELERATION OF RICHARDSON'S METHOD.
1SECTION : 5.2.1.2.2.1.2 (DECEMBER 1979) PAGE 2
SUBSECTION : RICHARDSON.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" RICHARDSON(U,LJ,UJ,LL,UL,INAP,RESIDUAL,A,B,N,DISCR,K,
RATECONV,DOMEIGVAL,OUT);
"VALUE" LJ,UJ,LL,UL,A,B;
"INTEGER" N,K,LJ,UJ,LL,UL;
"REAL" A,B,RATECONV,DOMIGVAL;
"BOOLEAN" INAP;
"ARRAY" U,DISCR;
"PROCEDURE" RESIDUAL, OUT;
"CODE" 33170;
THE MEANING OF THE FORMAL PARAMETERS IS:
U: <ARRAY IDENTIFIER>;
"ARRAY" U[LJ:UJ,LL:UL];
AFTER EACH ITERATION THE APPROXIMATE SOLUTION CALCULATED BY
THE PROCEDURE RICHARDSON IS STORED INTO U.
ENTRY: IF INAP IS CHOSEN TO BE "TRUE" THEN AN INITIAL
APPROXIMATION OF THE SOLUTION, OTHERWISE ARBITRARY;
EXIT: THE FINAL APPROXIMATION OF THE SOLUTION;
LJ,UJ: <ARITHMETIC EXPRESSION>;
LOWER AND UPPER BOUND FOR THE FIRST SUBSCRIPT OF U;
LL,UL: <ARITHMETIC EXPRESSION>;
LOWER AND UPPER BOUND FOR THE SECOND SUBSCRIPT OF U;
INAP: <BOOLEAN EXPRESSION>;
IF THE USER WISHES TO INTRODUCE AN INITIAL APPROXIMATION
INAP="TRUE" SHOULD BE CHOSEN; THE CHOICE INAP="FALSE" HAS
THE EFFECT THAT ALL COMPONENTS OF U ARE SET EQUAL TO 1
BEFORE THE FIRST ITERATION IS PERFORMED;
RESIDUAL: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS :
"PROCEDURE" RESIDUAL(U); "ARRAY" U;
SUPPOSE THAT THE SYSTEM OF EQUATIONS AT HAND IS AU= F;
FOR ANY ENTRY U THE PROCEDURE RESIDUAL SHOULD CALCULATE
THE RESIDUAL AU - F IN EACH POINT J,L, WHERE
LJ<=J<=UJ, LL<=L<=UL, AND SUBSTITUTE THESE VALUES IN THE
ARRAY U;
A,B: <ARITHMETIC EXPRESSION>;
IF ONE WISHES TO FIND THE SOLUTION OF THE BOUNDARY VALUE
PROBLEM, IN A AND B THE USER SHOULD GIVE A LOWER AND
UPPER BOUND FOR THE EIGENVALUES FOR WHICH THE CORRESPONDING
EIGENFUNCTIONS IN THE EIGENFUNCTION EXPANSION OF THE RESIDU
AL AU - F, WITH U = THE INITIAL APPROXIMATION, SHOULD BE
REDUCED; IF THE DOMINANT EIGENVALUE IS TO BE FOUND, ONE
SHOULD CHOOSE A GREATER THAN THIS EIGENVALUE
(SEE METHOD AND PERFORMANCE);
1SECTION : 5.2.1.2.2.1.2 (DECEMBER 1979) PAGE 3
N: <ARITHMETIC EXPRESSION>;
N GIVES THE TOTAL NUMBER OF ITERATIONS TO BE PERFORMED; THE
VALUE OF N SHOULD EITHER BE GIVEN, OR MADE DEPENDENT OF
SOME JENSEN PARAMETER; E.G. K AND RATECONV CAN SERVE
FOR THIS PURPOSE;
DISCR: <ARRAY IDENTIFIER>;
"ARRAY" DISCR[1:2];
AFTER EACH ITERATION THE PROCEDURE RICHARDSON DELIVERS
IN DISCR[1] THE EUCLIDEAN NORM OF THE RESIDUAL, AND
IN DISCR[2] THE MAXIMUM NORM OF THE RESIDUAL;
K: <VARIABLE>
K COUNTS THE NUMBER OF ITERATIONS RICHARDSON IS PERFORMING;
IT CAN SERVE AS A JENSEN PARAMETER FOR N AND OUT;
RATECONV: <VARIABLE>;
AFTER EACH ITERATION THE AVERAGE RATE OF CONVERGENCE IS
ASSIGNED TO RATECONV;
DOMEIGVAL: <VARIABLE>;
AFTER EACH ITERATION THE VALUE OF THE DOMINANT EIGENVALUE,
IF PRESENT, IS ASSIGNED TO DOMEIGVAL; IF THERE IS NO
DOMINANT EIGENVALUE, THE VALUE OF DOMEIGVAL IS MEANINGLESS,
WHICH MANIFESTS ITSELF BY SHOWING NO CONVERGENCE TO A
FIXED VALUE;
OUT: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE, TO BE WRITTEN BY THE USER,
READS :
"PROCEDURE" OUT(K); "VALUE" K; "INTEGER"K;
BY THIS PROCEDURE ONE HAS ACCESS TO THE FOLLOWING
QUANTITIES:
FOR 0<=K<=N THE K-TH ITERAND IN U,THE EUCLIDEAN AND
MAXIMUM NORM OF THE K-TH RESIDUAL IN DISCR[1] AND DISCR[2],
RESPECTIVELY;
FOR 0<K<=N ALSO THE AVERAGE RATE OF CONVERGENCE AND THE
APPROXIMATION TO THE DOMINANT EIGENVALUE, BOTH WITH RESPECT
TO THE K-TH ITERAND U, IN RATECONV AND DOMEIGVAL,
RESPECTIVELY;
MOREOVER, OUT CAN BE USED TO LET N BE DEPENDENT ON THE
ACCURACY REACHED IN APPROXIMATING THE DOMINANT EIGENVALUE.
DATA AND RESULTS: SEE REF[1],[2].
PROCEDURES USED: NONE.
REQUIRED CENTRAL MEMORY:
APPROXIMATELY 3*(UJ - LJ + 1) * (UL - LL + 1) WORDS.
1SECTION : 5.2.1.2.2.1.2 (DECEMBER 1979) PAGE 4
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
SUPPOSE THE SYSTEM OF EQUATIONS TO BE SOLVED READS AU = F, WHERE
A IS A MATRIX HAVING POSITIVE REAL EIGENVALUES. DENOTING THE
K-TH ITERATE BY U(K), U(K) BEING THE VECTOR U(K)[J,L], LJ<=J<=UJ,
LL<=L<=UL, THE SO-CALLED RESIDUAL WITH RESPECT TO THE K-TH ITERATE
IS DEFINED BY
R(K) = AU(K) - F.
A SECOND ORDER NON-STATIONARY ITERATIVE METHOD IS GIVEN BY
U(K+1) = BETA K * U(K) + (1 - BETA K) * U(K-1)
- OMEGA K * R(K),
OR, EQUIVALENTLY, IF U IS THE (UNKNOWN) EXACT SOLUTION OF AU = F,
U(K) - U = PK(A) (U(0) - U),
WHERE PK DENOTES A POLYNOMIAL OF DEGREE K. RICHARDSON'S METHOD
CONSISTS OF CHOOSING THIS POLYNOMIAL IN SUCH A WAY THAT AMONGST ALL
POLYNOMIALS PK(X) OF DEGREE K WITH PK(0)= 1 IT HAS MINIMAL MAXIMUM
NORM OVER THE INTERVAL [C,D], WHERE C > 0 SHOULD BE CHOSEN TO BE A
LOWER BOUND, AND D AN UPPER BOUND FOR THE EIGENVALUES OF A.
APPLICATION OF THIS POLYNOMIAL TO THE INITIAL ERROR U(0) - U HAS
THE EFFECT THAT EACH COMPONENT OF THE INITIAL ERROR IN ITS EIGEN-
FUNCTION EXPANSION IS REDUCED BY A FACTOR LESS OR EQUAL TO THE NORM
OF THE POLYNOMIAL.
THE POLYNOMIALS
PK(X) = CK((A+B-2*X)/(A-B)) / CK((A+B)/(A-B))
WHERE CK(Y) DENOTES THE K-TH CHEBYSHEV POLYNOMIAL, HAVE THE
DESIRED PROPERTIES. THUS, THE VALUES OF THE PARAMETERS BETA K
AND OMEGA K MAY BE DETERMINED FROM THE RECURRENCE RELATIONS FOR
CHEBESHEV POLYNOMIALS.
IN COMPUTATION U(K) - U IS NOT AVAILABLE, SO ONE USES R(K) AS
A MEASURE FOR THE ERROR.
THE ELEMENTS OF THE MATRIX A ARE NOT STORED, BUT
GENERATED WHEN NEEDED. MORE PRECISELY, THIS MEANS THAT THE
(UJ-LJ+1) * (UL-LL+1) COMPONENTS OF AU(K) - F ARE CALCULATED FOR
EACH PAIR (J,L) LJ<J<UJ, LL<L<UL. THE USER SHOULD INTRODUCE THE
EQUATION TO BE SOLVED IN THIS MANNER BY MEANS OF THE PROCEDURE
RESIDUAL.
CLEARLY, THE METHOD IS PARTICULARLY SUITABLE FOR SPARSE MATRICES,
FOR EXAMPLE MATRICES THAT ARE OBTAINED BY DISCRETIZING ELLIPTIC
PARTIAL DIFFERENTIAL EQUATIONS.
THE SHARPER THE BOUNDS C AND D FOR THE EIGENVALUES 0F A ARE,
THE BETTER APPROXIMATE SOLUTION ONE GETS FOR A GIVEN VALUE OF K,
SINCE THE ASYMPTOTIC RATE OF CONVERGENCE (K TO INFINITY) IS
2 * SQRT(C/D).
NOW LET ALPHA1 BE THE SMALLEST EIGENVALUE OF A. IF ONE CHOOSES
C > ALPHA1, THEN, STARTING WITH ANY INITIAL APPROXIMATION, FOR A
SUFFICIENTLY LARGE NUMBER OF ITERATIONS THE PROCEDURE RICHARDSON
WILL DELIVER AN APPROXIMATE VALUE FOR THIS EIGENVALUE.
1SECTION : 5.2.1.2.2.1.2 (OCTOBER 1974) PAGE 5
LET US EXPLAIN THIS FACT FOR THE CASE ALPHA1 < C < ALPHA2, WHERE
ALPHA2 IS THE SECOND SMALLEST EIGENVALUE OF A. THE POLYNOMIAL
PK(X) HAS SMALL MAXIMUM VALUE OVER THE INTERVAL [C,D] (WHICH, OF
COURSE, DEPENDS ON K), BUT BECOMES LARGE FOR X < A. SO, IF ONE
APPLIES PK(A) TO AN EIGENFUNCTION OF A, THIS EIGENFUNCTION WILL
ONLY BE REDUCED CONSIDERABLY IF IT CORRESPONDS TO AN EIGENVALUE
> C. CONSEQUENTLY, THE EIGENFUNCTION CORRESPONDING TO ALPHA1 WILL
BECOME DOMINANT IN THE EIGENFUNCTION EXPANSION OF
PK(A) (U(0) - U)
FOR SUFFICIENTLY LARGE K.
SEE REF[1],[2] FOR DETAILS.
REFERENCES:
[1].T.M.T.COOLEN, P.W.HEMKER, P.J.VAN DER HOUWEN AND
E.SLAGT.
ALGOL 60 PROCEDURES FOR INITIAL AND BOUNDARY VALUE PROBLEMS
(DUTCH).
MC-SYLLABUS 20, MATHEMATICAL CENTRE, 1973, AMSTERDAM.
[2].P.J.VAN DER HOUWEN.
FINITE DIFFERENCE METHODS FOR SOLVING PARTIAL DIFFERENTIAL
EQUATIONS.
MATHEMATICAL CENTRE TRACT NO. 20, 1968.
EXAMPLE OF USE:
THE APPROXIMATE SOLUTION OF THE BOUNDARY VALUE PROBLEM
- ((D/DX)**2 + (D/DY)**2) U(X,Y) = -2*(X*X+Y*Y), O<X,Y<PI,
U(X,0) = 0, U(X,PI) = PI*PI*X*X, 0 < X < PI,
U(0,Y) = 0, U(PI,Y) = PI*PI*X*X, 0 < Y < PI,
WHICH HAS THE ANALYTICAL SOLUTION X*X*Y*Y, MAY BE OBTAINED BY THE
FOLLOWING PROGRAM:
"BEGIN" "COMMENT" DIRICHLET PROBLEM FOR LAPLACE'S EQUATION;
"PROCEDURE" RICHARDSON(U,LJ,UJ,LL,UL,INAP,RESIDUAL,A,B,N,DISCR,K,
RATECONV,DOMEIGVAL,OUT); "CODE"33170;
"PROCEDURE" RESIDUAL(U); "ARRAY" U;
"BEGIN" "INTEGER" UJMIN1,ULMIN1,LJPLUS1;
"REAL" U2; "REAL" "ARRAY" U1[LJ:UJ];
UJMIN1:= UJ - 1; ULMIN1 := UL - 1; LJPLUS1:= LJ + 1;
"FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
"BEGIN" U1[J]:= U[J,LL]; U[J,LL]:= 0; "END";
1SECTION : 5.2.1.2.2.1.2 (DECEMBER 1979) PAGE 6
"FOR" L:= LL + 1 "STEP" 1 "UNTIL" ULMIN1 "DO"
"BEGIN" U1[LJ]:= U[LJ,L]; U[LJ,L]:= 0;
"FOR" J:= LJPLUS1"STEP" 1 "UNTIL" UJMIN1 "DO"
"BEGIN" U2:= U[J,L];
U[J,L]:=(4 * U2 - U1[J-1] - U1[J] - U[J+1,L] - U[J,L+1])
- F(J*H,L*H)*H2;
U1[J]:= U2
"END";
U[UJ,L]:= 0;
"END";
"FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO" U[J,UL]:= 0
"END" RESIDUAL;
"REAL" "PROCEDURE" F(X,Y); "VALUE" X,Y; "REAL" X,Y;
F:= -2*(X*X + Y*Y);
"REAL" "PROCEDURE" ANALSOL(X,Y); "VALUE" X,Y; "REAL" X,Y;
ANALSOL:= X*X*Y*Y;
"PROCEDURE" INITAPPR(U,J,L,G); "INTEGER" J,L; "ARRAY" U; "REAL" G;
"FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
"FOR" L:= LL "STEP" 1 "UNTIL" UL "DO"
U[J,L]:= "IF" J=LJ "OR" J=UJ "OR" L=LL "OR" L=UL "THEN"G "ELSE" 1;
"PROCEDURE"OUT1(K); "VALUE" K; "INTEGER" K;
"IF" K = N "THEN" OUTPUT(61,"("//"(" K DISCR[1] DISCR[2]
RATECONV")",//,+ZDB,3(+.7D"+ZDB)")",K,DISCR[1],DISCR[2],RATECONV);
"INTEGER" J,L,LJ,UJ,LL,UL,N,K;
"REAL" H,PI,D1,D2,H2,DOMEIGVAL,RATECONV,A,B;
"REAL" "ARRAY" DISCR[1:2];
OUTPUT(61,"("/"("GIVE LJ,UJ,LL,UL,N,A,B")"/")");
INPUT(60,"("")",LJ); INPUT(60,"("")",UJ);
INPUT(60,"("")",LL); INPUT(60,"("")",UL);
INPUT(60,"("")",N); INPUT(60,"("")",A); INPUT(60,"("")",B);
"BEGIN" "REAL" "ARRAY" U[LJ:UJ,LL:UL];
PI:=3.1415 92653 58979; H:= PI/(UJ - LJ); H2:= H * H;
INITAPPR(U,J,L,ANALSOL(J*H,L*H));
RICHARDSON(U,LJ,UJ,LL,UL,"TRUE",RESIDUAL,A,B,N,DISCR,K,
RATECONV ,DOMEIGVAL,OUT1);
"END"
"END"
IT DELIVERS WITH
LJ = 0, UJ = 11, LL = 0, UL = 11, N = 50, A = .163, B = 7.83
THE FOLLOWING RESULTS:
K DISCR[1] DISCR[2] RATECONV
+50 +.1401828" -3 +.4666866" -4 +.2921718" +0 .
1SECTION : 5.2.1.2.2.1.2 (DECEMBER 1979) PAGE 7
SUBSECTION : ELIMINATION.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" ELIMINATION(U,LJ,UJ,LL,UL,RESIDUAL,A,B,N,DISCR,K,
RATECONV,DOMEIGVAL,OUT);
"VALUE" LJ,UJ,LL,UL,A,B;
"INTEGER" N,K,LJ,UJ,LL,UL;
"REAL" A,B,RATECONV,DOMIGVAL;
"ARRAY" U,DISCR;
"PROCEDURE" RESIDUAL, OUT; "CODE" 33171;
THE MEANING OF THE FORMAL PARAMETERS IS:
U: <ARRAY IDENTIFIER>;
"ARRAY" U[LJ:UJ,LL:UL];
AFTER EACH ITERATION THE APPROXIMATE SOLUTION CALCULATED BY
THE PROCEDURE ELIMINATION IS STORED INTO U;
ENTRY: AN INITIAL APPROXIMATION OF THE SOLUTION, WHICH
IS OBTAINED BY USE OF RICHARDSON;
EXIT: THE FINAL APPROXIMATION OF THE SOLUTION;
LJ,UJ: <ARITHMETIC EXPRESSION>;
LOWER AND UPPER BOUND FOR THE FIRST SUBSCRIPT OF U;
LL,UL: <ARITHMETIC EXPRESSION>;
LOWER AND UPPER BOUND FOR THE SECOND SUBSCRIPT OF U;
RESIDUAL: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE READS :
"PROCEDURE" RESIDUAL(U); "ARRAY" U;
SUPPOSE THAT THE SYSTEM OF EQUATIONS AT HAND IS AU= F;
FOR ANY ENTRY U THE PROCEDURE RESIDUAL SHOULD CALCULATE
THE SO-CALLED RESIDUAL AU - F IN EACH POINT J,L, WHERE
LJ<=J<=UJ, LL<=L<=UL, AND SUBSTITUTE THESE VALUES IN THE
ARRAY U;
A,B: <ARITHMETIC EXPRESSION>;
A AND B SHOULD HAVE THE SAME VALUES AS IN THE PRECEDING
CALL OF RICHARDSON (SEE DESCRIPTION OF RICHARDSON);
N: <VARIABLE>;
THE NUMBER OF ITERATIONS THE PROCEDURE ELIMINATION NEEDS
TO ELIMINATE THE EIGENFUNCTION BELONGING TO THE DOMINANT
EIGENVALUE, IS ASSIGNED TO N;
DISCR: <ARRAY IDENTIFIER>;
"ARRAY" DISCR[1:2];
AFTER EACH ITERATION THE PROCEDURE ELIMINATION DELIVERS
IN DISCR[1] THE EUCLIDEAN NORM OF THE RESIDUAL, AND
IN DISCR[2] THE MAXIMUM NORM OF THE RESIDUAL;
K: <VARIABLE>
K COUNTS THE NUMBER OF ITERATIONS ELIMINATION IS PERFORMING
IT CAN SERVE AS A JENSEN PARAMETER FOR OUT;
RATECONV: <VARIABLE>;
AFTER EACH ITERATION THE AVERAGE RATE OF CONVERGENCE IS
ASSIGNED TO RATECONV;
1SECTION : 5.2.1.2.2.1.2 (DECEMBER 1979) PAGE 8
DOMEIGVAL: <ARITHMETIC EXPRESSION>;
BEFORE A CALL OF ELIMINATION THE VALUE OF THE EIGENVALUE
FOR WHICH THE CORRESPONDING EIGENFUNCTION HAS TO BE
ELIMINATED, SHOULD BE ASSIGNED TO DOMEIGVAL; IF AFTER
APPLICATION OF ELIMINATION THERE IS A NEW DOMINANT EIGEN-
FUNCTION, THEN DOMEIGVAL WILL BE EQUAL TO THE CORRESPOND-
ING EIGENVALUE; OTHERWISE, THE VALUE OF DOMEIGVAL BECOMES
MEANINGLESS;
OUT: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE, TO BE WRITTEN BY THE USER,
READS :
"PROCEDURE" OUT(K); "VALUE" K; "INTEGER"K;
BY THIS PROCEDURE ONE HAS ACCESS TO THE FOLLOWING
QUANTITIES:
FOR 0<=K<=N THE K-TH ITERAND IN U,THE EUCLIDEAN AND
MAXIMUM NORM OF THE K-TH RESIDUAL IN DISCR[1] AND DISCR[2],
RESPECTIVELY;
FOR 0<K<=N ALSO THE AVERAGE RATE OF CONVERGENCE WITH
RESPECT TO THE K-TH ITERAND U, IN RATECONV;
FOR K = N, POSSIBLY THE DOMINANT EIGENVALUE OF THE
COEFFICIENT MATRIX OF THE EQUATION AU= F, IN DOMEIGVAL.
DATA AND RESULTS: SEE REF[1],[2].
PROCEDURES USED:
RICHARDSON = CP33170,
TAN = CP35120,
TANH = CP35113,
ARCCOS = CP35122,
ZEROIN = CP34150.
REQUIRED CENTRAL MEMORY:
APPROXIMATELY 3*(UJ - LJ + 1) * (UL - LL + 1) WORDS.
1SECTION : 5.2.1.2.2.1.2 (DECEMBER 1979) PAGE 9
METHOD AND PERFORMANCE:
SEE THIS HEADING IN THE DESCRIPTION OF THE PROCEDURE RICHARDSON.
SOME ADDITIONAL REMARKS WILL BE MADE HERE.
IN ORDER TO USE ELIMINATION THE INITIAL APPROXIMATION OF THE
SOLUTION OF
AU = F
IS FIRST TREATED BY MEANS OF RICHARDSON'S METHOD, WHERE C IS
CHOSEN GREATER THAN THE SMALLEST EIGENVALUE. AFTER APPLICATION OF
RICHARDSON, THE EIGENFUNCTION CORRESPONDING TO THIS EIGENVALUE HAS
BECOME DOMINANT IN THE QUANTITY
PK(A) (U(0) - U),
WITH
PK(X) = CK((C+D-2*X)/(C-D)) / CK((C+D)/(C-D)),
WHEREAS THE CONTRIBUTION OF THE OTHER EIGENFUNCTIONS TO THE ERROR
U(K) - U AND TO R(K) HAS BEEN REDUCED CONSIDERABLY. CONSEQUENTLY
THE ERROR U(K) - U HAS VERY SMALL COMPONENTS IN THE SUBSPACE
SPANNED BY ALL EIGENVECTORS BUT THE "FIRST", IN WHICH DIRECTION IT
HAS A VERY LARGE COMPONENT.
THE CONTRIBUTION OF THE "FIRST" EIGENFUNCTION TO R(K) IS NOW
"ELIMINATED" BY APPLICATION OF A POLYNOMIAL OPERATOR E(A) SUCH
THAT E(X) HAS A ZERO IN THE FIRST EIGENVALUE.
THE POLYNOMIAL IS CHOSEN IN SUCH A WAY THAT A MAXIMAL RATE OF CON-
VERGENCE WITH RESPECT TO THE INITIAL APPROXIMATION USED IN
RICHARDSON IS OBTAINED.
FOR DETAILS SEE REF[1],[2].
REFERENCES:
[1].T.M.T.COOLEN, P.W.HEMKER, P.J.VAN DER HOUWEN AND
E.SLAGT.
ALGOL 60 PROCEDURES FOR INITIAL AND BOUNDARY VALUE PROBLEMS
(DUTCH).
MC-SYLLABUS 20, MATHEMATICAL CENTRE, 1976, AMSTERDAM.
[2].P.J.VAN DER HOUWEN.
FINITE DIFFERENCE METHODS FOR SOLVING PARTIAL DIFFERENTIAL
EQUATIONS.
MATHEMATICAL CENTRE TRACT NO. 20, 1968.
EXAMPLE OF USE:
THE APPROXIMATE SOLUTION OF THE BOUNDARY VALUE PROBLEM
- ((D/DX)**2 + (D/DY)**2) U(X,Y) = -2*(X*X+Y*Y), O<X,Y<PI,
U(X,0) = 0, U(X,PI) = PI*PI*X*X, 0 < X < PI,
U(0,Y) = 0, U(PI,Y) = PI*PI*X*X, 0 < Y < PI,
WHICH HAS THE ANALYTICAL SOLUTION X*X*Y*Y, MAY BE OBTAINED BY THE
FOLLOWING PROGRAM:
1SECTION : 5.2.1.2.2.1.2 (OCTOBER 1974) PAGE 10
"BEGIN" "COMMENT" DIRICHLET PROBLEM FOR LAPLACE'S EQUATION;
"PROCEDURE" RICHARDSON(U,LJ,UJ,LL,UL,INAP,RESIDUAL,A,B,DISCR,K,
RATECONV,DOMEIGVAL,OUT); "CODE"33170;
"PROCEDURE" ELIMINATION(U,LJ,UJ,LL,UL,RESIDUAL,A,B,DISCR,K,
RATECONV,DOMEIGVAL,OUT); "CODE"33171;
"PROCEDURE" RESIDUAL(U); "ARRAY" U;
"BEGIN" "INTEGER" UJMIN1,ULMIN1,LJPLUS1;
"REAL" U2; "REAL" "ARRAY" U1[LJ:UJ];
UJMIN1:= UJ - 1; ULMIN1 := UL - 1; LJPLUS1:= LJ + 1;
"FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
"BEGIN" U1[J]:= U[J,LL]; U[J,LL]:= 0; "END";
"FOR" L:= LL + 1 "STEP" 1 "UNTIL" ULMIN1 "DO"
"BEGIN" U1[LJ]:= U[LJ,L]; U[LJ,L]:= 0;
"FOR" J:= LJPLUS1"STEP" 1 "UNTIL" UJMIN1 "DO"
"BEGIN" U2:= U[J,L];
U[J,L]:=(4 * U2 - U1[J-1] - U1[J] - U[J+1,L] - U[J,L+1])
- F(J*H,L*H)*H2;
U1[J]:= U2
"END";
U[UJ,L]:= 0;
"END";
"FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO" U[J,UL]:= 0
"END" RESIDUAL;
"REAL" "PROCEDURE" F(X,Y); "VALUE" X,Y; "REAL" X,Y;
F:= -2*(X*X + Y*Y);
"REAL" "PROCEDURE" ANALSOL(X,Y); "VALUE" X,Y; "REAL" X,Y;
ANALSOL:= X*X*Y*Y;
"PROCEDURE" INITAPPR(U,J,L,G); "INTEGER" J,L; "ARRAY" U; "REAL" G;
"FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
"FOR" L:= LL "STEP" 1 "UNTIL" UL "DO"
U[J,L]:= "IF" J=LJ "OR" J=UJ "OR" L=LL "OR" L=UL "THEN"G "ELSE" 1;
"PROCEDURE"OUT3(K); "VALUE" K; "INTEGER" K;
"IF" K=P "THEN" OUTPUT(61,"("//,+ZDB,3(+.7D"+ZDB)")",K,DISCR[1],
DISCR[2],RATECONV);
"PROCEDURE"OUT1(K); "VALUE" K; "INTEGER" K;
"IF" K=N "THEN" OUTPUT(61,"("//"(" K DISCR[1] DISCR[2] R
TECONV")",//,+ZDB,3(+.7D"+ZDB)")",K,DISCR[1],DISCR[2],RATECONV);
"PROCEDURE" OUT2(K); "VALUE" K; "INTEGER" K;
"BEGIN"
"IF" K = 0 "THEN" D1:= D2:= 1 "ELSE"
"BEGIN" D2:= D1; D1:= DOMEIGVAL;
N:= "IF" ABS((D1 - D2)/D2) < 10**(-Q) "THEN" K "ELSE" NN;
OUT1(K)
"END"
"END" OUT2;
1SECTION : 5.2.1.2.2.1.2 (OCTOBER 1974) PAGE 11
"INTEGER" J,L,LJ,UJ,LL,UL,NN,N,P,K,Q;
"REAL" H,PI,D1,D2,H2,RATECONVR,RATECONVE,DOMEIGVAL,RATECONV,A,B,VAR;
"REAL" "ARRAY" DISCR[1:2];
OUTPUT(61,"("/"("GIVE LJ,UJ,LL,UL,N,Q,A,B")"/")");
INPUT(60,"("")",LJ); INPUT(60,"("")",UJ);
INPUT(60,"("")",LL); INPUT(60,"("")",UL);
INPUT(60,"("")", N); INPUT(60,"("")", Q);
INPUT(60,"("")", A); INPUT(60,"("")", B);
"BEGIN" "REAL" "ARRAY" U[LJ:UJ,LL:UL];
PI:=3.1415 92653 58979; H:= PI/(UJ - LJ); H2:= H * H;
INITAPPR(U,J,L,ANALSOL(J*H,L*H));
NN:= N;
RICHARDSON(U,LJ,UJ,LL,UL,"TRUE",RESIDUAL,A,B,N,DISCR,K,
RATECONV ,DOMEIGVAL,OUT2); RATECONVR:= RATECONV;
OUTPUT(61,"("//+.7D"+ZD4B"("DOMINANT EIGENVALUE")"")",DOMEIGVAL);
ELIMINATION(U,LJ,UJ,LL,UL,RESIDUAL,A ,B,P,DISCR,K,
RATECONV ,DOMEIGVAL,OUT3); RATECONVE:= RATECONV;
NN:= N + P; OUTPUT(61,"("//+Z2D13B"("TOTAL NUMBER OF ITERATIONS")"
")",NN);
OUTPUT(61,"("/+.7D"+ZD4B"("RATE OF CONVERGENCE WITH RESPECT TO")",
/17B"("THE ZEROTH ITERAND OF RICHARDSON")"")",
(N * RATECONVR + P * RATECONVE)/NN);
"END"
"END"
IT DELIVERS WITH
LJ = 0, UJ = 11, LL = 0, UL = 11, N = 50, Q = 4, A = .326, B = 7.83
THE FOLLOWING RESULTS:
K DISCR[1] DISCR[2] RATECONV
+45 +.4998463" -1 +.8903863" -2 +.2009943" +0
+.1620445" +0 DOMINANT EIGENVALUE
+7 +.3563865" -5 +.6714375" -6 +.1360086" +1
+52 TOTAL NUMBER OF ITERATIONS
+.3570259" +0 RATE OF CONVERGENCE WITH RESPECT TO
THE ZEROTH ITERAND OF RICHARDSON
1SECTION : 5.2.1.2.2.1.2 (OCTOBER 1974) PAGE 12
SOURCE TEXT(S):
"CODE"33170;
"PROCEDURE" RICHARDSON(U,LJ,UJ,LL,UL,INAP,RESIDUAL,A,B,N,DISCR,K,
RATECONV,DOMEIGVAL,OUT); "VALUE" LJ,UJ,LL,UL,A,B;
"INTEGER" N,K,LJ,UJ,LL,UL; "REAL" A,B,RATECONV,DOMEIGVAL; "BOOLEAN"
INAP; "ARRAY" U,DISCR; "PROCEDURE" RESIDUAL,OUT;
"BEGIN" "INTEGER" J,L; "REAL" X,Y,Z,Y0,C,D,ALFA,OMEGA,OMEGA0,
EIGMAX,EIGEUCL,EUCLRES,MAXRES,RCMAX,RCEUCL,MAXRES0,EUCLRES0;
"ARRAY" V,RES[LJ:UJ,LL:UL];
"PROCEDURE" CALPAR;
"COMMENT" CALPAR CALCULATES THE PARAMETERS ALFA AND OMEGA FOR
EACH ITERATION;
"BEGIN" ALFA:= Z/(Z - ALFA);
OMEGA:= 1/(X - OMEGA * Y)
"END" CALPAR;
"PROCEDURE" ITERATION;
"COMMENT" FIRST THE ITERATION FORMULA IS CONSTRUCTED;
"BEGIN" "REAL" AUXV,AUXU,AUXRES,EUCLUV,MAXUV;
EUCLUV:= EUCLRES:= MAXUV:= MAXRES:= 0;
"FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
"FOR" L:= LL "STEP" 1 "UNTIL" UL "DO" RES[J,L]:= V[J,L];
RESIDUAL(RES);
"FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
"FOR" L:= LL "STEP" 1 "UNTIL" UL "DO"
"BEGIN" AUXV:= U[J,L]; AUXU:= V[J,L]; AUXRES:= RES[J,L];
AUXV:= ALFA * AUXU - OMEGA * AUXRES + (1 - ALFA) * AUXV;
V[J,L]:= AUXV; U[J,L]:= AUXU;
"COMMENT" THE NORMS OF THE K-TH RESIDUAL AND THE DIFFERENCE
BETWEEN THE (K+1)-TH AND K-TH ITERAND ARE CALCULATED;
AUXU:= ABS(AUXU - AUXV); AUXRES:= ABS(AUXRES);
MAXUV:= "IF" MAXUV < AUXU "THEN" AUXU "ELSE" MAXUV;
MAXRES:= "IF" MAXRES < AUXRES "THEN" AUXRES "ELSE" MAXRES;
EUCLUV:= EUCLUV + AUXU * AUXU;
EUCLRES:= EUCLRES + AUXRES * AUXRES;
"END";
EUCLUV:= SQRT(EUCLUV); EUCLRES:= SQRT(EUCLRES);
DISCR[1]:= EUCLRES; DISCR[2]:= MAXRES;
"COMMENT" DOMEIGVAL IS EVALUATED;
MAXUV:= MAXRES/MAXUV; EUCLUV:= EUCLRES/EUCLUV;
EIGMAX:= MAXUV * (C - MAXUV)/(.25 * D - MAXUV);
EIGEUCL:= EUCLUV * (C - EUCLUV)/(.25 * D - EUCLUV);
DOMEIGVAL:= .5 * (EIGMAX + EIGEUCL);
"COMMENT" FINALLY THE RATE OF CONVERGENCE IS CALCULATED;
RCEUCL:= -LN(EUCLRES/EUCLRES0)/K;
RCMAX:= -LN(MAXRES/MAXRES0)/K;
RATECONV:= .5 * (RCEUCL + RCMAX)
"END" ITERATION; "COMMENT"
1SECTION : 5.2.1.2.2.1.2 (OCTOBER 1974) PAGE 13
;
"COMMENT" THE CONSTANTS FOR STARTING CALPAR ARE CALCULATED;
ALFA:= 2; OMEGA:= 4/(B + A); Y0:= (B + A)/(B - A);
X:= .5 * (B + A); Y:= (B - A) * (B - A)/16; Z:= 4 * Y0 * Y0;
"COMMENT" THE CONSTANTS NEEDED FOR DOMEIGVAL ARE CALCULATED;
C:= A * B; C:= SQRT(C); D:= SQRT(A) + SQRT(B); D:= D * D;
"COMMENT" THE INITIAL APPROXIMATION IS PUT INTO ARRAY U;
"IF" ^INAP "THEN"
"BEGIN" "FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
"FOR" L:= LL "STEP" 1 "UNTIL" UL "DO" U[J,L]:= 1
"END";
"COMMENT" THE ZEROTH ITERATION IS NOW PERFORMED;
K:= 0;
"FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
"FOR" L:= LL "STEP" 1 "UNTIL" UL "DO" RES[J,L]:= U[J,L];
RESIDUAL(RES);
OMEGA0:= 2/(B+A);
"BEGIN" "REAL" AUXRES0;
MAXRES0:= EUCLRES0:= 0;
"FOR" J:= LJ "STEP" 1 "UNTIL" UJ "DO"
"FOR" L:= LL "STEP" 1 "UNTIL" UL "DO"
"BEGIN" AUXRES0:= RES[J,L];
V[J,L]:= U[J,L] - OMEGA0 * AUXRES0;
AUXRES0:= ABS(AUXRES0);
MAXRES0:= "IF" MAXRES0 < AUXRES0 "THEN" AUXRES0 "ELSE" MAXRES0;
EUCLRES0:= EUCLRES0 + AUXRES0 * AUXRES0
"END";
EUCLRES0:= SQRT(EUCLRES0)
"END";
DISCR[1]:= EUCLRES0; DISCR[2]:= MAXRES0;
OUT(K);
"IF" K >= N "THEN" "GOTO" FINALLY;
NEXT STEP:
K:= K + 1; CALPAR; ITERATION; OUT(K);
"IF" K < N "THEN" "GOTO" NEXT STEP;
FINALLY:
"END" RICHARDSON;
"EOP"
1SECTION : 5.2.1.2.2.1.2 (OCTOBER 1974) PAGE 14
"CODE"33171;
"PROCEDURE" ELIMINATION(U,LJ,UJ,LL,UL,RESIDUAL,A,B,N,DISCR,K,RATECONV
DOMEIGVAL,OUT); "VALUE" LJ,UJ,LL,UL,A,B; "INTEGER" LJ,UJ,LL,UL,N,K;
"REAL" A,B,RATECONV,DOMEIGVAL; "ARRAY" U,DISCR;
"PROCEDURE" RESIDUAL,OUT;
"BEGIN" "REAL" PI,AUXCOS,C,D;
"REAL" "PROCEDURE" ARCCOS(X); "CODE" 35122;
"REAL" "PROCEDURE" TAN(X); "CODE" 35120;
"REAL" "PROCEDURE" TANH(X); "CODE" 35113;
"PROCEDURE" RICHARDSON(U,LJ,UJ,LL,UL,INAP,RESIDUAL,A,B,N,DISCR,
K,RATECONV,DOMEIGVAL,OUT); "CODE"33170;
"BOOLEAN" "PROCEDURE" ZEROIN(X,Y,FX,TOLX); "CODE"34150;
"REAL" "PROCEDURE" OPTPOL(X); "VALUE" X; "REAL" X;
"BEGIN" "REAL" W,Y;
W:= (B * COS(.5*PI/X) + DOMEIGVAL) / (B - DOMEIGVAL);
"IF" W < -1 "THEN" W:= -1;
"IF" ABS(W) <= 1 "THEN"
"BEGIN" Y:= ARCCOS(W);
OPTPOL:= 2 * SQRT(A/B) + TAN(X*Y) *
(Y - B*PI*SIN(.5*PI/X)*.5 / (X * (B-DOMEIGVAL) *
SQRT(ABS(1-W*W))))
"END" "ELSE"
"BEGIN" Y:= LN(W + SQRT(ABS(W*W-1)));
OPTPOL:= 2 * SQRT(A/B) - TANH(X*Y) * (Y + B*PI*SIN(.5*PI/X)*
.5/(X*(B-DOMEIGVAL)*SQRT(ABS(W*W-1))))
"END"
"END" OPTPOL;
PI:= 3.1415 92653 58979;
C:= 1;
"IF" OPTPOL(C) < 0 "THEN"
"BEGIN" D:= .5 * PI * SQRT(ABS(B/DOMEIGVAL));
M: D:= D + D;
"IF" ZEROIN(C,D,OPTPOL(C),C*"-3) "THEN" N:= ENTIER(C+.5)
"ELSE" "GOTO" M;
"END" "ELSE" N:= 1;
AUXCOS:= COS(.5*PI/N);
RICHARDSON(U,LJ,UJ,LL,UL,"TRUE",RESIDUAL,
(2*DOMEIGVAL + B*(AUXCOS-1))/(AUXCOS+1),B,N,DISCR,K,RATECONV,
DOMEIGVAL,OUT)
"END" ELIMINATION;
"EOP"
1SECTION : 6.7 (OCTOBER 1974) PAGE 1
AUTHOR: S.P.N. VAN KAMPEN.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 740410.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS FIVE PROCEDURES:
A) THE PROCEDURE ERRORFUNCTION COMPUTES THE ERROR FUNCTION AND
COMPLEMENTARY ERROR FUNCTION FOR A REAL ARGUMENT, I.E.
ERF(X) = 2 / SQRT(PI) * INTEGRAL FROM 0 TO X OF EXP(-T ** 2)DT
AND
ERFC(X) = 2 / SQRT(PI) * INTEGRAL FROM X TO INFINITY OF
EXP(-T ** 2)DT
= 1 - ERF(X),
(SEE E.G. [1] EQ. 7.1.1 AND 7.1.2);
THESE FORMULAS ARE RELATED TO THE NORMAL OR GAUSSIAN PROBABILITY
FUNCTION:
P(X) = 1 / SQRT(2 * PI) * INTEGRAL FROM - INFINITY TO X OF
EXP(-T ** 2 / 2)DT
= (1 + ERF(X / SQRT(2))) / 2
AND
Q(X) = 1 / SQRT(2 * PI) * INTEGRAL FROM X TO INFINITY OF
EXP(-T ** 2 / 2)DT
= ERFC(X / SQRT(2)) / 2,
(SEE E.G. [1] EQ. 26.2.2, 26.2.3 AND 26.2.29).
B) THE AUXILIARY PROCEDURE NONEXPERFC COMPUTES
EXP(X * X) * ERFC(X).
C) THE PROCEDURE INVERSE ERROR FUNCTION CALCULATES THE INVERSE OF
THE ERROR FUNCTION DEFINED BY:
Y = INVERF(X),
WHERE
X = ERF(Y) =
= 2 / SQRT(PI) * INTEGRAL FROM 0 TO Y OF EXP(-T ** 2) DT,
(SEE THE PROCEDURE ERRORFUNCTION (THIS SECTION) ).
D) THE PROCEDURE FRESNEL CALCULATES THE FRESNEL INTEGRALS C(X) AND
S(X) DEFINED BY
C(X) = INTEGRAL FROM 0 TO X OF COS(PI / 2 * T * T)DT
AND
S(X) = INTEGRAL FROM 0 TO X OF SIN(PI / 2 * T * T)DT
(SEE [1] EQ. 7.3.1 AND 7.3.2);
1SECTION : 6.7 (OCTOBER 1974) PAGE 2
E) THE AUXILIARY PROCEDURE FG CALCULATES F(X) AND G(X) DEFINED BY
F(X) = (0.5 - S(X))COS(PI / 2 * X * X) -
(0.5 - C(X))SIN(PI / 2 * X * X)
AND
G(X) = (0.5 - C(X))COS(PI / 2 * X * X) +
(0.5 - S(X))SIN(PI / 2 * X * X)
(SEE [1] EQ. 7.3.5 AND 7.3.6).
KEYWORDS:
ERROR FUNCTION,
COMPLEMENTARY ERROR FUNCTION,
NORMAL PROBABILITY FUNCTION,
GAUSSIAN PROBABILITY FUNCTION,
FRESNEL INTEGRALS,
INVERSE ERROR FUNCTION.
1SECTION : 6.7 (OCTOBER 1974) PAGE 3
SUBSECTION: ERRORFUNCTION.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS :
"PROCEDURE" ERRORFUNCTION(X, ERF, ERFC);
"VALUE" X; "REAL" X, ERF, ERFC;
THE MEANING OF THE FORMAL PARAMETERS IS :
X: <ARITHMETIC EXPRESSION>;
ENTRY: THE (REAL) ARGUMENT OF ERF(X) AND ERFC(X);
ERF: <VARIABLE>;
EXIT: THE VALUE OF ERF(X),
ERFC: <VARIABLE>;
EXIT: THE VALUE OF ERFC(X).
PROCEDURES USED: NONEXPERFC = CP35022.
RUNNING TIME: ABOUT 0.001 100 SEC.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
SEE METHOD AND PERFORMANCE OF NONEXPERFC (THIS SECTION).
1SECTION : 6.7 (OCTOBER 1974) PAGE 4
SUBSECTION: NONEXPERFC.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS :
"REAL" "PROCEDURE" NONEXPERFC(X); "VALUE" X; "REAL" X;
NONEXPERFC DELIVERS THE VALUE OF EXP(X * X) * ERFC(X);
THE MEANING OF THE FORMAL PARAMETERS IS :
X: <ARITHMETIC EXPRESSION>;
ENTRY: THE (REAL) ARGUMENT OF NONEXPERFC.
PROCEDURES USED: ERRORFUNCTION = CP35021.
RUNNING TIME: ABOUT 0.000 900 SEC.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
IF ABS(X) <= 0.5 THE VALUES OF ERF(X) AND ERFC(X) ARE COMPUTED IN
THE PROCEDURE ERRORFUNCTION BY MEANS OF RATIONAL CHEBYSHEV
APPROXIMATION AS GIVEN IN [2]. ON THIS INTERVAL THE VALUE OF
NONEXPERFC(X) = EXP(X * X) * ERFC(X) IS COMPUTED BY CALLING THE
PROCEDURE ERRORFUNCTION.
IF ABS(X) > 0.5 THE VALUES OF ERF(X) AND ERFC(X) ARE COMPUTED BY
CALLING THE PROCEDURE NONEXPERFC, WHILE THE VALUE OF NONEXPERFC(X)
IS COMPUTED BY MEANS OF RATIONAL CHEBYSHEV APPROXIMATIONS AS GIVEN
IN [2].
THE COMPUTED VALUES OF ERF(X) AND ERFC(X) ARE COMPARED WITH HIGHER
PRECISION VALUES USING 4000 PSEUDO-RANDOM ARGUMENTS. IT APPEARED
THAT ERF(X) IS COMPUTED WITH AN AVERAGE RELATIVE ERROR 1.93"-15
AND A MAXIMUM RELATIVE ERROR 1.35"-14.
IF X < 6 ERFC(X) IS COMPUTED WITH AN AVERAGE RELATIVE ERROR
8.87"-15 AND A MAXIMUM RELATIVE ERROR 1.55"-13.
IF X <= 26 ERFC(X) IS COMPUTED WITH AN AVERAGE RELATIVE ERROR
5.71"-14 AND A MAXIMUM RELATIVE ERROR 2.70"-12.
IF X > 26 ERFC(X)=0, BECAUSE IN THIS CASE ERFC(X) IS LESS THAN THE
SMALLEST REPRESENTABLE POSITIVE NUMBER ON THE CD CYBER 73-28.
FOR THIS REASON IT IS ADVISABLE TO COMPUTE FOR X > 26 NONEXPERFC(X)
INSTEAD OF ERFC(X).
IF X < -26.2 THE PROCEDURE NONEXPERFC WILL BE TERMINATED ABNORMALLY
BY CAUSE OF OVERFLOW.
REFERENCES: SEE REFERENCES [1] AND [2] OF THE PROCEDURE FG (THIS
SECTION).
1SECTION : 6.7 (OCTOBER 1974) PAGE 5
EXAMPLE OF USE:
WE COMPUTE THE VALUES OF
ERF(1) = 0.84270 07929 49714 8693,
ERFC(1) = 0.15729 92070 50285 1307
AND NONEXPERFC(100) =
EXP(100 * 100) * ERFC(100) = 0.5416 13782 89843 2905"-2;
"BEGIN"
"PROCEDURE" ERRORFUNCTION(X, ERF, ERFC); "CODE" 35021;
"REAL" "PROCEDURE" NONEXPERFC(X); "CODE" 35022;
"REAL" ERF, ERFC, P;
ERRORFUNCTION(1, ERF, ERFC);
P:= NONEXPERFC(100);
OUTPUT(61, "(""(" ERF(1) = ")", +D.5DB5DB5D, /,
"(" ERFC(1) = ")", +D.5DB5DB5D, /,
"(" NONEXPERFC(100) = ")", +.5DB5DB5D"+D")",
ERF, ERFC, P);
"END"
THIS PROGRAM DELIVERS:
ERF(1) = +0.84270 07929 49710
ERFC(1) = +0.15729 92070 50280
NONEXPERFC(100) = +.56416 13782 98940"-2.
1SECTION : 6.7 (OCTOBER 1974) PAGE 6
SUBSECTION : INVERSE ERROR FUNCTION.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" INVERSE ERROR FUNCTION(X, ONEMINX, INVERF);
"VALUE" X, ONEMINX; "REAL" X, ONEMINX, INVERF;
THE MEANING OF THE FORMAL PARAMETERS IS :
X: <ARITHMETIC EXPRESSION>;
ENTRY:
THE ARGUMENT OF THE FUNCTION INVERF;
IT IS NECESSARY THAT -1 < X < 1;
IF ABS(X) > 0.8 THE VALUE OF X IS NOT USED IN THE
PROCEDURE;
ONEMINX: <ARITHMETIC EXPRESSION>;
ENTRY:
IF ABS(X) <= 0.8 THE VALUE OF ONEMINX IS NOT USED IN
THE PROCEDURE; IF ABS(X) > 0.8 ONEMINX HAS TO CONTAIN
THE VALUE OF 1 - ABS(X); IN THE CASE THAT ABS(X) IS IN
THE NEIGHBOURHOOD OF 1, CANCELLATION OF DIGITS TAKE
PLACE IN THE CALCULATION OF 1 - ABS(X); IF THE VALUE
1-ABS(X) IS KNOWN EXACTLY FROM ANOTHER SOURCE, ONEMINX
HAS TO CONTAIN THIS VALUE, WHICH WILL GIVE BETTER
RESULTS;
INVERF: <VARIABLE>;
EXIT: THE RESULT OF THE PROCEDURE.
PROCEDURES USED: CHEPOLSER = CP31046.
RUNNING TIME: ABOUT 0.003 800 SEC.
LANGUAGE: ALGOL 60.
1SECTION : 6.7 (OCTOBER 1974) PAGE 7
METHOD AND PERFORMANCE:
THE FUNCTION VALUE INVERF IS CALCULATED ON DIFFERENT INTERVALS BY
MEANS OF CHEBYSHEV POLYNOMIALS, OF WHICH THE COEFFICIENTS ARE GIVEN
IN [1].
ON THE COMPUTED RESULTS WE USED THE TESTS:
EPS1:= ABS(ERF(INVERF(X)) / X - 1),
EPS2:= ABS(INVERF(ERF(Y)) / Y - 1),
EPS3:= ABS((1 - ERF(INVERF(1 - X))) / X - 1).
IF ABS(X) < 0.9 UPPER BOUNDS FOR EPS1 AND EPS2 ARE 7.1"-15 AND
4.1"-14 RESP.
IF 0.9 < ABS(X) < 1 CANCELLATION OF DIGITS TAKE PLACE IN THE
CALCULATION OF 1 - ABS(X). THIS CANCELLED DIGITS ARE ALSO LOST IN
THE RESULT. IF THE VALUE OF 1 - ABS(X) IS KNOWN EXACTLY AND GIVEN
IN ONEMINX , EPS1 AND EPS2 HAVE THE SAME UPPER BOUND AS BEFORE.
IF ABS(X) <= 0.99 AND THE VALUE OF 1 - ABS(X) IS KNOWN EXACTLY
EPS3 <= 3.6"-14.
FOR "-300 <= 1 - ABS(X) < "-2 WE FOUND EPS3 <= 2.2"-11.
REFERENCES:
1. ANTHONY J. STRECOK.
ON THE CALCULATION OF THE INVERSE OF THE ERROR FUNCTION.
MATH. OF COMP., V. 22, 1968, PP144 - 158.
EXAMPLE OF USE:
IN THE FOLLOWING PROGRAM WE COMPUTE THE VALUES OF INVERF(0.6) AND
INVERF(1 - "-150):
"BEGIN"
"PROCEDURE" INVERSE ERROR FUNCTION(X, ONEMINX, INVERF);
"CODE" 35023;
"REAL" INVERF1, INVERF2;
INVERSE ERROR FUNCTION(0.6, 0, INVERF1);
INVERSE ERROR FUNCTION(1, "-150, INVERF2);
OUTPUT(61,"(""(" X = ")", +D.D, "(" 1 - X = ")", +D.3D"+2ZD,
"(" INVERF = ")", +.5DB5DB5D"+D, /")",
0.6, 0.4, INVERF1);
OUTPUT(61,"(""(" X = ")", +D.D, "(" 1 - X = ")", +D.3D"+2ZD,
"(" INVERF = ")", +.5DB5DB5D"+D, /")",
1 - "-150, "-150, INVERF2)
"END"
THIS PROGRAM DELIVERS:
X = +0.6 1 - X = +4.000" -1 INVERF = +.59511 60814 50000"+0
X = +1.0 1 - X = +1.000"-150 INVERF = +.18490 44855 00090"+2
1SECTION : 6.7 (OCTOBER 1974) PAGE 8
SUBSECTION: FRESNEL.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS :
"PROCEDURE" FRESNEL(X, C, S); "VALUE" X; "REAL" X, C, S;
THE MEANING OF THE FORMAL PARAMETERS IS :
X: <ARITHMETIC EXPRESSION>;
ENTRY: THE (REAL) ARGUMENT OF C(X) AND S(X);
C: <VARIABLE>;
EXIT: THE VALUE OF C(X);
S: <VARIABLE>;
EXIT: THE VALUE OF S(X).
PROCEDURES USED: FG = CP35028.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
SEE METHOD AND PERFORMANCE OF THE PROCEDURE FG (THIS SECTION).
REFERENCES :
SEE REF. [1] AND [3] OF THE PROCEDURE FG (THIS SECTION).
1SECTION : 6.7 (OCTOBER 1974) PAGE 9
SUBSECTION: FG.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS :
"PROCEDURE" FG(X, F, G); "VALUE" X; "REAL" X, F, G;
THE MEANING OF THE FORMAL PARAMETERS IS :
X: <ARITHMETIC EXPRESSION>;
ENTRY: THE (REAL) ARGUMENT OF F(X) AND G(X);
F: <VARIABLE>;
EXIT: THE VALUE OF F(X);
G: <VARIABLE>;
EXIT: THE VALUE OF G(X).
PROCEDURES USED: FRESNEL = CP35027.
RUNNING TIME: ABOUT 0.001 400 SEC.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
IF ABS(X) <= 1.6 THE FRESNEL INTEGRALS ARE COMPUTED WITH RATIONAL
CHEBYSHEV APPROXIMATIONS AS GIVEN IN [3]. ON THIS INTERVAL THE
FUNCTIONS F AND G ARE CALCULATED BY MEANS OF THE EQUATIONS GIVEN IN
THE BRIEF DESCRIPTION.
IF ABS(X) > 1.6 THE FUNCTIONS F AND G ARE COMPUTED WITH RATIONAL
CHEBYSHEV APPROXIMATIONS AS GIVEN IN [3]. IN THIS CASE THE FRESNEL
INTEGRALS ARE COMPUTED BY MEANS OF
C(X) = 0.5 + F(X)SIN(PI / 2 * X * X) - G(X)COS(PI / 2 * X * X)
AND
S(X) = 0.5 - F(X)COS(PI / 2 * X * X) - G(X)SIN(PI / 2 * X * X).
IF X < 0 WE USE THE RELATIONS
C(-X) = -C(X), S(-X) = -S(X), F(-X) = -F(X) AND G(-X) = -G(X).
THE FUNCTION VALUES ARE COMPUTED WITH A RELATIVE PRECISION OF
ABOUT "-14.
1SECTION : 6.7 (OCTOBER 1974) PAGE 10
REFERENCES:
[1].M.ABRAMOWITZ AND I.A.STEGUN (ED.).
HANDBOOK OF MATHEMATICAL FUNCTIONS.
DOVER PUBLICATIONS, INC., NEW YORK, 1965.
[2].W.J.CODY.
RATIONAL CHEBYSHEV APPROXIMATIONS FOR THE ERROR FUNCTION.
MATH. COMP. V. 23, 1969, PP631-637.
[3].W.J.CODY.
CHEBYSHEV APPROXIMATIONS FOR THE FRESNEL INTEGRALS.
MATH. COMP. V. 22, 1968, PP450-453.
EXAMPLE OF USE:
IN THE FOLLOWING PROGRAM WE COMPUTE THE VALUES OF C(X), S(X), F(X)
AND G(X) FOR X = 1;
"BEGIN"
"PROCEDURE" FRESNEL(X, C, S); "CODE" 35027;
"PROCEDURE" FG(X, F, G); "CODE" 35028;
"REAL" C, S, F, G;
FRESNEL(1, C, S);
FG(1, F, G);
OUTPUT(61, "(""(" C(1) = ")", +.5DB5D,
"(" S(1) = ")", +.5DB5D, /")", C, S);
OUTPUT(61, "(""(" F(1) = ")", +.5DB5D,
"(" G(1) = ")", +.5DB5D")", F, G)
"END"
THIS PROGRAM DELIVERS:
C(1) = +.77989 34004 S(1) = +.43825 91474
F(1) = +.27989 34004 G(1) = +.06174 08526
1SECTION : 6.7 (OCTOBER 1974) PAGE 11
SOURCE TEXT(S) :
0"CODE" 35021;
"PROCEDURE" ERRORFUNCTION(X, ERF, ERFC);
"VALUE" X; "REAL" X, ERF, ERFC;
"IF" X > 26 "THEN" "BEGIN" ERF:= 1; ERFC:= 0 "END" "ELSE"
"IF" X < -5.5 "THEN" "BEGIN" ERF:= -1; ERFC:= 2 "END" "ELSE"
"BEGIN" "REAL" ABSX, C, P, Q;
"REAL" "PROCEDURE" NONEXPERFC(X); "CODE" 35022;
ABSX:= ABS(X);
"IF" ABSX <= 0.5 "THEN"
"BEGIN" C:= X * X; P:= ((-0.35609 84370 18154"-1 * C +
0.69963 83488 61914"+1) * C + 0.21979 26161 82942"+2) * C +
0.24266 79552 30532"+3;
Q:= ((C +
0.15082 79763 04078"+2) * C + 0.91164 90540 45149"+2) * C +
0.21505 88758 69861"+3;
ERF:= X * P / Q; ERFC:= 1 - ERF
"END" "ELSE"
"BEGIN" ERFC:= EXP(-X * X) * NONEXPERFC(ABSX);
ERF:= 1 - ERFC;
"IF" X < 0 "THEN"
"BEGIN" ERF:= -ERF; ERFC:= 2 - ERFC "END"
"END"
"END" ERRORFUNCTION;
"EOP"
"CODE" 35023;
"PROCEDURE" INVERSE ERROR FUNCTION(X, ONEMINX, INVERF);
"VALUE" X, ONEMINX; "REAL" X, ONEMINX, INVERF;
"BEGIN" "REAL" ABSX, P, BETAX;
"REAL" "ARRAY" A[0 : 23];
"REAL" "PROCEDURE" CHEPOLSER(N, X, A); "CODE" 31046;
ABSX:= ABS(X);
"IF" ABSX > 0.8 "AND" ONEMINX > 0.2 "THEN" ONEMINX:= 0;
"IF" ABSX <= 0.8 "THEN"
"BEGIN"
A[ 0]:= 0.99288 53766 18941; A[ 1]:= 0.12046 75161 43104;
A[ 2]:= 0.01607 81993 42100; A[ 3]:= 0.00268 67044 37162;
A[ 4]:= 0.00049 96347 30236; A[ 5]:= 0.00009 88982 18599;
A[ 6]:= 0.00002 03918 12764; A[ 7]:= 0.00000 43272 71618;
A[ 8]:= 0.00000 09380 81413; A[ 9]:= 0.00000 02067 34720;
A[10]:= 0.00000 00461 59699; A[11]:= 0.00000 00104 16680;
A[12]:= 0.00000 00023 71501; A[13]:= 0.00000 00005 43928;
A[14]:= 0.00000 00001 25549; A[15]:= 0.00000 00000 29138;
A[16]:= 0.00000 00000 06795; A[17]:= 0.00000 00000 01591;
A[18]:= 0.00000 00000 00374; A[19]:= 0.00000 00000 00088;
A[20]:= 0.00000 00000 00021; A[21]:= 0.00000 00000 00005;
INVERF:= CHEPOLSER(21, X * X / 0.32 - 1, A) * X
"END" "ELSE"
"IF" ONEMINX >= 25"-4 "THEN"
"BEGIN" "COMMENT"
1SECTION : 6.7 (MARCH 1977) PAGE 12
;
A[ 0]:= 0.91215 88034 17554; A[ 1]:= -0.01626 62818 67664;
A[ 2]:= 0.00043 35564 72949; A[ 3]:= 0.00021 44385 70074;
A[ 4]:= 0.00000 26257 51076; A[ 5]:= -0.00000 30210 91050;
A[ 6]:= -0.00000 00124 06062; A[ 7]:= 0.00000 00624 06609;
A[ 8]:= -0.00000 00005 40125; A[ 9]:= -0.00000 00014 23208;
A[10]:= 0.00000 00000 34384; A[11]:= 0.00000 00000 33584;
A[12]:= -0.00000 00000 01458; A[13]:= -0.00000 00000 00810;
A[14]:= 0.00000 00000 00053; A[15]:= 0.00000 00000 00020;
BETAX:= SQRT(- LN((1 + ABSX) * ONEMINX));
P:= -1.54881 30423 7326 * BETAX + 2.56549 01231 4782;
P:= CHEPOLSER(15, P, A);
INVERF:= "IF" X < 0 "THEN" - BETAX * P "ELSE" BETAX * P
"END" "ELSE"
"IF" ONEMINX >= 5"-16 "THEN"
"BEGIN"
A[ 0]:= 0.95667 97090 20493; A[ 1]:= -0.02310 70043 09065;
A[ 2]:= -0.00437 42360 97508; A[ 3]:= -0.00057 65034 22651;
A[ 4]:= -0.00001 09610 22307; A[ 5]:= 0.00002 51085 47025;
A[ 6]:= 0.00001 05623 36068; A[ 7]:= 0.00000 27544 12330;
A[ 8]:= 0.00000 04324 84498; A[ 9]:= -0.00000 00205 30337;
A[10]:= -0.00000 00438 91537; A[11]:= -0.00000 00176 84010;
A[12]:= -0.00000 00039 91289; A[13]:= -0.00000 00001 86932;
A[14]:= 0.00000 00002 72923; A[15]:= 0.00000 00001 32817;
A[16]:= 0.00000 00000 31834; A[17]:= 0.00000 00000 01670;
A[18]:= -0.00000 00000 02036; A[19]:= -0.00000 00000 00965;
A[20]:= -0.00000 00000 00220; A[21]:= -0.00000 00000 00010;
A[22]:= 0.00000 00000 00014; A[23]:= 0.00000 00000 00006;
BETAX:= SQRT(- LN((1 + ABSX) * ONEMINX));
P:= -0.55945 76313 29832 * BETAX + 2.28791 57162 6336;
P:= CHEPOLSER(23, P, A);
INVERF:= "IF" X < 0 "THEN" - BETAX * P "ELSE" BETAX * P
"END" "ELSE" "IF" ONEMINX >= 3.13 "-294 "THEN"
"BEGIN"
A[ 0]:= 0.98857 50640 66189; A[ 1]:= 0.01085 77051 84599;
A[ 2]:= -0.00175 11651 02763; A[ 3]:= 0.00002 11969 93207;
A[ 4]:= 0.00001 56648 71404; A[ 5]:= -0.00000 05190 41687;
A[ 6]:= -0.00000 00371 35790; A[ 7]:= 0.00000 00012 17431;
A[ 8]:= -0.00000 00001 76812; A[ 9]:= -0.00000 00000 11937;
A[10]:= 0.00000 00000 00380; A[11]:= -0.00000 00000 00066;
A[12]:= -0.00000 00000 00009;
BETAX:= SQRT(- LN((1 + ABSX) * ONEMINX));
P:= -9.19999 23588 3015 / SQRT(BETAX) + 2.79499 08201 2460;
P:= CHEPOLSER(12, P, A);
INVERF:= "IF" X < 0 "THEN" - BETAX * P "ELSE" BETAX * P
"END" "ELSE" INVERF:= SIGN(X) * 26
"END" INVERSE ERROR FUNCTION;
"EOP"
1SECTION : 6.7 (OCTOBER 1974) PAGE 13
0"CODE" 35022;
"REAL" "PROCEDURE" NONEXPERFC(X); "VALUE" X; "REAL" X;
"BEGIN" "REAL" ABSX, ERF, ERFC, C, P, Q;
"PROCEDURE" ERRORFUNCTION(X, ERF, ERFC); "CODE" 35021;
ABSX:= ABS(X);
"IF" ABSX <= 0.5 "THEN"
"BEGIN" ERRORFUNCTION(X, ERF, ERFC);
NONEXPERFC:= EXP(X * X) * ERFC
"END" "ELSE"
"IF" ABSX < 4 "THEN"
"BEGIN" C:= ABSX; P:= ((((((-0.13686 48573 82717"-6 * C +
0.56419 55174 78974"+0) * C + 0.72117 58250 88309"+1) * C +
0.43162 22722 20567"+2) * C + 0.15298 92850 46940"+3) * C +
0.33932 08167 34344"+3) * C + 0.45191 89537 11873"+3) * C +
0.30045 92610 20162"+3;
Q:= ((((((C +
0.12782 72731 96294"+2) * C + 0.77000 15293 52295"+2) * C +
0.27758 54447 43988"+3) * C + 0.63898 02644 65631"+3) * C +
0.93135 40948 50610"+3) * C + 0.79095 09253 27898"+3) * C +
0.30045 92609 56983"+3;
NONEXPERFC:= "IF" X > 0 "THEN" P / Q "ELSE"
EXP(X * X) * 2 - P / Q
"END" "ELSE"
"BEGIN" C:= 1 / X / X; P:= (((0.22319 24597 34185"-1 * C +
0.27866 13086 09648"-0) * C + 0.22695 65935 39687"-0) * C +
0.49473 09106 23251"-1) * C + 0.29961 07077 03542"-2;
Q:= (((C +
0.19873 32018 17135"+1) * C + 0.10516 75107 06793"+1) * C +
0.19130 89261 07830"+0) * C + 0.10620 92305 28468"-1;
C:= (C * (-P) / Q + 0.56418 95835 47756) / ABSX;
NONEXPERFC:= "IF" X > 0 "THEN" C "ELSE" EXP(X * X) * 2 - C
"END"
"END" NONEXPERFC;
"EOP"
0"CODE" 35027;
"PROCEDURE" FRESNEL(X, C, S); "VALUE" X; "REAL" X, C, S;
"BEGIN" "REAL" ABSX, X3, X4, A, P, Q, F, G, C1, S1;
"PROCEDURE" FG(X, F, G); "CODE" 35028;
ABSX:= ABS(X);
"IF" ABSX <= 1.2 "THEN"
"BEGIN" A:= X * X; X3:= A * X; X4:= A * A;
P:= (((5.47711 38568 2687"-6 * X4 - 5.28079 65137 2623"-4)
* X4 + 1.76193 95254 3491"-2) * X4 - 1.99460 89882 6184"-1)
* X4 + 1;
Q:= (((1.18938 90142 2876"-7 * X4 + 1.55237 88527 6994"-5)
* X4 + 1.09957 21502 5642"-3) * X4 + 4.72792 11201 0453"-2)
* X4 + 1;
C:= X * P / Q;
P:= (((6.71748 46662 5141"-7 * X4 - 8.45557 28435 2777"-5)
* X4 + 3.87782 12346 3683"-3) * X4 - 7.07489 91514 4523"-2)
* X4 + 5.23598 77559 8299"-1; "COMMENT"
1SECTION : 6.7 (OCTOBER 1974) PAGE 14
;
Q:= (((5.95281 22767 8410"-8 * X4 + 9.62690 87593 9034"-6)
* X4 + 8.17091 94215 2134"-4) * X4 + 4.11223 15114 2384"-2)
* X4 + 1;
S:= X3 * P / Q
"END" "ELSE"
"IF" ABSX <= 1.6 "THEN"
"BEGIN" A:= X * X; X3:= A * X; X4:= A * A;
P:=((((-5.68293 31012 1871"-8 * X4 + 1.02365 43505 6106"-5)
* X4 - 6.71376 03469 4922"-4) * X4 + 1.91870 27943 1747"-2)
* X4 - 2.07073 36033 5324"-1) * X4 + 1.00000 00000 0111"+0;
Q:=((((4.41701 37406 5010"-10 * X4 + 8.77945 37789 2369"-8)
* X4 + 1.01344 63086 6749"-5) * X4 + 7.88905 24505 2360"-4)
* X4 + 3.96667 49695 2323"-2) * X4 + 1;
C:= X * P / Q;
P:=((((-5.76765 81559 3089"-9 * X4 + 1.28531 04374 2725"-6)
* X4 - 1.09540 02391 1435"-4) * X4 + 4.30730 52650 4367"-3)
* X4 - 7.37766 91401 0191"-2) * X4 + 5.23598 77559 8344"-1;
Q:=((((2.05539 12445 8580"-10 * X4 + 5.03090 58124 6612"-8)
* X4 + 6.87086 26571 8620"-6) * X4 + 6.18224 62019 5473"-4)
* X4 + 3.53398 34276 7472"-2) * X4 + 1;
S:= X3 * P / Q
"END" "ELSE"
"IF" ABSX < "15 "THEN"
"BEGIN" FG(X, F, G);
A:= X * X;
A:= (A - ENTIER(A / 4) * 4) * 1.57079 63267 9490;
C1:= COS(A); S1:= SIN(A);
A:= "IF" X < 0 "THEN" -0.5 "ELSE" 0.5;
C:= F * S1 - G * C1 + A;
S:= -F * C1 - G * S1 + A
"END" "ELSE" C:= S:= SIGN(X) * 0.5
"END" FRESNEL;
"EOP"
0"CODE" 35028;
"PROCEDURE" FG(X, F, G); "VALUE" X; "REAL" X, F, G;
"BEGIN" "REAL" ABSX, C, S, C1, S1, A, XINV, X3INV, C4, P, Q;
"PROCEDURE" FRESNEL(X, C, S); "CODE" 35027;
ABSX:= ABS(X);
"IF" ABSX <= 1.6 "THEN"
"BEGIN" FRESNEL(X, C, S);
A:= X * X * 1.57079 63267 9490; C1:= COS(A); S1:= SIN(A);
A:= "IF" X < 0 "THEN" -0.5 "ELSE" 0.5;
P:= A - C; Q:= A - S;
F:= Q * C1 - P * S1;
G:= P * C1 + Q * S1
"END" "ELSE"
"IF" ABSX <= 1.9 "THEN"
"BEGIN" XINV:= 1 / X; A:= XINV * XINV;
X3INV:= A * XINV; C4:= A * A; "COMMENT"
1SECTION : 6.7 (OCTOBER 1974) PAGE 15
;
P:= (((1.35304 23554 0388"+1 * C4 + 6.98534 26160 1021"+1)
* C4 + 4.80340 65557 7925"+1) * C4 + 8.03588 12280 3942"+0)
* C4 + 3.18309 26850 4906"-1;
Q:= (((6.55630 64008 3916"+1 * C4 + 2.49561 99380 5172"+2)
* C4 + 1.57611 00558 0123"+2) * C4 + 2.55491 61843 5795"+1)
* C4 + 1;
F:= XINV * P / Q;
P:=((((2.05421 43249 8501"+1 * C4 + 1.96232 03797 1663"+2)
* C4 + 1.99182 81867 8903"+2) * C4 + 5.31122 81348 0989"+1)
* C4 + 4.44533 82755 0512"+0) * C4 + 1.01320 61881 0275"-1;
Q:=((((1.01379 48339 6003"+3 * C4 + 3.48112 14785 6545"+3)
* C4 + 2.54473 13318 1822"+3) * C4 + 5.83590 57571 6429"+2)
* C4 + 4.53925 01967 3689"+1) * C4 + 1;
G:= X3INV * P / Q
"END" "ELSE"
"IF" ABSX <= 2.4 "THEN"
"BEGIN" XINV:= 1 / X; A:= XINV * XINV;
X3INV:= A * XINV; C4:= A * A;
P:=((((7.17703 24936 5140"+2 * C4 + 3.09145 16157 4430"+3)
* C4 + 1.93007 64078 6716"+3) * C4 + 3.39837 13492 6984"+2)
* C4 + 1.95883 94102 1969"+1) * C4 + 3.18309 88182 2017"-1;
Q:=((((3.36121 69918 0551"+3 * C4 + 1.09334 24898 8809"+4)
* C4 + 6.33747 15585 1144"+3) * C4 + 1.08535 06750 0650"+3)
* C4 + 6.18427 13817 2887"+1) * C4 + 1;
F:= XINV * P / Q;
P:=((((3.13330 16306 8756"+2 * C4 + 1.59268 00608 5354"+3)
* C4 + 9.08311 74952 9594"+2) * C4 + 1.40959 61791 1316"+2)
* C4 + 7.11205 00178 9783"+0) * C4 + 1.01321 16176 1805"-1;
Q:=((((1.15149 83237 6261"+4 * C4 + 2.41315 56721 3370"+4)
* C4 + 1.06729 67803 0581"+4) * C4 + 1.49051 92279 7329"+3)
* C4 + 7.17128 59693 9302"+1) * C4 + 1;
G:= X3INV * P / Q
"END" "ELSE"
"BEGIN" XINV:= 1 / X; A:= XINV * XINV;
X3INV:= A * XINV; C4:= A * A;
P:=((((2.61294 75322 5142"+4 * C4 + 6.13547 11361 4700"+4)
* C4 + 1.34922 02817 1857"+4) * C4 + 8.16343 40178 4375"+2)
* C4 + 1.64797 71284 1246"+1) * C4 + 9.67546 03296 7090"-2;
Q:=((((1.37012 36481 7226"+6 * C4 + 1.00105 47890 0791"+6)
* C4 + 1.65946 46262 1853"+5) * C4 + 9.01827 59623 1524"+3)
* C4 + 1.73871 69067 3649"+2) * C4 + 1;
F:= (C4 * (-P) / Q + 0.31830 98861 83791) * XINV;
P:=(((((1.72590 22465 4837"+6 * C4 + 6.66907 06166 8636"+6)
* C4 + 1.77758 95083 8030"+6) * C4 + 1.35678 86781 3756"+5)
* C4 + 3.87754 14174 6378"+3) * C4 + 4.31710 15782 3358"+1)
* C4 + 1.53989 73381 9769"-1;
Q:=(((((1.40622 44112 3580"+8 * C4 + 9.38695 86253 1635"+7)
* C4 + 1.62095 60050 0232"+7) * C4 + 1.02878 69305 6688"+6)
* C4 + 2.69183 18039 6243"+4) * C4 + 2.86733 19497 5899"+2)
* C4 + 1;
G:= (C4 * (-P) / Q + 0.10132 11836 42338) * X3INV
"END"
"END" FG;
"EOP"
1SECTION : 2.2.2.2 (DECEMBER 1978) PAGE 1
AUTHOR: C.G. VAN DER LAAN.
CONTRIBUTORS: C.G. VAN DER LAAN, M. VOORINTHOLT.
INSTITUTE: REKENCENTRUM RIJKSUNIVERSITEIT GRONINGEN.
RECEIVED: 780601.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS FOUR PROCEDURES ABOUT CHEBYSHEV POLYNOMIALS
OF THE FIRST KIND:
CHEPOLSUM: EVALUATES A (FINITE) SUM OF CHEBYSHEV POLYNOMIALS,
ODDCHEPOLSUM: EVALUATES A (FINITE) SUM OF CHEBYSHEV POLYNOMIALS
OF ODD DEGREE,
CHEPOL: EVALUATES A CHEBYSHEV POLYNOMIAL,
ALLCHEPOL: EVALUATES ALL CHEBYSHEV POLYNOMIALS WITH DEGREE LESS
THAN A GIVEN POSITIVE INTEGER.
KEYWORDS:
(FINITE) SUM OF (SHIFTED) CHEBYSHEV POLYNOMIALS OF THE FIRST KIND,
GOERTZEL,WATT,CLENSHAW,GENERALIZED HORNER ALGORITHM,
LINEAR THREE TERM (INHOMOGENEOUS) RECURRENCE RELATION.
REFERENCES:
CLENSHAW, C.W. (1962):
CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS.
MATH. TAB. NAT. PHYS. LAB. 5, LONDON. HMO.
FOX, L. & I.B. PARKER (1972):
CHEBYSHEV POLYNOMIALS IN NUMERICAL ANALYSIS.
OXFORD UNIVERSITY PRESS.
RIVLIN, T.J. (1974):
THE CHEBYSHEV POLYNOMIALS.
WILEY.
STOER,J.(1972):
EINFUEHRUNG IN DIE NUMERISCHE MATHEMATIK 1.
HEIDELBERGER TASCHENBUECHER 105,SPRINGER-VERLAG.
1SECTION : 2.2.2.2 (DECEMBER 1978) PAGE 2
SUBSECTION: CHEPOLSUM.
CALLING SEQUENCE:
THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
"REAL""PROCEDURE"CHEPOLSUM(N,X,A);
"VALUE"N,X;"INTEGER"N;"REAL"X;"ARRAY"A;
"CODE"31046;
CHEPOLSUM:=THE VALUE OF THE CHEBYSHEV SUM
A[0] + A[1]*T[1](X) + .... + A[N]*T[N](X),
WHERE T[1](X),....,T[N](X) ARE CHEBYSHEV POLYNOMIALS
OF THE FIRST KIND, OF DEGREE 1,....,N, RESPECTIVELY.
THE MEANING OF THE FORMAL PARAMETERS IS :
N : <ARITHMETIC EXPRESSION>;
ENTRY: THE DEGREE OF THE POLYNOMIAL REPRESENTED BY THE
CHEBYSHEV SUM (N>=0);
X : <ARITHMETIC EXPRESSION>;
ENTRY: THE ARGUMENT OF THE CHEBYSHEV POLYNOMIALS , ABS(X)<=1;
A : <ARRAY IDENTIFIER>;
"ARRAY" A[0:N];
ENTRY: THE COEFFICIENTS OF THE CHEBYSHEV SUM MUST BE GIVEN IN
ARRAY A, WHERE A[K] IS THE COEFFICIENT OF THE CHEBYSHEV
POLYNOMIAL OF DEGREE K, 0<=K<=N.
PROCEDURES USED: NONE.
RUNNING TIME: PROPORTIONAL TO N.
METHOD AND PERFORMANCE:
N N / 2*X 1 \ K / A[K] \
SUM A[K]*T[K](X) = (1,X) * SUM [ ] * [ ]
K=0 K=0 \ -1 0 / \ 0 /
WE USE THE CLENSHAW OR GENERALIZED HORNER ALGORITHM:
N
SUM A[K]*T[K](X) = (1,X) *
K=0
/ / A[0] \ / 2*X 1 \ / / A[1] \ / 2*X 1 \ / A[N] \ \ \
[ [ ] + [ ] * [ [ ] +..+ [ ] * [ ] ]..].
\ \ 0 / \ -1 0 / \ \ 0 / \ -1 0 / \ 0 / / /
1SECTION 2.2.2.2 (DECEMBER 1978) PAGE 3
THIS PROCEDURE MAY BE USED:
- TO EVALUATE THE SUM OF CHEBYSHEV POLYNOMIALS OF EVEN DEGREE
N
SUM A[K]*T[2*K](X) ,
K=0
BY THE CALL OF
CHEPOLSUM(N,"IF" ABS(X)<ARREB "THEN" -1 "ELSE" 2*X*X-1,A)
BECAUSE OF
T[2*K](X) = T[K](T[2](X));
(ARREB DENOTES THE MACHINE PRECISION)
- TO EVALUATE THE SUM OF SHIFTED CHEBYSHEV POLYNOMIALS FOR 0<=X<=1
N
SUM A[K]*T'[K](X) ,
K=0
BY THE CALL OF
CHEPOLSUM(N,2*X-1,A)
BECAUSE OF
T'[K](X) = T[K](2*X-1).
EXAMPLE OF USE :
THE POLYNOMIAL : 1 + 1/2*T[1](X) + 1/4*T[2](X) IS EVALUATED FOR
X = -1,0,1, WHERE T[1](X) AND T[2](X) ARE THE CHEBYSHEV POLYNOMIALS
OF FIRST AND SECOND DEGREE, RESPECTIVELY.
"BEGIN""ARRAY"A[0:2];
"REAL""PROCEDURE"CHEPOLSUM(N,X,A);
"VALUE"N,X;"INTEGER"N;"REAL"X;"ARRAY"A;
"CODE"31046;
A[2]:=.25;A[1]:=.5;A[0]:=1;
OUTPUT(61,"("3(BZ.DD)")",CHEPOLSUM(2,-1,A),CHEPOLSUM(2,0,A),
CHEPOLSUM(2,1,A))
"END"
.75 .75 1.75
SUBSECTION: ODDCHEPOLSUM.
CALLING SEQUENCE:
THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
"REAL""PROCEDURE" ODDCHEPOLSUM(N,X,A);
"VALUE"N,X;"INTEGER"N;"REAL"X;"ARRAY"A;
"CODE"31059;
1SECTION : 2.2.2.2 (DECEMBER 1978) PAGE 4
ODDCHEPOLSUM:= THE VALUE OF THE CHEBYSHEV SUM
A[1]*T[1](X) + .... + A[N]*T[2*N+1](X),
WHERE T[1](X),....,T[2*N+1](X) ARE CHEBYSHEV POLYNOMIALS
OF THE FIRST KIND, OF (ODD) DEGREE 1,....,2*N+1.
THE MEANING OF THE FORMAL PARAMETERS IS:
N : <ARITHMETIC EXPRESSION>;
ENTRY: THE DEGREE OF THE POLYNOMIAL REPRESENTED BY
THE CHEBYSHEV SUM IS 2*N+1 (N>=0);
X : <ARITHMETIC EXPRESSION>;
ENTRY: THE ARGUMENT OF THE CHEBYSHEV POLYNOMIALS, ABS(X)<=1;
A : <ARRAY IDENTIFIER>;
"ARRAY" A[0:N];
ENTRY: THE COEFFICIENTS OF THE CHEBYSHEV SUM MUST BE GIVEN IN
ARRAY A, WHERE A[K] IS THE COEFFICIENT OF THE CHEBYSHEV
POLYNOMIAL OF DEGREE 2*K+1, 0<=K<=N.
PROCEDURES USED: NONE.
RUNNING TIME: PROPORTIONAL TO N.
METHOD AND PERFORMANCE:
FROM THE REPRESENTATION, FOR ABS(X)<=1,
N N / 2*T[2](X) -1 \ K /A[K] \
SUM A[K]*T[2*K+1](X) = X * (1,-1) * SUM [ ] * [ ]
K=0 K=0 \ 1 0 / \ 0 /
WE USE THE CLENSHAW OR GENERALIZED HORNER ALGORITHM:
N / / A[0] \
SUM A[K]*T[2*K+1](X) = X * (1,-1) * [ [ ] +
K=0 \ \ 0 /
/ 2*T[2](X) -1 \ / / A[1] \ / 2*T[2](X) -1 \ / A[N] \ \ \
[ ] * [ [ ]+...+[ ] * [ ] ]..].
\ 1 0 / \ \ 0 / \ 1 0 / \ 0 / / /
THIS PROCEDURE MAY BE USED TO EVALUATE THE SUM OF SHIFTED CHEBYSHEV
POLYNOMIALS OF ODD DEGREE FOR 0<=X<=1,
N
SUM A[K]*T'[2*K+1](X) ,
K=0
BY THE CALL OF
ODDCHEPOLSUM(N,2*X-1,A)
BECAUSE OF
T'[K](X) = T[K](2*X-1).
1SECTION : 2.2.2.2 (DECEMBER 1978) PAGE 5
EXAMPLE OF USE:
THE POLYNOMIAL 1/2*T[1](X) + 1/5*T[3](X) IS EVALUATED FOR X=-1,0,1,
WHERE T[1](X) AND T[3](X) ARE CHEBYSHEV POLYNOMIALS OF THE FIRST
AND THIRD DEGREE, RESPECTIVELY.
"BEGIN"
"ARRAY"A[0:1];
"REAL""PROCEDURE"ODDCHEPOLSUM(N,X,A);
"VALUE"N,X;"INTEGER"N;"REAL"X;"ARRAY"A;
"CODE"31049;
A[1]:=.2;A[0]:=.5;
OUTPUT(61,"("/,3(B,-Z.DD)")",ODDCHEPOLSUM(1,-1,A),
ODDCHEPOLSUM(1,0,A),
ODDCHEPOLSUM(1,1,A));
"END"
-.70 .00 .70
SUBSECTION: CHEPOL.
CALLING SEQUENCE:
THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
"REAL""PROCEDURE"CHEPOL(N,X);
"VALUE"N,X;"INTEGER"N;"REAL"X;
"CODE"31042;
CHEPOL:=THE VALUE OF THE CHEBYSHEV POLYNOMIAL OF THE FIRST KIND OF
DEGREE N FOR THE ARGUMENT X.
THE MEANING OF THE FORMAL PARAMETERS IS:
N : <ARITHMETIC EXPRESSION>;
ENTRY: THE DEGREE OF THE POLYNOMIAL (N>=0);
X : <ARITHMETIC EXPRESSION>;
ENTRY: THE ARGUMENT OF THE CHEBYSHEV POLYNOMIAL, ABS(X)<=1.
PROCEDURES USED: NONE.
RUNNING TIME: PROPORTIONAL TO N.
METHOD AND PERFORMANCE: SEE ALLCHEPOL (NEXT SUBSECTION).
EXAMPLE OF USE: SEE NEXT SUBSECTION.
1SECTION : 2.2.2.2 (DECEMBER 1978) PAGE 6
SUBSECTION: ALLCHEPOL.
CALLING SEQUENCE:
THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:
"PROCEDURE"ALLCHEPOL(N,X,T);
"VALUE"N,X;"INTEGER"N;"REAL"X;"REAL""ARRAY"T;
"CODE"31043;
THE MEANING OF THE FORMAL PARAMETERS IS:
N : <ARITHMETIC EXPRESSION>;
ENTRY: THE DEGREE OF THE LAST POLYNOMIAL (N>=0);
X : <ARITHMETIC EXPRESSION>;
ENTRY: THE ARGUMENT OF THE CHEBYSHEV POLYNOMIALS, ABS(X)<=1;
T : <ARRAY IDENTIFIER>;
"ARRAY" T[0:N];
EXIT: THE VALUES OF THE CHEBYSHEV POLYNOMIALS OF THE FIRST
KIND OF DEGREES 0,1,...,N , FOR THE ARGUMENT X, ARE
DELIVERED IN T[0],T[1],...,T[N], RESPECTIVELY.
PROCEDURES USED: NONE.
RUNNING TIME: PROPORTIONAL TO N.
METHOD AND PERFORMANCE:
FOR A DESCRIPTION OF THE ALGORITHM SEE STOER,1972,P.21.
THE MAXIMUM (ABSOLUTE) VALUE OF THE CHEBYSHEV POLYNOMIAL
EQUALS 1 AS A NORMALIZATION.
AN UPPER BOUND FOR THE (ABSOLUTE) ERROR IS A QUADRATIC FUNCTION
OF THE DEGREE OF THE CHEBYSHEV POLYNOMIAL.THIS UPPER BOUND IS A
ROUGH OVER-ESTIMATE FOR THE SPECIAL CASE ABS(X)<.5 (STOER,1972,
P. 21-24).
EXAMPLE OF USE :
BY THE PROCEDURE (ALL)CHEPOL THE CHEBYSHEV POLYNOMIALS OF THE FIRST
KIND OF DEGREES 0,1,2 ARE EVALUATED AT -1,0,1.
1SECTION : 2.2.2.2 (DECEMBER 1978) PAGE 7
"BEGIN"
"ARRAY"T[0:2];
"REAL""PROCEDURE"CHEPOL(N,X);
"VALUE"N,X;"INTEGER"N;"REAL"X;
"CODE"31042;
"PROCEDURE"ALLCHEPOL(N,X,T);
"VALUE"N,X;"INTEGER"N;"REAL"X;"REAL""ARRAY"T;
"CODE"31043;
ALLCHEPOL(2,-1,T);OUTPUT(61,"("/,3(-DB)")",T[0],T[1],T[2]);
ALLCHEPOL(2, 0,T);OUTPUT(61,"("/,3(-DB)")",T[0],T[1],T[2]);
ALLCHEPOL(2, 1,T);OUTPUT(61,"("/,3(-DB)")",T[0],T[1],T[2]);
OUTPUT(61,"("/,3(/,-D)")",CHEPOL(2,-1),CHEPOL(2,0),CHEPOL(2,1))
"END"
1 -1 1
1 0 -1
1 1 1
1
-1
1
SOURCE TEXT(S):
"CODE"31046;
"REAL" "PROCEDURE" CHEPOLSUM(N,X,A);
"VALUE" N,X;"INTEGER" N;"REAL" X;"ARRAY" A;
"IF" N=0 "THEN" CHEPOLSUM:=A[0] "ELSE"
"IF" N=1 "THEN" CHEPOLSUM:=A[0]+A[1]*X "ELSE"
"BEGIN" "INTEGER" K;"REAL" H,R,S,TX;
TX:=X+X;R:=A[N];
H:=A[N-1]+R*TX;
"FOR" K:=N-2 "STEP" -1 "UNTIL" 1 "DO"
"BEGIN" S:=R;R:=H;
H:=A[K]+R*TX-S
"END"K;
CHEPOLSUM:=A[0]-R+H*X
"END" CHEPOLSUM;
"EOP"
1SECTION : 2.2.2.2 (DECEMBER 1978) PAGE 8
"CODE"31059;
"REAL""PROCEDURE"ODDCHEPOLSUM(N,X,A);
"VALUE"N,X;"INTEGER"N;"REAL"X;"ARRAY"A;
"COMMENT" ODDCHEPOLSUM:=A[0]T[1](X)+A[1]T[3](X)+....+A[N]T[2N+1](X);
"IF" N=0 "THEN" ODDCHEPOLSUM:=X*A[0] "ELSE"
"IF" N=1 "THEN" ODDCHEPOLSUM:=X*(A[0]+A[1]*(4*X*X-3)) "ELSE"
"BEGIN"
"INTEGER" K;
"REAL" H,R,S,Y;
Y:=4*X*X-2;
R:=A[N];
H:=A[N-1]+R*Y;
"FOR" K:=N-2 "STEP" -1 "UNTIL" 0 "DO"
"BEGIN"
S:=R;
R:=H;
H:=A[K]+R*Y-S;
"END" K;
ODDCHEPOLSUM:=X*(H-R);
"END" ODDCHEPOLSUM;
"EOP"
"CODE"31042;
"REAL""PROCEDURE"CHEPOL(N,X);
"VALUE"N,X;"INTEGER"N;"REAL"X;
"IF" N = 0 "THEN" CHEPOL :=1 "ELSE"
"IF" N = 1 "THEN" CHEPOL :=X "ELSE"
"BEGIN""INTEGER"I;"REAL"T1,T2,H,X2;
T2:=X;T1:=1;X2:=X+X;
"FOR"I:=2"STEP"1"UNTIL"N"DO"
"BEGIN"H:=X2*T2-T1;T1:=T2;T2:=H"END";
CHEPOL:=H
"END"CHEPOL;
"EOP"
"CODE"31043;
"PROCEDURE"ALLCHEPOL(N,X,T);
"VALUE"N,X;"INTEGER"N;"REAL"X;"REAL""ARRAY"T;
"IF" N = 0 "THEN" T[0] :=1 "ELSE"
"IF" N = 1 "THEN" "BEGIN" T[0] := 1; T[1] := X "END" "ELSE"
"BEGIN""INTEGER"I;"REAL"T1,T2,H,X2;
T[0]:=T1:=1;T[1]:=T2:=X;X2:=X+X;
"FOR"I:=2"STEP"1"UNTIL"N"DO"
"BEGIN"T[I]:=H:=X2*T2-T1;T1:=T2;T2:=H"END"
"END"ALLCHEPOL;
"EOP"
1SECTION : 5.2.1.1.1.2.E (OCTOBER 1974) PAGE 1
AUTHOR: J.G. VERWER.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 740809.
BRIEF DESCRIPTION:
GMS SOLVES AN INITIAL VALUE PROBLEM, GIVEN AS AN AUTONOMOUS
SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS DY/DX = F(Y), BY
MEANS OF A THIRD ORDER GENERALIZED LINEAR MULTISTEP METHOD;
IN PARTICULAR THIS PROCEDURE IS SUITABLE FOR THE INTEGRATION
OF STIFF EQUATIONS.
KEYWORDS:
DIFFERENTIAL EQUATIONS,
INITIAL VALUE PROBLEM,
AUTONOMOUS SYSTEM,
STIFF EQUATIONS,
GENERALIZED LINEAR MULTISTEP METHOD.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" GMS(X, XE, R, Y, H, HMIN, HMAX, DELTA, DERIVATIVE,
JACOBIAN, AETA, RETA, N, JEV, LU, NSJEV,
LINEAR, OUT);
"VALUE" R;
"REAL" X, XE, H, HMIN, HMAX, DELTA, AETA, RETA;
"INTEGER" R, N, JEV, NSJEV, LU;
"BOOLEAN" LINEAR;
"ARRAY" Y;
"PROCEDURE" DERIVATIVE, JACOBIAN, OUT;
GMS INTEGRATES THE SYSTEM OF DIFFERENTIAL EQUATIONS DY/DX = F(Y)
FROM X = X0 TO X = XE;
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <VARIABLE>;
THE INDEPENDENT VARIABLE X;
ENTRY: THE INITIAL VALUE OF X;
EXIT : THE END VALUE OF X;
XE: <ARITHMETIC EXPRESSION>;
ENTRY: THE END VALUE OF X;
R: <ARITHMETIC EXPRESSION>;
ENTRY: THE NUMBER OF DIFFERENTIAL EQUATIONS;
1SECTION : 5.2.1.1.1.2.E (MARCH 1977) PAGE 2
Y: <ARRAY IDENTIFIER>;
"ARRAY" Y[1:R];
THE DEPENDENT VARIABLE;
ENTRY: THE INITIAL VALUE OF Y;
EXIT : THE SOLUTION Y AT THE POINT X AFTER EACH
INTEGRATION STEP;
H: <ARITHMETIC EXPRESSION>;
ENTRY: THE STEPLENGTH WHEN THE INTEGRATION HAS TO BE
PERFORMED WITHOUT THE STEPSIZE MECHANISM, OTHER-
WISE THE INITIAL STEPLENGTH (SEE THE PARAMETERS
HMIN AND HMAX);
HMIN, HMAX: <ARITHMETIC EXPRESSION>;
ENTRY: MINIMAL AND MAXIMAL STEPLENGTH BY WHICH THE INTE-
GRATION IS ALLOWED TO BE PERFORMED;
BY PUTTING HMIN = HMAX THE STEPSIZE MECHANISM IS
ELIMINATED; IN THIS CASE THE GIVEN VALUES FOR HMIN AND
HMAX ARE IRRELEVANT, WHILE THE INTEGRATION IS PERFORMED
WITH THE STEPLENGTH GIVEN BY H;
DELTA: <ARITHMETIC EXPRESSION>;
ENTRY: THE REAL PART OF THE POINT AT WHICH EXPONENTIAL
FITTING IS DESIRED;
(SEE METHOD AND PERFORMANCE);
ALTERNATIVES:
DELTA = (AN ESTIMATE OF) THE REAL PART OF THE LARGEST
EIGENVALUE IN MODULUS OF THE JACOBIAN MATRIX OF THE
SYSTEM;
DELTA <= -10**15, IN ORDER TO OBTAIN ASYMPTOTIC STABILITY;
DELTA = 0, IN ORDER TO OBTAIN A HIGHER ORDER OF ACCURACY
IN CASE OF LINEAR EQUATIONS;
DERIVATIVE: <PROCEDURE IDENTIFIER>;
"PROCEDURE" DERIVATIVE(Y); "ARRAY" Y;
<REPLACEMENT OF THE I-TH COMPONENT OF THE SOLUTION Y BY
THE I-TH COMPONENT OF THE DERIVATIVE F(Y), I = 1,..., R>;
JACOBIAN: <PROCEDURE IDENTIFIER>;
"PROCEDURE" JACOBIAN(J, Y); "ARRAY" J, Y;
WHEN IN GMS JACOBIAN IS CALLED THE ARRAY Y CONTAINS
THE VALUES OF THE DEPENDENT VARIABLE;
UPON COMPLETION OF A CALL OF JACOBIAN THE ARRAY J SHOULD
CONTAIN THE VALUES OF THE JACOBIAN MATRIX OF F(Y);
AETA, RETA: <ARITHMETIC EXPRESSION>;
ENTRY: MEASURE OF THE ABSOLUTE AND RELATIVE LOCAL
ACCURACY REQUIRED;
THESE VALUES ARE IRRELEVANT WHEN THE INTEGRATION IS PER-
FORMED WITHOUT THE STEPSIZE MECHANISM;
N: <VARIABLE>;
EXIT : THE NUMBER OF INTEGRATION STEPS;
JEV: <VARIABLE>;
EXIT: THE NUMBER OF JACOBIAN EVALUATIONS;
LU: <VARIABLE>;
EXIT: THE NUMBER OF LU-DECOMPOSITIONS;
1SECTION : 5.2.1.1.1.2.E (OCTOBER 1974) PAGE 3
NSJEV: <VARIABLE>;
ENTRY: NUMBER OF INTEGRATION STEPS PER
JACOBIAN EVALUATION;
THE VALUE OF NSJEV IS RELEVANT ONLY WHEN THE INTEGRATION
IS PERFORMED WITHOUT THE STEPSIZE MECHANISM AND THE
SYSTEM TO BE SOLVED IS NON-LINEAR;
LINEAR: <BOOLEAN EXPRESSION>;
ENTRY: TRUE WHEN THE SYSTEM TO BE INTEGRATED IS LINEAR,
OTHERWISE FALSE;
IF LINEAR IS TRUE THE STEPSIZE MECHANISM IS AUTOMATICALLY
ELIMINATED;
OUT: <PROCEDURE IDENTIFIER>;
"PROCEDURE" OUT;
<BY MEANS OF OUT ONE MAY PRINT THE VALUES OF THE RELEVANT
PARAMETERS OCCURRING IN THE PARAMETERLIST; OUT IS CALLED
AFTER EACH INTEGRATION STEP>;
DATA AND RESULTS: SEE REF[2].
PROCEDURES USED:
VECVEC = CP34010,
MATVEC = CP34011,
MATMAT = CP34013,
ELMROW = CP34024,
ELMVEC = CP34020,
DUPVEC = CP31030,
GSSELM = CP34231,
SOLELM = CP34061,
COLCST = CP31131,
MULVEC = CP31020.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: 8 * R + 3 * R * R;
RUNNING TIME:
DEPENDS STRONGLY ON THE DIFFERENTIAL EQUATION TO BE SOLVED.
LANGUAGE: ALGOL 60.
1SECTION : 5.2.1.1.1.2.E (MARCH 1977) PAGE 4
METHOD AND PERFORMANCE:
THE PROCEDURE GMS DESCRIBES AN IMPLEMENTATION OF A THIRD ORDER
THREE-STEP GENERALIZED LINEAR MULTISTEP METHOD WITH QUASI-ZERO
PARASITIC ROOTS AND QUASI-ADAPTIVE STABILITY FUNCTION. IN PARTI-
CULAR THE ALGORITHM IS DEVELOPED FOR THE INTEGRATION OF STIFF
SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS. THE PROCEDURE SUPPLIES
THE ADDITIONAL STARTING VALUES AND PERFORMS A STEPSIZE CONTROL
WHICH IS BASED ON THE NON-LINEARITY OF THE DIFFERENTIAL EQUATION.
BY THIS CONTROL THE JACOBIAN MATRIX IS INCIDENTALLY EVALUATED. IT
IS POSSIBLE TO ELIMINATE THE STEPSIZE CONTROL. THEN, ONE HAS
TO GIVE THE NUMBER OF INTEGRATION STEPS PER JACOBIAN EVALUATION.
FOR LINEAR EQUATIONS THE STEPSIZE CONTROL IS AUTOMATICALLY ELIMIN-
ATED, WHILE THE PROCEDURE PERFORMS ONE EVALUATION OF THE JACOBIAN.
MOREOVER, IN THIS CASE THE THREE-STEP SCHEME IS REDUCED TO A ONE-
STEP SCHEME. THE PROCEDURE USES ONE FUNCTION EVALUATION PER INTE-
GRATION STEP AND IT DOES NOT REJECT INTEGRATION STEPS. EACH CHANGE
IN THE STEPLENGTH OR EACH REEVALUATION OF THE JACOBIAN COSTS ONE
LU-DECOMPOSITION. IT IS POSSIBLE TO FIT EXPONENTIALLY, THIS FIT-
TING IS EQUIVALENT TO FITTING IN THE SENSE OF LINIGER AND
WILLOUGHBY, ONLY WHEN THE JACOBIAN MATRIX IS EVALUATED AT EACH IN-
TEGRATION STEP. WHEN THE SYSTEM TO BE INTEGRATED IS NON-LINEAR
AND THE JACOBIAN MATRIX IS NOT EVALUATED AT EACH INTEGRATION STEP,
IT IS RECOMMENDED TO FIT AT INFINITY (DELTA <= -10**15).
DETAILS ARE GIVEN IN REFERENCE 2.
REFERENCES:
[1] HOUWEN, P. J. VAN DER AND VERWER, J. G.,
GENERALIZED LINEAR MULTISTEP METHODS 1, DEVELOPMENT OF ALGO-
RITHMS WITH ZERO-PARASITIC ROOTS,
REPORT NW 10/74, MATHEMATISCH CENTRUM, AMSTERDAM 1974.
[2] VERWER, J. G.,
GENERALIZED LINEAR MULTISTEP METHODS 2, NUMERICAL
APPLICATIONS, REPORT NW 12/74, MATHEMATISCH CENTRUM,
AMSTERDAM, 1974.
EXAMPLE OF USE:
WE CONSIDER THE DIFFERENTIAL EQUATION
DY1/DX = -1000 * Y1 * (Y1 + Y2 -1.999987),
DY2/DX = -2500 * Y2 * (Y1 + Y2 - 2),
ON THE INTERVAL [0,50], WITH INITIAL VALUE Y1(0) = Y2(0) = 1.
THE REFERENCE SOLUTION AT X = 50 IS GIVEN BY:
Y1(50) = .5976546988,
Y2(50) = 1.4023434075.
1SECTION : 5.2.1.1.1.2.E (MARCH 1977) PAGE 5
"BEGIN"
"PROCEDURE" DER(Y); "ARRAY" Y;
"BEGIN" "REAL" Y1, Y2;
Y1:= Y[1]; Y2:= Y[2];
Y[1]:= -1000 * Y1 * (Y1 + Y2 - 1.999987);
Y[2]:= -2500 * Y2 * (Y1 + Y2 - 2)
"END" DER;
"PROCEDURE" JAC(J, Y); "ARRAY" J, Y;
"BEGIN" "REAL" Y1, Y2; Y1:= Y[1]; Y2:= Y[2];
J[1,1]:= 1999.987 - 1000 * (2 * Y1 + Y2);
J[1,2]:= -1000 * Y1;
J[2,1]:= -2500 * Y2;
J[2,2]:= 2500 * (2 - Y1 - 2 * Y2)
"END" JAC;
"PROCEDURE" OUTP;
"IF" X = 50 "THEN"
"BEGIN" "REAL" YE1, YE2;
YE1:= .5976546988; YE2:= 1.4023434075;
OUTPUT(61, "("
"("X = ")", 2D2B,
"("N = ")", 3ZD2B,
"("JEV = ")", 3ZD2B,
"("LU = ")", 3ZD, 2/,
"("Y1 = ")", +.13D"+2D2B,
"("REL. ERR. = ")", .2D"+2D, /,
"("Y2 = ")", +.13D"+2D2B,
"("REL. ERR. = ")", .2D"+2D")",
X, N, JEV, LU, Y[1], ABS((Y[1] - YE1) / YE1),
Y[2], ABS((Y[2] - YE2) / YE2))
"END" OUTP;
"INTEGER" N, JEV, LU; "REAL" X;
"ARRAY" Y[1:2];
"PROCEDURE" GMS(X, XE, R, Y, H, HMIN, HMAX, DELTA,
DERIVATIVE, JACOBIAN, AETA, RETA, N,
JEV, LU, NSJEV, LINEAR, OUT); "CODE" 33191;
GMS(X, 50, 2, Y, .01, .001, .5, -"15,
DER, JAC, "-5, "-5, N, JEV,
LU, 0, "FALSE", OUTP)
"END"
THIS PROGRAM DELIVERS:
X = 50 N = 109 JEV = 3 LU = 12
Y1 = +.5976547958004"+00 REL. ERR. = .16"-06
Y2 = +.1402343310813"+01 REL. ERR. = .69"-07
1SECTION : 5.2.1.1.1.2.E (MARCH 1977) PAGE 6
SOURCE TEXT:
"CODE" 33191;
"PROCEDURE" GMS(X, XE, R, Y, H, HMIN, HMAX, DELTA, DERIVATIVE,
JACOBIAN, AETA, RETA, N, JEV, LU, NSJEV,
LINEAR, OUT);
"VALUE" R;
"REAL" X, XE, H, HMIN, HMAX, DELTA, AETA, RETA;
"INTEGER" R, N, JEV, NSJEV, LU;
"BOOLEAN" LINEAR;
"ARRAY" Y;
"PROCEDURE" DERIVATIVE, JACOBIAN, OUT;
"BEGIN"
"INTEGER" I, J, K, L, NSJEV1, COUNT, COUNT1, KCHANGE;
"REAL" A, A1, ALFA, E, S1, S2, Z1, X0, XL0, XL1,
XL2, ETA, H0, H1, Q, Q1, Q2, Q12, Q22, Q1Q2, DISCR;
"BOOLEAN" UPDATE, CHANGE, REEVAL, STRATEGY;
"INTEGER" "ARRAY" RI, CI[1:R];
"ARRAY" AUX[1:9], BD1, BD2[1:3,1:3], Y1,
Y0[1:R], HJAC, H2JAC2, RQZ[1:R,1:R], YL, FL[1:3 * R];
"REAL" "PROCEDURE" VECVEC(L, U, SHIFT, A, B); "CODE" 34010;
"REAL" "PROCEDURE" MATVEC(L, U, I, A, B); "CODE" 34011;
"REAL" "PROCEDURE" MATMAT(L, U, I, J, A, B); "CODE" 34013;
"PROCEDURE" ELMROW(L, U, I, J, A, B, X); "CODE" 34024;
"PROCEDURE" ELMVEC(L, U, SHIFT, A, B, X); "CODE" 34020;
"PROCEDURE" DUPVEC(L, U, SHIFT, A, B); "CODE" 31030;
"PROCEDURE" GSSELM(A, N, AUX, RI, CI); "CODE" 34231;
"PROCEDURE" SOLELM(A, N, RI, CI, B); "CODE" 34061;
"PROCEDURE" COLCST(L, U, J, A, X); "CODE" 31131;
"PROCEDURE" MULVEC(L, U, SHIFT, A, B, X); "CODE" 31020;
"PROCEDURE" INITIALIZATION;
"BEGIN" LU:= JEV:= N:= NSJEV1:= KCHANGE:= 0; X0:= X; DISCR:= 0;
K:=1; H1:= H0:= H; COUNT:= -2; AUX[2]:= "-14; AUX[4]:= 8;
DUPVEC(1, R, 0, YL, Y); REEVAL:= CHANGE:= "TRUE";
STRATEGY:= HMIN ^= HMAX "AND" ^LINEAR; Q1:= -1; Q2:= -2;
COUNT1:= 0; XL0:= XL1:= XL2:= 0
"END" INITIALIZATION; "COMMENT"
1SECTION : 5.2.1.1.1.2.E (OCTOBER 1974) PAGE 7
;
"PROCEDURE" COEFFICIENT;
"BEGIN" XL2:= XL1; XL1:= XL0; XL0:= X0;
"IF" CHANGE "THEN"
"BEGIN" "IF" N > 2 "THEN"
"BEGIN" Q1:= (XL1 - XL0) / H1;
Q2:= (XL2 - XL0) / H1
"END";
Q12:= Q1 * Q1; Q22:= Q2 * Q2; Q1Q2:= Q1 * Q2;
A:= -(3 * ALFA + 1) / 12;
BD1[1,3]:= 1 + (1 / 3 - (Q1 + Q2) * .5) / Q1Q2;
BD1[2,3]:= (1 / 3 - Q2 * .5) / (Q12 - Q1Q2);
BD1[3,3]:= (1 / 3 - Q1 * .5) / (Q22 - Q1Q2);
BD2[1,3]:= -ALFA * .5 + A * (1 - Q1 - Q2) / Q1Q2;
BD2[2,3]:= A * (1 - Q2) / (Q12 - Q1Q2);
BD2[3,3]:= A * (1 - Q1) / (Q22 - Q1Q2);
"IF" STRATEGY "OR" N <= 2 "THEN"
"BEGIN" BD1[2,2]:= 1 / (2 * Q1);
BD1[1,2]:= 1 - BD1[2,2];
BD2[2,2]:= -(3 * ALFA + 1) / (12 * Q1);
BD2[1,2]:= -BD2[2,2] - ALFA * .5
"END"
"END"
"END" COEFFICIENT;
"PROCEDURE" OPERATOR CONSTRUCTION;
"BEGIN" "IF" REEVAL "THEN"
"BEGIN" JACOBIAN(HJAC, Y);
JEV:= JEV + 1; NSJEV1:= 0;
"IF" DELTA <= -"15 "THEN" ALFA:= 1 / 3 "ELSE"
"BEGIN" Z1:= H1 * DELTA;
A:= Z1 * Z1 + 12; A1:= 6 * Z1;
"IF" ABS(Z1) < .1 "THEN"
ALFA:= (Z1 * Z1 / 140 - 1) * Z1 / 30 "ELSE"
"IF" Z1 < -33 "THEN"
ALFA:= (A + A1) / (3 * Z1 * (2 + Z1)) "ELSE"
"BEGIN" E:= EXP(Z1); ALFA:= ((A - A1) *
E - A - A1) / (((2 - Z1) * E - 2 - Z1) *
Z1 * 3)
"END"
"END";
S1:= -(1 + ALFA) * .5; S2:= (ALFA * 3 + 1) / 12
"END"; "COMMENT"
1SECTION : 5.2.1.1.1.2.E (OCTOBER 1974) PAGE 8
;
A:= H1 / H0; A1:= A * A;
"IF" REEVAL "THEN" A:= H1;
"IF" A ^= 1 "THEN"
"FOR" J:= 1 "STEP" 1 "UNTIL" R "DO"
COLCST(1, R, J, HJAC, A);
"FOR" I:= 1 "STEP" 1 "UNTIL" R "DO"
"BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" R "DO"
"BEGIN" Q:= H2JAC2[I,J]:= "IF" REEVAL "THEN"
MATMAT(1, R, I, J, HJAC, HJAC)
"ELSE" H2JAC2[I,J] * A1;
RQZ[I,J]:= S2 * Q
"END";
RQZ[I,I]:= RQZ[I,I] + 1;
ELMROW(1, R, I, I, RQZ, HJAC, S1)
"END";
GSSELM(RQZ, R, AUX, RI, CI); LU:= LU + 1;
REEVAL:= UPDATE:= "FALSE"
"END" OPERATOR CONSTRUCTION;
"PROCEDURE" DIFFERENCE SCHEME;
"BEGIN" "IF" COUNT ^= 1 "THEN"
"BEGIN" DUPVEC(1, R, 0, FL, YL);
DERIVATIVE(FL); N:= N + 1; NSJEV1:= NSJEV1 + 1
"END";
MULVEC(1, R, 0, Y0, YL, (1 - ALFA) / 2 - BD1[1,K]);
"FOR" L:= 2 "STEP" 1 "UNTIL" K "DO"
ELMVEC(1, R, R * (L - 1), Y0, YL, -BD1[L,K]);
"FOR" L:= 1 "STEP" 1 "UNTIL" K "DO"
ELMVEC(1, R, R * (L - 1), Y0, FL, H1 * BD2[L,K]);
"FOR" I:= 1 "STEP" 1 "UNTIL" R "DO"
Y[I]:= MATVEC(1, R, I, HJAC, Y0);
MULVEC(1, R, 0, Y0, YL, (1 - 3 * ALFA) / 12 - BD2[1,K]);
"FOR" L:= 2 "STEP" 1 "UNTIL" K "DO"
ELMVEC(1, R, R * (L - 1), Y0, YL, -BD2[L,K]);
"FOR" I:= 1 "STEP" 1 "UNTIL" R "DO"
Y[I]:= Y[I] + MATVEC(1, R, I, H2JAC2, Y0);
DUPVEC(1, R, 0, Y0, YL);
"FOR" L:= 1 "STEP" 1 "UNTIL" K "DO"
ELMVEC(1, R, R * (L - 1), Y0, FL, H1 * BD1[L,K]);
ELMVEC(1, R, 0, Y, Y0, 1); SOLELM(RQZ, R, RI, CI, Y)
"END" DIFFERENCE SCHEME;
"PROCEDURE" NEXT INTEGRATION STEP;
"BEGIN" "FOR" L:= 2, 1 "DO"
"BEGIN" DUPVEC(L * R + 1, (L + 1) * R, -R, YL, YL);
DUPVEC(L * R + 1, (L + 1) * R, -R, FL, FL)
"END";
DUPVEC(1, R, 0, YL, Y)
"END" NEXT INTEGRATION STEP; "COMMENT"
1SECTION : 5.2.1.1.1.2.E (OCTOBER 1974) PAGE 9
;
"PROCEDURE" TEST ACCURACY;
"BEGIN" K:= 2;
DUPVEC(1, R, 0, Y1, Y); DIFFERENCE SCHEME; K:= 3;
ETA:= AETA + RETA * SQRT(VECVEC(1, R, 0, Y1, Y1));
ELMVEC(1, R, 0, Y, Y1, -1);
DISCR:= SQRT(VECVEC(1, R, 0, Y, Y));
DUPVEC(1, R, 0, Y, Y1)
"END" TEST ACCURACY;
"PROCEDURE" STEPSIZE;
"BEGIN" X0:= X; H0:= H1;
"IF" N <= 2 "AND" ^LINEAR "THEN" K:= K + 1;
"IF" COUNT = 1 "THEN"
"BEGIN" A:= ETA / (.75 * (ETA + DISCR)) + .33;
H1:= "IF" A <= .9 "OR" A >= 1.1 "THEN" A * H0
"ELSE" H0; COUNT:= 0;
REEVAL:= A <= .9 "AND" NSJEV1 ^= 1;
COUNT1:= "IF" A >= 1 "OR" REEVAL "THEN" 0 "ELSE"
COUNT1 + 1; "IF" COUNT1 = 10 "THEN"
"BEGIN" COUNT1:= 0; REEVAL:= "TRUE";
H1:= A * H0
"END"
"END" "ELSE"
"BEGIN" H1:= H; REEVAL:= NSJEV = NSJEV1 "AND"
^STRATEGY "AND" ^LINEAR
"END";
"IF" STRATEGY "THEN" H1:= "IF" H1 > HMAX
"THEN" HMAX "ELSE" "IF" H1 < HMIN "THEN" HMIN "ELSE" H1;
X:= X + H1; "IF" X >= XE "THEN"
"BEGIN" H1:= XE - X0; X:= XE "END";
"IF" N <= 2 "AND" ^LINEAR "THEN" REEVAL:= "TRUE";
"IF" H1 ^= H0 "THEN"
"BEGIN" UPDATE:= "TRUE"; KCHANGE:= 3 "END";
"IF" REEVAL "THEN" UPDATE:= "TRUE";
CHANGE:= KCHANGE > 0 "AND" ^LINEAR;
KCHANGE:= KCHANGE - 1
"END" STEPSIZE;
INITIALIZATION; OUT; X:= X + H1;
OPERATOR CONSTRUCTION;
BD1[1,1]:= 1; BD2[1,1]:= -ALFA * .5;
"IF" ^LINEAR "THEN" COEFFICIENT;
NEXT STEP: DIFFERENCE SCHEME;
"IF" STRATEGY "THEN" COUNT:= COUNT + 1;
"IF" COUNT = 1 "THEN" TEST ACCURACY;
OUT; "IF" X >= XE "THEN" "GOTO" END;
STEPSIZE; "IF" UPDATE "THEN" OPERATOR CONSTRUCTION;
"IF" ^LINEAR "THEN" COEFFICIENT;
NEXT INTEGRATION STEP; "GOTO" NEXT STEP;
END:
"END" GMS;
"EOP"
1SECTION : 5.1.1.1.2 � (OCTOBER 1975) PAGE 1
AUTHOR: T.J. DEKKER.
CONTRIBUTORS: T.J. DEKKER AND T.H.P. REYMER.
INSTITUTE: UNIVERSITY OF AMSTERDAM.
RECEIVED: 750615.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS ONE PROCEDURE FOR FINDING A ZERO OF A GIVEN
DIFFERENTIABLE FUNCTION IN A GIVEN INTERVAL;
ZEROINDER APPROXIMATES A ZERO MAINLY BY MEANS OF CONFLUENT 3-POINT
RATIONAL INTERPOLATION USING NOT ONLY VALUES OF THE
GIVEN FUNCTION BUT ALSO OF ITS DERIVATIVE.
ZEROINDER IS TO PREFER TO ZEROIN OR ZEROINRAT (SECTION 5.1.1.1.1),
IF THE DERIVATIVE IS (MUCH) CHEAPER TO EVALUATE THAN THE FUNCTION.
KEYWORDS:
ZERO SEARCHING,
ANALYTIC EQUATIONS,
SINGLE NONLINEAR EQUATION,
DERIVATIVE AVAILABLE.
CALLING SEQUENCE:
THE HEADING OF THIS PROCEDURE READS:
"BOOLEAN" "PROCEDURE" ZEROINDER(X, Y, FX, DFX, TOLX);
"REAL" X, Y, FX, DFX, TOLX;
ZEROINDER SEARCHES FOR A ZERO OF A DIFFERENTIABLE REAL FUNCTION F
DEFINED ON A CERTAIN INTERVAL J;
ZEROINDER := "TRUE" WHEN A (SUFFICIENTLY SMALL) SUBINTERVAL
OF J CONTAINING A ZERO OF F HAS BEEN FOUND; OTHERWISE,
ZEROINDER := "FALSE";
LET DF AND T DENOTE REAL FUNCTIONS DEFINED ON J, WHERE DF
IS THE DERIVATIVE OF F AND T A TOLERANCE FUNCTION DEFINING
THE REQUIRED PRECISION OF THE ZERO; ( FOR INSTANCE
T(X) = ABS(X) * RE + AE,
WHERE RE AND AE ARE THE REQUIRED RELATIVE AND ABSOLUTE
PRECISION RESPECTIVELY);
1SECTION : 5.1.1.1.2 � (OCTOBER 1975) PAGE 2
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <REAL VARIABLE>;
A JENSEN VARIABLE; THE ACTUAL PARAMETERS FOR FX, DFX AND
TOLX (MAY) DEPEND ON THE ACTUAL PARAMETER FOR X;
ENTRY: ONE ENDPOINT OF INTERVAL J IN WHICH A ZERO IS
SEARCHED FOR;
EXIT: A VALUE APPROXIMATING THE ZERO WITHIN THE TOLERANCE
2 * T(X) WHEN ZEROINDER HAS THE VALUE "TRUE", AND A
PRESUMABLY WORTHLESS ARGUMENT VALUE OTHERWISE;
Y: <REAL VARIABLE>;
ENTRY: THE OTHER ENDPOINT OF INTERVAL J IN WHICH A ZERO IS
SEARCHED FOR; UPON ENTRY X < Y AS WELL AS Y < X IS
ALLOWED;
EXIT: THE OTHER STRADDLING APPROXIMATION OF THE ZERO,
I.E. UPON EXIT THE VALUES OF Y AND X SATISFY
1. F(X) * F(Y) <= 0, 2. ABS(X - Y) <= 2 * T(X) AND
3. ABS(F(X)) <= ABS(F(Y)) WHEN ZEROINDER HAS THE
VALUE "TRUE", AND A PRESUMABLY WORTHLESS ARGUMENT
VALUE SATISFYING CONDITIONS 2 AND 3 BUT NOT 1
OTHERWISE;
FX: <ARITHMETIC EXPRESSION>;
DEFINES FUNCTION F AS A FUNCTION DEPENDING ON THE ACTUAL
PARAMETER FOR X;
DFX: <ARITHMETIC EXPRESSION>;
DEFINES THE DERIVATIVE DF OF F AS A FUNCTION DEPENDING ON
THE ACTUAL PARAMETER FOR X;
TOLX: <ARITHMETIC EXPRESSION>;
DEFINES THE TOLERANCE FUNCTION T WHICH MAY DEPEND ON THE
ACTUAL PARAMETER FOR X;
ONE SHOULD CHOOSE TOLX POSITIVE AND NEVER SMALLER THAN THE
PRECISION OF THE MACHINE'S ARITHMETIC AT X, I.E. IN THIS
ARITHMETIC X + TOLX AND X - TOLX SHOULD ALWAYS YIELD
VALUES DISTINCT FROM X; OTHERWISE THE PROCEDURE MAY GET
INTO A LOOP.
PROCEDURES USED: DWARF = CP30003;
REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED IN ZEROINDER.
RUNNING TIME: SEE METHOD AND PERFORMANCE.
LANGUAGE: ALGOL 60.
1SECTION : 5.1.1.1.2 � (OCTOBER 1975) PAGE 3
METHOD AND PERFORMANCE:
THE METHOD USED IS CONFLUENT 3-POINT RATIONAL INTERPOLATION, I.E.
THE INTERPOLATION FUNCTION OF THE FORM (X - A) / (BX + C) IS FITTED
EXACTLY AT 3 POINTS 2 OF WHICH ARE COINCIDING; MOREOVER, IN ORDER
TO IMPROVE THE GLOBAL BEHAVIOUR OF THE PROCESS, LINEAR
INTERPOLATION ON THE FUNCTION F / DF OR BISECTION ARE INCLUDED IN
SOME STEPS;
THE PERFORMANCE IS AS FOLLOWS:
THE NUMBER OF EVALUATIONS OF FX, DFX AND TOLX IS AT MOST
4 LOG(ABS(X - Y)) / TAU,
WHERE X AND Y ARE THE ARGUMENT VALUES GIVEN UPON ENTRY, LOG DENOTES
THE BASE 2 LOGARITHM, AND TAU IS THE MINIMUM OF THE TOLERANCE
FUNCTION ON J (I.E. ZEROINDER REQUIRES AT MOST 4 TIMES THE NUMBER
OF STEPS REQUIRED FOR BISECTION); IF UPON ENTRY X AND Y SATISFY
F(X) * F(Y) <= 0, THEN CONVERGENCE TO A ZERO IS GUARANTEED, SO THAT
UPON EXIT X AND Y SATISFY CONDITION 1 TO 3 MENTIONED ABOVE (SEE
PARAMETER Y) AND ZEROINDER HAS THE VALUE "TRUE";
FOR A SIMPLE ZERO OF F, THE ASYMPTOTIC ORDER OF CONVERGENCE
APPROXIMATELY EQUALS 2.414.
EXAMPLE OF USE:
THE ZERO OF THE FUNCTION EXP(-X * 3) * (X - 1) + X ** 3, IN THE
INTERVAL [0, 1], MAY BE CALCULATED BY THE FOLLOWING PROGRAM:
"BEGIN" "REAL" X, Y;
"BOOLEAN" "PROCEDURE" ZEROINDER(X,Y,FX,DFX,TOLX); "CODE" 34453;
"REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
F:= EXP(-X * 3) * ( X - 1) + X ** 3;
"REAL" "PROCEDURE" DF(X); "VALUE" X; "REAL" X;
DF:= EXP(-X * 3) * (-3 * X + 4) + 3 * X ** 2;
X:= 0; Y:= 1;
"IF" ZEROINDER(X, Y, F(X), DF(X), ABS(X) * "-14 + "-14) "THEN"
OUTPUT(71, "("/4B,"("CALCULATED ZERO AND FUNCTION VALUE:")",
/8B, 2(B+.15D"+3D4B), /4B,
"("OTHER STRADDLING APPROXIMATION AND FUNCTION VALUE:")",
/8B, 2(B+.15D"+3D4B)")", X, F(X), Y, F(Y))
"ELSE" OUTPUT(71, "("/4B, "("NO ZERO FOUND")"")")
"END"
RESULT:
CALCULATED ZERO AND FUNCTION VALUE:
+.489702748548240"+000 -.444089209850060"-015
OTHER STRADDLING APPROXIMATION AND FUNCTION VALUE:
+.489702748548250"+000 +.177635683940030"-013
1SECTION : 5.1.1.1.2 � (OCTOBER 1975) PAGE 4
REFERENCES:
[1]: BUS, J.C.P. AND DEKKER, T.J.,
TWO EFFICIENT ALGORITHMS WITH GUARANTEED CONVERGENCE FOR
FINDING A ZERO OF A FUNCTION.
MATHEMATICAL CENTRE, REPORT NW 13/74, AMSTERDAM (1974),
(ALSO TO APPEAR IN TOMS 1975).
[2]: OSTROWSKI, A.M.,
SOLUTION OF EQUATIONS AND SYSTEMS OF EQUATIONS.
ACADEMIC PRESS 1966. CHAPTERS 3 AND 11.
SOURCE TEXT(S):
0"CODE" 34453;
"BOOLEAN" "PROCEDURE" ZEROINDER(X, Y, FX, DFX, TOLX);
"REAL" X, Y, FX, DFX, TOLX;
"BEGIN" "INTEGER" EXT;
"REAL" B, FB, DFB, A, FA, DFA, C, FC, DFC, D, W, MB,
TOL, M, P, Q, DW;
"REAL" "PROCEDURE" DWARF; "CODE" 30003;
DW:= DWARF;
B:= X; FB:= FX; DFB:= DFX; A:= X:= Y; FA:= FX; DFA:= DFX;
INTERPOLATE: C:= A; FC:= FA; DFC:= DFA; EXT:= 0;
EXTRAPOLATE: "IF" ABS(FC) < ABS(FB) "THEN"
"BEGIN" A:= B; FA:= FB; DFA:= DFB; B:= X:= C; FB:= FC;
DFB:= DFC; C:= A; FC:= FA; DFC:= DFA
"END" INTERCHANGE;
TOL:= TOLX; M:= (C + B) * 0.5; MB:= M - B;
"IF" ABS(MB) > TOL "THEN"
"BEGIN" "IF" EXT > 2 "THEN" W:= MB "ELSE"
"BEGIN" TOL:= TOL * SIGN(MB);
D:= "IF" EXT = 2 "THEN" DFA "ELSE" (FB - FA) / (B - A);
P:= FB * D * (B - A);
Q:= FA * DFB - FB * D;
"IF" P < 0 "THEN" "BEGIN" P:= -P; Q:= -Q "END";
W:= "IF" P < DW "OR" P <= Q * TOL "THEN" TOL "ELSE"
"IF" P < MB * Q "THEN" P / Q "ELSE" MB;
"END"; A:= B; FA:= FB; DFA:= DFB;
X:= B:= B + W; FB:= FX; DFB:= DFX;
"IF" ("IF" FC >= 0 "THEN" FB >= 0 "ELSE" FB <= 0) "THEN"
"GOTO" INTERPOLATE "ELSE"
"BEGIN" EXT:= "IF" W = MB "THEN" 0 "ELSE" EXT + 1;
"GOTO" EXTRAPOLATE
"END"
"END"; Y:= C;
ZEROINDER:= "IF" FC >= 0 "THEN" FB <= 0 "ELSE" FB >= 0
"END" ZEROINDER;
"EOP"
1SECTION : 5.1.2.1.1 ( DECEMBER 1978 ) PAGE 1
AUTHOR : J. C. P. BUS
INSTITUTE : MATHEMATICAL CENTRE.
RECEIVED : 741101.
BRIEF DESCRIPTION :
THIS SECTION CONTAINS THE PROCEDURE MININ, FOR MINIMIZING
A FUNCTION OF ONE VARIABLE IN A GIVEN INTERVAL;
KEYWORDS :
MINIMIZATION,
FUNCTIONS OF ONE VARIABLE.
CALLING SEQUENCE :
THE HEADING OF THIS PROCEDURE IS :
"REAL" "PROCEDURE" MININ(X, A, B, FX, TOLX);
"REAL" X, A, B, FX, TOLX;
"CODE" 34433;
MININ DELIVERS THE CALCULATED MINIMUM VALUE OF THE FUNCTION,
DEFINED BY FX, ON THE INTERVAL WITH ENDPOINTS A AND B.
THE MEANING OF THE FORMAL PARAMETERS IS :
X : <REAL VARIABLE>;
A JENSEN VARIABLE; THE ACTUAL PARAMETERS FOR FX AND TOLX
DEPEND ON X;
EXIT : THE CALCULATED APPROXIMATION OF THE POSITION OF THE
MINIMUM;
A, B : <REAL VARIABLE>;
ENTRY : THE ENDPOINTS OF THE INTERVAL ON WHICH A MINIMUM
IS SEARCHED FOR;
EXIT : THE ENDPOINTS OF THE INTERVAL WITH LENGTH LESS THAN
4 * TOL(X) SUCH THAT A < X < B;
FX : <ARITHMATIC EXPRESSION>;
THE FUNCTION ISGIVEN BY THE ACTUAL PARAMETER FX, WHICH
DEPENDS ON X;
TOLX : <ARITHMETIC EXPRESSION>;
THE TOLERANCE IS GIVEN BY THE ACTUAL PARAMETER TOLX, WHICH MAY
DEPEND ON X; A SUITABLE TOLERANCE FUNCTION IS : ABS(X)*RE + AE,
WHERE RE IS THE RELATIVE PRECISION DESIRED AND AE IS AN ABSOLUTE
PRECISION WHICH SHOULD NOT BE CHOSEN EQUAL TO ZERO.
1SECTION : 5.1.2.1.1 ( DECEMBER 1978 ) PAGE 2
DATA AND RESULTS :
THE USER SHOULD BE AWARE OF THE FACT THAT THE CHOICE OF TOLX MAY
HIGHLY AFFECT THE BEHAVIOUR OF THE ALGORITHM, ALTHOUGH CONVERGENCE
TO A POINT FOR WHICH THE GIVEN FUNCTION IS MINIMAL ON THE INTERVAL
IS ASSURED; THE ASYMPTOTIC BEHAVIOUR WILL USUALLY BE FINE AS LONG AS
THE NUMERICAL FUNCTION IS STRICTLY DELTA-UNIMODAL ON THE GIVEN
INTERVAL (SEE [1]) AND THE TOLERANCE FUNCTION SATISFIES TOL(X)>=DELTA,
FOR ALL X IN THE GIVEN INTERVAL.
PROCEDURES USED : NONE.
REQUIRED CENTRAL MEMORY : NO AUXILIARY ARRAYS ARE DECLARED IN MININ.
METHOD AND PERFORMANCE :
MININ IS A SLIGHTLY MODIFIED VERSION OF THE ALGORITHM GIVEN IN [1].
EXAMPLE OF USE:
THE FOLLOWING PROGRAM MAY BE USED TO CALCULATE THE MINIMUM OF THE
FUNCTION F(X) = SUM(((I * 2 - 5)/(X - I ** 2)) ** 2; I = 1 (1) 20)
ON THE INTERVAL [1 + TOL, 4 - TOL] (SEE [1]).
"BEGIN"
"REAL" "PROCEDURE" MININ(X, A, B, FX, TOLX); "CODE" 34433;
"REAL" M, X, A, B; "INTEGER" CNT;
"REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
"BEGIN" "INTEGER" I; "REAL" S;
S:= 0; "FOR" I:= 1 "STEP" 1 "UNTIL" 20 "DO"
S:= S + ((I * 2 - 5) / (X - I ** 2)) ** 2;
CNT:= CNT + 1; F:= S
"END" F;
"REAL" "PROCEDURE" TOL(X); "VALUE" X; "REAL" X;
TOL:= ABS(X) * "-7 + "-7;
A:= 1 + TOL(1); B:= 4 - TOL(4); CNT:= 0;
M:= MININ(X, A, B, F(X), TOL(X));
OUTPUT(61, "("4B,"("MINIMUM IS ")",N,/4B,
"("FOR X IS ")",N,/4B,
"("IN THE INTERVAL WITH ENDPOINTS ")",/8B,2(N),/4B,
"("THE NUMBER OF FUNCTION EVALUATIONS NEEDED IS ")",2ZD,/")",
M, X, A, B, CNT)
"END"
1SECTION : 5.1.2.1.1 ( DECEMBER 1978 ) PAGE 3
RESULTS :
MINIMUM IS +3.6766990169019"+000
FOR X IS +3.0229153387991"+000
IN THE INTERVAL WITH ENDPOINTS
+3.0229149365075"+000 +3.0229157410906"+000
THE NUMBER OF FUNCTION EVALUATIONS NEEDED IS 13
SOURCE TEXT:
"CODE" 34433;
"REAL" "PROCEDURE" MININ(X, A, B, FX, TOLX);
"REAL" X, A, B, FX, TOLX;
"BEGIN" "COMMENT" SEE BRENT, 1973, P79;
"REAL" Z, C, D, E, M, P, Q, R, TOL, T, U, V, W, FU, FV, FW, FZ;
C:= (3 - SQRT(5)) / 2; "IF" A > B "THEN"
"BEGIN" Z:= A; A:= B; B:= Z "END";
W:= X:= A; FW:= FX; Z:= X:= B; FZ:= FX; "IF" FZ > FW "THEN"
"BEGIN" Z:= W; W:= X; V:= FZ; FZ:= FW; FW:= V "END";
V:= W; FV:= FW; E:= 0;
LOOP: M:= (A + B) * 0.5; TOL:= TOLX; T:= TOL * 2;
"IF" ABS(Z - M) > T - (B - A) * 0.5 "THEN"
"BEGIN" P:= Q:= R:= 0; "IF" ABS(E) > TOL "THEN"
"BEGIN" R:= (Z - W) * (FZ - FV);
Q:= (Z - V) * (FZ - FW); P:= (Z - V) * Q - (Z - W) * R;
Q:= (Q - R) * 2; "IF" Q>0 "THEN" P:= -P "ELSE" Q:= -Q;
R:= E; E:= D
"END";
"IF" ABS(P) < ABS(Q * R * 0.5) "AND" P > (A - Z) * Q
"AND" P < (B - Z) * Q "THEN"
"BEGIN" D:= P / Q; U:= Z + D;
"IF" U - A < T "OR" B - U < T "THEN"
D:= "IF" Z < M "THEN" TOL "ELSE" -TOL
"END" "ELSE"
"BEGIN" E:= ("IF" Z < M "THEN" B "ELSE" A) - Z; D:= C * E
"END";
U:= X:= Z + ("IF" ABS(D) >= TOL "THEN" D "ELSE" "IF" D > 0
"THEN" TOL "ELSE" -TOL); FU:= FX;
"IF" FU <= FZ "THEN"
"BEGIN" "IF" U < Z "THEN" B:= Z "ELSE" A:= Z;
V:= W; FV:= FW; W:= Z; FW:= FZ; Z:= U; FZ:= FU
"END" "ELSE"
"BEGIN" "IF" U < Z "THEN" A:= U "ELSE" B:= U;
"IF" FU <= FW "THEN"
"BEGIN" V:= W; FV:= FW; W:= U; FW:= FU "END" "ELSE"
"IF" FU <= FV "OR" V = W "THEN"
"BEGIN" V:= U; FV:= FU "END"
"END"; "GOTO" LOOP
"END"; X:= Z; MININ:= FZ
"END" MININ;
1SECTION : 5.1.2.1.2 � (OCTOBER 1975) PAGE 1
AUTHOR: J. C. P. BUS.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 741101.
BRIEF DESCRIPTION:�
THIS SECTION CONTAINS THE PROCEDURE MININDER, FOR MINIMIZING A
FUNCTION OF ONE VARIABLE IN A GIVEN INTERVAL, WHEN THE ANALYTICAL
DERIVATIVE OF THE FUNCTION IS AVAILABLE.
KEYWORDS :
MINIMIZATION,
FUNCTIONS OF ONE VARIABLE.
CALLING SEQUENCE:
THE HEADING OF THIS PROCEDURE IS:
"REAL" "PROCEDURE" MININDER(X, Y, FX, DFX, TOLX);
"REAL" X, Y, FX, DFX, TOLX;
MININDER DELIVERS THE CALCULATED MINIMUM VALUE OF THE FUNCTION,
DEFINED BY FX, ON THE INTERVAL WITH END POINTS A AND B.
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <REAL VARIABLE>;�
A JENSEN VARIABLE; THE ACTUAL PARAMETERS FOR FX, DFX AND�
TOLX DEPEND ON X;
ENTRY: ONE OF THE END POINTS OF THE INTERVAL ON WHICH THE
FUNCTION HAS TO BE MINIMIZED;
EXIT: THE CALCULATED APPROXIMATION OF THE POSITION OF THE
MINIMUM;�
Y: <REAL VARIABLE>;�
ENTRY: THE OTHER END POINT OF THE INTERVAL ON WHICH THE�
FUNCTION HAS TO BE MINIMIZED;
EXIT: A VALUE SUCH THAT ABS(X - Y) <= TOL(X) * 3;
FX: <ARITHMETIC EXPRESSION>;
THE FUNCTION IS GIVEN BY THE ACTUAL PARAMETER FX WHICH
DEPENDS ON X;
DFX: <ARITHMETIC EXPRESSION>;
THE DERIVATIVE OF THE FUNCTION IS GIVEN BY THE ACTUAL
PARAMETER DFX WHICH DEPENDS ON X; FX AND DFX ARE EVALUATED
SUCCESSIVELY FOR A CERTAIN VALUE OF X;
TOLX: <ARITHMETIC EXPRESSION>;
THE TOLERANCE IS GIVEN BY THE ACTUAL PARAMETER TOLX, WHICH
MAY DEPEND ON X; A SUITABLE TOLERANCE FUNCTION IS:
ABS(X) * RE + AE, WHERE RE IS THE RELATIVE PRECISION
DESIRED AND AE IS AN ABSOLUTE PRECISION WHICH SHOULD NOT
BE CHOSEN EQUAL TO ZERO.
1SECTION : 5.1.2.1.2 � (OCTOBER 1975) PAGE 2
DATA AND RESULTS:
THE USER SHOULD BE AWARE OF THE FACT THAT THE CHOICE OF TOLX
MAY HIGHLY AFFECT THE BEHAVIOUR OF THE ALGORITHM, ALTHOUGH
CONVERGENCE TO A POINT FOR WHICH THE GIVEN FUNCTION IS
MINIMAL ON THE GIVEN INTERVAL IS ASSURED; THE ASYMPTOTIC
BEHAVIOUR WILL USUALLY BE FINE AS LONG AS THE NUMERICAL FUNCTION IS
STRICTLY DELTA-UNIMODAL ON THE GIVEN INTERVAL (SEE [1]) AND THE
TOLERANCE FUNCTION SATISFIES TOL(X) >= DELTA, FOR ALL X IN THE
GIVEN INTERVAL; LET THE VALUE OF DFX AT THE BEGIN AND END POINT OF
THE INITIAL INTERVAL BE DENOTED BY DFA AND DFB, RESPECTIVELY, THEN,
FINDING A GLOBAL MINIMUM IS ONLY GUARANTEED IF THE FUNCTION IS
CONVEX AND DFA <= 0 AND DFB >= 0; IF THESE CONDITIONS ARE NOT
SATISFIED, THEN A LOCAL MINIMUM OR A MINIMUM AT ONE OF THE END
POINTS MIGHT BE FOUND.
PROCEDURES USED: NONE.
REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED IN MININDER.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
MININDER HAS ALMOST THE SAME STRUCTURE AS THE PROCEDURE GIVEN IN�
[1]; HOWEVER, CUBIC INTERPOLATION (SEE [2]) IS USED INSTEAD OF
QUADRATIC INTERPOLATION TO APPROXIMATE THE MINIMIUM.
REFERENCES:
[1]: BRENT, R.P.
ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES. CH.5.
PRENTICE HALL, 1973.
[2]: DAVIDON, W.C.
VARIABLE METRIC METHODS FOR MINIMIZATION.
REP. A.N.L. 5990, 1959.
1SECTION : 5.1.2.1.2 � (OCTOBER 1975) PAGE 3
EXAMPLE OF USE:
THE FOLLOWING PROGRAM MAY BE USED TO CALCULATE THE MINIMUM OF THE
FUNCTION F(X) = SUM(((I * 2 - 5)/(X - I ** 2)) ** 2; I = 1 (1) 20)
ON THE INTERVAL [1.01,3.99] (SEE [1]).
"BEGIN"
"REAL" "PROCEDURE" MININDER(X, Y, FX, DFX, TOLX); "CODE" 34435;
"REAL" M, X, Y; "INTEGER" CNT;�
"REAL" "PROCEDURE" F(X); "VALUE" X; "REAL" X;
"BEGIN" "INTEGER" I; "REAL" S;�
S:= 0; "FOR" I:= 1 "STEP" 1 "UNTIL" 20 "DO"
S:= S + ((I * 2 - 5) / (X - I ** 2)) ** 2;
CNT:= CNT + 1; F:= S
"END" F;
"REAL" "PROCEDURE" DF(X); "VALUE" X; "REAL" X;
"BEGIN" "INTEGER" I; "REAL" S;�
S:= 0; "FOR" I:= 1 "STEP" 1 "UNTIL" 20 "DO"
S:= S + (I * 2 - 5) ** 2 / (X - I ** 2) ** 3;
DF:= -S * 2
"END" DF;
"REAL" "PROCEDURE" TOL(X); "VALUE" X; "REAL" X;
TOL:= ABS(X) * "-7 + "-7;
X:= 1.01; Y:= 3.99; CNT:= 0;
M:= MININDER(X, Y, F(X), DF(X),TOL(X));
OUTPUT(61 ,"("4B,"("MINIMUM IS ")",N,/4B,
"("FOR X IS ")",N,/4B,
"("AND Y IS ")",N,/4B,
"("THE NUMBER OF FUNCTION EVALUATIONS NEEDED IS ")",2ZD,/")",
M, X, Y, CNT)
"END"
RESULTS:
MINIMUM IS +3.6766990169021"+000
FOR X IS +3.0229155250302"+000
AND Y IS +3.0229151227386"+000
THE NUMBER OF FUNCTION EVALUATIONS NEEDED IS 9
1SECTION : 5.1.2.1.2 (NOVEMBER 1976) PAGE 4
SOURCE TEXT(S):
"CODE"34435;
"REAL" "PROCEDURE" MININDER(X, Y, FX, DFX, TOLX);
"REAL" X, Y, FX, DFX, TOLX;
"BEGIN" "COMMENT" THE FUNCTION IS APPROXIMATED BY A CUBIC AS
NIVEN BY DAVIDON, 1958, THE STRUCTURE IS SIMILAR TO THE
STRUCTURE OF THE PROGRAM GIVEN BY BRENT, 1973, THIS IS
A REVISION OF 760407;
"INTEGER" SGN;
"REAL" A, B, C, FA, FB, FU, DFA, DFB, DFU, E, D, TOL, BA,
Z, P, Q, S;
"IF" X <= Y "THEN"
"BEGIN" A:= X; FA:= FX; DFA:= DFX;
B:= X:= Y; FB:= FX; DFB:= DFX
"END" "ELSE"
"BEGIN" B:= X; FB:= FX; DFB:= DFX;
A:= X:= Y; FA:= FX; DFA:= DFX
"END";
C:= (3 - SQRT(5)) / 2; D:= B - A; E:= D * 2; Z:= E * 2;
LOOP: BA:= B - A; TOL:= TOLX; "IF" BA >= TOL * 3 "THEN"
"BEGIN" "IF" ABS(DFA) <= ABS(DFB) "THEN"
"BEGIN" X:=A; SGN:= 1 "END" "ELSE"
"BEGIN" X:= B; SGN:= -1 "END";
"IF" DFA <= 0 "AND" DFB >= 0 "THEN"
"BEGIN" Z:= (FA - FB) * 3 / BA + DFA + DFB;
S:= SQRT(Z ** 2 - DFA * DFB);
P:= "IF" SGN = 1 "THEN" DFA - S - Z "ELSE"
DFB + S - Z; P:= P * BA;
Q:= DFB - DFA + S * 2; Z:= E; E:= D;
D:= "IF" ABS(P) <= ABS(Q) * TOL "THEN" TOL * SGN
"ELSE" -P / Q
"END" "ELSE" D:= BA;
"IF" ABS(D) >= ABS(Z * 0.5) "OR" ABS(D) > BA * 0.5 "THEN"
"BEGIN" E:= BA; D:= C * BA * SGN "END";
X:= X + D; FU:= FX; DFU:= DFX;
"IF" DFU >= 0 "OR" (FU >= FA "AND" DFA <= 0) "THEN"
"BEGIN" B:= X; FB:= FU; DFB:= DFU "END" "ELSE"
"BEGIN" A:= X; FA:= FU; DFA:= DFU "END";
"GOTO" LOOP
"END"; "IF" FA <= FB "THEN"
"BEGIN" X:= A; Y:= B; MININDER:= FA "END" "ELSE"
"BEGIN" X:= B; Y:= A; MININDER:= FB "END"
"END" MININDER;
"EOP"
1SECTION : 5.1.2.2.2 (DECEMBER 1979) PAGE 1
AUTHOR: J.C.P. BUS.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 741101.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS THE PROCEDURE PRAXIS;
PRAXIS MINIMIZES A FUNCTION OF SEVERAL VARIABLES.
KEYWORDS:
MINIMIZATION,
FUNCTION OF SEVERAL VARIABLES.
CALLING SEQUENCE:
THE HEADING OF THIS PROCEDURE IS:
"PROCEDURE" PRAXIS(N, X, FUNCT, IN, OUT); "VALUE" N;
"INTEGER" N; "ARRAY" X, IN, OUT; "REAL" "PROCEDURE" FUNCT;"CODE" 34432;
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE NUMBER OF VARIABLES OF THE FUNCTION TO BE MINIMIZED;
X: <ARRAY IDENTIFIER>;
"ARRAY" X[1 : N];
THE VARIABLES OF THE FUNCTION;
ENTRY: AN APPROXIMATION OF THE POSITION OF THE MINIMUM;
EXIT: THE CALCULATED POSITION OF THE MINIMUM;
FUNCT: <PROCEDURE IDENTIFIER>;
THE HEADING OF THIS PROCEDURE SHOULD BE:
"REAL" "PROCEDURE" FUNCT(N, X); "VALUE" N;
"INTEGER" N; "ARRAY" X;
FUNCT SHOULD DELIVER THE VALUE OF THE FUNCTION TO BE
MINIMIZED, AT THE POINT GIVEN BY X[1:N];
THE MEANING OF THE FORMAL PARAMETERS IS:
N: <ARITHMETIC EXPRESSION>;
THE NUMBER OF VARIABLES;
X: <ARRAY IDENTIFIER>; "ARRAY" X[1:N];
THE VALUES OF THE VARIABLES FOR WHICH THE FUNCTION HAS
TO BE EVALUATED;
IN: <ARRAY IDENTIFIER>;
"ARRAY" IN[0:9];
ENTRY:
IN[0]: THE MACHINE PRECISION; FOR THE CYBER 73 A
SUITABLE VALUE IS "-14;
1SECTION : 5.1.2.2.2 � (OCTOBER 1975) PAGE 2
IN[1], IN[2]: RELATIVE AND ABSOLUTE TOLERANCE,
RESPECTIVELY, FOR THE STEPVECTOR
(RELATIVE TO THE CURRENT ESTIMATES OF
THE VARIABLES); THE PROCESS IS TERMINATED WHEN
IN IN[8] + 1 SUCCESSIVE ITERATION STEPS THE
EUCLIDEAN NORM OF THE STEP VECTOR IS LESS THAN
(IN[1] * NORM(X) + IN[2]) * 0.5;
IN[1] SHOULD BE CHOSEN IN AGREEMENT WITH THE
PRECISION IN WHICH THE FUNCTION IS CALCULATED;
USUALLY IN[1] SHOULD BE CHOSEN SUCH THAT
IN[1] >= SQRT(IN[0]); IN[0] SHOULD BE CHOSEN
DIFFERENT FROM ZERO.
IN[3], IN[4] ARE NEITHER USED NOR CHANGED;
IN[5]: THE MAXIMUM NUMBER OF FUNCTION EVALUATIONS
ALLOWED (I.E. CALLS OF FUNCT);
IN[6]: THE MAXIMUM STEP SIZE; IN[6] SHOULD BE EQUAL TO
THE MAXIMUM EXPECTED DISTANCE BETWEEN THE GUESS
AND THE MINIMUM; IF IN[6] IS TOO SMALL OR TOO
LARGE, THEN THE INITIAL RATE OF CONVERGENCE
WILL BE SLOW;
IN[7]: THE MAXIMUM SCALING FACTOR; THE VALUE OF IN[7]
MAY BE USED TO OBTAIN AUTOMATIC SCALING OF THE
VARIABLES; HOWEVER, THIS SCALING IS WORTHWHILE
BUT MAY BE UNRELIABLE; THEREFORE, THE USER
SHOULD TRY TO SCALE HIS PROBLEM HIMSELF AS WELL
AS POSSIBLE AND SET IN[7]:= 1; IN EITHER CASE,
IN[7] SHOULD NOT BE CHOSEN GREATER THAN 10;
IN[8]: THE PROCESS TERMINATES IF NO SUBSTANTIAL
IMPROVEMENT OF THE VALUES OF THE VARIABLES IS
OBTAINED IN IN[8] + 1 SUCCESSIVE ITERATION
STEPS (SEE IN[1], IN[2]); IN[8] = 4 IS VERY
CAUTIOUS; USUALLY, IN[8] = 1 IS SATISFACTORY;
IN[9]: IF THE PROBLEM IS KNOWN TO BE ILL-CONDITIONED
(SEE [1]), THEN THE VALUE OF IN[9] SHOULD BE
NEGATIVE, OTHERWISE IN[9] >= 0;
OUT: <ARRAY IDENTIFIER>;
"ARRAY" OUT[1:6];
EXIT:
OUT[1]: THIS VALUE GIVES INFORMATION ABOUT THE
TERMINATION OF THE PROCESS;
OUT[1] = 0: NORMAL TERMINATION;
OUT[1] = 1: THE PROCESS IS BROKEN OFF, BECAUSE,
AT THE END OF AN ITERATION STEP,
THE NUMBER OF CALLS OF FUNCT
EXCEEDED THE VALUE GIVEN IN IN[5];
OUT[1] = 2: THE PROCESS IS BROKEN OFF, BECAUSE
THE CONDITION OF THE PROBLEM IS TOO
BAD;
OUT[2]: THE CALCULATED MINIMUM OF THE FUNCTION;
OUT[3]: THE VALUE OF THE FUNCTION AT THE INITIAL GUESS;
OUT[4]: THE NUMBER OF FUNCTION EVALUATIONS NEEDED TO
OBTAIN THIS RESULT;
OUT[5]: THE NUMBER OF LINE SEARCHES (SEE [1]);
OUT[6]: THE STEP SIZE IN THE LAST ITERATION STEP.
1SECTION : 5.1.2.2.2 (DECEMBER 1979) PAGE 3
PROCEDURES USED:
INIVEC = CP31010,
INIMAT = CP31011,
DUPVEC = CP31030,
DUPMAT = CP31035,
DUPCOLVEC = CP31034,
MULROW = CP31021,
MULCOL = CP31022,
VECVEC = CP34010,
TAMMAT = CP34014,
MATTAM = CP34015,
ICHROWCOL = CP34033,
ELMVECCOL = CP34021,
QRISNGVALDEC = CP34273,
SETRANDOM = CP11014,
RANDOM = CP11015,
DWARF = CP30003.
REQUIRED CENTRAL MEMORY:
EXECUTION FIELD LENGTH: ONE ARRAY OF LENGTH N SQUARED AND FIVE
ARRAYS OF LENGTH N ARE DECLARED;
RUNNING TIME:
THE NUMBER OF ITERATION STEPS DEPENDS STRONGLY ON THE PROBLEM TO BE
SOLVED.
METHOD AND PERFORMANCE:
THIS PROCEDURE IS ADOPTED FROM [1].
REFERENCES:
[1] R. P. BRENT,
ALGORITHMS FOR MINIMIZATION WITHOUT DERIVATIVES, CH. 7.
PRENTICE HALL, 1973.
1SECTION : 5.1.2.2.2 (DECEMBER 1979) PAGE 4
EXAMPLE OF USE:
THE FOLLOWING PROGRAM MAY BE USED TO CALCULATE THE MINIMUM OF THE
FUNCTION F(X) = 100 * (X[2] - X[1] ** 2) ** 2 + (1 - X[1]) ** 2,
USING (-1.2, 1) AS AN INITIAL ESTIMATE.
"BEGIN"
"PROCEDURE" PRAXIS(N, X, FUNCT, IN, OUT); "CODE" 34432;
"ARRAY" X[1:2], IN[0:9], OUT[1:6];
"REAL" "PROCEDURE" F(N, X); "VALUE" N; "INTEGER" N; "ARRAY" X;
F:= (X[2] - X[1] ** 2) ** 2 * 100 + (1 - X[1]) ** 2;
IN[0]:= "-14; IN[1]:= IN[2]:= "-6; IN[5]:= 250;
IN[6]:= 1; IN[7]:= 1; IN[8]:= 1; IN[9]:= 1;
X[1]:= -1.2; X[2]:= 1;
PRAXIS(2, X, F, IN, OUT);
"IF" OUT[1] = 0 "THEN" OUTPUT(61,"(""(" NORMAL TERMINATION")"
,//")");
OUTPUT(61, "("4B,"("MINIMUM IS ")",N,/,4B,
"("FOR X IS ")", 2(N),/,4B,
"("THE INITIAL FUNCTION VALUE WAS ")" ,N,/,4B,
"("THE NUMBER OF FUNCTION EVALUATIONS NEEDED WAS ")",3ZD,/,4B,
"("THE NUMBER OF LINE SEARCHES WAS ")",3ZD,/,4B,
"("THE STEP SIZE IN THE LAST ITERATION STEP WAS ")", N,/")",
OUT[2], X[1], X[2], OUT[3], OUT[4], OUT[5], OUT[6])
"END"
RESULTS:
NORMAL TERMINATION
MINIMUM IS +1.5694986738789"-021
FOR X IS +1.0000000000389"+000 +1.0000000000785"+000
THE INITIAL FUNCTION VALUE WAS +2.4200000000001"+001
THE NUMBER OF FUNCTION EVALUATIONS NEEDED WAS 189
THE NUMBER OF LINE SEARCHES WAS 72
THE STEP SIZE IN THE LAST ITERATION STEP WAS +5.3830998470105"-009
1SECTION : 5.1.2.2.2 (DECEMBER 1979) PAGE 5
SOURCE TEXT(S):
0"CODE" 34432;
"PROCEDURE" PRAXIS(N, X, FUNCT, IN, OUT);
"VALUE" N; "INTEGER" N;
"ARRAY" X, IN, OUT;
"REAL" "PROCEDURE" FUNCT;
"BEGIN"
"COMMENT"THIS PROCEDURE MINIMIZES FUNCT(N,X),WITH THE
PRINCIPAL AXIS METHOD (SEE BRENT,R.P, 1973, ALGORITHMS
FOR MINIMIZATION WITHOUT DERIVATIVES,CH.7);
"PROCEDURE" INIVEC(L, U, A, X); "CODE" 31010;
"PROCEDURE" INIMAT(L, U, K, V, A, X); "CODE" 31011;
"PROCEDURE" DUPVEC(L, U, K, A, X); "CODE" 31030;
"PROCEDURE" DUPMAT(L, U, K, V, A, B); "CODE" 31035;
"PROCEDURE" DUPCOLVEC(L, U, K, A, B); "CODE" 31034;
"PROCEDURE" MULROW(L, U, I, J, A, B, X); "CODE" 31021;
"PROCEDURE" MULCOL(L, U, I, J, A, B, X); "CODE" 31022;
"REAL" "PROCEDURE" VECVEC(L, U, S, A, B); "CODE" 34010;
"REAL" "PROCEDURE" TAMMAT(L, U, I, J, A, B); "CODE" 34014;
"REAL" "PROCEDURE" MATTAM(L, U, I, J, A, B); "CODE" 34015;
"PROCEDURE" ICHROWCOL(L, U, I, J, A); "CODE" 34033;
"PROCEDURE" ELMVECCOL(L, U, I, A, B, X); "CODE" 34021;
"INTEGER" "PROCEDURE" QRISNGVALDEC(A,M,N,VAL,V,EM); "CODE" 34273;
"PROCEDURE" SETRANDOM(X); "CODE" 11014;
"REAL" "PROCEDURE" RANDOM; "CODE" 11015;
"REAL" "PROCEDURE" DWARF; "CODE" 30003;
"PROCEDURE" SORT;
"BEGIN" "INTEGER" I, J, K; "REAL" S;
"FOR" I:= 1 "STEP" 1 "UNTIL" N - 1 "DO"
"BEGIN" K:= I; S:= D[I];
"FOR" J:= I+1 "STEP" 1 "UNTIL" N "DO" "IF" D[J]>S "THEN"
"BEGIN" K:= J; S:= D[J] "END";
"IF" K>I "THEN"
"BEGIN" D[K]:= D[I]; D[I]:= S;
"FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"
"BEGIN" S:=V[J,I]; V[J,I]:= V[J,K]; V[J,K]:= S
"END"
"END"
"END"
"END" SORT;
1SECTION : 5.1.2.2.2 � (OCTOBER 1975) PAGE 6
;
"PROCEDURE" MIN(J, NITS, D2, X1, F1, FK); "VALUE" J, NITS, FK;
"INTEGER" J, NITS; "REAL" D2, X1, F1; "BOOLEAN" FK;
"BEGIN"
"REAL" "PROCEDURE" FLIN(L); "VALUE" L; "REAL" L;
"BEGIN" "INTEGER" I; "ARRAY" T[1:N];
"IF" J > 0 "THEN"
"BEGIN" "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
T[I]:= X[I] + L * V[I,J]
"END" "ELSE"
"BEGIN" "COMMENT" SEARCH ALONG PARABOLIC SPACE CURVE;
QA:= L * (L - QD1) / (QD0 * (QD0 + QD1));
QB:= (L + QD0) * (QD1 - L) /(QD0 * QD1);
QC:= L * (L + QD0) / (QD1 * (QD0 + QD1));
"FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
T[I]:= QA * Q0[I] +QB * X[I] + QC * Q1[I]
"END";
NF:= NF + 1; FLIN:= FUNCT(N, T)
"END" FLIN;
"INTEGER" K; "BOOLEAN" DZ;
"REAL" X2, XM, F0, F2, FM, D1, T2, S, SF1, SX1;
SF1:= F1; SX1:= X1;
K:= 0; XM:= 0; F0:= FM:= FX; DZ:= D2 < RELTOL;
S:= SQRT(VECVEC(1,N,0,X,X));
T2:= M4 * SQRT(ABS(FX) / ("IF" DZ "THEN" DMIN "ELSE" D2)
+ S * LDT) + M2 * LDT; S:= S * M4 + ABSTOL;
"IF" DZ "AND" T2 > S "THEN" T2:= S;
"IF"T2<SMALL"THEN"T2:= SMALL;
"IF"T2>0.01*H "THEN"T2:= 0.01*H;
"IF"FK"AND"F1<=FM "THEN"
"BEGIN"XM:=X1; FM:= F1 "END";
"IF" ^ FK"OR"ABS(X1)<T2"THEN"
"BEGIN"X1:="IF"X1>0 "THEN"T2"ELSE"-T2;
F1:= FLIN(X1)
"END";
"IF"F1<= FM"THEN"
"BEGIN"XM:= X1; FM:= F1 "END";
L0: "IF" DZ "THEN"
"BEGIN" "COMMENT"EVALUATE FLIN AT ANOTHER POINT
AND ESTIMATE THE SECOND DERIVATIVE;
X2:= "IF" F0 < F1 "THEN" -X1 "ELSE" X1 * 2;
F2:= FLIN(X2); "IF"F2 <= FM "THEN"
"BEGIN" XM:= X2; FM:= F2 "END";
D2:=(X2*(F1-F0)-X1*(F2-F0))/(X1*X2*(X1-X2))
"END";
"COMMENT"ESTIMATE FIRST DERIVATIVE AT 0;
D1:=(F1-F0)/X1-X1*D2; DZ:="TRUE";
X2:= "IF"D2<=SMALL"THEN"
("IF"D1<0"THEN"H"ELSE"-H)
"ELSE"-0.5*D1/D2;
"IF"ABS(X2)>H"THEN"X2:="IF"X2>0"THEN"H"ELSE"-H;
"COMMENT"
1SECTION : 5.1.2.2.2 (DECEMBER 1979) PAGE 7
;
L1: F2:=FLIN(X2);
"IF"K<NITS"AND"F2>F0"THEN"
"BEGIN"K:=K+1;
"IF"F0<F1"AND"X1*X2>0"THEN" "GOTO"L0;
X2:= 0.5*X2; "GOTO"L1
"END";
NL:= NL+1;
"IF"F2>FM"THEN"X2:=XM"ELSE"FM:=F2;
D2:="IF"ABS(X2*(X2-X1))>SMALL"THEN"
(X2*(F1-F0)-X1*(FM-F0))/(X1*X2*(X1-X2))
"ELSE" "IF"K>0"THEN"0"ELSE"D2;
"IF"D2<=SMALL"THEN"D2:=SMALL;
X1:=X2; FX:=FM;
"IF"SF1<FX"THEN"
"BEGIN" FX:=SF1; X1:=SX1 "END";
"IF"J>0"THEN"ELMVECCOL(1,N,J,X,V,X1)
"END" MIN;
"PROCEDURE"QUAD;
"BEGIN" "INTEGER" I; "REAL" L, S;
S:= FX; FX:= QF1; QF1:= S; QD1:= 0;
"FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
"BEGIN"S:=X[I]; X[I]:= L:= Q1[I]; Q1[I]:= S;
QD1:= QD1 + (S - L) ** 2
"END";
L:=QD1:=SQRT(QD1); S:= 0;
"IF"(QD0*QD1>DWARF)"AND"NL>=3*N*N"THEN"
"BEGIN"MIN(0,2,S,L,QF1,"TRUE");
QA:= L*(L-QD1)/(QD0*(QD0+QD1));
QB:=(L+QD0)*(QD1-L)/(QD0*QD1);
QC:= L*(L+QD0)/(QD1*(QD0+QD1))
"END" "ELSE"
"BEGIN" FX:= QF1; QA:= QB:= 0; QC:= 1 "END";
QD0:= QD1;"FOR"I:= 1"STEP"1"UNTIL"N"DO"
"BEGIN"S:=Q0[I]; Q0[I]:=X[I];
X[I]:= QA*S + QB*X[I]+QC*Q1[I]
"END"
"END" QUAD;
"BOOLEAN" ILLC;
"INTEGER" I, J, K, K2, NL, MAXF, NF, KL, KT, KTM;
"REAL" S, SL, DN, DMIN, FX, F1, LDS, LDT, SF, DF, QF1, QD0,
QD1, QA, QB, QC, M2, M4, SMALL, VSMALL, LARGE, VLARGE, SCBD,
LDFAC,T2, MACHEPS, RELTOL, ABSTOL, H;
"ARRAY" V[1:N,1:N], D, Y, Z, Q0, Q1[1:N];
MACHEPS:= IN[0]; RELTOL:= IN[1]; ABSTOL:= IN[2]; MAXF:= IN[5];
H:= IN[6]; SCBD:= IN[7]; KTM:= IN[8]; ILLC:= IN[9] < 0;
SMALL:= MACHEPS ** 2; VSMALL:= SMALL ** 2;
LARGE:= 1/SMALL; VLARGE:= 1/VSMALL;
M2:= RELTOL; M4:= SQRT(M2); SETRANDOM(0.5);
LDFAC:= "IF" ILLC "THEN" 0.1 "ELSE" 0.01;
KT:=NL:=0; NF:=1; OUT[3]:= QF1:=FX:=FUNCT(N,X);
"COMMENT"
1SECTION : 5.1.2.2.2 (DECEMBER 1979) PAGE 8
;
ABSTOL:=T2:= SMALL+ABS(ABSTOL); DMIN:= SMALL;
"IF" H<ABSTOL*100"THEN"H:=ABSTOL*100; LDT:=H;
INIMAT(1,N,1,N,V,0);
"FOR"I:=1"STEP"1"UNTIL"N"DO"V[I,I]:= 1;
D[1]:= QD0:= 0; DUPVEC(1,N,0,Q1,X);
INIVEC(1,N,Q0,0);
"COMMENT"MAIN LOOP;
L0: SF:=D[1]; D[1]:= S:= 0;
MIN(1,2,D[1],S,FX,"FALSE");
"IF" S <= 0 "THEN" MULCOL(1, N, 1, 1, V, V, -1);
"IF" SF <= 0.9 * D[1] "OR" 0.9 * SF >= D[1] "THEN"
INIVEC(2,N,D,0);
"FOR" K:= 2"STEP"1"UNTIL"N"DO"
"BEGIN" DUPVEC(1,N,0,Y,X); SF:=FX;
ILLC:= ILLC "OR" KT>0;
L1: KL:=K; DF:= 0; "IF" ILLC "THEN"
"BEGIN" "COMMENT"RANDOM STOP TO GET OFF
RESULTION VALLEY;
"FOR"I:= 1 "STEP"1"UNTIL"N"DO"
"BEGIN"S:=Z[I]:=(0.1*LDT+T2*10**KT)
*(RANDOM-0.5);
ELMVECCOL(1,N,I,X,V,S)
"END";
FX:= FUNCT(N,X); NF:= NF+1
"END";
"FOR"K2:= K "STEP" 1 "UNTIL" N "DO"
"BEGIN" SL:=FX; S:= 0;
MIN (K2, 2, D[K2], S, FX, "FALSE");
S:="IF" ILLC "THEN" D[K2] * (S + Z[K2]) ** 2
"ELSE"SL-FX;"IF"DF<S"THEN"
"BEGIN"DF:=S;KL:= K2"END";
"END";
"IF" ^ILLC "AND" DF < ABS(100 * MACHEPS * FX) "THEN"
"BEGIN" ILLC:= "TRUE"; "GOTO" L1 "END";
"FOR" K2:= 1"STEP" 1"UNTIL"K-1"DO"
"BEGIN" S:= 0; MIN(K2, 2, D[K2], S, FX, "FALSE") "END";
F1:= FX; FX:= SF; LDS:= 0;
"FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
"BEGIN" SL:= X[I]; X[I]:= Y[I]; SL:= Y[I]:= SL - Y[I];
LDS:= LDS + SL * SL
"END"; LDS:= SQRT(LDS);
"IF" LDS > SMALL "THEN"
"BEGIN" "FOR" I:= KL - 1 "STEP" -1 "UNTIL" K "DO"
"BEGIN" "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO"
V[J, I + 1]:= V[J,I]; D[I + 1]:= D[I]
"END";
D[K]:= 0; DUPCOLVEC(1, N, K, V, Y);
MULCOL(1, N, K, K, V, V, 1 / LDS);
MIN(K, 4, D[K], LDS, F1, "TRUE"); "IF" LDS <= 0 "THEN"
"BEGIN" LDS:= LDS; MULCOL(1, N, K, K, V, V, -1) "END"
"END";
LDT:= LDFAC * LDT; "IF" LDT < LDS "THEN" LDT:= LDS;
T2:= M2 * SQRT(VECVEC(1, N, 0, X, X)) + ABSTOL;
"COMMENT"
1SECTION : 5.1.2.2.2 � (OCTOBER 1975) PAGE 9
;
KT:= "IF" LDT > 0.5 * T2 "THEN" 0 "ELSE" KT + 1;
"IF" KT > KTM "THEN" "BEGIN" OUT[1]:= 0; "GOTO" L2 "END"
"END";
QUAD;
DN:= 0;"FOR"I:= 1"STEP"1"UNTIL"N"DO"
"BEGIN"D[I]:= 1/SQRT(D[I]);
"IF"DN<D[I]"THEN"DN:=D[I]
"END";
"FOR"J:= 1"STEP"1"UNTIL"N"DO"
"BEGIN"S:= D[J]/DN; MULCOL(1,N,J,J,V,V,S)"END";
"IF"SCBD>1"THEN"
"BEGIN"S:=VLARGE; "FOR"I:=1 "STEP"1"UNTIL"N"DO"
"BEGIN" SL:= Z[I]:= SQRT(MATTAM(1, N, I, I, V, V));
"IF"SL<M4"THEN"Z[I]:= M4;
"IF" S>SL "THEN" S:= SL
"END";
"FOR"I:=1"STEP"1"UNTIL"N"DO"
"BEGIN"SL:=S/Z[I];Z[I]:= 1/SL;
"IF"Z[I]>SCBD"THEN"
"BEGIN"SL:=1/SCBD; Z[I]:= SCBD"END";
MULROW(1, N, I, I, V, V, SL)
"END"
"END";
"FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
ICHROWCOL(I + 1, N, I, I, V);
"BEGIN" "ARRAY" A[1:N,1:N], EM[0:7];
EM[0]:= EM[2]:= MACHEPS;
EM[4]:= 10 * N; EM[6]:= VSMALL;
DUPMAT(1, N, 1, N, A, V);
"IF" QRISNGVALDEC(A, N, N, D, V, EM) ^= 0 "THEN"
"BEGIN" OUT[1]:= 2; "GOTO" L2 "END";
"END";
"IF"SCBD>1"THEN"
"BEGIN" "FOR"I:=1"STEP"1"UNTIL"N"DO"
MULROW(1,N,I,I,V,V,Z[I]);
"FOR"I:= 1"STEP"1"UNTIL"N"DO"
"BEGIN"S:= SQRT(TAMMAT(1,N,I,I,V,V));
D[I]:= S*D[I]; S:= 1/S;
MULCOL(1,N,I,I,V,V,S)
"END"
"END";
"FOR" I:= 1 "STEP" 1 "UNTIL" N "DO"
"BEGIN" S:= DN * D[I];
D[I]:= "IF" S > LARGE "THEN" VSMALL "ELSE"
"IF" S < SMALL "THEN" VLARGE "ELSE" S ** (-2)
"END";
SORT;
DMIN:= D[N]; "IF" DMIN < SMALL "THEN" DMIN:= SMALL;
ILLC:= (M2 * D[1]) > DMIN;
"IF" NF < MAXF "THEN" "GOTO" L0 "ELSE" OUT[1]:= 1;
L2: OUT[2]:= FX;
OUT[4]:= NF; OUT[5]:= NL; OUT[6]:= LDT
"END"PRAXIS;
"EOP"