1SECTION : 1.1.8 (OCTOBER 1975) PAGE 1 AUTHORS: C.G. VAN DER LAAN AND J.C.P. BUS. CONTRIBUTOR: J.C.P. BUS. INSTITUTE: MATHEMATICAL CENTRE. RECEIVED: 740921. BRIEF DESCRIPTION: INFNRMVEC CALCULATES THE INFINITY-NORM OF A VECTOR; INFNRMROW CALCULATES THE INFINITY-NORM OF A ROW VECTOR; INFNRMCOL CALCULATES THE INFINITY-NORM OF A COLUMN VECTOR; INFNRMMAT CALCULATES THE INFINITY-NORM OF A MATRIX; ONENRMVEC CALCULATES THE ONE-NORM OF A VECTOR; ONENRMROW CALCULATES THE ONE-NORM OF A ROW VECTOR; ONENRMCOL CALCULATES THE ONE-NORM OF A COLUMN VECTOR; ONENRMMAT CALCULATES THE ONE-NORM OF A MATRIX; ABSMAXMAT CALCULATES FOR A GIVEN MATRIX THE MODULUS OF AN ELEMENT WHICH IS OF MAXIMUM ABSOLUTE VALUE; KEYWORDS: VECTOR NORMS, MATRIX NORMS. 1SECTION : 1.1.8 (DECEMBER 1979) PAGE 2 SUBSECTION: INFNRMVEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" INFNRMVEC(L, U, K, A); "VALUE" L, U; "INTEGER" L, U, K; "ARRAY" A; INFNRMVEC := MAX( ABS(A[I]), I= L, ..., U ); THE MEANING OF THE FORMAL PARAMETERS IS: L, U: ; ENTRY:THE LOWER BOUND AND UPPER BOUND OF THE INDEX OF THE VECTOR A, RESPECTIVELY; K: ; EXIT:THE FIRST INDEX FOR WHICH ABS(A[I]), I = L, ..., U, IS MAXIMAL; A: ; "ARRAY" A[L:U]. PROCEDURES USED: NONE. SUBSECTION: INFNRMROW. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" INFNRMROW(L, U, I, K, A); "VALUE" L, U, I; "INTEGER" L, U, I, K; "ARRAY" A; INFNRMROW := MAX( ABS(A[I,J]), J= L, ..., U ); THE MEANING OF THE FORMAL PARAMETERS IS: L, U: ; ENTRY:THE LOWER BOUND AND UPPER BOUND OF THE COLUMN INDEX OF THE ROW VECTOR A, RESPECTIVELY; I: ; ENTRY:THE ROW INDEX; K: ; EXIT:THE FIRST INDEX FOR WHICH ABS(A[I,J]), J = L, ..., U, IS MAXIMAL; A: ; "ARRAY" A[I:I,L:U]. PROCEDURES USED: NONE. 1SECTION : 1.1.8 (DECEMBER 1979) PAGE 3 SUBSECTION: INFNRMCOL. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" INFNRMCOL(L, U, J, K, A); "VALUE" L, U, J; "INTEGER" L, U, J, K; "ARRAY" A; INFNRMCOL := MAX( ABS(A[I,J]), I= L, ..., U ); THE MEANING OF THE FORMAL PARAMETERS IS: L, U: ; ENTRY:THE LOWER BOUND AND UPPER BOUND OF THE ROW INDEX OF THE COLUMN VECTOR A, RESPECTIVELY; J: ; ENTRY:THE COLUMN INDEX; K: ; EXIT:THE FIRST INDEX FOR WHICH ABS(A[I,J]), I = L, ..., U, IS MAXIMAL; A: ; "ARRAY" A[L:U,J:J]. PROCEDURES USED: NONE. SUBSECTION: INFNRMMAT. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" INFNRMMAT(LR, UR, LC, UC, KR, A); "VALUE" LR, UR, LC, UC; "INTEGER" LR, UR, LC, UC, KR; "ARRAY" A; INFNRMMAT := MAX( ONENRMROW(LC, UC, I, A), I=LR, ..., UR ); THE MEANING OF THE FORMAL PARAMETERS IS: LR, UR: ; ENTRY:THE LOWER BOUND AND UPPER BOUND OF THE ROW INDEX, RESPECTIVELY; LC, UC: ; ENTRY:THE LOWER BOUND AND UPPER BOUND OF THE COLUMN INDEX, RESPECTIVELY; KR: ; EXIT:THE FIRST ROW INDEX FOR WHICH THE ONE-NORM IS MAXIMAL; A: ; "ARRAY" A[LR:UR,LC:UC]. PROCEDURES USED: ONENRMROW. 1SECTION : 1.1.8 (DECEMBER 1979) PAGE 4 SUBSECTION: ONENRMVEC. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" ONENRMVEC(L, U, A); "VALUE" L, U; "INTEGER" L, U; "ARRAY" A; ONENRMVEC := SUM( ABS(A[I]), I= L, ..., U ); THE MEANING OF THE FORMAL PARAMETERS IS: L, U: ; ENTRY:THE LOWER BOUND AND UPPER BOUND OF THE INDEX OF THE VECTOR A, RESPECTIVELY; A: ; "ARRAY" A[L:U]. PROCEDURES USED: NONE. SUBSECTION: ONENRMROW. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" ONENRMROW(L, U, I, A); "VALUE" L, U, I; "INTEGER" L, U, I; "ARRAY" A; ONENRMROW := SUM( ABS(A[I,J]), J= L, ..., U ); THE MEANING OF THE FORMAL PARAMETERS IS: L, U: ; ENTRY:THE LOWER BOUND AND UPPER BOUND OF THE COLUMN INDEX OF THE ROW VECTOR A, RESPECTIVELY; I: ; ENTRY: THE ROW INDEX; A: ; "ARRAY" A[I:I,L:U]. PROCEDURES USED: NONE. 1SECTION : 1.1.8 (DECEMBER 1979) PAGE 5 SUBSECTION: ONENRMCOL. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" ONENRMCOL(L, U, J, A); "VALUE" L, U, J; "INTEGER" L, U, J; "ARRAY" A; ONENRMCOL := SUM( ABS(A[I,J]), I= L, ..., U ); THE MEANING OF THE FORMAL PARAMETERS IS: L, U: ; ENTRY:THE LOWER BOUND AND UPPER BOUND OF THE ROW INDEX OF THE COLUMN VECTOR A, RESPECTIVELY; J: ; ENTRY: THE COLUMN INDEX; A: ; "ARRAY" A[L:U,J:J]. PROCEDURES USED: NONE. SUBSECTION: ONENRMMAT. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" ONENRMMAT(LR, UR, LC, UC, KC, A); "VALUE" LR, UR, LC, UC; "INTEGER" LR, UR, LC, UC, KC; "ARRAY" A; ONENRMMAT := MAX(ONENRMCOL(LR, UR, J, A), J=LC, ..., UC); THE MEANING OF THE FORMAL PARAMETERS IS: LR, UR: ; ENTRY:THE LOWER BOUND AND UPPER BOUND OF THE ROW INDEX, RESPECTIVELY; LC, UC: ; ENTRY:THE LOWER BOUND AND UPPER BOUND OF THE COLUMN INDEX, RESPECTIVELY; KC: ; EXIT:THE FIRST COLUMN INDEX FOR WHICH THE ONE-NORM IS MAXIMAL; A: ; "ARRAY" A[LR:UR,LC:UC]. PROCEDURES USED: ONENRMCOL. 1SECTION : 1.1.8 (OCTOBER 1975) PAGE 6 SUBSECTION: ABSMAXMAT. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" ABSMAXMAT(LR, UR, LC, UC, KR, KC, A); "VALUE" LR, UR, LC, UC; "INTEGER" LR, UR, LC, UC, KR, KC; "ARRAY" A; ABSMAXMAT := MAX( ABS(A[I,J]), I= LR, ..., UR, J= LC, ..., UC ); THE MEANING OF THE FORMAL PARAMETERS IS: LR, UR: ; ENTRY:THE LOWER BOUND AND UPPER BOUND OF THE ROW INDEX, RESPECTIVELY; LC, UC: ; ENTRY:THE LOWER BOUND AND UPPER BOUND OF THE COLUMN INDEX, RESPECTIVELY; KR, KC: ; EXIT:THE ROW AND COLUMN INDEX OF AN ELEMENT FOR WHICH THE MODULUS IS MAXIMAL; A: ; "ARRAY" A[LR:UR,LC:UC]. PROCEDURES USED: INFNRMCOL. LANGUAGE: COMPASS. METHOD AND PERFORMANCE: ELEMENTARY. 1SECTION : 1.1.8 (OCTOBER 1975) PAGE 7 SOURCE TEXT(S): THE FOLLOWING PROCEDURES ARE WRITTEN IN COMPASS, AN EQUIVALENT ALGOL 60 TEXT OF THESE COMPASS ROUTINES IS GIVEN. 0"CODE" 31061; "REAL" "PROCEDURE" INFNRMVEC(L, U, K, A); "VALUE" L, U; "INTEGER" L, U, K; "ARRAY" A; "BEGIN" "REAL" R, MAX; MAX:= 0; K:= L; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" "BEGIN" R:= ABS(A[L]); "IF" R > MAX "THEN" "BEGIN" MAX:= R; K:= L "END" "END"; INFNRMVEC:= MAX "END" INFNRMVEC; "EOP" 0"CODE" 31062; "REAL" "PROCEDURE" INFNRMROW(L, U, I, K, A); "VALUE" L, U, I; "INTEGER" L, U, I, K; "ARRAY" A; "BEGIN" "REAL" R, MAX; MAX:= 0; K:= L; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" "BEGIN" R:= ABS(A[I,L]); "IF" R > MAX "THEN" "BEGIN" MAX:= R; K:= L "END" "END"; INFNRMROW:= MAX "END" INFNRMROW; "EOP" 0"CODE" 31063; "REAL" "PROCEDURE" INFNRMCOL(L, U, J, K, A); "VALUE" L, U, J; "INTEGER" L, U, J, K; "ARRAY" A; "BEGIN" "REAL" R, MAX; MAX:= 0; K:= L; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" "BEGIN" R:= ABS(A[L,J]); "IF" R > MAX "THEN" "BEGIN" MAX:= R; K:= L "END" "END"; INFNRMCOL:= MAX "END" INFNRMCOL; "EOP" 1SECTION : 1.1.8 (OCTOBER 1975) PAGE 8 0"CODE" 31064; "REAL" "PROCEDURE" INFNRMMAT(LR, UR, LC, UC, KR, A); "VALUE" LR, UR, LC, UC; "INTEGER" LR, UR, LC, UC, KR; "ARRAY" A; "BEGIN" "REAL" R, MAX; "REAL" "PROCEDURE" ONENRMROW(L, U, I, A); "CODE" 31066; MAX:= 0; KR:= LR; "FOR" LR:= LR "STEP" 1 "UNTIL" UR "DO" "BEGIN" R:= ONENRMROW(LC, UC, LR, A); "IF" R > MAX "THEN" "BEGIN" MAX:= R; KR:= LR "END" "END"; INFNRMMAT:= MAX "END" INFNRMMAT; "EOP" 0"CODE" 31065; "REAL" "PROCEDURE" ONENRMVEC(L, U, A); "VALUE" L, U; "INTEGER" L, U; "ARRAY" A; "BEGIN" "REAL" SUM; SUM:= 0; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" SUM:= SUM + ABS(A[L]); ONENRMVEC:= SUM "END" ONENRMVEC; "EOP" 0"CODE" 31066; "REAL" "PROCEDURE" ONENRMROW(L, U, I, A); "VALUE" L, U, I; "INTEGER" L, U, I; "ARRAY" A; "BEGIN" "REAL" SUM; SUM:= 0; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" SUM:= SUM + ABS(A[I,L]); ONENRMROW:= SUM "END" ONENRMROW; "EOP" 1SECTION : 1.1.8 (OCTOBER 1975) PAGE 9 0"CODE" 31067; "REAL" "PROCEDURE" ONENRMCOL(L, U, J, A); "VALUE" L, U, J; "INTEGER" L, U, J; "ARRAY" A; "BEGIN" "REAL" SUM; SUM:= 0; "FOR" L:= L "STEP" 1 "UNTIL" U "DO" SUM:= SUM + ABS(A[L,J]); ONENRMCOL:= SUM "END" ONENRMCOL; "EOP" 0"CODE" 31068; "REAL" "PROCEDURE" ONENRMMAT(LR, UR, LC, UC, KC, A); "VALUE" LR, UR, LC, UC; "INTEGER" LR, UR, LC, UC, KC; "ARRAY" A; "BEGIN" "REAL" MAX, R; "REAL" "PROCEDURE" ONENRMCOL(L, U, J, A); "CODE" 31067; MAX:= 0; KC:= LC; "FOR" LC:= LC "STEP" 1 "UNTIL" UC "DO" "BEGIN" R:= ONENRMCOL(LR, UR, LC, A); "IF" R > MAX "THEN" "BEGIN" MAX:= R; KC:= LC "END" "END"; ONENRMMAT:= MAX "END" ONENRMMAT; "EOP" 0"CODE" 31069; "REAL" "PROCEDURE" ABSMAXMAT(LR, UR, LC, UC, I, J, A); "VALUE" LR, UR, LC, UC; "INTEGER" LR, UR, LC, UC, I, J; "ARRAY" A; "BEGIN" "INTEGER" II; "REAL" MAX, R; "REAL" "PROCEDURE" INFNRMCOL(L, U, I, K, A); "CODE" 31063; MAX:= 0; I:= LR; J:= LC; "FOR" LC:= LC "STEP" 1 "UNTIL" UC "DO" "BEGIN" R:= INFNRMCOL(LR, UR, LC, II, A); "IF" R > MAX "THEN" "BEGIN" MAX:= R; I:= II; J:= LC "END" "END"; ABSMAXMAT:= MAX "END" ABSMAXMAT; "EOP" 1SECTION : 6.10.4 (OCTOBER 1975) PAGE 1 AUTHOR : P.W.HEMKER. CONTRIBUTOR : F.GROEN. INSTITUTE : MATHEMATICAL CENTRE. RECEIVED : 740620. BRIEF DESCRIPTION : THIS SECTION CONTAINS TWO PROCEDURES FOR THE EVALUATION OF AIRY FUNCTIONS AND COMPUTING THEIR ZEROS. FOR THE DEFINITION OF THESE FUNCTIONS SEE REF[1]. AIRY EVALUATES THE AIRY FUNCTIONS AI(Z) AND BI(Z) AND THEIR DERIVATIVES. AIRYZEROS COMPUTES THE ZEROS AND ASSOCIATED VALUES OF THE AIRY FUNCTIONS AI(Z) AND BI(Z) AND THEIR DERIVATIVES. KEYWORDS : AIRY FUNCTION, DERIVATIVE AIRY FUNCTION, ZERO OF AIRY FUNCTION. 1SECTION : 6.10.4 (OCTOBER 1975) PAGE 2 SUBSECTION : AIRY. CALLING SEQUENCE : THE HEADING OF THE PROCEDURE READS : "PROCEDURE" AIRY(X,AI,AID,BI,BID,EXPON,FIRST); "VALUE" X,FIRST; "BOOLEAN" FIRST; "REAL" X,AI,AID,BI,BID,EXPON; THE MEANING OF THE FORMAL PARAMETERS IS : X: ; ENTRY : THE REAL ARGUMENT OF THE AIRY FUNCTIONS. AI: ; EXIT : THE VALUE OF THE AIRY FUNCTION AI IS GIVEN BY : EXP( -EXPON ) * AI. NOTE : IF X < 9 THEN EXPON = 0. AID: ; EXIT : THE VALUE OF THE DERIVATIVE OF THE AIRY FUNCTION AI IS GIVEN BY : EXP( -EXPON ) * AID. NOTE : IF X < 9 THEN EXPON = 0. BI: ; EXIT : THE VALUE OF THE AIRY FUNCTION BI IS GIVEN BY : EXP( EXPON ) * BI. NOTE : IF X < 9 THEN EXPON = 0. BID: ; EXIT : THE VALUE OF THE DERIVATIVE OF THE AIRY FUNCTION BI IS GIVEN BY : EXP( EXPON ) * BID. NOTE : IF X < 9 THEN EXPON = 0. EXPON: ; EXIT : IF X < 9 THEN 0 ELSE 2/3 * X ** (3/2). FIRST: ; FIRST SHOULD BE "FALSE" UNLESS THE PROCEDURE IS CALLED FOR THE FIRST TIME. IF FIRST IS "TRUE" THEN TWO OWN ARRAYS OF COEFFICIENTS ARE BUILT UP. PROCEDURES USED : NONE. 1SECTION : 6.10.4 (OCTOBER 1975) PAGE 3 REQUIRED CENTRAL MEMORY : TWO OWN ARRAYS OF ORDER 10 ARE DECLARED. RUNNING TIME : IF 2.5 <= X <= 8 THEN ABOUT 8"-3 SEC., ELSE BETWEEN 3"-3 AND 4"-3 SEC. ON THE CYBER 73/28. LANGUAGE : ALGOL 60. METHOD AND PERFORMANCE : SEE REF[2] OF THE SUBSECTION AIRYZEROS (THIS SECTION). REFERENCES : SEE REFERENCES OF THE SUBSECTION AIRYZEROS (THIS SECTION). EXAMPLE OF USE : "BEGIN" "REAL" A,B,C,D,E; "PROCEDURE" AIRY(X,AI,AID,BI,BID,EXPON,FIRST);"CODE" 35140; AIRY(9.654894,A,B,C,D,E,"TRUE"); OUTPUT(61,"("/,"("AI (9.654894) = ")",N")",A*EXP(-E)); OUTPUT(61,"("/,"("AID(9.654894) = ")",N")",B*EXP(-E)); OUTPUT(61,"("/,"("BI (9.654894) = ")",N")",C*EXP( E)); OUTPUT(61,"("/,"("BID(9.654894) = ")",N")",D*EXP( E)); "END" RESULTS : AI (9.654894) = +3.2873525549165"-010 AID(9.654894) = -1.0297999323482"-009 BI (9.654894) = +1.5583887049670"+008 BID(9.654894) = +4.8010374682654"+008 1SECTION : 6.10.4 (OCTOBER 1975) PAGE 4 SUBSECTION : AIRYZEROS. CALLING SEQUENCE : THE HEADING OF THE PROCEDURE READS : "REAL" "PROCEDURE" AIRYZEROS(N,D,ZAI,VAI); "VALUE" N,D; "INTEGER" N,D; "ARRAY" ZAI,VAI; AIRYZEROS := THE N-TH ZERO OF THE SELECTED AIRY-FUNCTION. THE MEANING OF THE FORMAL PARAMETERS IS : N : ; ENTRY : THE NUMBER OF ZEROS TO BE CALCULATED; D : ; ENTRY : AN INTEGER WHICH SELECTS THE REQUIRED AIRY FUNCTION. D = 0, 1, 2 OR 3. ZAI : ; "ARRAY" ZAI[1 : N]; EXIT : ZAI[J] CONTAINS THE J-TH ZERO OF THE SELECTED AIRY-FUNCTION : IF D = 0 THEN AI(Z), IF D = 1 THEN (D/DX) AI(X), IF D = 2 THEN BI(X), IF D = 3 THEN (D/DX) BI(X); VAI : ; "ARRAY" VAI[1 : N]; EXIT: VAI[J] CONTAINS THE VALUE AT X = ZAI[J] OF THE FOLLOWING FUNCTION : IF D = 0 THEN (D/DX) AI(X), IF D = 1 THEN AI(X), IF D = 2 THEN (D/DX) BI(X), IF D = 3 THEN BI(X); PROCEDURES USED : AIRY = CP35140; REQUIRED CENTRAL MEMORY : NO AUXILIARY ARRAYS ARE DECLARED. RUNNING TIME : DEPENDENT ON THE VALUES OF N AND D. IN MOST CASES THE RUNNING TIME IS LESS THAN N * 0.01 SEC. ON THE CYBER 73/28. LANGUAGE : ALGOL 60. 1SECTION : 6.10.4 (OCTOBER 1975) PAGE 5 METHOD AND PERFORMANCE : A FIRST APPROXIMATION OF THE ZEROS OF THE SELECTED AIRY-FUNCTION IS CALCULATED BY MEANS OF THE ASYMPTOTIC EXPANSION ( SEE THE FORMULAS 10.4.94 - 10.4.105 OF REF[1] ); THIS VALUE IS CORRECTED BY THE (REPEATED) USE OF A QUADRATIC INTERPOLATION RULE. THE COMPUTED ZEROS WILL SATISFY AT LEAST ONE OF THE FOLLOWING CONDITIONS : 1: THE ABSOLUTE VALUE OF THE SELECTED AIRY-FUNCTION AT A COMPUTED ZERO IS LESS THAN "-12. NOTE: THE VALUES OF THE AIRY-FUNCTIONS ARE CALCULATED BY MEANS OF THE PROCEDURE AIRY (THIS SECTION). 2: THE RELATIVE PRECISION OF THE COMPUTED ZERO IS "-14. THE ASSOCIATED VALUES ( DELIVERED IN THE ARRAY VAI ) ARE ALSO CALCULATED BY MEANS OF THE PROCEDURE AIRY (THIS SECTION). REFERENCES : [1] : M.ABRAMOWITZ AND I.A.STEGUN, HANDBOOK OF MATHMATICAL FUNCTIONS, DOVER PUBLICATIONS, INC. NEW YORK, 1965. [2] : R.G.GORDON, EVALUATION OF AIRY FUNCTIONS, THE JOURNAL OF CHEMICAL PHYSICS, VOLUME 51, 1969, PP. 23-24. EXAMPLE OF USE : "BEGIN" "ARRAY" ZBI,VBID[1 : 3]; "REAL" "PROCEDURE" AIRYZEROS(N,D,ZAI,VAI); "CODE"35145; OUTPUT(61,"("/"("THE THIRD ZERO OF BI(X) IS")"/,N, /"("THE VALUE OF (D/DX)BI(X) IN THIS POINT IS")"/,N")" ,AIRYZEROS(3,2,ZBI,VBID),VBID[3]) "END" RESULTS : THE THIRD ZERO OF BI(X) IS -4.8307378416626"+000 THE VALUE OF (D/DX)BI(X) IN THIS POINT IS +8.3699101261986"-001 1SECTION : 6.10.4 (OCTOBER 1975) PAGE 6 SOURCE TEXT(S): 0"CODE" 35140; "PROCEDURE" AIRY(Z,AI,AID,BI,BID,EXPON,FIRST); "VALUE" Z,FIRST; "BOOLEAN" FIRST; "REAL" Z,AI,AID,BI,BID,EXPON; "BEGIN" "REAL" S,T,U,V,SC,TC,UC,VC,X,K1,K2,K3,K4, C,ZT,SI,CO,EXPZT,SQRTZ,WWL,PL,PL1,PL2,PL3; "OWN" "REAL" C1,C2,SQRT3,SQRT1OPI,PIO4; "OWN" "REAL" "ARRAY" XX,WW[1:10]; "INTEGER" N,L; "IF" FIRST "THEN" "BEGIN" SQRT3:= 1.73205080756887729; SQRT1OPI:= 0.56418958354775629; PIO4:= 0.78539816339744831; C1:= 0.35502 80538 87817; C2:= 0.25881 94037 92807; XX[ 1]:= 1.40830 81072 180964 "+1; XX[ 2]:= 1.02148 85479 197331 "+1; XX[ 3]:= 7.44160 18450 450930 ; XX[ 4]:= 5.30709 43061 781927 ; XX[ 5]:= 3.63401 35029 132462 ; XX[ 6]:= 2.33106 52303 052450 ; XX[ 7]:= 1.34479 70824 609268 ; XX[ 8]:= 6.41888 58369 567296 "-1; XX[ 9]:= 2.01003 45998 121046 "-1; XX[10]:= 8.05943 59172 052833 "-3; WW[ 1]:= 3.15425 15762 964787"-14; WW[ 2]:= 6.63942 10819 584921"-11; WW[ 3]:= 1.75838 89061 345669"- 8; WW[ 4]:= 1.37123 92370 435815"- 6; WW[ 5]:= 4.43509 66639 284350"- 5; WW[ 6]:= 7.15550 10917 718255"- 4; WW[ 7]:= 6.48895 66103 335381"- 3; WW[ 8]:= 3.64404 15875 773282"- 2; WW[ 9]:= 1.43997 92418 590999"- 1; WW[10]:= 8.12311 41336 261486"- 1; "END"; EXPON:= 0; "IF" Z >= -5.0 "AND" Z <= 8 "THEN" "BEGIN" U:= V:= T:= UC:= VC:= TC:= 1; S:= SC:= 0.5; N:= 0; X:= Z*Z*Z; "FOR" N:= N+3 "WHILE" ABS(U)+ABS(V)+ABS(S)+ABS(T) > "-18 "DO" "BEGIN" U:=U*X/(N*(N-1)); V:= V*X/(N*(N+1)); S:=S*X/(N*(N+2)); T:= T*X/(N*(N-2)); UC:= UC+U; VC:= VC+V; SC:= SC+S; TC:= TC+T "END"; "COMMENT" 1SECTION : 6.10.4 (OCTOBER 1975) PAGE 7 ; BI:= SQRT3 * (C1*UC + C2*Z*VC); BID:=SQRT3 * (C1*Z*Z*SC +C2*TC); "IF" Z<2.5 "THEN" "BEGIN" AI:= C1*UC - C2*Z*VC; AID:= C1*SC*Z*Z - C2*TC; "GOTO" END "END" "END"; K1:= K2:= K3:= K4:= 0; SQRTZ:= SQRT(ABS(Z)); ZT:= 0.66666 66666 66667 * ABS(Z)*SQRTZ; C:= SQRT1OPI/SQRT(SQRTZ); "IF" Z<0 "THEN" "BEGIN" Z:= -Z; CO:= COS(ZT-PIO4); SI:= SIN(ZT-PIO4); "FOR" L:= 1 "STEP" 1 "UNTIL" 10 "DO" "BEGIN" WWL:= WW[L]; PL:= XX[L]/ZT; PL2:=PL*PL; PL1:= 1+PL2; PL3:= PL1*PL1; K1:= K1 + WWL/PL1; K2:= K2 + WWL*PL/PL1; K3:= K3 + WWL*PL*(1+PL*(2/ZT+PL))/PL3; K4:= K4 + WWL*(-1-PL*(1+PL*(ZT-PL))/ZT)/PL3; "END"; AI:= C*(CO*K1+SI*K2); AID:= 0.25*AI/Z - C*SQRTZ*(CO*K3+SI*K4); BI:= C*(CO*K2-SI*K1); BID:= 0.25*BI/Z - C*SQRTZ*(CO*K4-SI*K3); "END" "ELSE" "BEGIN" "IF" Z < 9 "THEN" EXPZT:= EXP(ZT) "ELSE" "BEGIN" EXPZT:= 1; EXPON:= ZT "END"; "FOR" L:= 1 "STEP" 1 "UNTIL" 10 "DO" "BEGIN" WWL:= WW[L]; PL:= XX[L]/ZT; PL1:= 1+PL; PL2:= 1-PL; K1:= K1 + WWL/PL1; K2:= K2 + WWL*PL/(ZT*PL1*PL1); K3:= K3 + WWL/PL2; K4:= K4 + WWL*PL/(ZT*PL2*PL2); "END"; AI:= 0.5*C*K1/EXPZT; AID:= AI*(-.25/Z-SQRTZ) + 0.5*C*SQRTZ*K2/EXPZT; "IF" Z >= 8 "THEN" "BEGIN" BI:= C*K3*EXPZT; BID:= BI*(SQRTZ-0.25/Z) - C*K4*SQRTZ*EXPZT; "END"; "END"; END: "END" AIRY; "EOP" 1SECTION : 6.10.4 (OCTOBER 1975) PAGE 8 0"CODE" 35145; "REAL" "PROCEDURE" AIRYZEROS(N,D,ZAI,VAI); "VALUE" N,D; "INTEGER" N,D; "ARRAY" ZAI,VAI; "BEGIN" "BOOLEAN" A, FOUND; "INTEGER" I; "REAL" C,E,R,ZAJ,ZAK,VAJ,DAJ,KAJ,ZZ; "PROCEDURE" AIRY(A,B,C,D,E,F,G); "CODE" 35140; A := D = 0 "OR" D = 2; R := "IF" D = 0 "OR" D = 3 "THEN" -1.1780 97245 09617 "ELSE" -3.5342 91735 28852; "COMMENT" R := "IF" D = 0 "OR" D = 3 "THEN" -3 * PI / 8 "ELSE" -9 * PI / 8; AIRY(0,ZAJ,VAJ,DAJ,KAJ,ZZ,"TRUE"); "FOR" I := 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" R := R + 4.7123 88980 38469; "COMMENT" R := R + 3 * PI / 2; ZZ := R * R; ZAJ := "IF" I = 1 "AND" D = 1 "THEN" -1.01879 297 "ELSE" "IF" I = 1 "AND" D = 2 "THEN" -1.17371 322 "ELSE" R ** 0.66666 66666 66667 * ( "IF" A "THEN" - ( 1 + ( 5/48 - ( 5/36 - ( 77125/82944 - ( 1080 56875 / 69 67296 - (16 23755 96875 / 3344 30208) /ZZ)/ZZ)/ZZ)/ZZ)/ZZ) "ELSE" - ( 1 - ( 7/48 - ( 35/288 - ( 1 81223 / 2 07360 - ( 186 83371 / 12 44160 - ( 9 11458 84361 / 1911 02976 ) /ZZ)/ZZ)/ZZ)/ZZ)/ZZ)); "IF" D <= 1 "THEN" AIRY(ZAJ,VAJ,DAJ,C,E,ZZ,"FALSE") "ELSE" AIRY(ZAJ,C,E,VAJ,DAJ,ZZ,"FALSE"); FOUND := ABS( "IF" A "THEN" VAJ "ELSE" DAJ ) < "-12; "FOR" C := C "WHILE" "NOT" FOUND "DO" "BEGIN" "IF" A "THEN" "BEGIN" KAJ := VAJ / DAJ; ZAK := ZAJ - KAJ * (1 + ZAJ * KAJ * KAJ) "END" "ELSE" "BEGIN" KAJ := DAJ / (ZAJ * VAJ); ZAK := ZAJ - KAJ * (1 + KAJ * (KAJ * ZAJ + 1 / ZAJ)) "END"; "IF" D <= 1 "THEN" AIRY(ZAK,VAJ,DAJ,C,E,ZZ,"FALSE") "ELSE" AIRY(ZAK,C,E,VAJ,DAJ,ZZ,"FALSE"); FOUND := ABS(ZAK - ZAJ) < "-14 * ABS(ZAK) "OR" ABS("IF" A "THEN" VAJ "ELSE" DAJ) < "-12; ZAJ := ZAK "END"; VAI[I] := "IF" A "THEN" DAJ "ELSE" VAJ; ZAI[I] := ZAJ; "END"; AIRYZEROS := ZAI[N]; "END" AIRYZEROS; "EOP" 1SECTION : 2.2.1.1 (OCTOBER 1975) PAGE 1 AUTHOR: C.G. VAN DER LAAN. INSTITUTE: REKENCENTRUM RIJKSUNIVERSITEIT GRONINGEN. RECEIVED: 741114. BRIEF DESCRIPTION: THIS SECTION CONTAINS THE PROCEDURES POL, TAYPOL, NORDERPOL AND DERPOL. POL EVALUATES A POLYNOMIAL; DERPOL EVALUATES THE FIRST K DERIVATIVES OF A POLYNOMIAL; NORDERPOL EVALUATES THE FIRST K NORMALIZED DERIVATIVES OF A POLYNOMIAL (I.E. J-TH DERIVATIVE/(J FACTORIAL),J=0,1,...,K<=DEGREE; TAYPOL EVALUATES X**J*(J-TH DERIVATIVE)/(J FACTORIAL), J=0,1,...,K<=DEGREE. KEYWORDS: POLYNOMIAL EVALUATION, (NORMALIZED) DERIVATIVES, DERIVATIVES OF A POLYNOMIAL. SUBSECTION: POL. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL""PROCEDURE" POL(N,X,A); "VALUE"N,X;"INTEGER"N;"REAL"X;"ARRAY"A; POL:=A[0]+A[1]*X+...+A[N]*X**N; THE MEANING OF THE FORMAL PARAMETERS IS: N: ; ENTRY: THE DEGREE OF THE POLYNOMIAL; X: ; ENTRY: THE ARGUMENT OF THE POLYNOMIAL; A: ; "ARRAY"A[0:N]; ENTRY: THE COEFFICIENTS OF THE POLYNOMIAL A[0]+A[1]*X+...+A[N]*X**N. PROCEDURES USED: NONE. 1SECTION : 2.2.1.1 (OCTOBER 1975) PAGE 2 RUNNING TIME: PROPORTIONAL TO N (N MULTIPLICATIONS AND ADDITIONS). LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: THE METHOD USED FOR EVALUATION IS HORNER'S RULE (SYNTHETIC DIVISION). THE ERROR GROWTH IS GIVEN BY A LINEAR FUNCTION OF THE DEGREE OF THE POLYNOMIAL (SEE VAN DER LAAN, STOER(1972) P. 29 (EX. 11) OR WILKINSON(1963) P. 36,37). SUBSECTION: DERPOL. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" DERPOL(N,K,X,A); "VALUE"N,K,X;"INTEGER"N,K;"REAL"X;"ARRAY"A; THE MEANING OF THE FORMAL PARAMETERS IS: N: ; ENTRY: THE DEGREE OF THE POLYNOMIAL; K: ; ENTRY: THE FIRST K DERIVATIVES ARE TO BE CALCULATED; X: ; ENTRY: THE ARGUMENT OF THE POLYNOMIAL; A: ; "ARRAY"A[0:N]; ENTRY: THE COEFFICIENTS OF THE POLYNOMIAL A[0]+A[1]*X+...+A[N]*X**N; EXIT: THE J-TH DERIVATIVE IS DELIVERED IN A[J], J=0,1,...,K<= DEGREE; THE OTHER ELEMENTS ARE GENERALLY ALTERED. PROCEDURES USED : NORDERPOL = CP31242 RUNNING TIME: THE NUMBER OF ADDITIONS IS (K+1)*(N-K/2) AND THE NUMBER OF MULTIPLICATIONS IS N, IN FIRST ORDER. LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: SEE TAYPOL(THIS SECTION). 1SECTION : 2.2.1.1 (OCTOBER 1975) PAGE 3 SUBSECTION: NORDERPOL. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" NORDERPOL(N,K,X,A); "VALUE"N,K,X;"INTEGER"N,K;"REAL"X;"ARRAY"A; THE MEANING OF THE FORMAL PARAMETERS IS: N: ; ENTRY: THE DEGREE OF THE POLYNOMIAL; K: ; THE FIRST K NORMALIZED DERIVATIVES ARE TO BE CALCULATED (I.E. (J-TH DERIVATIVE)/(J FACTORIAL), J=0,1,...,K<=DEGREE). X: ; ENTRY: THE ARGUMENT OF THE POLYNOMIAL; A: ; "ARRAY"A[0:N]; ENTRY: THE COEFFICIENTS OF THE POLYNOMIAL A[0]+A[1]*X+...+A[N]*X**N; EXIT: THE J-TH NORMALIZED DERIVATIVE IS DELIVERED IN A[J] J=0,1,...,K<=DEGREE; THE OTHER ELEMENTS ARE GENERALLY ALTERED. PROCEDURES USED: NONE. REQUIRED CENTRAL MEMORY: AN AUXILIARY ARRAY OF ORDER N + 1 IS DECLARED. RUNNING TIME: THE NUMBER OF ADDITIONS IS (K+1)*(N-K/2) AND THE NUMBER OF MULTIPLICATIONS/DIVISIONS IS 2 * N + K. LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: SEE TAYPOL(THIS SECTION). 1SECTION : 2.2.1.1 (OCTOBER 1975) PAGE 4 SUBSECTION: TAYPOL. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" TAYPOL(N,K,X,A); "VALUE"N,K,X;"INTEGER"N,K;"REAL"X;"ARRAY"A; THE MEANING OF THE FORMAL PARAMETERS IS: N: ; ENTRY: THE DEGREE OF THE POLYNOMIAL; K: ; ENTRY: THE FIRST K TERMS X**J*(J-TH DERIVATIVE)/(J FACTORIAL), J=0,1,...,K<=DEGREE, ARE TO BE CALCULATED; X: ; ENTRY: THE ARGUMENT OF THE POLYNOMIAL; A: ; "ARRAY"A[0:N]; ENTRY: THE COEFFICIENTS OF THE POLYNOMIAL A[0]+A[1]*X+...+A[N]*X**N; EXIT: THE J-TH TERM X**J*(J-TH DERIVATIVE)/(J FACTORIAL), IS DELIVERED IN A[J], J=0,1,...,K<=DEGREE; THE OTHER ELEMENTS ARE GENERALLY ALTERED. PROCEDURES USED: NONE. RUNNING TIME: THE NUMBER OF ADDITIONS IS (K+1)*(N-K/2) AND THE NUMBER OF MULTIPLICATIONS IS 2 * N. LANGUAGE: ALGOL 60. METHOD AND PERFORMANCE: THE METHOD OF EVALUATION IS GIVEN BY TRAUB AND SHAW(1972,1974). LET X**J*(J-TH DERIVATIVE OF THE POLYNOMIAL)/(J FACTORIAL)= (J OVER J)*A[J]*X**J+(J+1 OVER J)*A[J+1]*X**(J+1)+...+(N OVER J)* A[N]*X**N,THEN THE J-TH DERIVATIVE (UP TO A FACTOR) CAN BE OBTAINED FROM THE BINOMIAL COEFFICIENTS FOLLOWED BY EVALUATION OF THE ABOVE INPRODUCT. THE SHAW AND TRAUB ALGORITHM PERFORMS THE BUILDING UP OF THE BINOMIAL COEFFICIENTS IMPLICITLY. WE HAVE NOT IMPLEMENTED THE MORE SOPHISTICATED ALGORITHM, BASED ON DIVISORS OF N+1, BECAUSE OF THE MORE COMPLEX APPEARANCE OF THE IMPLEMENTATION AND BECAUSE OF THE DIFFICULTY IN CHOSING THE MOST EFFICIENT DIVISOR. OUR (RESTRICTED) IMPLEMENTATION OF THE ONE-PARAMETER FAMILY OF ALGORITHMS PRESERVES THE LINEAR NUMBER OF MULTIPLICATIONS (2*N (NORDERPOL, TAYPOL) AND 3*N (DERPOL)). THE ABSOLUTE ERROR IS OF ORDER MAX((N OVER N)*A[N]*X**(N-K),..., (N OVER K)*A[K]), FOR THE K-TH NORMALIZED DERIVATIVE (SEE VAN DER LAAN OR WOZNIAKOWSKI). 1SECTION : 2.2.1.1 (DECEMBER 1979) PAGE 5 REFERENCES: [1].SHAW,M. AND J. TRAUB: ON THE NUMBER OF MULTIPLICATIONS FOR THE EVALUATION OF A POLYNOMIAL AND SOME OF ITS DERIVATIVES (21 P.). JOURN. ACM, 1974, VOL. 21, NO. 1, P. 161-167. [2].STOER,J. : EINFUEHRUNG IN DIE NUMERISCHE MATHEMATIK 1. SPRINGER, 1972. [3].VAN DER LAAN C.G.: ORTHOGONAL POLYNOMIALS IN NUMERICAL ANALYSIS. 1. ERROR ANALYSIS OF RECURRENCE RELATIONS IN FLOATING-POINT COMPUTATION, (TO APPEAR). [4].WILKINSON,J.H. : ROUNDING ERRORS IN ALGEBRAIC PROCESSES. HSO, NOTES ON APPLIED SCIENCES NO. 32, 1963. [5].WOZNIAKOWSKI,H.: ROUNDING ERROR ANALYSIS FOR THE EVALUATION OF A POLYNOMIAL AND SOME OF ITS DERIVATIVES. SIAM J. OF NUM. AN. VOL 11, NO. 4, P. 780-787. EXAMPLE OF USE: AS A FORMAL TEST OF THE PROCEDURE DERPOL THE DERIVATIVES OF THE POLYNOMIAL 3*X**3-2*X**2+X-1 ARE CALCULATED AT X=1. "BEGIN""ARRAY"A[0:3]; "PROCEDURE"DERPOL(N,K,X,A);"CODE"31243; A[3]:=3;A[2]:=-2;A[1]:=1;A[0]:=-1; DERPOL(3,3,1,A); OUTPUT(61,"(" "("THE 0-TH UNTIL AND INCLUDING THE 3-TH DERIVATIVES :")", 4(BZDB)")",A[0],A[1],A[2],A[3]); "END" EXAMPLE OF USE; THE 0-TH UNTIL AND INCLUDING THE 3-TH DERIVATIVES : 1 6 14 18 SOURCE TEXT(S): 0"CODE" 31040; "REAL" "PROCEDURE" POL(N,X,A); "VALUE" N,X;"INTEGER" N;"REAL" X;"ARRAY" A; "BEGIN" "REAL" R; R:= 0; "FOR" N:= N "STEP" -1 "UNTIL" 0 "DO" R:=R*X + A[N]; POL:= R "END" POL; "EOP" 1SECTION : 2.2.1.1 (OCTOBER 1975) PAGE 6 0"CODE" 31241; "PROCEDURE" TAYPOL(N,K,X,A); "VALUE" N,K,X; "INTEGER" N,K;"REAL" X;"ARRAY" A; "IF" X^= 0 "THEN" "BEGIN" "INTEGER" I,J,NM1; "REAL" XJ,AA,H; XJ:=1; "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" XJ:=XJ*X;A[J]:=A[J]*XJ "END"; AA:=A[N];NM1:=N-1; "FOR" J:= 0 "STEP" 1 "UNTIL" K "DO" "BEGIN" H:=AA; "FOR" I:= NM1 "STEP" -1 "UNTIL" J "DO" H:= A[ I]:=A[I]+H "END" "END" "ELSE" "FOR" K:= K "STEP" -1 "UNTIL" 1 "DO" A[K]:=0; "EOP" 0"CODE" 31242; "PROCEDURE" NORDERPOL (N,K,X,A); "VALUE" N,K,X; "INTEGER" N,K;"REAL" X;"ARRAY" A; "IF" X^= 0 "THEN" "BEGIN" "INTEGER" I,J,NM1; "REAL" XJ,AA,H; "ARRAY" XX[0:N]; XJ:=1; "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" XJ:=XX[J]:=XJ*X;A[J]:=A[J]*XJ "END"; H:=AA:=A[N];NM1:=N-1; "FOR" I:= NM1 "STEP" -1 "UNTIL" 0 "DO" H:= A[I]:=A[I]+H; "FOR" J:= 1 "STEP" 1 "UNTIL" K "DO" "BEGIN" H:=AA; "FOR" I:= NM1 "STEP" -1 "UNTIL" J "DO" H:= A[ I]:=A[I]+H; A[J]:=H/XX[J] "END" "END" NORDERPOL ; "EOP" 0"CODE" 31243; "PROCEDURE" DERPOL (N,K,X,A); "VALUE" N,K,X; "INTEGER" N,K;"REAL" X;"ARRAY" A; "BEGIN" "INTEGER" J; "REAL" FAC; "PROCEDURE"NORDERPOL(N,K,X,A);"CODE" 31242; FAC:=1; NORDERPOL (N,K,X,A); "FOR" J:= 2 "STEP" 1 "UNTIL" K "DO" "BEGIN" FAC:=FAC*J;A[J]:=A[J]*FAC "END" "END" DERPOL ; "EOP" 1SECTION : 6.10.3 (DECEMBER 1978) PAGE 1 AUTHOR: M. BAKKER. INSTITUTE: MATHEMATICAL CENTRE. BRIEF DESCRIPTION: THIS SECTION CONTAINS THE PROCEDURES SPHER BESS J: THIS PROCEDURE CALCULATES THE SPHERICAL BESSEL FUNCTIONS J[K+.5](X)*SQRT(PI/(2*X)),K=0, ..., N, WHERE J[K+.5](X) DENOTES THE BESSEL FUNCTION OF THE FIRST KIND OF ORDER K+.5; X>= 0; SPHER BESS Y: THIS PROCEDURE CALCULATES THE SPHERICAL BESSEL FUNCTIONS Y[K+.5](X)*SQRT(PI/(2*X)), K=0, ..., N, WHERE Y[K+.5](X) DENOTES THE BESSEL FUNCTION OF THE THIRD KIND OF ORDER K+.5; X SHOULD BE POSITIVE; SPHER BESS I: THIS PROCEDURE CALCULATES THE MODIFIED SPHERICAL BESSEL FUNCTIONS I[K+.5](X)*SQRT(PI/(2*X))), K=0, ..., N, WHERE I[K+.5](X) DENOTES THE MODIFIED BESSEL FUNCTION OF THE FIRST KIND OF ORDER K+.5; X>=0; NONEXP SPHER BESS I: THIS PROCEDURE CALCULATES THE MODIFIED SPHERICAL BESSEL FUNCTIONS MULTIPIED BY EXP(-X) EXP(-X)*I[K+.5](X)*SQRT(PI/(2*X)), K=0, ...,N, WHERE I[K+.5](X) DENOTES THE MODIFIED BESSEL FUNCTION OF THE FIRST KIND OF ORDER K+.5; X>= 0; SPHER BESS K: THIS PROCEDURE CALCULATES THE MODIFIED SPHERICAL BESSEL FUNCTIONS K[I+.5](X)*SQRT(PI/(2*X)), I=0, ...,N, WHERE K[I+.5](X) DENOTES THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER I+.5; X>0; NONEXP SPHER BESS K: THIS PROCEDURE CALCULATES THE MODIFIED SPHERICAL BESSEL FUNCTIONS MULTIPLIED BY EXP(+X) EXP(+X)*K[I+.5](X)*SQRT(PI/(2*X)), I=0, ..., N, WHERE K[I+.5](X) DENOTES THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER I+.5; X>0; 1SECTION : 6.10.3 (DECEMBER 1978) PAGE 2 KEYWORDS: BESSEL FUNCTIONS, SPHERICAL BESSEL FUNCTIONS, MODIFIED SPHERICAL BESSEL FUNCTIONS. REFERENCES: [1]. ABRAMOWITZ, M., AND STEGUN, I. (EDS), HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICAL TABLES. APPL. MATH. SER. 55, U.S. GOVT. PRINTING OFFICE, WASHINGTON, D.C. , 1974. [2]. GAUTSCHI, W., COMPUTATIONAL ASPECTS OF THREE TERM RECURRENCE RELATIONS. SIAM REVIEW, VOLUME 9(1967), NUMBER 1, P.24 FF. SUBSECTION: SPHER BESS J. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" SPHER BESS J (X, N, J); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" J; "CODE" 35150; THE MEANING OF THE FORMAL PARAMETERS IS: X: < ARITHMETIC EXPRESSION >; THE VALUE OF THE ARGUMENT; X > = 0; N: < ARITHMETIC EXPRESSION >; THE UPPER BOUND OF THE INDICES OF THE ARRAY J; N > = 0; J: < ARRAY IDENTIFIER >; "ARRAY" J[0:N]; EXIT: J[K] HAS THE VALUE OF THE SPHERICAL BESSEL FUNCTION J[K+.5](X) * SQRT(PI/(2*X)), 0< = K < = N; PROCEDURES USED: START = CP 35185. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. 1SECTION : 6.10.3 (DECEMBER 1978) PAGE 3 METHOD AND PERFORMANCE: AT FIRST THE RATIO OF TWO CONSEQUENT ARRAY ELEMENTS IS COMPUTED BY MEANS OF A BACKWARD RECURRENCE FORMULA USING MILLER 'S METHOD (SEE[2, P.46-52]) AND HENCE ALL THE ARRAY ELEMENTS ARE COMPUTED SINCE THE ZEROTH ELEMENT IS KNOWN TO BE SIN(X)/X. THE STARTING VALUE IS COMPUTED BY START. RUNNING TIME: ROUGHLY PROPERTIONAL TO THE MAXIMUM OF X AND N. EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X ; "ARRAY" J[0:2]; "INTEGER" N; "PROCEDURE" SPHER BESS J (X, N, J); "CODE" 35150; X:= 1.5; N:= 2; SPHER BESS J(X, N, J); OUTPUT(61, "("/, "("X=")" D.D, B"("N=")"D, 3(3B-.14D"-ZD)")", X, N, J[0], J[1], J[2]) "END" PRINTS THE FOLLOWING RESULTS: X=1.5 N=2 .66499665773603"0 .3961729707122"0 .12734928368841"0 SUBSECTION: SPHER BESS Y. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" SPHER BESS Y(X, N, Y); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" Y; "CODE" 35151; THE MEANING OF THE FORMAL PARAMETERS IS : X: < ARITHMETIC EXPRESSION >; THE ARGUMENT OF THE BESSEL FUNCTIONS; X > 0; N: < ARITHMETIC EXPRESSION >; THE UPPER BOUND OF THE INDICES OF THE ARRAY Y; N > = 0; Y: < ARRAY IDENTIFIER >; "ARRAY" Y[0:N]; EXIT: Y[K] HAS THE VALUE OF THE K-TH SPHERICAL BESSEL FUNCTION OF THE SECOND KIND; 1SECTION : 6.10.3 (DECEMBER 1978) PAGE 4 PROCEDURES USED: NONE. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. METHOD AND PERFORMANCE: Y[0] AND Y[1] ARE GIVEN IN [1, P.438, FORMULA 10.1.12] AND Y[2], ..., Y[N] ARE COMPUTED BY USING THE RECURRENCE FORMULA Y[K]:= ((2*K-1)/X) * Y[K-1] - Y[K-2], K > = 2. EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X; "INTEGER" N; "ARRAY" Y[0:2]; "PROCEDURE" SPHER BESS Y(X, N, Y); "CODE" 35151; X:= 1.5707 96326 79489; "COMMENT" X= PI/2; N:= 2; SPHER BESS Y(X, N, Y); OUTPUT(61, "("5(4B-.10D"-ZD)")", X, N, Y) "END" PRINTS THE FOLLOWING RESULTS: .15707963271"1 .2000000000"1 -.6223649549"-13 -.6366197724"0 -.1215854200"0 SUBSECTION: SPHER BESS I. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" SPHER BESS I(X, N, I); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" I; "CODE" 35152; THE MEANING OF THE FORMAL PARAMETERS IS: X: < ARITHMETIC EXPRESSION >; THE ARGUMENT OF THE BESSEL FUNCTIONS; X > = 0; N: < ARITHMETIC EXPRESSION >; THE UPPER BOUND OF THE INDICES OF THE ARRAY I; N > = 0; I: < ARRAY IDENTIFIER >; "ARRAY" I[0:N]; EXIT: I[K] HAS THE VALUE OF THE MODIFIED SPHERICAL BESSEL FUNCTION AS DESCRIBED IN [1, CH.10.2]. 1SECTION : 6.10.3 (DECEMBER 1978) PAGE 5 METHOD AND PERFORMANCE: AT FIRST THE NONEXPONENTIAL MODIFIED SPHERICAL BESSEL FUNCTIONS ARE COMPUTED BY USING THE PROCEDURE NONEXP SPHER BESS I; AFTERWARDS THEY ARE MULTIPLIED BY EXP(X). REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. PROCEDURES USED: NONEXP SPHER BESS I = CP 35154. EXAMPLE OF USE: THE PROGRAM SHOWS THAT THE RESULTS OF SPHER BESS I AND NONEXP SPHER BESS I DIFFER ONLY BY A FACTOR EXP(X): "BEGIN" "REAL" X, EXPX; "INTEGER" N; "ARRAY" I1, I2[0:3]; "PROCEDURE" SPHER BESS I(X, N, I); "CODE" 35152; "PROCEDURE" NE SPHER BESS I(X, N, NEI); "CODE" 35154; X:=1; EXPX:= EXP(X); N:= 3; SPHER BESS I(X, N,I1); NESPHER BESS I(X, N, I2);"FOR" N:=0, 1, 2, 3 "DO" OUTPUT(61, "("/ZD, 2(5B-.14D"-ZD)")", N, I1[N], I2[N]*EXPX) "END" RESULTS: 0 .11752011936438" 1 .11752011936438" 1 1 .36787944117144" 0 .36787944117144" 0 2 .71562870129474"-1 .71562870129474"-1 3 .10065090524070"-1 .10065090524070"-1 SUBSECTION: NONEXP SPHER BESS I. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" NONEXP SPHER BESS I(X, N, I); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" I; "CODE" 35154; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTIONS; X >= 0; N: ; THE UPPER BOUND OF THE INDICES OF THE ARRAY I; N >= 0; I: ; "ARRAY" I[0:N]; EXIT: I[K] HAS THE VALUE OF THE FUNCTION I[K+.5](X)*EXP(-X)*SQRT(PI/(2*X)), K=0, ..., N, N >=0. 1SECTION : 6.10.3 (DECEMBER 1978) PAGE 6 PROCEDURES USED: SINH = CP 35111, START = CP 35185. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. METHOD AND PERFORMANCE: THE RATIO OF TWO SUBSEQUENT ELEMENTS IS COMPUTED USING A BACKWARD RECURRENCE FORMULA ACCORDING MILLER'S METHOD (SEE[2]); SINCE THE ZEROETH ELEMENT IS KNOWN TO BE (1-EXP(-2*X))/(2*X), THE OTHER ELEMENTS FOLLOW IMMEDIATELY.THE STARTING VALUE IS COMPUTED BY START. EXAMPLE OF USE: SEE SPHER BESS I. SUBSECTION: SPHER BESS K. CALLING SEQUENCE : THE HEADING OF THE PROCEDURE READS: "PROCEDURE" SPHER BESS K(X, N, K); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" K; "CODE" 35153; THE MEANING OF THE FORMAL PARAMETERS IS: X: < ARITHMETIC EXPRESSION >; THE ARGUMENT VALUE; X > 0; N: < ARITHMETIC EXPRESSION >; THE UPPER BOUND OF THE INDICES OF THE ARRAY K; N > = 0; K: < ARRAY IDENTIFIER >; "ARRAY" K[0:N]; EXIT: K[J] HAS THE VALUE OF THE J-TH MODIFIED SPHERICAL BESSEL FUNCTION OF THE THIRD KIND, 0 < = J < = N. PROCEDURES USED: NON EXP SPHER BESS K = CP 35155. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. METHOD AND PERFORMANCE: AT FIRST THE NONEXPONENTIAL MODIFIED SPHERICAL BESSEL FUNCTIONS OF THE THIRD KIND ARE COMPUTED BY THE PROCEDURE NONEXP SPHER BESS K; AFTERWARDS THEY ARE MULTIPLIED BY EXP(-X). 1SECTION : 6.10.3 (DECEMBER 1978) PAGE 7 EXAMPLE OF USE: THE FOLLOWING PROGRAM SHOWS THAT THE RESULTS OF THE PROCEDURES SPHER BESS K EN NONEXP SPHER BESS K DIFFER ONLY BY A FACTOR EXP(X); "BEGIN" "REAL" X, EXPX; "INTEGER" N; "ARRAY" K1, K2[0:3]; "PROCEDURE" SPHER BESS K (X, N, K); "CODE" 35153; "PROCEDURE" NESPHER BESS K (X, N, K); "CODE" 35155; X:= 2; EXPX:= EXP(-X); N:= 3; SPHER BESS K (X, N, K1); NESPHER BESS K (X, N, K2); "FOR" K:= 0, 1, 2, 3 "DO" OUTPUT(61, "("/D, 2(5B-.14D"-ZD)")", N, K1[N], K2[N]*EXPX) "END" RESULTS: 0 .10629208289691"0 .10629208289691"0 1 .15943812434536"0 .15943812434536"0 2 .34544926941495"0 .34544926941494"0 3 .10230612978828"1 .10230612978828"1 SUBSECTION: NONEXP SPHER BESS K. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" NON EXP SPHER BESS K(X, N, K); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" K; "CODE" 35155; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTIONS; X > 0; N: ; THE UPPER BOUND OF THE INDICES OF THE ARRAY K; N >= 0; K: ; "ARRAY" K[0:N]; EXIT: K[J] HAS THE VALUE OF THE FUNCTION K[J+.5](X)*EXP(X)*SQRT(PI/(2*X)), J=0,...,N. PROCEDURES USED: NONE. REQUIRED CENTRAL MEMORY : NO AUXILIARY ARRAYS ARE DECLARED. METHOD AND PERFORMANCE: THE FUNCTIONS ARE COMPUTED BY USING THE (NUMERICALLY STABLE) RECURRENCE FORMULA : K[J]=((2*J-1)/X)*K[J-1]+K[J-2], J >=2, K[0]=PI/(2*X), K[1]=K[0]*(1+1/X) . EXAMPLE OF USE: SEE SPHER BESS K. 1SECTION : 6.10.3 (DECEMBER 1978) PAGE 8 SOURCE TEXT(S): "CODE" 35150; "COMMENT" SPHERICAL BESSEL FUNCTIONS J[.5](X), , J[N+.5](X); "PROCEDURE" SPHER BESS J(X, N, J); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" J; "IF" X = 0 "THEN" "BEGIN" J[0]:= 1; "FOR" N:= N "STEP" -1 "UNTIL" 1 "DO" J[N]:=0 "END" "ELSE" "IF" N = 0 "THEN" "BEGIN" "REAL" X2; "IF" ABS(X) < .015 "THEN" "BEGIN" X2:= X * X / 6; J[0]:= 1 + X2 * (X2 * .3 - 1) "END" "ELSE" J[0]:= SIN(X)/X "END" "ELSE" "BEGIN" "INTEGER" M; "REAL" R, S; "INTEGER" "PROCEDURE" START(X,N,T); "CODE" 35185; R:= 0; M:= START(X,N,0); "FOR" M:= M "STEP" - 1 "UNTIL" 1 "DO" "BEGIN" R:= 1 / ((M + M + 1) / X - R); "IF" M <= N "THEN" J[M]:= R "END"; "IF" X < .015 "THEN" "BEGIN" S:= X * X / 6; J[0]:= R:= S * (S * .3 - 1) + 1 "END" "ELSE" J[0]:= R:= SIN(X) / X; "FOR" M:= 1 "STEP" 1 "UNTIL" N "DO" J[M]:= R:= J[M] * R; "END" SPHER BESS J; "EOP" "CODE" 35151; "COMMENT" SPHERICAL BESSEL FUNCTIONS Y[.5](X), , Y[N+.5](X); "PROCEDURE" SPHER BESS Y(X, N, Y); "VALUE" X, N; "INTEGER" N; "REAL" X; "ARRAY" Y; "IF" N=0 "THEN" Y[0]:= - COS(X)/X "ELSE" "BEGIN" "REAL" YI, YI1, YI2; "INTEGER" I; YI2:= Y[0]:= -COS(X)/X; YI1:= Y[1]:= (YI2 - SIN(X))/X; "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO" "BEGIN" Y[I]:= YI:= -YI2 + (I+I-1) * YI1/X; YI2:= YI1; YI1:= YI "END" "END"; "EOP" 1SECTION : 6.10.3 (DECEMBER 1978) PAGE 9 "CODE" 35152; "COMMENT" SPHERICAL BESSEL FUNCTIONS I[.5](X), , I[N+.5](X); "PROCEDURE" SPHER BESS I(X, N, I); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" I; "IF" X= 0 "THEN" "BEGIN" I[0]:=1; "FOR" N:= N "STEP" -1 "UNTIL" 1 "DO" I[N]:= 0 "END" "ELSE" "BEGIN" "REAL" EXPX; "PROCEDURE" NONEXP SPHER BESS I(X, N, I); "CODE" 35154; EXPX:= EXP(X); NONEXP SPHER BESS I(X, N, I); "FOR" N:= N "STEP" - 1 "UNTIL" 0 "DO" I [N]:= I [N] * EXPX "END" SPHER BESS I; "EOP" "CODE" 35153; "COMMENT" MODIFIED SPHERICAL BESSEL FUNCTIONS K[.5](X), , K[N+.5](X); "PROCEDURE" SPHER BESS K(X, N, K); "VALUE" X, N; "INTEGER" N; "REAL" X; "ARRAY" K; "BEGIN" "REAL" EXPX; "PROCEDURE" NONEXP SPHER BESS K(X, N, K); "CODE" 35155; EXPX:= EXP(-X); NONEXP SPHER BESS K(X, N, K); "FOR" N:= N "STEP" -1 "UNTIL" 0 "DO" K[N]:= K[N] * EXPX "END"; "EOP" 1SECTION : 6.10.3 (DECEMBER 1978) PAGE 10 "CODE" 35154; "PROCEDURE" NONEXP SPHER BESS I(X, N, I); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" I; "IF" X= 0 "THEN" "BEGIN" I[0]:=1; "FOR" N:= N "STEP" -1 "UNTIL" 1 "DO" I[N]:= 0 "END" "ELSE" "BEGIN" "REAL" X2, R, S; "INTEGER" M; "REAL" "PROCEDURE" SINH(X); "CODE" 35111; "INTEGER" "PROCEDURE" START(X,N,T); "CODE" 35185; X2:= X+X; I[0]:= X2:= "IF" X = 0 "THEN" 1 "ELSE" "IF" X2 < 0.7 "THEN" SINH(X) / (X * EXP(X)) "ELSE" (1-EXP(-X2))/X2; "IF" N= 0 "THEN" "GO TO" EXIT; R:= 0; M:= START(X,N,1); "FOR" M:= M "STEP" -1 "UNTIL" 1 "DO" "BEGIN" R:= 1/((M+M+1)/X+R); "IF" M <= N "THEN" I[M]:= R "END"; "FOR" M:= 1 "STEP" 1 "UNTIL" N "DO" I[M]:= X2:= X2 * I[M]; EXIT: "END"; "EOP" "CODE" 35155; "PROCEDURE" NONEXP SPHER BESS K(X, N, K); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" K; "BEGIN" "INTEGER" I; "REAL" KI, KI1, KI2; X:= 1/X; K[0]:= KI2:= X*1.5707963267949; "IF" N=0 "THEN" "GO TO" EXIT; K[1]:= KI1:= KI2 * (1+X); "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO" "BEGIN" K[I]:= KI:= KI2 + (I+I-1) * X * KI1; KI2:= KI1; KI1:= KI "END"; EXIT: "END"; "EOP" 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 1 AUTHORS: M.BAKKER AND N.M.TEMME. CONTRIBUTOR: R.MONTIJN. INSTITUTE: MATHEMATICAL CENTRE. RECEIVED: 781101. BRIEF DESCRIPTION: THIS SECTION CONTAINS THE PROCEDURES: BESS JAPLUSN: THIS PROCEDURE CALCULATES THE BESSEL FUNCTIONS OF THE FIRST KIND OF ORDER A+K (0<=K<=N, 0<=A<1) AND ASSIGNS THEM TO AN ARRAY. THE ARGUMENT MUST BE NON-NEGATIVE. BESS YA01: THIS PROCEDURE CALCULATES THE BESSEL FUNCTIONS OF THE SECOND KIND (ALSO CALLED NEUMANN'S FUNCTIONS) OF ORDER A AND A+1 AND ARGUMENT X>0. BESS YAPLUSN: THIS PROCEDURE GENERATES AN ARRAY OF BESSEL FUNCTIONS OF THE SECOND KIND OF ORDER A+N, N=0, 1, 2, ..., NMAX, AND ARGUMENT X>0. THE BESSEL FUNCTIONS OF THE SECOND KIND CORRESPOND TO THE FUNCTION DEFINED IN FORMULA 9.1.2 OF REFERENCE [1]. BESS PQA01: THIS PROCEDURE IS AN AUXILIARY PROCEDURE FOR THE COMPUTATION OF THE BESSEL FUNCTIONS FOR LARGE VALUES OF THEIR ARGUMENT. BESS ZEROS: THIS PROCEDURE CALCULATES THE FIRST N ZEROS OF A BESSEL FUNCTION OF THE FIRST OR THE SECOND KIND OR ITS DERRIVATIVE. START: THIS IS AN AUXILIARY PROCEDURE WHICH COMPUTES A STARTING VALUE OF AN ALGORITHM USED IN SEVERAL BESSEL FUNCTION PROCEDURES. 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 2 KEYWORDS: BESSEL FUNCTION, BESSEL FUNCTION OF THE SECOND KIND, NEUMANN'S FUNCTION, ZEROS OF BESSEL FUNCTIONS. REFERENCES: [1]. ABRAMOWITZ, M., AND STEGUN, I. (EDS), HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICAL TABLES. APPL. MATH. SER. 55, U.S. GOVT. PRINTING OFFICE, WASHINGTON, D.C. , 1974. [2]. GAUTSCHI, W., COMPUTATIONAL ASPECTS OF THREE TERM RECURRENCE RELATIONS. SIAM REVIEW, VOLUME 9(1967), NUMBER 1, P.24 FF. [3]. TEMME, N.M. ON THE NUMERICAL EVALUATION OF THE ORDINARY BESSEL FUNCTION OF THE SECOND KIND. J. COMP. PHYS., 21, P. 343 FF, 1976. [4]. WATSON, G.N. A TREATISE ON THE THEORY OF BESSEL FUNCTIONS. CAMBRIDGE UNIV. PRESS, LONDON AND NEW YORK, 1945. [5]. TEMME, N.M., SPECIALE FUNCTIES, IN: COLLOQUIUM NUMERIEKE PROGRAMMATUUR, J.C.P. BUS (RED.), MC SYLLABUS 29.1B, MATHEMATICAL CENTRE, AMSTERDAM, 1976. [6]. TEMME, N.M., AN ALGOLRITHM WITH ALGOL 60 IMPLEMENTATION FOR THE CALCULATION OF THE ZEROS OF A BESSEL FUNCTION, REPORT TW 179 MATHEMATICAL CENTRE, AMSTERDAM, 1978. SUBSECTION: BESS JAPLUSN. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS JAPLUSN(A, X, N, JA); "VALUE" A, X, N; "INTEGER" N; "REAL" A, X; "ARRAY" JA; "CODE" 35180; 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 3 THE MEANING OF THE FORMAL PARAMETERS IS: A: < ARITHMETIC EXPRESSION > ; THE NONINTEGER PART OF THE ORDER; 0 <= A < 1; X: < ARITHMETIC EXPRESSION >; THE ARGUMENT VALUE; X > = 0; N: < ARITHMETIC EXPRESSION >; THE UPPER BOUND OF THE INDICES OF THE ARRAY JA; JA: < ARRAY IDENTIFIER >; "ARRAY" JA[0:N]; EXIT: JA[K] IS ASSIGNED THE VALUE OF THE BESSEL FUNCTION OF THE FIRST KIND J[K+A](X), 0 < = K < = N. PROCEDURES USED: BESS J = CP 35162, SPHER BESS J = CP 35150, GAMMA = CP 35061, START = CP 35185. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. METHOD AND PERFORMANCE: IN ALL THE CASES THE BESSEL FUNCTIONS ARE COMPUTED ACCORDING TO THE MILLER METHOD DISCRIBED IN [2, P.46-52]. THE STARTING VALUE IS COMPUTED BY THE PROCEDURE START. RUNNING TIME: ROUGHLY PROPORTIONAL TO THE MAXIMUM OF X AND N. EXAMPLE OF USE: "BEGIN" "INTEGER" N; "REAL" A, X; "ARRAY" JA[0:2]; "PROCEDURE" BESS JAPLUSN(A, X, N, JA); "CODE" 35180; X:= 2; A:= .78; N:= 2; BESS JAPLUSN(A, X, N, JA); OUTPUT(61, "("/, "("X=")"D, 3B"("A=")".DD, 3B"("N=")"D, /, 3(3B-.14D"-ZD)")", X, A, N, JA[0], JA[1], JA[2]) "END" RESULTS: X=2 A= .78 N=2 .57306126928364"0 .41529475124424" 0 .16616338793111" 0 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 4 SUBSECTION: BESS YA01. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS YA01(A, X, YA, YA1); "VALUE" A, X; "REAL" A, X, YA, YA1; "CODE" 35181; THE MEANING OF THE FORMAL PARAMETERS IS: A: ; THE ORDER; X: ; THE ARGUMENT. THIS ARGUMENT SHOULD SATISFY X>0; YA: ; EXIT: THE NEUMANN FUNCTION OF ORDER A AND ARGUMENT X; YA1: ; EXIT: THE NEUMANN FUNCTION OF ORDER A+1. PROCEDURES USED: RECIP GAMMA = CP 35060; BESS PQA01 = CP 35183; SINH = CP 35111. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. METHOD AND PERFORMANCE: FOR 0=3 THE PROCEDURE CALLS FOR THE PROCEDURE BESS PQA01 (SEE SUBSECTION BESS PQA01). THE RELATIVE ACCURACY IS ABOUT "-13, EXCEPT FOR LARGE VALUES OF X; IN THAT CASE THE ACCURACY MAINLY DEPENDS ON THE ACCURACY OF THE FUNCTIONS SIN(X) AND COS(X). 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 5 EXAMPLE OF USE: THE PROGRAM: "BEGIN" "REAL" P, Q; "PROCEDURE" BESS YA01(A, X, YA, YA1); "CODE" 35181; BESS YA01(0, 1, P, Q); OUTPUT(61, "("2(N)")", P, Q) "END" YIELDS THE FOLLOWING RESULTS +8.8256964215677"-002 -7.8121282130028"-001. SUBSECTION: BESS YAPLUSN. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS YAPLUSN(A, X, NMAX, YAN); "VALUE" A, X, NMAX; "REAL" A, X; "INTEGER" NMAX; "ARRAY" YAN; "CODE" 35182; THE MEANING OF THE FORMAL PARAMETERS IS: A: ; THE ORDER; X: ; THE ARGUMENT; THIS ARGUMENT SHOULD SATISFY X>0; NMAX: ; THE UPPER BOUND OF THE INDICES OF THE ARRAY YAN; YAN: ; "ARRAY" YAN[0:NMAX]; NMAX>=0; EXIT: THE VALUES OF THE BESSEL FUNCTIONS OF THE SECOND KIND OF ORDER A+K, FOR THE ARGUMENT X ARE ASSIGNED TO YAN[K],0<=K<=NMAX. PROCEDURES USED: BESS YA01 = CP 35181. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. 1SECTION: 6.10.1 (DECEMBER 1978) PAGE 6 METHOD AND PERFORMANCE: THE RECURRENCE RELATION YAN[N+1]= -YAN[N-1] + 2*(N+A)*YAN[N]/X IS USED. THE INITIAL VALUES ARE OBTAINED FROM THE PROCEDURE BESS YA01. THE RECURRENCE RELATION IS NUMERICALLY STABLE IN THE FORWARD DIRECTION (IF A >= 0). EXAMPLE OF USE: THE PROGRAM: "BEGIN" "ARRAY" YAN[0:2]; "PROCEDURE" BESS YAPLUSN(A, X, NMAX, YAN); "CODE" 35182; BESS YAPLUSN(0, 1, 2, YAN); OUTPUT(61, "("3(N)")", YAN[0], YAN[1], YAN[2]) "END" YIELDS THE FOLLOWING RESULTS +8.8256964215677"-002 -7.8121282130028"-001 -1.6506826068163"+000. SUBSECTION: BESS PQA01. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS PQA01(A, X, PA, QA, PA1, QA1); "VALUE" X, A; "REAL" X, A, PA, QA, PA1, QA1; "CODE" 35183; THE MEANING OF THE FORMAL PARAMETERS IS: A: ; THE ORDER; X: ; THE ARGUMENT. THIS ARGUMENT SHOULD SATISFY X>0; PA: ; EXIT: THIS FUNCTION CORRESPONDS TO THE FUNCTION P(X, A) DEFINED ON P. 205 OF REFERENCE [4]. SEE ALSO REFERENCE [1], FORMULA 9.2.6; QA: ; EXIT: THIS FUNCTION CORRESPONDS TO THE FUNCTION Q(X, A) DEFINED ON P.205 OF REFERENCE [4]. SEE ALSO REFERENCE [1], FORMULA 9.2.6; 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 7 PA1: ; EXIT: THE FUNCTION P(X, A+1); QA1: ; EXIT: THE FUNCTION Q(X, A+1). PROCEDURES USED: BESS JAPLUSN = CP35180, BESS YA01 = CP35181. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. METHOD AND PERFORMANCE: X < 3 : PA, QA, PA1, QA1 ARE COMPUTED FROM THE RELATIONS PA = B * (YA0 * S + JA0 * C), QA = B * (YA0 * C - JA0 * S), PA1 = B * (JA1 * S - YA1 * C), QA1 = B * (JA1 * C + YA1 * S), WHERE B = SQRT(HALFPI * X), C = COS(X - HALFPI * (A + .5)), S = SIN(X - HALFPI * (A + .5)), HALFPI = 1.57079 63267 9489, YA0 = Y[A](X), YA1 = Y[A + 1](X), JA0 = J[A](X), JA1 = J[A + 1](X); X >= 3: THE METHOD IS DESCRIBED IN REFERENCE [3]. IT DEPENDS ON USING A MILLER ALGORITHM FOR CONFLUENT HYPERGEOMETRIC FUNCTIONS. THE ACCURACY IS ABOUT "-13 AND IS BETTER FOR LARGE X. THE FUNCTIONS PA AND QA CAN ALSO BE USED FOR THE COMPUTATION OF THE BESSEL FUNCTION J OF THE FIRST KIND. SEE REFERENCE[1], FORMULA 9.2.5. 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 8 EXAMPLE OF USE: FROM SOME PROPERTIES OF THE BESSEL FUNCTIONS IT CAN BE PROVED THAT PA*PA1+QA*QA1=1, WHATEVER X AND A. IN THE FOLLOWING PROGRAM WE VERIFY THIS RELATION. "BEGIN" "REAL" A, X, P, Q, R, S; "PROCEDURE" BESS PQA01(A, X, PA, QA, PA1, QA1); "CODE" 35183; "FOR" X:= 1, 3, 5, 10, 15, 20, 50 "DO" "BEGIN" BESS PQA01(0, X, P, Q, R, S); OUTPUT(61, "("BB, D.2D"+3D")", ABS(P*R+Q*S-1)) "END" "END" THIS PROGRAM GIVES THE FOLLOWING RESULTS: 1.42"-014 7.11"-015 7.11"-015 7.11"-015 1.42"-014 0.00"+000 2.13"-014. SUBSECTION: BESS ZEROS. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS ZEROS(A,N,Z,D); "VALUE" A,N,D; "REAL" A; "INTEGER" N,D; "ARRAY" Z; "CODE" 35184; THE MEANING OF THE FORMAL PARAMETERS IS: A: ; THE ORDER OF THE BESSEL FUNCTION, A>=0. N: ; THE NUMBER OF ZEROS TO BE EVALUATED, N>=1. Z: ; "ARRAY" Z[1:N]; EXIT: Z[J] IS THE J-TH ZERO OF THE SELECTED BESSEL FUNCTON; D: ; THE CHOICE OF D DETERMINES THE TYPE OF THE BESSEL FUNCTION OF WHICH THE ZEROS ARE COMPUTED: IF D=1 THEN JA , IF D=2 THEN YA , IF D=3 THEN JA-PRIME, IF D=4 THEN YA-PRIME. 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 9 PROCEDURES USED: BESS PQA01 = CP 35183. REQUIRED CENTRAL MEMMORY: NO AUXILIARY ARRAYS ARE USED. RUNNING TIME: DEPENDS ON THE VALUES OF A AND N AND ON THE MUMBER OF ITERATIONS IN THE ALGORITHM. FROM TESTS IT FOLLOWS THAT FOR EACH ZERO AT MOST 3 EVALUATIONS OF THE PROCEDURE BESS PQA01 ARE NEEDED. METHOD AND PERFORMANCE: A FIRST APPROXIMATION OF THE ZEROS OF THE SELECTED BESSEL FUNCTION IS CALCULATED BY MEANS OF THE ASYMPTOTIC EXPANTIONS ( SEE THE FORMULAS 9.5.12, 9.5.13 ( FOR A < 3 ) AND 9.5.22, 9.5.24( FOR A >= 3 ) OF REF [1] ). THIS VALUE IS CORRECTED BY THE USE OF A FOURTH ORDER NEWTON-RAPHSON METHOD AS DISCRIBED ON P. 179 OF REF [6]. MORE DETAILS CAN BE FOUND IN REF [7]. A RELATIVE PRECISION OF 13 DIGITS IS PERSUED. THE COMPUTATION OF A ZERO IS TERMINATED IF THIS ACCURRACY IS ACHIEVED OR IF MORE THAN 5 ITERATIONS ARE NEEDED. THE PROCEDURE DOES NOT CHECK ON THE RANGE OF THE PARAMETERS A,N AND D. EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" A; "INTEGER" N,D; "ARRAY" Z[1:2]; "PROCEDURE"BESS ZEROS(A,N,Z,D); "CODE" 35184; A:=3.14; N:= 2; D:= 2; BESS ZEROS(A,N,Z,D); OUTPUT(61,"M"N,/,N")",Z[1],Z[2]) "END" PRINTS THE FIRST TWO ZEROS OF THE BESSEL FUNCTION Y OF THE ORDER 3.14; THE RESULT IS: +4.6847847078799"+000 +8.2765898338392"+000 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 10 SUBSECTION: START. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "INTEGER" "PROCEDURE" START(X,N,T); "VALUE" X,N,T; "REAL" X; "INTEGER" N,T; "CODE" 35185; START:= A STARTING VALUE FOR THE MILLER ALGORITHM FOR COMPUTING AN ARRAY OF BESSEL FUNCTIONS; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTIONS, X > 0; N: ; THE NUMBER OF BESSEL FUNCTIONS TO BE COMPUTED, N >= 0; T: ; THE TYPE OF BESSEL FUNCTION IN QUESTION, T = 0 CORRESPONDS TO ORDINARY BESSEL FUNCTIONS, T = 1 CORRESPONDS TO MODIFIED BESSEL FUNCTIONS. PROCEDURES USED: NONE. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. METHOD AND PERFORMANCE: THE PROCEDURE IS CALLED IN THE FOLLOWING PROCEDURES: BESS J CODE 35162 NON EXP BESS I CODE 35177 BESS JAPLUSN CODE 35180 BESS KAPLUSN CODE 35192 NON EXP BESS IAPLUSN CODE 35193 SPHER BESS J CODE 35150 NON EXP SPHER BESS I CODE 35154. 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 11 IN THESE PROCEDURES AN ARRAY OF BESSEL FUNCTIONS IS GENERATED BY USING MILLER 'S ALGORITHM (SEE REF[5]). FOR STARTING THIS ALGORITHM ONE NEEDS AN INTEGER NU WHICH CAN BE COMPUTED BY USING GAUTSCHI 'S ESTIMATES OF THE ERROR ( SEE REF[5,FORMULA (5.11)] ). WE COMPUTE THIS STARTING VALUE NU BY USING ASYMPTOTIC APPROXIMA- TIONS OF THE BESSEL FUNCTIONS, AS GIVEN IN REF[1, FORMULA 9.3.7, 9.3.8, 9.7.7, AND 9.7.8]. GAUTSCHI USED DIFFERENT FORMULAS, BUT THOSE USED HERE GIVE FOR LARGE X AND N MORE REALISTIC ESTIMATES. THE PERSUED ACCURACY IN THE ABOVE MENTIONED PROCEDURES IS ABOUT "-14 . FOR OBTAINING AN ACCURACY OF "-D THE NUMBERS 36 AND 18 APPEARING IN THE FOURTH AND SIXTH LINE OF THE SOURCE TEXT OF START SHOULD BE REPLACED BY (D+1)* LN(10) AND .5*(D+1)* LN(10), RESPECTIVELY. FOR MODIFIED BESSEL FUNCTIONS THE ACCURRACY IS IN A RELATIVE SENSE; FOR ORDINARY BESSEL FUNCTIONS THE ACCURRACY IS ABSOLUTE IF THE ORDER OF THE BESSEL FUNCTION IS SMALLER THAN X, OTHERWISE IT IS RELATIVE. RUNNING TIME: NEGLECTABLE IF COMPARED WITH THE TIME NEEDED FOR THE BESSEL FUNCTION PROCEDURES. EXAMPLE OF USE: SEE THE ABOVE MENTIONED PROCEDURES. SOURCE TEXT(S): "CODE" 35180; "PROCEDURE" BESS JAPLUSN(A, X, N, JA); "VALUE" A, X, N; "INTEGER" N; "REAL" X, A; "ARRAY" JA; "IF" X = 0 "THEN" "BEGIN" JA[0]:= "IF" A = 0 "THEN" 1 "ELSE" 0; "FOR" N:= N "STEP" -1 "UNTIL" 1 "DO" JA[N]:= 0 "END" "ELSE" "IF" A = 0 "THEN" "BEGIN" "PROCEDURE" BESS J(X, N, J); "CODE" 35162; BESS J(X, N, JA) "END" "ELSE" "IF" A = .5 "THEN" "BEGIN" "REAL" S; "PROCEDURE" SPHER BESS J(X, N, J); "CODE" 35150; S:= SQRT(X) * .797 884 560 802 865; "COMMENT" S = SQRT(2X / PI); SPHER BESS J(X, N, JA); "FOR" N:= N "STEP" - 1 "UNTIL" 0 "DO" JA[N]:= JA[N] * S "END" "ELSE" "COMMENT" 1SECTION: 6.10.1 (DECEMBER 1978) PAGE 12 ; "BEGIN" "REAL" A2, X2, R, S, L, LABDA; "INTEGER" K, M, NU; "REAL" "PROCEDURE" GAMMA(Y); "CODE" 35061; "INTEGER" "PROCEDURE" START(X,N,T); "CODE" 35185; L:= 1; NU:= START(X,N,0); "FOR" M:= 1 "STEP" 1 "UNTIL" NU "DO" L:= L * (M+A) / (M+1); R:= S:= 0; X2:= 2 / X; K:= -1; A2:= A + A; "FOR" M:= NU+NU "STEP" - 1 "UNTIL" 1 "DO" "BEGIN" R:= 1 / (X2 * (A + M) - R); "IF" K = 1 "THEN" LABDA:= 0 "ELSE" "BEGIN" L:= L * (M + 2) / (M + A2); LABDA:= L * (M + A) "END"; S:= R * (LABDA + S); K:= -K; "IF" M<= N "THEN" JA[M]:= R "END"; JA[0]:= R:= 1 / GAMMA(1 + A) / (1 + S) / X2 ** A; "FOR" M:= 1 "STEP" 1 "UNTIL" N "DO" JA[M]:= R:= R * JA[M]; "END" BESS JAPLUSN; "EOP" "CODE" 35181; "PROCEDURE" BESS YA01(A,X,YA,YA1);"VALUE" A,X; "REAL" A,X,YA,YA1; "IF" A = 0 "THEN" "BEGIN" "PROCEDURE" BESS Y01(X,Y0,Y1); "CODE" 35163; BESS Y01(X,YA,YA1) "END" "ELSE" "BEGIN" "REAL" B,C,D,E,F,G,H,P,PI,Q,R,S;"INTEGER" N,NA; "BOOLEAN" REC,REV; PI:=4*ARCTAN(1);NA:=ENTIER(A+.5);REC:=A>=.5; REV:=A<-.5;"IF" REV "OR" REC "THEN" A:=A-NA; "IF" A=-.5 "THEN" "BEGIN" P:=SQRT(2/PI/X);F:=P*SIN(X);G:=-P*COS(X) "END" "ELSE" "IF" X<3 "THEN" "BEGIN" "REAL" "PROCEDURE" RECIP GAMMA(X,ODD,EVEN); "CODE" 35060; "REAL" "PROCEDURE" SINH(X); "CODE" 35111; B:=X/2;D:=-LN(B);E:=A*D; C:="IF" ABS(A)<"-8 "THEN" 1/PI "ELSE" A/SIN(A*PI); S:="IF" ABS(E)<"-8 "THEN" 1 "ELSE" SINH(E)/E; E:=EXP(E);G:=RECIP GAMMA(A,P,Q)*E;E:=(E+1/E)/2; F:=2*C*(P*E+Q*S*D);E:=A*A; P:=G*C;Q:=1/G/PI;C:=A*PI/2; R:="IF" ABS(C)<"-8 "THEN" 1 "ELSE" SIN(C)/C;R:=PI*C*R*R; C:=1;D:=-B*B;YA:=F+R*Q;YA1:=P; "FOR" N:=1,N+1 "WHILE" ABS(G/(1+ABS(YA)))+ABS(H/(1+ABS(YA1)))>"-15 "DO" "BEGIN" F:=(F*N+P+Q)/(N*N-E);C:=C*D/N; P:=P/(N-A);Q:=Q/(N+A); G:=C*(F+R*Q);H:=C*P-N*G; YA:=YA+G;YA1:=YA1+H; "END"; F:=-YA;G:=-YA1/B "END" "ELSE" "COMMENT" 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 13 ; "BEGIN" "PROCEDURE" BESS PQA01(A,X,PA,QA,PA1,QA1); "CODE" 35183; B:=X-PI*(A+.5)/2;C:=COS(B);S:=SIN(B); D:=SQRT(2/X/PI); BESS PQA01(A,X,P,Q,B,H); F:=D*(P*S+Q*C);G:=D*(H*S-B*C) "END"; "IF" REV "THEN" "BEGIN" X:=2/X;NA:=-NA-1; "FOR" N:=0 "STEP" 1 "UNTIL" NA "DO" "BEGIN" H:=X*(A-N)*F-G;G:=F;F:=H "END" "END" "ELSE" "IF" REC "THEN" "BEGIN" X:=2/X; "FOR" N:=1 "STEP" 1 "UNTIL" NA "DO" "BEGIN" H:=X*(A+N)*G-F;F:=G;G:=H "END" "END"; YA:=F;YA1:=G "END" BESS YA01; "EOP" "CODE" 35182; "PROCEDURE" BESS YAPLUSN(A, X, NMAX, YAN); "VALUE" A, X, NMAX; "REAL" A, X; "INTEGER" NMAX; "ARRAY" YAN; "BEGIN" "INTEGER" N; "REAL" Y1; "PROCEDURE" BESS YA01(A, X, YA, YA1); "CODE" 35181; BESS YA01(A, X, YAN[0], Y1); A:= A-1; X:= 2/X; "IF" NMAX > 0 "THEN" YAN[1]:= Y1; "FOR" N:= 2 "STEP" 1 "UNTIL" NMAX "DO" YAN[N]:= -YAN[N-2] + (A+N)*X*YAN[N-1] "END" BESS YAPLUSN; "EOP" "CODE" 35183; "PROCEDURE" BESS PQA01(A,X,PA,QA,PA1,QA1);"VALUE" A,X; "REAL" A,X,PA,PA1,QA,QA1; "IF" A = 0 "THEN" "BEGIN" "PROCEDURE" BESS PQ0(X,P0,Q0); "CODE" 35165; "PROCEDURE" BESS PQ1(X,P1,Q1); "CODE" 35166; BESS PQ0(X,PA,QA); BESS PQ1(X,PA1,QA1) "END" "ELSE" "BEGIN" "INTEGER" N,NA; "REAL" B, PI, P0, Q0 ; "BOOLEAN" REC, REV; PI:= 4 * ARCTAN(1); REV:=A<-.5;"IF" REV "THEN" A:=-A-1; REC:=A>=.5;"IF" REC "THEN" "BEGIN" NA:=ENTIER(A+.5);A:=A-NA "END"; "IF" A=-.5 "THEN" "BEGIN" PA:=PA1:=1;QA:=QA1:=0 "END" "ELSE" "IF" X >= 3 "THEN" "COMMENT" 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 14 ; "BEGIN" "REAL" C,D,E,F,G,H,P,Q,R,S; C:=.25 - A*A; B:= X + X; F:= R:= 1; G:= -X; S:= 0; E:=(X*COS(A*PI)/PI*"15)**2; "FOR" N:=2,N+1 "WHILE" (P*P + Q*Q)*N*N0 "DO" "BEGIN" R:=(N+1)*(2-P)-2;S:=B+(N+1)*Q;D:=(N-1+C/N)/ (R*R+S*S);P:=D*R;Q:=D*S;E:=F; F:=P*(E+1)-G*Q;G:=Q*(E+1)+P*G "END"; F:=1+F; D:=F*F + G*G; PA:=F/D;QA:=-G/D;D:=A+.5-P;Q:=Q+X; PA1:=(PA*Q-QA*D)/X; QA1:=(QA*Q+PA*D)/X "END" "ELSE" "BEGIN" "REAL" C, S, CHI, YA, YA1; "ARRAY" JA[0:1]; "PROCEDURE" BESS JAPLUSN(A, X, N, JA); "CODE" 35180; "PROCEDURE" BESS YA01(A, X, YA, YA1); "CODE" 35181; B:= SQRT(PI * X / 2); CHI:= X - PI * (A / 2 + .25); C:= COS(CHI); S:= SIN(CHI); BESS YA01(A, X, YA, YA1); BESS JAPLUSN(A, X, 1, JA); PA:= B * (YA * S + C * JA[0]); QA:= B * (C * YA - S * JA[0]); PA1:= B * (S * JA[1] - C * YA1); QA1:= B * (JA[1] * C + YA1 * S) "END"; "IF" REC "THEN" "BEGIN" X:=2/X;B:=(A+1)*X; "FOR" N:=1 "STEP" 1 "UNTIL" NA "DO" "BEGIN" P0:=PA-QA1*B; Q0:=QA+PA1*B; PA:=PA1;PA1:=P0; QA:=QA1; QA1:=Q0; B:=B+X "END" "END"; "IF" REV "THEN" "BEGIN" P0:=PA1;PA1:=PA;PA:=P0;Q0:=QA1;QA1:=QA;QA:=Q0 "END" "END" BESS PQA01; "EOP" 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 15 "CODE" 35184; "PROCEDURE" BESS ZEROS(A,N,Z,D); "VALUE" A,N,D; "REAL" A;"ARRAY" Z; "INTEGER" N,D; "COMMENT" COMPUTES Z[1],...Z[N],THE FIRST N ZEROS OF A BESSEL FUNCTION. THE CHOICE OF D DETERMINES THE TYPE OF THE BESSEL FUNCTION : IF D=1 THEN JA ELSE IF D=2 THEN YA ELSE IF D=3 THEN JA-PRIME ELSE IF D=4 THEN YA-PRIME. A IS THE ORDER OF THE BESSEL FUNCTION, IT MUST BE NON-NEGATIVE.; "BEGIN" "REAL" AA,A2,B,BB,C,CHI,CO,MU,MU2,MU3,MU4,P,PI,PA,PA1,P0,P1,PP1, Q,QA,QA1,Q1,QQ1,RO,SI,T,TT,U,V,W,X,XX,X4,Y; "INTEGER" J,S; "REAL" "PROCEDURE" FI(Y); "VALUE" Y; "REAL" Y; "COMMENT" COMPUTES FI FROM THE EQUATION TAN(FI)-FI=Y , WHERE Y>=0. THE RELATIVE ACCURACY IS AT LEAST 5 DIGITS; "IF" Y=0 "THEN" FI:=0 "ELSE" "IF" Y>"5 "THEN" FI:=1.570796 "ELSE" "BEGIN" "REAL" R,P,PP; "IF" Y<1 "THEN" "BEGIN" P:=(3*Y)**(1/3); PP:=P*P; P:=P*(1+PP*(-210+PP*(27-2*PP))/1575) "END" "ELSE" "BEGIN" P:=1/(Y+1.570796); PP:=P*P; P:= 1.570796-P*(1+PP*(2310+PP*(3003+PP*(4818+PP* (8591+PP*16328))))/3465) "END"; PP:=(Y+P)*(Y+P); R:=(P-ARCTAN(P+Y))/PP; FI:=P-(1+PP)*R*(1+R/(P+Y)) "END" FI; "PROCEDURE" BESS PQA01(A,X,PA,QA,PA1,QA1); "CODE" 35183; "REAL" "PROCEDURE" R; "BEGIN" BESS PQA01(A,X,PA,QA,PA1,QA1); CHI:=X-PI*(A/2+0.25); SI :=SIN(CHI); CO:=COS(CHI); R:= "IF" D=1 "THEN" (PA*CO-QA*SI)/(PA1*SI+QA1*CO) "ELSE" "IF" D=2 "THEN" (PA*SI+QA*CO)/(QA1*SI-PA1*CO) "ELSE" "IF" D=3 "THEN" A/X-(PA1*SI+QA1*CO)/(PA*CO-QA*SI) "ELSE" A/X-(QA1*SI-PA1*CO)/(PA*SI+QA*CO) "END" R; "COMMENT" 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 16 ; PI:=4*ARCTAN(1); AA:=A*A; MU:=4*AA; MU2:=MU*MU; MU3:=MU*MU2; MU4:=MU2*MU2; "IF" D<3 "THEN" "BEGIN" P:=7*MU-31; P0:=MU-1; P1:=4*(253*MU2-3722*MU+17869)/15/P*P0; Q1:=8*( 83*MU2- 982*MU+ 3779)/ 5/P "END" "ELSE" "BEGIN" P:=7*MU2+82*MU-9; P0:=MU+3; P1:=(4048*MU4+131264*MU3-221984*MU2-417600*MU+1012176)/60/P; Q1:=1.6*(83*MU3+2075*MU2-3039*MU+3537)/P "END"; T:="IF" D=1"OR"D=4 "THEN" 0.25 "ELSE" 0.75; TT:=4*T; "IF" D<3 "THEN" "BEGIN" PP1:= 5/48; QQ1:= -5/36 "END" "ELSE" "BEGIN" PP1:=-7/48; QQ1:= 35/288 "END"; Y:= 3*PI/8; BB:= "IF" A>=3 "THEN" A **(-2/3) "ELSE" 0.0 ; "FOR" S:=1 "STEP" 1 "UNTIL" N "DO" "BEGIN" "IF" A=0"AND"S=1"AND"D=3 "THEN" "BEGIN" X:=0; J:=0 "END" "ELSE" "BEGIN" "IF" S >= 3*A -8 "THEN" "BEGIN" B:=(S+A/2-T)*PI; C:=1/B/B/64; X:=B-1/B/8*(P0-P1*C)/(1-Q1*C) "END" "ELSE" "BEGIN" "IF" S=1 "THEN" "BEGIN" X:= "IF" D=1 "THEN" -2.33811 "ELSE" "IF" D=2 "THEN" -1.17371 "ELSE" "IF" D=3 "THEN" -1.01879 "ELSE" -2.29444 "END" "ELSE" "BEGIN" X:= Y*(4*S-TT); V:= 1/X/X; X:= -X**(2/3)*(1+V*(PP1+QQ1*V)) "END"; U:=X*BB; V:=FI(2/3*(-U)**1.5); W:=1/COS(V); XX:=1-W*W; C:=SQRT(U/XX); X:=W*(A+C/A/U* ("IF" D<3 "THEN" -5/48/U-C*(-5/24/XX+1/8) "ELSE" 7/48/U+C*(-7/24/XX+3/8))) "END"; J:=0; L1: XX:=X*X; X4:=XX*XX; A2:=AA-XX; RO:=R; J:=J+1; "IF" D<3 "THEN" "BEGIN" U:=RO; P:=(1-4*A2)/6/X/(2*A+1); Q:=(2*(XX-MU)-1-6*A)/3/X/(2*A+1) "END" "ELSE" "BEGIN" U:=-XX*RO/A2; V:=2*X*A2/(AA+XX)/3; W:=A2*A2*A2; Q:=V*(1+( MU2+32*MU*XX+48*X4)/32/W); P:=V*(1+(-MU2+40*MU*XX+48*X4)/64/W) "END"; W:=U*(1+P*RO)/(1+Q*RO); X:=X+W; "IF" ABS(W/X)>"-13 "AND"J<5 "THEN" "GOTO" L1 "END"; Z[S]:=X "END" "END" BESS ZEROS; "EOP" 1SECTION : 6.10.1 (DECEMBER 1978) PAGE 17 "CODE" 35185; "INTEGER" "PROCEDURE" START(X,N,T); "VALUE" X,N,T; "REAL" X; "INTEGER" N,T; "BEGIN" "REAL"P,Q,R,Y; "INTEGER" S; S:= 2*T-1; P:= 36/X-T; R:= N/X; "IF" R>1 "OR" T=1 "THEN" "BEGIN" Q:= SQRT(R*R+S); R:= R*LN(Q+R)-Q "END" "ELSE" R:= 0; Q:= 18/X+R; R:= "IF" P>Q "THEN" P "ELSE" Q; P:= SQRT(2*(T+R)); P:= X*((1+R)+P)/(1+P); Y:= 0; "FOR" Q:= Y, Y "WHILE" P>Q "OR" P = 0. NONEXP BESS IAPLUSN: THIS PROCEDURE GENERATES AN ARRAY OF MODIFIED BESSEL FUNCTIONS OF THE FIRST KIND OF ORDER A + N, N = 0, ..., NMAX, 0<=A <1 AND ARGUMENT X > = 0 MULTIPLIED BY THE FACTOR EXP(-X). THUS, APART FROM THE EXPONENTIAL FACTOR THE ARRAY ENTRIES ARE THE SAME AS THOSE COMPUTED BY BESS IAPLUSN. BESS KA01: THIS PROCEDURE CALCULATES THE MODIFIED BESSEL FUNCTIONS OF THE THIRD KIND OF ORDER A AND A+1, AND ARGUMENT X, X>0; BESS KAPLUSN: THIS PROCEDURE GENERATES AN ARRAY OF MODIFIED BESSEL FUNCTIONS OF THE THIRD KIND OF ORDER A+N, N=0, 1, ..., NMAX, AND ARGUMENT X>0. THE MODIFIED BESSEL FUNCTIONS CORRESPOND TO THE FUNCTION DEFINED IN FORMULA 9.6.2 OF REFERENCE[1]; NONEXP BESS KA01: THIS PROCEDURE CALCULATES THE MODIFIED BESSEL FUNCTIONS OF THE THIRD KIND OF ORDER A AND A + 1, AND ARGUMENT X, X > 0, MULTIPLIED BY THE FACTOR EXP(X). THUS, APART FROM THE EXPONENTIAL FACTOR, THE FUNCTIONS ARE THE SAME AS THOSE COMPUTED BY BESS KA01; 1SECTION : 6.10.2 (DECEMBER 1978) PAGE 2 NONEXP BESS KAPLUSN: THIS PROCEDURE GENERATES AN ARRAY OF MODIFIED BESSEL FUNCTIONS OF THE THIRD KIND OF ORDER A + N, N = 0, 1,..., NMAX, AND ARGUMENT X>0 MULTIPLIED BY THE FACTOR EXP(X). THUS, APART FROM THE EXPONENTIAL FACTOR, THE FUNCTIONS ARE THE SAME AS THOSE COMPUTED BY THE PROCEDURE BESS KAPLUSN. KEYWORDS: BESSEL FUNCTION, MODIFIED BESSEL FUNCTION, MODIFIED BESSEL FUNCTION OF THE THIRD KIND. REFERENCES: [1]. ABRAMOWITZ, M., AND STEGUN, I. (EDS.), HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICAL TABLES. APPL. MATH. SER. 55, U.S. GOVT. PRINTING OFFICE, WASHINGTON, D.C. (1964). [2]. GAUTSCHI, W., COMPUTATIONAL ASPECTS OF THREE TERM RECURRENCE RELATIONS. SIAM REVIEW, VOLUME 9, (1967), NUMBER 1, P.24. [3]. TEMME, N.M., ON THE NUMERICAL EVALUATION OF THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND. J. COMP. PHYSICS, VOL. 19, (1975), NUMBER 3, P. 324. SUBSECTION: BESS IAPLUSN. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS IAPLUSN(A, X, N, IA); "VALUE" X, N, A; "REAL" X, A; "INTEGER" N; "ARRAY" IA; "CODE" 35190; 1SCETION : 6.10.2 (DECEMBER 1979) PAGE 3 THE MEANING OF THE FORMAL PARAMETERS IS: A: < ARITHMETIC EXPRESSION >; THE NONINTEGER PART OF THE ORDER OF THE BESSEL FUNCTIONS; 0 < = A < 1; X: < ARITHMETIC EXPRESSION >; THE ARGUMENT OF THE BESSEL FUNCTIONS; X > = 0; N: < ARITHMETIC EXPRESSION >; THE UPPER BOUND OF THE INDICES OF THE ARRAY IA; N>= 0; IA: < ARRAY IDENTIFIER >; "ARRAY" IA[0:N]; N > = 0; EXIT: THE VALUES OF THE MODIFIED BESSEL FUNCTIONS OF THE FIRST KIND , OF ORDER A+K AND ARGUMENT X, I[A+K](X) ARE ASSIGNED TO THE ARRAY IA. PROCEDURES USED: NONEXP BESS IAPLUSN = CP 35193, BESS I = CP 35172, NONEXP SPHER BESS I = CP 35154. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. METHOD AND PERFORMANCE: SEE SUBSECTION NONEXP BESS IAPLUSN. RUNNING TIME: ROUGHLY PROPORTIONAL TO THE MAXIMUM OF X AND N. EXAMPLE OF USE: "BEGIN" "REAL" X, A; "ARRAY" IA[0:2] ; "PROCEDURE" BESS IAPLUSN(A, X, N, IA); "CODE" 35190; A:= .25; X:= 2; BESS IAPLUSN(A, X, 2, IA); OUTPUT(61,"("2(4BD.DD),/,3(4B-.14D"-ZD)")", A, X, IA[0], IA[1], IA[2]) "END" PRINTS THE FOLLOWING RESULTS: 0.25 2.00 .22033544516736" 1 .13401967589829" 1 .52810850294501" 0 1SECTION : 6.10.2 (DECEMBER 1978) PAGE 4 SUBSECTION: NONEXP BESS IAPLUSN. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" NONEXP BESS IAPLUSN(A, X, N, IA); "VALUE" A, X, N; "REAL" A, X; "INTEGER" N; "ARRAY" IA; "CODE" 35193; THE MEANING OF THE FORMAL PARAMETERS IS: A: < ARITHMETIC EXPRESSION >; THE NONINTEGER PART OF THE ORDER A+N, 0 < = A < 1; X: < ARITHMETIC EXPRESSION >; THE ARGUMENT OF THE BESSEL FUNCTIONS; X >= 0; N: < ARITHMETIC EXPRESSION >; THE UPPER BOUND OF THE INDICES OF THE ARRAY IA; N>= 0; IA: < ARRAY IDENTIFIER >; "ARRAY" IA[0:N]; N > = 0; EXIT: IA[K] HAS THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE FIRST KIND OF ORDER A + K AND ARGUMENT X MULTIPLIED BY EXP (-X), 0 < = K < = N. PROCEDURES USED: NONEXP BESS I = CP 35177 NONEXP SPHER BESS I = CP 35154 GAMMA = CP 35061 START = CP 35185 REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. METHOD AND PERFORMANCE: IN ALL THE CASES THE BESSEL FUNCTIONS ARE COMPUTED ACCORDING TO THE MILLER METHOD DESCRIBED IN [2, P.46-52]. THE STARTING VALUE IS COMPUTED BY THE PROCEDURE START (SECTION 6.10.1). RUNNING TIME: ROUGHLY PROPORTIONAL TO THE MAXIMUM OF X AND N. 1SECTION : 6.10.2 (DECEMBER 1979) PAGE 5 EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X, A; "ARRAY" IA[0:2]; "PROCEDURE" NON EXPBESS IAPLUSN(A, X, N, IA); "CODE" 35193; A:= .25; X:= 2; NON EXPBESS IAPLUSN(A, X, 2, IA); OUTPUT(61,"("2(4BD.DD),/,3(4B-.14D"-ZD)")", A, X, IA[0], IA[1], IA[2]) "END" PRINTS THE FOLLOWING RESULTS: 0.25 2.00 .29819159878790" 0 .18137590796974" 0 .71471713825726" -1 SUBSECTION: BESS KA01. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS KA01(A, X, KA, KA1); "VALUE" A, X; "REAL" A, X, KA, KA1; "CODE" 35191; THE MEANING OF THE FORMAL PARAMETERS IS: A: ; THE ORDER; X: ; THE ARGUMENT. THIS ARGUMENT SHOULD SATISFY X>0; KA: ; EXIT: THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER A AND ARGUMENT X; KA1: ; EXIT: THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER A+1 AND ARGUMENT X. PROCEDURES USED: RECIP GAMMA = CP 35060; NONEXP BESS KA01 = CP 35194; SINH = CP 35111. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRYS ARE DECLARED. 1SECTION : 6.10.2 (DECEMBER 1978) PAGE 6 METHOD AND PERFORMANCE: FOR 0=1 THE PROCEDURE CALLS FOR THE PROCEDURE NONEXP BESS KA ( SEE SUBSECTION NONEXP BESS KA). THE RELATIVE ACCURACY IS ABOUT "-13, EXCEPT FOR LARGE VALUES OF X; IN THAT CASE THE ACCURACY ALSO DEPENDS ON THE RELATIVE ACCURACY OF THE EXPONENTIAL FUNCTION. IF ONE IS INTERESTED IN THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND TIMES THE FACTOR EXP(X), THE PROCEDURE NONEXP BESS KA SHOULD BE USED. EXAMPLE OF USE: THE PROGRAM: "BEGIN" "REAL" P, Q; "PROCEDURE" BESS KA01(A, X, KA, KA1); "CODE" 35191; BESS KA01(0, 1, P, Q); OUTPUT(61, "("2(N)")", P, Q) "END" YIELDS THE FOLLOWING RESULTS +4.2102443824071"-001 +6.0190723019724"-001. SUBSECTION: BESS KAPLUSN. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS KAPLUSN(A, X, NMAX, KAN); "VALUE" A, X, NMAX; "INTEGER" NMAX; "REAL" A, X; "ARRAY" KAN; "CODE" 35192; THE MEANING OF THE FORMAL PARAMETERS IS: A: ; THE ORDER. IT IS ADVISED TO TAKE A >= 0; X: ; THE ARGUMENT. THIS ARGUMENT SHOULD SATISFY X>0; NMAX: ; THE UPPER BOUND OF THE INDICES OF THE ARRAY KAN; KAN: ; "ARRAY" KAN[0:NMAX]; NMAX>=0; EXIT: THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER N+A IS ASSIGNED TO KAN[N], 0 <= N <= NMAX. 1SECTION : 6.10.2 (DECEMBER 1978) PAGE 7 PROCEDURES USED: BESS KA01 = CP 35191. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. METHOD AND PERFORMANCE: THE RECURRENCE RELATION KAN[N+1]=KAN[N-1]+2*(N+A)*KAN[N]/X IS USED. THE STARTING VALUES ARE OBTAINED FROM THE PROCEDURE BESS KA01. IF A>=0, RECURSION IS NUMERICALLY STABLE IN THE FORWARD DIRECTION. IF ONE IS INTERESTED IN THE MODIFIED BESSEL FUNCTIONS OF THE THIRD KIND TIMES THE FACTOR EXP(X), THE PROCEDURE NONEXP BESS KAPLUSN SHOULD BE USED. EXAMPLE OF USE: THE PROGRAM: "BEGIN" "ARRAY" KAN[0:2]; "PROCEDURE" BESS KAPLUSN(A, X, NMAX, KAN); "CODE" 35193; BESS KAPLUSN(0, 1, 2, KAN); OUTPUT(61, "("3(N)")", KAN[0], KAN[1], KAN[2]) "END" YIELDS THE FOLLOWING RESULTS +4.2102443824071"-001 +6.0190723019724"-001 +1.6248388986352"+000. SUBSECTION: NONEXP BESS KA01. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" NONEXP BESS KA01(A, X, KA, KA1); "VALUE" A, X; "REAL" A, X, KA, KA1; "CODE" 35194; 1SECTION : 6.10.2 (DECEMBER 1978) PAGE 8 THE MEANING OF THE FORMAL PARAMETERS IS: A: ; THE ORDER; X: ; THE ARGUMENT. THIS ARGUMENT SHOULD SATISFY X>0; KA: ; EXIT: KA HAS THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER A AND ARGUMENT X TIMES THE FACTOR EXP(X); KA1: ; EXIT: THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER A+1 AND ARGUMENT X TIMES THE FACTOR EXP(X). PROCEDURES USED: BESS KA01 = CP 35191. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. METHOD AND PERFORMANCE: FOR 0=1 THE BESSEL FUNCTIONS ARE COMPUTED WITH A MILLER ALGORITHM FOR CONFLUENT HYPERGEOMETRIC FUNCTIONS. THE METHOD IS DESCRIBED IN REFERENCE [3]. FOR ALL VALUES OF X CONSIDERED (X>0) THE FUNCTIONS DELIVERED ARE EQUAL TO THE VALUES COMPUTED BY THE PROCEDURE BESS KA01, APART FROM AN EXPONENTIAL FACTOR. THE RELATION BETWEEN THE TWO PROCEDURES WILL BE DESCRIBED BY THE PROGRAM: "BEGIN" "REAL" A, X, KA, NEKA, KA1, NEKA1; "PROCEDURE" BESS KA01(A, X, KA, KA1); "CODE" 35191; "PROCEDURE" NONEXP BESS KA(A, X, KA, KA1); "CODE" 35194; A:= .3; X:= 3.14; BESS KA01(A, X, KA, KA1); NONEXP KA 01(A, X, NEKA, NEKA1) "END" THEN WE HAVE KA = EXP(-X)*NEKA, KA1 = EXP(-X)*NEKA1. THE RELATIVE ACCURACY IS ABOUT "-13. 1SECTION : 6.10.2 (DECEMBER 1978) PAGE 9 EXAMPLE OF USE: THE PROGRAM: "BEGIN" "REAL" P, Q; "PROCEDURE" NONEXP BESS KA(A, X, KA, KA1); "CODE" 35194; NONEXP BESS KA(0, 2, P, Q); OUTPUT(61, "("2(N)")", P, Q) "END" YIELDS THE FOLLOWING RESULTS: 8.4156821507078"-001 +1.0334768470687"+000. SUBSECTION: NONEXP BESS KAPLUSN. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" NONEXP BESS KAPLUSN(A, X, NMAX, KAN); "VALUE" A, X, NMAX; "REAL" A, X; "INTEGER" NMAX; "ARRAY" KAN; "CODE" 35195; NONEXP BESS KAPLUSN GENERATES AN ARRAY OF MODIFIED BESSEL FUNCTIONS THE THIRD KIND OF ARGUMENT X AND ORDERS A+N, N=0, 1,..., NMAX TIMES THE FACTOR EXP(X). THE MEANING OF THE FORMAL PARAMETERS IS: A: ; THE ORDER. IT IS ADVISED TO TAKE A >= 0; X: ; THE ARGUMENT. THIS ARGUMENT SHOULD SATISFY X>0; NMAX: ; THIS PARAMETER SHOULD SATISFY NMAX>=0; NMAX INDICATES THE MAXIMUM NUMBER OF FUNCTION VALUES TO BE GENERATED; KAN: ; "ARRAY" KAN[0:NMAX]; NMAX>=0; EXIT: KAN[N] IS THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER N+A AND OF ARGUMENT X (N=0(1)NMAX) TIMES THE FACTOR EXP(X). PROCEDURES USED: NONEXP BESS KA = CP 35194. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. 1SECTION : 6.10.2 (DECEMBER 1978) PAGE 10 METHOD AND PERFORMANCE: THE RECURRENCE RELATION KAN[N+1]=KAN[N-1]+2*(N+A)*KAN[N]/X IS USED. THE STARTING VALUES ARE OBTAINED FROM THE PROCEDURE NONEXP BESS KA. IF A>=0, RECURSION IS NUMERICALLY STABLE IN THE FORWARD DIRECTION. FOR ALL VALUES OF X AND NMAX CONSIDERED (X>0) THE FUNCTIONS DELIVERED ARE EQUAL TO THE VALUES COMPUTED BY THE PROCEDURE BESS KAPLUSN,APART FROM AN EXPONENTIAL FACTOR. THE RELATION BETWEEN THE TWO PROCEDURES WILL BE DESCRIBED BY THE PROGRAM: "BEGIN" "REAL" X, A; "ARRAY" KA, NEKA[0:10]; "PROCEDURE" BESS KAPLUSN(A, X, NMAX, KA); "CODE" 35193; "PROCEDURE" NONEXP BESS KAPLUSN(A, X, NMAX, KAN); "CODE" 35195; X:= 2.78; A:= .96; BESS KAPLUSN(A, X, 10, KA); NONEXP BESS KAPLUSN(A, X, 10, NEKA) "END" THEN WE HAVE KA[N] = EXP(-X)*NEKA[N], N=0, 1, ..., 10. EXAMPLE OF USE: THE PROGRAM: "BEGIN" "ARRAY" KAN[0:2]; "PROCEDURE" NONEXP BESS KAPLUSN(A, X, NMAX, KAN); "CODE" 35195; NONEXP BESS KAPLUSN(0, 5, 2, KAN); OUTPUT(61, "("3(N)")", KAN[0], KAN[1], KAN[2]) "END" YIELDS THE FOLLOWING RESULTS: +5.4780756431352"-001 +6.0027385878831"-001 +7.8791710782884"-001. 1SECTION : 6.10.2 (DECEMBER 1979) PAGE 11 SOURCE TEXT(S): "CODE" 35190; "COMMENT" COMPUTATION OF I[A](X), , I[N+A](X); "PROCEDURE" BESS IAPLUSN(A, X, N, IA); "VALUE" A, X, N; "INTEGER" N; "REAL" X, A; "ARRAY" IA; "IF" X= 0 "THEN" "BEGIN" IA[0]:= "IF" A= 0 "THEN" 1 "ELSE" 0; "FOR" N:= N "STEP" -1 "UNTIL" 1 "DO" IA[N]:= 0 "END" "ELSE" "IF" A= 0 "THEN" "BEGIN" "PROCEDURE" BESS I(X, N, I); "CODE" 35172; BESS I(X, N, IA); "END" "ELSE" "IF" A= .5 "THEN" "BEGIN" "REAL" C; "PROCEDURE" NONEXP SPHER BESSI(X, N, I); "CODE" 35154; C:= .797 884 560 802 865 * SQRT(ABS(X)) * EXP (ABS (X)); NONEXP SPHER BESSI(X, N, IA); "FOR" N:= N "STEP" -1 "UNTIL" 0 "DO" IA[N]:= C*IA[N] "END" "ELSE" "BEGIN" "REAL" EXPX; "PROCEDURE" NONEXP BESS IAPLUSN(A, X, N, IA); "CODE" 35193; EXPX:= EXP(ABS(X)); NONEXP BESS IAPLUSN(A, X, N, IA); "FOR" N:= N "STEP" -1 "UNTIL" 0 "DO" IA[N]:= EXPX * IA[N] "END" BESS IAPLUSN; "EOP" "CODE" 35191; "PROCEDURE" BESS KA01(A, X, KA, KA1); "VALUE" A, X; "REAL" A, X, KA, KA1; "IF" A = 0 "THEN" "BEGIN" "PROCEDURE" BESS K01(X,K0,K1); "CODE" 35173; BESS K01(X,KA,KA1) "END" "ELSE" "BEGIN" "REAL" F, G, H, PI; "INTEGER" N, NA; "BOOLEAN" REC, REV; PI:= 4 * ARCTAN(1); REV:= A < -.5; "IF" REV "THEN" A:= -A-1; REC:= A >= .5; "IF" REC "THEN" "BEGIN" NA:=ENTIER(A+.5); A:= A - NA "END"; "IF" A = .5 "THEN" F:= G:= SQRT(PI / X / 2) * EXP (-X) "ELSE" "IF" X < 1 "THEN" "BEGIN" "COMMENT" 1SECTION : 6.10.2 (DECEMBER 1978) PAGE 12 ; "REAL" A1, B, C, D, E, P, Q, S; "REAL" "PROCEDURE" RECIP GAMMA(X, ODD, EVEN); "CODE" 35060; "REAL" "PROCEDURE" SINH(X); "CODE" 35111; B:=X/2;D:=-LN(B);E:=A*D;C:=A*PI; C:="IF" ABS(C)<"-15 "THEN" 1 "ELSE" C/SIN(C); S:="IF" ABS(E)<"-15 "THEN" 1 "ELSE" SINH(E)/E; E:=EXP(E);A1:=(E+1/E)/2;G:=RECIP GAMMA(A,P,Q)*E; KA:=F:=C*(P*A1+Q*S*D);E:=A*A; P:=.5*G*C;Q:=.5/G;C:=1;D:=B*B;KA1:=P; "FOR" N:=1,N+1 "WHILE" H/KA+ABS(G)/KA1>"-15 "DO" "BEGIN" F:=(F*N+P+Q)/(N*N-E);C:=C*D/N; P:=P/(N-A);Q:=Q/(N+A);G:=C*(P-N*F); H:=C*F;KA:=KA+H;KA1:=KA1+G "END"; F:=KA;G:=KA1/B "END" "ELSE" "BEGIN" "REAL" EXPON; "PROCEDURE" NONEXP BESS KA01(A, X, KA, KA1); "CODE" 35194; EXPON:= EXP(-X); NONEXP BESS KA01(A, X, KA, KA1); F:= EXPON * KA; G:= EXPON * KA1 "END"; "IF" REC "THEN" "BEGIN" X:= 2 / X; "FOR" N:= 1 "STEP" 1 "UNTIL" NA "DO" "BEGIN" H:= F + (A + N) * X * G; F:= G; G:= H "END" "END"; "IF" REV "THEN" "BEGIN" KA1:= F; KA:= G "END" "ELSE" "BEGIN" KA:= F; KA1:= G "END" "END" BESS KA01; "EOP" "CODE" 35192; "PROCEDURE" BESS KAPLUSN(A, X, NMAX, KAN); "VALUE" A, X, NMAX; "REAL" A, X; "INTEGER" NMAX; "ARRAY" KAN; "BEGIN" "INTEGER" N; "REAL" K1; "PROCEDURE" BESS KA01(A, X, KA, KA1);"CODE" 35191; BESS KA01(A, X, KAN[0], K1); A:= A-1; X:= 2/X; "IF" NMAX > 0 "THEN" KAN[1]:= K1; "FOR" N:= 2 "STEP" 1 "UNTIL" NMAX "DO" KAN[N]:= KAN[N-2] + (A+N) * X * KAN[N-1] "END" BESS KAPLUSN; "EOP" 1SECTION : 6.10.2 (DECEMBER 1979) PAGE 13 "CODE" 35193; "COMMENT" COMPUTATION OF NONEXPONENTIAL MODIFIED BESSEL FUNCTIONS OF FRACTIONAL ORDER; "PROCEDURE" NONEXP BESS IAPLUSN(A, X, N, IA); "VALUE" A, X, N; "REAL" X, A; "INTEGER" N; "ARRAY" IA; "IF" X= 0 "THEN" "BEGIN" IA[0]:= "IF" A= 0 "THEN" 1 "ELSE" 0; "FOR" N:= N "STEP" -1 "UNTIL" 1 "DO" IA[N]:= 0 "END" "ELSE" "IF" A= 0 "THEN" "BEGIN" "PROCEDURE" NONEXP BESSI(X, N, I); "CODE" 35177; NONEXP BESSI(X, N, IA) "END" "ELSE" "IF" A= .5 "THEN" "BEGIN" "REAL" C; "PROCEDURE" NONEXP SPHER BESSI(X, N, I); "CODE" 35154; C:= .797 884 560 802 865 * SQRT(X); NONEXP SPHER BESSI(X, N, IA); "FOR" N:= N "STEP" -1 "UNTIL" 0 "DO" IA[N]:= C * IA[N] "END" "ELSE" "BEGIN" "INTEGER" M, NU; "REAL" R, S, LABDA, L, A2, X2; "REAL" "PROCEDURE" GAMMA(X); "CODE" 35061; "INTEGER" "PROCEDURE" START(X,N,T); "CODE" 35185; A2:= A+A; X2:= 2/X; L:=1; NU:= START(X,N,1) ; R:= S:= 0; "FOR" M:= 1 "STEP" 1 "UNTIL" NU "DO" L:= L * (M+A2)/(M+1); "FOR" M:= NU "STEP" -1 "UNTIL" 1 "DO" "BEGIN" R:= 1/(X2 *(A+M)+R); L:= L*(M+1)/(M+A2); LABDA:= L*(M+A) * 2; S:= R * (LABDA + S); "IF" M <= N "THEN" IA[M]:= R "END"; IA[0]:= R:= 1/(1+S)/GAMMA(1+A)/X2 **A; "FOR" M:= 1 "STEP" 1 "UNTIL" N "DO" IA[M]:= R:= IA[M] * R; "END"; "EOP" "CODE" 35194; "PROCEDURE" NONEXP BESS KA01(A, X, KA, KA1); "VALUE" A, X; "REAL" A, X, KA, KA1; "IF" A = 0 "THEN" "BEGIN" "PROCEDURE" NONEXP BESS K01(X,K0,K1); "CODE" 35178; NONEXP BESS K01(X,KA,KA1) "END" "ELSE" "BEGIN" "REAL" F, G, H, PI; "INTEGER" N, NA; "BOOLEAN" REC, REV; PI:= 4 * ARCTAN(1); REV:= A < -.5; "IF" REV "THEN" A:= -A-1; REC:= A >= .5; "IF" REC "THEN" "BEGIN" NA:=ENTIER(A+.5); A:= A - NA "END"; "IF" A = -.5 "THEN" F:= G:= SQRT(PI / X / 2) "ELSE" "IF" X < 1 "THEN" "BEGIN" "COMMENT" 1SECTION : 6.10.2 (DECEMBER 1978) PAGE 14 ; "REAL" EXPON; "PROCEDURE" BESS KA01(A, X, KA, KA1);"CODE" 35191; EXPON:= EXP(X); BESS KA01(A, X, KA, KA1); F:= EXPON * KA; G:= EXPON * KA1 "END" "ELSE" "BEGIN" "REAL" B, C, E, P, Q; C:=.25-A*A;B:=X+X;G:=1;F:=0;E:=COS(A*PI)/PI*X*"15; "FOR" N:=1,N+1 "WHILE" H*N0 "DO" "BEGIN" P:=(N-1+C/N)/(E+(N+1)*(2-P));Q:=P*(1+Q) "END"; F:=SQRT(PI/B)/(1+Q);G:=F*(A+X+.5-P)/X "END"; "IF" REC "THEN" "BEGIN" X:= 2 / X; "FOR" N:= 1 "STEP" 1 "UNTIL" NA "DO" "BEGIN" H:= F + (A + N) * X * G; F:= G; G:= H "END" "END"; "IF" REV "THEN" "BEGIN" KA1:= F; KA:= G "END" "ELSE" "BEGIN" KA:= F; KA1:= G "END" "END" NONEXP BESS KA01; "EOP" "CODE" 35195; "PROCEDURE" NONEXP BESS KAPLUSN(A, X, NMAX, KAN); "VALUE" A, X, NMAX; "REAL" A, X; "INTEGER" NMAX; "ARRAY" KAN; "BEGIN" "INTEGER" N; "REAL" K1; "PROCEDURE" NONEXP BESS KA(A, X, KA, KA1); "VALUE" A, X; "REAL" A, X, KA, KA1; "CODE" 35194; NONEXP BESS KA(A, X, KAN[0], K1); A:= A-1; X:= 2/X; "IF" NMAX > 0 "THEN" KAN[1]:= K1; "FOR" N:= 2 "STEP" 1 "UNTIL" NMAX "DO" KAN[N]:= KAN[N-2] + (A+N)*X*KAN[N-1]; "END" NONEXP BESS KAPLUSN; "EOP" 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 1 AUTHORS: M. BAKKER AND N.M. TEMME. INSTITUTE: MATHEMATICAL CENTRE. RECEIVED: 780601. BRIEF DESCRIPTION: THIS SECTION CONTAINS THE FOLLOWING PROCEDURES: BESS J0; COMPUTES THE ORDINARY BESSEL FUNCTION OF THE FIRST KIND OF ORDER ZERO WITH ARGUMENT X; BESS J1; COMPUTES THE ORDINARY BESSEL FUNCTION OF THE FIRST KIND OF ORDER ONE WITH ARGUMENT X; BESS J; GENERATES AN ARRAY OF ORDINARY BESSEL FUNCTIONS OF THE FIRST KIND OF ORDER L (L = 0,...,N) WITH ARGUMENT X; BESS Y01; COMPUTES THE ORDINARY BESSEL FUNCTIONS OF THE SECOND KIND OF ORDERS ZERO AND ONE WITH ARGUMENT X; X > 0; BESS Y; GENERATES AN ARRAY OF ORDINARY BESSEL FUNCTIONS OF THE SECOND KIND OF ORDER L ( L = 0,...N) WITH ARGUMENT X; X> 0; BESS PQ0; THIS PROCEDURE IS AN AUXILIARY PROCEDURE FOR THE COMPUTATION OF THE ORDINARY BESSEL FUNCTIONS OF ORDER ZERO FOR LARGE VALUES OF THEIR ARGUMENT; BESS PQ1; THIS PROCEDURE IS AN AUXILIARY PROCEDURE FOR THE COMPUTATION OF THE ORDINARY BESSEL FUNCTIONS OF ORDER ONE FOR LARGE VALUES OF THEIR ARGUMENT. 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 2 KEYWORDS: BESSEL FUNCTION, ORDINARY BESSEL FUNCTION OF THE FIRST KIND, ORDINARY BESSEL FUNCTION OF THE SECOND KIND. REFERENCES: [1] ABRAMOWITZ, M., AND STEGUN, I. (EDS), HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND MATHEMATICAL TABLES. APPL. MATH. SER. 55, U.S. GOVT. PRINTING OFFICE, WASHINGTON, D.C. (1964). [2] C.W. CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS, NAT. PHYS. LAB. MATH. TABLES, VOL. 5, HER MAJESTY'S STATIONARY OFFICE, LONDON (1962). [3] W. GAUTSCHI, COMPUTATIONAL ASPECTS OF THREE TERM RECURRENCE RELATIONS, SIAM REVIEW, VOL. 9, 24-82 (1967). SUBSECTION: BESS J0. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS : "REAL" "PROCEDURE" BESS J0(X); "VALUE" X; "REAL" X; "CODE" 35160; BESS J0 DELIVERS THE ORDINARY BESSEL FUNCTION OF THE FIRST KIND OF ORDER ZERO WITH ARGUMENT X; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTION. PROCEDURES USED: BESS PQ0 = CP 35165. REQUIRED CENTRAL MEMORY: NO ARRAYS ARE USED. 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 3 RUNNING TIME: FOR ABS(X) < 8: LESS THAN 3 MS, FOR ABS(X) >= 8: LESS THAN 5 MS, ON THE CYBER 73/28. METHOD AND PERFORMANCE: CHEBYSHEV SERIES FROM [2]. EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X; "REAL" "PROCEDURE" BESS J0(X); "CODE" 35160; X:= 1; OUTPUT(61,"("/,D,6B-.14D "-ZD")", X, BESS J0(X)) "END" PRINTS THE FOLLOWING RESULTS: 1 .76519768655794" 0 SUBSECTION: BESS J1. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" BESS J1(X); "VALUE" X; "REAL" X; "CODE" 35161; BESS J1 DELIVERS THE ORDINARY BESSEL FUNCTION OF THE FIRST KIND OF ORDER ONE WITH ARGUMENT X; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTION. PROCEDURES USED: BESS PQ1 = CP 35166. REQUIRED CENTRAL MEMORY: NO ARRAYS ARE USED. 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 4 RUNNING TIME: FOR ABS(X) < 8: LESS THAN 3 MS, FOR ABS(X) >= 8: LESS THAN 5 MS, ON THE CYBER 73/28. METHOD AND PERFORMANCE: CHEBYSHEV SERIES FROM [2]. EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X; "REAL" "PROCEDURE" BESS J1(X); "CODE" 35161; X:= 1; OUTPUT(61,"("/,D,6B-.14D "-ZD")", X, BESS J1(X))) "END" DELIVERS THE FOLLOWING RESULTS: 1 .44005058574492" 0 SUBSECTION: BESS J. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS : "PROCEDURE" BESS J(X,N,J); "VALUE" X,N; "INTEGER" N; "REAL" X; "ARRAY" J; "CODE" 35162; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTIONS; N: ; THE UPPER BOUND OF THE INDICES OF ARRAY J; N >= 0; J: ; "ARRAY" J[0:N]; EXIT: J[L] IS THE ORDINARY BESSEL FUNCTION OF THE FIRST KIND OF ORDER L AND ARGUMENT X. PROCEDURES USED: START = CP 35185; 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 5 REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. RUNNING TIME: ROUGHLY PROPORTIONAL TO THE MAXIMUM OF 1.359 * X + 72 AND N + 18. METHOD AND PERFORMANCE: MILLER'S ALGORITHM, SEE [3]. EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X; "ARRAY" J[0:1]; "REAL" "PROCEDURE" BESS J0(X); "CODE" 35160; "REAL" "PROCEDURE" BESS J1(X); "CODE" 35161; "PROCEDURE" BESS J(X,N,J,); "CODE" 35162; "FOR" X:= 1,5,10,25 "DO" "BEGIN" BESS J(X,1,J); OUTPUT(61,"("ZZ.D, 2(BB-.D"-ZD),/")", X, J[0] - BESS J0(X),J[1] - BESS J1(X)) "END" "END" DELIVERS THE FOLLOWING RESULTS: 1.0 .2"-13 .2"-13 5.0 -.8"-14 -.4"-14 10.0 -.4"-14 .4"-14 25.0 -.1"-14 -.9"-15 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 6 SUBSECTION: BESS Y01. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS : "PROCEDURE" BESS Y01(X,Y0,Y1); "VALUE" X; "REAL" X,Y0,Y1; "CODE" 35163; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTIONS; X > 0; Y0: ; EXIT: Y0 HAS THE VALUE OF THE ORDINARY BESSEL FUNCTION OF THE SECOND KIND OF ORDER 0 AND ARGUMENT X; Y1: ; EXIT: Y1 HAS THE VALUE OF THE ORDINARY BESSEL FUNCTION OF THE SECOND KIND OF ORDER 1 AND ARGUMENT X. PROCEDURES USED: BESS J0 = CP 35160, BESS J1 = CP 35161, BESS PQ0 = CP 35165, BESS PQ1 = CP 35166. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. RUNNING TIME: ABOUT 15 MS, ON THE CYBER 73/28. METHOD AND PERFORMANCE: CHEBYSHEV SERIES FROM [2]. EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X,Y0,Y1; "PROCEDURE" BESS Y01(X,Y0,Y1); "CODE" 35163; X:= 1; BESS Y01(X,Y0,Y1); OUTPUT(61,"("/,4BD.D,2(B-.14D"-ZD)")",X,Y0,Y1) "END" 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 7 DELIVERS THE FOLLOWING RESULTS: 1.0 988256964215676" -1 -.78121282130028" 0 SUBSECTION: BESS Y. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS Y(X,N,Y); "VALUE" X,N; "INTEGER" N; "REAL" X; "ARRAY" Y; "CODE" 35164; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT; THIS ARGUMENT SHOULD SATISFY X > 0; N: ; THE UPPER BOUND OF THE INDICES OF THE ARRAY Y; N >= 0; Y: ; "ARRAY" Y[0:N]; EXIT: Y[I] IS THE VALUE OF THE ORDINARY BESSEL FUNCTION OF THE SECOND KIND OF ORDER I (I = 0,...,N) AND ARGUMENT X. PROCEDURES USED: BESS Y01 = CP 35163. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED. RUNNING TIME: DEPENDS ON N; SEE BESS Y01. METHOD AND PERFORMANCE: Y[0] AND Y[1] ARE COMPUTED BY USING BESS Y01 (CP 35163); THE REMAINING Y[I] ARE COMPUTED BY USING THE RECURRENCE RELATION Y[I+1]:= Y[I] * 2 * I/X - Y[I-1], I >= 1. 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 8 EXAMPLE OF USE: THE PROGRAM "BEGIN" "ARRAY" Y[0:2]; "PROCEDURE" BESS Y(X,N,Y); "CODE" 35164; BESS Y(1,2,Y); OUTPUT(61,"("3(-.14D"-ZD)")", Y[0], Y[1], Y[2]) "END" PRINTS THE FOLLOWING RESULTS: 8.8256964215676"- 2 -7.8121282130028"- 1 -1.6506826068162" 0 SUBSECTION: BESS PQ0. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS PQ0(X,P0,Q0); "VALUE" X; "REAL" X,P0,Q0; "CODE" 35165; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT; THIS ARGUMENT SHOULD SATISFY X > 0; P0: ; EXIT: P0 CORRESPONDS WITH THE FUNCTION P(X,0) DEFINED IN [1,FORMULAS 9.2.5 AND 9.2.6]; Q0: ; EXIT: Q0 CORRESPONDS WITH THE FUNCTION Q(X,0) DEFINED IN [1,FORMULAS 9.2.5 AND 9.2.6]. PROCEDURES USED: BESS J0 = CP 35160, BESS Y01 = CP 35163. REQUIRED CENTRAL MEMORY: NO ARRAYS ARE USED. RUNNING TIME: ABOUT 15 MS, ON THE CYBER 73/28. 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 9 METHOD AND PERFORMANCE: FOR X >= 8 CHEBYSHEV SERIES FROM [2], FOR X < 8 WITH BESS J0 AND BESS Y01. EXAMPLE OF USE: SEE SUBSECTION BESS PQ1. SUBSECTION: BESS PQ1. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS PQ1(X,P1,Q1); "VALUE" X; "REAL" X,P1,Q1; "CODE" 35166; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT; THIS ARGUMENT SHOULD SATISFY X > 0; P1: ; EXIT: P1 CORRESPONDS WITH THE FUNCTION P(X,1) DEFINED IN [1,FORMULAS 9.2.5 AND 9.2.6]; Q1: ; EXIT: Q1 CORRESPONDS WITH THE FUNCTION Q(X,1) DEFINED IN [1,FORMULAS 9.2.5 AND 9.2.6]. PROCEDURES USED: BESS J1 = CP 35161, BESS Y01 = CP 35163. REQUIRED CENTRAL MEMORY: NO ARRAYS ARE USED. RUNNING TIME: ABOUT 15 MS, ON THE CYBER 73/28. METHOD AND PERFORMANCE: FOR X >= 8 CHEBYSHEV SERIES FROM [2], FOR X < 8 WITH BESS J1 AND BESS Y01. 1SECTION 6.9.1 (DECEMBER 1978) PAGE 10 EXAMPLE OF USE: FROM THE WRONSKIAN RELATION [1,9.1.16] IT CAN BE SHOWN THAT P0 * P1 + Q0 * Q1 = 1, WHATEVER X. IN THE FOLLOWING PROGRAM WE VERIFY THIS RELATION. "BEGIN" "REAL" X,P,Q,R,S; "PROCEDURE" BESS PQ0(X,P0,Q0); "CODE" 35165; "PROCEDURE" BESS PQ1(X,P1,Q1); "CODE" 35166; "FOR" X:= 1,3,5,10 "DO" "BEGIN" BESSPQ0(X,P,Q); BESSPQ1(X,R,S); OUTPUT(61,"("BB,D.2D"+3D")", ABS(P*R + Q*S -1)) "END" "END" THE RESULTS ARE: 4.97"-014 4.26"-014 5.68"-014 7.11"-015 SOURCE TEXT(S): "CODE" 35160; "REAL" "PROCEDURE" BESS J0(X); "VALUE" X; "REAL" X; "IF" X=0 "THEN" BESS J0:= 1 "ELSE" "IF" ABS(X) < 8 "THEN" "BEGIN" "REAL" Z, Z2, AR, B0, B1, B2; X:= X/8; Z:= 2*X*X - 1; Z2:= Z + Z; B1:= B2:= 0; "FOR" AR:=-.75885"-15, +.4125321 "-13, -.194383469 "-11, +.7848696314 "-10, -.267925353056 "- 8, +.7608163592419 "- 7, -.176194690776215"- 5, +.324603288210051"- 4, -.46062616620628 "- 3, +.48191800694676 "- 2, -.34893769411409 "- 1, +.158067102332097 , -.37009499387265 "- 0, +.265178613203337 , -.872344235285222"- 2 "DO" "BEGIN" B0:= Z2*B1-B2+AR; B2:= B1; B1:= B0 "END"; BESS J0:= Z*B1 - B2 + .15772 79714 7489 "END" "ELSE" "BEGIN" "REAL" C, COSX, SINX, P0, Q0; "PROCEDURE" BESS PQ0(X, P0, Q0); "CODE" 35165; X:= ABS(X); C:= .79788 45608 02865 / SQRT(X); COSX:= COS(X-.70685 83470 57703" 1); SINX:= SIN(X-.70685 83470 57703" 1); BESS PQ0(X, P0, Q0); BESSJ0:= C * (P0 * COSX - Q0 * SINX) "END" BESS J0; "EOP" 1SECTION 6.9.1 (DECEMBER 1978) PAGE 11 "CODE" 35161; "REAL" "PROCEDURE" BESS J1(X); "VALUE" X; "REAL" X; "IF" X=0 "THEN" BESS J1:= 0 "ELSE" "IF" ABS(X) < 8 "THEN" "BEGIN" "REAL" Z, Z2, AR, B0, B1, B2; X:= X/8; Z:= 2*X*X - 1; Z2:= Z + Z; "COMMENT" COMPUTATION OF J1; B1:= B2:= 0; "FOR" AR:= -.19554 "-15, +.1138572 "-13, -.57774042 "-12, +.2528123664 "-10, -.94242129816 "- 9, +.2949707007278 "- 7, -.76175878054003 "- 6, +.158870192399321"- 4, -.260444389348581"- 3, +.324027018268386"- 2, -.291755248061542"- 1, +.177709117239728"- 0, -.661443934134543"- 0, +.128799409885768"+ 1, -.119180116054122"+ 1 "DO" "BEGIN" B0:= Z2*B1-B2+AR; B2:= B1; B1:= B0 "END"; BESS J1:= X * (Z * B1 - B2 + .64835 87706 05265) "END" "ELSE" "BEGIN" "REAL" C, COSX, SINX, P1, Q1; "INTEGER" SGNX; "PROCEDURE" BESS PQ1(X, P1, Q1); "CODE" 35166; SGNX:= SIGN(X); X:= ABS(X); C:= .79788 45608 02865 / SQRT(X); COSX:= COS(X-.70685 83470 57703"+1); SINX:= SIN(X-.70685 83470 57703"+1); BESS PQ1(X, P1, Q1); BESS J1:= SGNX * C * (P1*SINX + Q1*COSX) "END" BESS J1; "EOP" "CODE" 35162; "PROCEDURE" BESS J(X, N, J); "VALUE" X, N; "REAL"X; "INTEGER" N; "ARRAY" J; "IF" X=0 "THEN" "BEGIN" J[0]:= 1; "FOR" N:= N "STEP" -1 "UNTIL" 1 "DO" J[N]:= 0 "END" "ELSE" "COMMENT" 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 12 ; "BEGIN""REAL" X2, R, S; "INTEGER" L, M, NU, SIGNX; "INTEGER" "PROCEDURE" START(X,N,T); "CODE" 35185; SIGNX:= SIGN(X); X:= ABS(X); R:= S:= 0; X2:= 2/X; L:= 0; NU:= START(X,N,0); "FOR" M:= NU "STEP" -1 "UNTIL" 1 "DO" "BEGIN" R:= 1/(X2*M-R); L:= 2-L; S:= R*(L+S); "IF" M<=N "THEN" J[M]:= R "END"; J[0]:= R:= 1/(1+S); "FOR" M:= 1 "STEP" 1 "UNTIL" N "DO" J[M]:= R:= R*J[M]; "IF" SIGNX < 0 "THEN" "FOR" M:= 1 "STEP" 2 "UNTIL" N "DO" J[M]:= -J[M]; "END" BESSELJ; "EOP" "CODE" 35163; "PROCEDURE" BESS Y01(X, Y0, Y1); "VALUE" X; "REAL" X, Y0, Y1; "IF" X< 8 "THEN" "BEGIN" "REAL" Z, Z2, C, LNX, AR, B0, B1, B2; "REAL" "PROCEDURE" J0(X); "CODE" 35160; "REAL" "PROCEDURE" J1(X); "CODE" 35161; C:= .63661 97723 67581; LNX:= C * LN(X); C:= C/X; X:= X/8; Z:= 2*X*X - 1; Z2:= Z + Z; "COMMENT" COMPUTATION OF Y0; B1:= B2:= 0; "FOR" AR:= +.164349 "-14, -.8747341 "-13, +.402633082 "-11, -.15837552542 "- 9, +.524879478733 "- 8, -.14407233274019 "- 6, +.32065325376548 "- 5, -.563207914105699"- 4, +.753113593257774"- 3, -.72879624795521 "- 2, +.471966895957634"- 1, -.177302012781143"- 0, +.261567346255047"- 0, +.179034314077182"- 0, -.274474305529745"DO" "BEGIN" B0:= Z2*B1-B2+AR; B2:= B1; B1:= B0 "END"; Y0:= LNX * J0(8*X)+Z*B1-B2-.33146 11320 3285"-1; "COMMENT" COMPUTATION OF Y1; B1:= B2:= 0; "COMMENT" 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 13 ; "FOR" AR:= +.42773 "-15, -.2440949 "-13, +.121143321 "-11, -.5172121473 "-10, +.187547032473 "- 8, -.5688440039919 "- 7, +.141662436449235"- 5, -.283046401495148"- 4, +.440478629867099"- 3, -.51316411610611 "- 2, +.423191803533369"- 1, -.226624991556755"- 0, +.675615780772188"- 0, -.767296362886646"- 0, -.128697384381350"- 0"DO" "BEGIN" B0:= Z2*B1-B2+AR; B2:= B1; B1:= B0 "END"; Y1:= LNX * J1(X*8)-C + X * (Z*B1-B2+.20304 10588 593425"-1) "END" "ELSE" "BEGIN" "REAL" C, COSX, SINX, P0, Q0, P1, Q1; "PROCEDURE" BESS PQ0(X, P0, Q0); "CODE" 35165; "PROCEDURE" BESS PQ1(X, P1, Q1); "CODE" 35166; C:= .79788 45608 02865 / SQRT(X); BESS PQ0(X, P0, Q0); BESS PQ1(X, P1, Q1); X:= X-.70685 83470 57703"1; COSX:= COS(X); SINX:= SIN(X); Y0:= C * (P0*SINX + Q0*COSX); Y1:= C * (Q1*SINX - P1*COSX) "END" BESS Y01; "EOP" "CODE" 35164; "PROCEDURE" BESS Y(X, N, Y); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" Y; "BEGIN" "INTEGER" I; "REAL" Y0, Y1, Y2; "PROCEDURE" BESS Y01(X, Y0, Y1); "CODE" 35163; BESS Y01(X, Y0, Y1); Y[0]:= Y0; "IF" N > 0 "THEN" Y[1]:= Y1 ; X:= 2/X; "FOR" I:=2 "STEP" 1 "UNTIL" N "DO" "BEGIN" Y[I]:= Y2:= (I-1)*X*Y1 - Y0; Y0:= Y1; Y1:= Y2 "END" "END" BESS Y; "EOP" 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 14 "CODE" 35165; "PROCEDURE" BESS PQ0(X, P0, Q0); "VALUE" X; "REAL" X, P0, Q0; "IF" X < 8 "THEN" "BEGIN" "REAL" B, COSX, SINX, J0X, Y0; "REAL" "PROCEDURE" J0(X); "CODE" 35160; "PROCEDURE" BESS Y01(X, Y0, Y1); "CODE" 35163; B:= SQRT(X) * 1.2533 14137 31550; BESS Y01(X, Y0, J0X); J0X:= J0(X); X:= X-.78539 81633 97448; COSX:= COS(X); SINX:= SIN(X); P0:= B * (Y0 * SINX + J0X * COSX); Q0:= B * (Y0 * COSX - J0X * SINX) "END" "ELSE" "BEGIN" "REAL" X2, AR, B0, B1, B2, Y; Y:= 8/X; X:= 2*Y*Y-1; X2:= X+X; B1:= B2:= 0; "FOR" AR:= -.10012 "-15, +.67481 "-15, -.506903 "-14, +.4326596 "-13, -.43045789 "-12, +.516826239 "-11, -.7864091377 "-10, +.163064646352 "- 8, -.5170594537606 "- 7, +.30751847875195 "- 5, -.536522046813212"- 3 "DO" "BEGIN" B0:= X2 * B1 - B2 + AR; B2:= B1; B1:= B0 "END"; P0:= X * B1 - B2 + .99946034934752; "COMMENT" COMPUTATION OF Q0; B1:= B2:= 0; "FOR" AR:= -.60999 "-15, +.425523 "-14, -.3336328 "-13, +.30061451 "-12, -.320674742 "-11, +.4220121905 "-10, -.72719159369 "- 9, +.1797245724797 "- 7, -.74144984110606 "- 6, +.683851994261165"- 4 "DO" "BEGIN" B0:= X2 * B1 - B2 + AR; B2:= B1; B1:= B0 "END"; Q0:=(X * B1 - B2 -.015555854605337) * Y "END" BESS PQ0; "EOP" 1SECTION : 6.9.1 (DECEMBER 1978) PAGE 15 "CODE" 35166; "PROCEDURE" BESS PQ1(X, P1, Q1); "VALUE" X; "REAL" X, P1, Q1; "IF" X < 8 "THEN" "BEGIN" "REAL" B, COSX, SINX, J1X, Y1; "REAL" "PROCEDURE" J1(X); "CODE" 35161; "PROCEDURE" BESS Y01(X, Y0, Y1); "CODE" 35163; B:= SQRT(X) * 1.253314137 31550; BESS Y01(X, J1X, Y1); J1X:= J1(X); X:= X-.78539 81633 97448; COSX:= COS(X); SINX:= SIN(X); P1:= B * (J1X * SINX - Y1 * COSX); Q1:= B * (J1X * COSX + Y1 * SINX) "END" "ELSE" "BEGIN" "REAL" X2, AR, B0, B1, B2, Y; Y:= 8 / X; X:= 2 * Y * Y - 1; X2 := X + X; "COMMENT" COMPUTATION OF P1; B1:= B2:= 0; "FOR" AR:= +.10668"-15, -.72212 "-15, +.545267 "-14, -.4684224 "-13, +.46991955 "-12, -.570486364 "-11, +.881689866 "-10, -.187189074911 "- 8, +.6177633960644 "- 7, -.39872843004889 "- 5, +.89898983308594 "- 3 "DO" "BEGIN" B0:= B1 * X2 - B2 + AR; B2:= B1; B1:= B0 "END"; P1:= X * B1 - B2 + 1.0009030408600137; "COMMENT" COMPUTATION OF Q1; B1:= B2:= 0; "FOR" AR:= -.10269 "-15, +.65083 "-15, -.456125 "-14, +.3596777 "-13, -.32643157 "-12, +.351521879 "-11, -.4686363688 "-10, +.82291933277 "- 9, -.2095978138408 "- 7, +.91386152579555 "- 6, -.96277235491571 "- 4 "DO" "BEGIN" B0:= X2 * B1 - B2 + AR; B2:= B1; B1:= B0 "END"; Q1:=(X * B1 - B2 + .46777787069535" -1) * Y "END" BESS PQ1; "EOP" 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 1 AUTHORS: M. BAKKER AND N.M. TEMME. INSTITUTE: MATHEMATICAL CENTRE. RECEIVED: 750201. BRIEF DESCRIPTION: THIS SECTION CONTAINS THE FOLLOWING PROCEDURES: BESS I0; COMPUTES THE MODIFIED BESSEL FUNCTION OF THE FIRST KIND OF ORDER ZERO WITH ARGUMENT X; BESS I1; COMPUTES THE MODIFIED BESSEL FUNCTION OF THE FIRST KIND OF ORDER ONE WITH ARGUMENT X; BESS I; GENERATES AN ARRAY OF MODIFIED BESSEL FUNCTIONS OF THE FIRST KIND OF ORDER L (L = 0, ..., N) WITH ARGUMENT X; BESS K01; COMPUTES THE MODIFIED BESSEL FUNCTIONS OF THE THIRD KIND OF ORDERS ZERO AND ONE WITH ARGUMENT X; X > 0; BESS K; GENERATES AN ARRAY OF MODIFIED BESSEL FUNCTIONS OF THE THIRD KIND OF ORDER L ( L = 0, ..., N) WITH ARGUMENT X; X > 0; NONEXP BESS I0; DOES THE SAME AS BESS I0, BUT THE RESULT IS MULTIPLIED BY EXP(-ABS(X)); NONEXP BESS I1; DOES THE SAME AS BESS I1, BUT THE RESULT IS MULTIPLIED BY EXP(-ABS(X)); NONEXP BESS I; DOES THE SAME AS BESS I, BUT THE ARRAY ELEMENTS ARE MULTIPLIED BY EXP(-ABS(X)); NONEXP BESS K01; DOES THE SAME AS BESS K01, BUT THE RESULTS ARE MULTIPLIED BY EXP(X); 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 2 NONEXP BESS K; DOES THE SAME AS BESS K, BUT THE ARRAY ELEMENTS ARE MULTIPLIED BY EXP(X). KEYWORDS: BESSEL FUNCTIONS, MODIFIED BESSEL FUNCTIONS, INTEGER ORDER. REFERENCES: [1] M.ABRAMOWITZ AND I.A. STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, DOVER PUBLICATIONS, INC., NEW YORK, 1968. [2] D.B.HUNTER, THE CALCULATION OF SOME BESSEL FUNCTIONS, MATHEMATICS OF COMPUTATION (1964), P. 123. [3] YUDELL LUKE, THE SPECIAL FUNCTIONS AND THEIR APPROXIMATIONS, VOLUME 2, ACADEMIC PRESS, NEW YORK AND LONDON (1969). [4] C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS, NAT. PHYS. LAB. MATH. TABLES, VOLUME 5, HER MAJESTY,S STATIONARY OFFICE, LONDON (1962). [5] W.GAUTSCHI, COMPUTATIONAL ASPECTS OF THREE TERM RECURRENCE RELATIONS, SIAM REVIEWS, VOLUME 9 (1967), P. 24. [6] J.M.BLAIR, RATIONAL CHEBYSHEV APPROXIMATIONS FOR THE MODIFIED BESSEL FUNCTIONS I0(X) AND I1(X); MATHEMATICS OF COMPUTATIONS,VOLUME 28, NR 126, APRIL 1974, P. 581-583. 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 3 SUBSECTION: BESS I0. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" BESS I0(X); "VALUE" X; "REAL" X; "CODE" 35170; BESS I0 DELIVERS THE MODIFIED BESSEL FUNCTION OF THE FIRST KIND OF ORDER ZERO WITH ARGUMENT X; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTION. PROCEDURES USED: NONEXP BESS I0 = CP35175. REQUIRED CENTRAL MEMORY: NO ARRAYS ARE USED. RUNNING TIME: FOR X = 0 BESS I0 IS ASSIGNED ITS VALUE IMMEDIATELY; FOR 0 < ABS(X) <= 15.0 17 MULTIPLICATIONS AND ONE DIVISION ARE REQUIRED; FOR ABS(X) > 15.0 11 MULTIPLICATIONS, 3 DIVISIONS, ONE EVALUATION OF THE SQUARE ROOT AND ONE EVALUATION OF THE EXPONENNTIAL FUNCTION ARE REQUIRED. METHOD AND PERFORMANCE: RATIONAL APPROXIMATION, SEE [6]. EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X; "REAL" "PROCEDURE" BESS I0(X); "CODE" 35170; X:= 1; OUTPUT(61,"("/,D,6B-.14D"-ZD")", X, BESS I0(X)) "END" 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 4 PRINTS THE FOLLOWING RESULTS: 1 .12660658777520" 1 SUBSECTION: BESS I1. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" BESS I1(X); "VALUE" X; "REAL" X; "CODE" 35171; BESS I1 DELIVERS THE MODIFIED BESSEL FUNCTION OF THE FIRST KIND OF ORDER ONE WITH ARGUMENT X; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTION. PROCEDURES USED: NONEXP BESS I1 = CP35176. REQUIRED CENTRAL MEMORY: NO ARRAYS ARE USED. RUNNING TIME: FOR X = 0 BESS I1 IS ASSIGNED ITS VALUE IMMEDIATELY; FOR 0 < ABS(X) <= 15.0 17 MULTIPLICATIONS AND ONE DIVISION ARE REQUIRED; FOR ABS(X) > 15.0 12 MULTIPLICATIONS, 3 DIVISIONS, ONE EVALUATION OF THE SQUARE ROOT AND ONE EVALUATION OF THE EXPONENTIAL FUNCTION ARE REQUIRED. METHOD AND PERFORMANCE: RATIONAL APPROXIMATION, SEE [6]. 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 5 EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X; "REAL" "PROCEDURE" BESS I1(X); "CODE" 35171; X:= 1; OUTPUT(61,"("/,D,6B-.14D"-ZD")", X, BESS I1(X)) "END" PRINTS THE FOLLOWING RESULTS: 1 .56515910399252" 0 SUBSECTION: BESS I. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS I(X, N, I); "VALUE" X, N; "INTEGER" N; "REAL" X; "ARRAY" I; "CODE" 35172; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTIONS; N: ; THE UPPER BOUND OF THE INDICES OF THE ARRAY I; I: ; "ARRAY" I[0 : N]; EXIT: I[L] POSSESSES THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE FIRST KIND OF ORDER L (0 <= L <= N). METHOD AND PERFORMANCE: SEE NON EXP BESS I (THIS SECTION). PROCEDURES USED : NONEXP BESS I = CP 35177. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE USED. 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 6 RUNNING TIME: ROUGHLY PROPORTIONAL TO THE MAXIMUM OF 1.359 * X + 72 AND N + 18. EXAMPLE OF USE : THE FOLLOWING PROGRAM CHECKS FOR X = 1 (1) 20 THE WRONSKIAN RELATION X * (I[N - 1] * K[N] + I[N] * K[N - 1]) - 1 = 0 FOR N = 1 (1) 5; THE PROGRAM READS: "BEGIN" "REAL" X; "INTEGER" N; "ARRAY" I, K[0:5]; "PROCEDURE" BESS I(X, N, I); "CODE" 35172; "PROCEDURE" BESS K(X, N, K); "CODE" 35174; "FOR" X:= 1 "STEP" 1 "UNTIL" 20 "DO" "BEGIN" OUTPUT(61,"("/ZD")", X); BESS I(X, 5, I); BESS K(X, 5, K); "FOR" N:= 1, 2, 3, 4, 5 "DO" OUTPUT(61,"("BB-.D"-ZD")", X * (I[N] * K[N - 1] + I[N - 1] * K[N]) - 1) "END" "END" THE RESULTS ARE: 1 .0" 0 .0" 0 -.7"-14 -.7"-14 -.7"-14 2 .0" 0 .0" 0 .0" 0 .0" 0 .0" 0 3 .7"-14 .7"-14 .0" 0 .0" 0 .0" 0 4 .7"-14 .0" 0 .0" 0 .0" 0 .0" 0 5 .0" 0 .7"-14 .7"-14 .0" 0 .0" 0 6 .0" 0 .0" 0 .0" 0 .0" 0 -.7"-14 7 .0" 0 .0" 0 .0" 0 .0" 0 .0" 0 8 -.1"-13 -.1"-13 -.1"-13 -.1"-13 -.1"-13 9 .0" 0 .0" 0 .0" 0 -.7"-14 -.7"-14 10 .0" 0 .0" 0 .0" 0 .0" 0 .0" 0 11 .0" 0 .0" 0 .0" 0 .0" 0 .0" 0 12 .0" 0 .0" 0 .0" 0 .0" 0 .0" 0 13 .7"-14 .7"-14 .0" 0 .7"-14 .0" 0 14 .0" 0 .7"-14 .0" 0 .0" 0 .0" 0 15 .0" 0 .0" 0 .0" 0 .0" 0 .0" 0 16 .0" 0 .0" 0 .0" 0 .0" 0 -.7"-14 17 .7"-14 .0" 0 .0" 0 .0" 0 .0" 0 18 .7"-14 .0" 0 .0" 0 .0" 0 -.7"-14 19 .7"-14 .0" 0 .0" 0 .0" 0 .0" 0 20 .0" 0 .0" 0 .0" 0 .0" 0 -.7"-14 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 7 SUBSECTION: BESS K01. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS K01(X, K0, K1); "VALUE" X; "REAL" X, K0, K1; "CODE" 35173; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTIONS; X > 0; K0: ; EXIT: K0 HAS THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER 0 WITH ARGUMENT X; K1: ; EXIT: K1 HAS THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER ONE. PROCEDURES USED: NONEXP BESS K01 = CP35178 REQUIRED CENTRAL MEMORY: NO ARRAYS ARE USED. RUNNING TIME: DEPENDS ON THE VALUE OF X; THE GLOBAL VALUES IN MILLISECONDS ARE: 0 < X <= 1.5 : 2.2 MS, 1.5 < X <= 5.0 : 5.5 MS, 5.0 < X : 2.3 MS, ON THE CYBER 73/28. METHOD AND PERFORMANCE: FOR THE COMPUTATION OF K0 AND K1 THREE DIFFERENT METHODS ARE USED DEPENDING ON THE VALUE OF X: FOR 0 < X <= 1.5 K0 AND K1 ARE EVALUATED BY MEANS OF TAYLOR SERIES EXPANSIONS (SEE [1], P. 375, FORMULA 9.6.13); FOR X > 1.5 K0 AND K1 ARE COMPUTED BY MEANS OF A CALL OF THE CODE PROCEDURE NONEXP BESS K01 (SEE DESCRIPTION AHEAD) AND MULTIPLICATION BY EXP(- X). 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 8 EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X, K0, K1; "PROCEDURE" BESS K01(X, K0, K1); "CODE" 35173; "FOR" X:= .5, 1.5, 2.5 "DO" "BEGIN" BESS K01(X, K0, K1); OUTPUT(61,"("/,4BD.D,2(B-.14D"-ZD)")",X,K0,K1) "END" "END" PRINTS THE FOLLOWING RESULTS: 0.5 .92441907122766" 0 .16564411200033" 1 1.5 .21380556264754" 0 .27738780045683" 0 2.5 .62347553200366" -1 .73890816347746" -1 SUBSECTION: BESS K. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" BESS K(X, N, K); "VALUE" X, N; "INTEGER" N; "REAL" X; "ARRAY" K; "CODE" 35174; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTIONS; X > 0; N: ; THE UPPER BOUND OF THE INDICES OF THE ARRAY K; N >= 0; K: ; "ARRAY" K[0 : N]; EXIT: K[I] POSSESSES THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER I (0 <= I <= N). PROCEDURES USED: K01 = CP 35173. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE USED. RUNNING TIME : DEPENDS ON THE VALUE OF X (SEE TABLE BELONGING TO BESS K01) AND N. 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 9 METHOD AND PERFORMANCE: K[0], ..., K[N] ARE COMPUTED ACCORDING TO THE RECURRENCE RELATION K[I + 1] = K[I - 1] + (2 * I / X) * K[I], I = 2, ..., N, (SEE [1], P. 376, FORMULA 9.6.26). EXAMPLE OF USE: THE PROGRAM "BEGIN" "ARRAY" K[0 : 2]; "REAL" X; "PROCEDURE" BESS K(X, N, K); "CODE" 35174; "FOR" X:= .5, 1.0, 1.5, 2.0 "DO" "BEGIN" BESS K(X, 2, K); OUTPUT(61,"("/D.D,3(BB.12D"-D)")",X,K) "END" "END" PRINTS THE FOLLOWING RESULTS: 0.5 .924419071228"0 .165644112000"1 .755018355124"1 1.0 .421024438241"0 .601907230197"0 .162483889864"1 1.5 .213805562648"0 .277387800457"0 .583655963257"0 2.0 .113893872750"0 .139865881817"0 .253759754566"0 SUBSECTION: NONEXP BESS I0. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" NONEXP BESS I0(X); "VALUE" X; "REAL" X; "CODE" 35175; NONEXP BESS I0 DELIVERS THE MODIFIED BESSEL FUNCTION OF THE FIRST KIND OF ORDER 0 WITH ARGUMENT X MULTIPLIED BY EXP(-ABS(X)). THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTION. PROCEDURES USED: BESS I0 = CP35170. 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 10 REQUIRED CENTRAL MEMORY: NO ARRAYS ARE USED. RUNNING TIME: FOR X = 0 NONEXP BESS I0 IS ASSIGNED ITS VALUE IMMEDIATELY; FOR 0 < ABS(X) <= 15.0 18 MULTIPLICATIONS, ONE DIVISION AND ONE EVALUATION OF THE EXPONENTIAL FUNCTION ARE REQUIRED; FOR ABS(X) > 15.0 10 MULTIPLICATIONS, 3 DIVISIONS AND ONE EVALUATION OF THE SQUARE ROOT ARE REQUIRED. METHOD AND PERFORMANCE: SEE [6]. EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X; "REAL" "PROCEDURE" NONEXP BESS I0(X); "CODE" 35175; X:= 1; OUTPUT(61,"("/,D,6B-.14D"-ZD")", X, NONEXP BESS I0(X)) "END" PRINTS THE FOLLOWING RESULTS: 1 .46575960759364" 0 SUBSECTION: NONEXP BESS I1. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "REAL" "PROCEDURE" NONEXP BESS I1(X); "VALUE" X; "REAL" X; "CODE" 35176; NONEXP BESS I1 DELIVERS THE MODIFIED BESSEL FUNCTION OF THE FIRST KIND OF ORDER 1 WITH ARGUMENT X MULTIPLIED BY EXP(-ABS(X)). THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTION. 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 11 PROCEDURES USED: BESS I1 = CP35171. REQUIRED CENTRAL MEMORY: NO ARRAYS ARE USED. RUNNING TIME: FOR X = 0 NONEXP BESS I1 IS ASSIGNED ITS VALUE IMMEDIATELY; FOR 0 < ABS(X) <= 15.0 18 MULTIPLICATIONS, ONE DIVISION AND ONE EVALUATION OF THE EXPONENTIAL FUNCTION ARE REQUIRED; FOR X > 15.0 11 MULTIPLICATIONS, 3 DIVISIONS AND ONE EVALUATION OF THE SQUARE ROOT ARE REQUIRED. METHOD AND PERFORMANCE: SEE [6]. EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X; "REAL" "PROCEDURE" NONEXP BESS I1(X); "CODE" 35176; X:= 1; OUTPUT(61,"("/,D,6B-.14D"-ZD")", X, NONEXP BESS I1(X)) "END" DELIVERS THE FOLLOWING RESULTS: 1 .20791041534972" 0 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 12 SUBSECTION: NONEXP BESS I. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" NONEXP BESS I(X, N, I); "VALUE" X, N; "INTEGER" N; "REAL" X; "ARRAY" I; "CODE" 35177; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTIONS; N: ; THE UPPER BOUND OF THE INDICES OF THE ARRAY I; N >= 0; I: ; "ARRAY" I[0:N]; EXIT: I[L] POSSESSES THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE FIRST KIND OF ORDER L (L=0,..,N) MULTIPLIED BY EXP (- ABS(X)). PROCEDURES USED: START = CP 35185; REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE USED. RUNNING TIME: ROUGHLY PROPORTIONAL TO THE MAXIMUM OF 1.359*X + 72 AND N+18. METHOD AND PERFORMANCE: SEE [5]. EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X; "ARRAY" I[0:2]; "PROCEDURE" NONEXP BESS I(X, N, I); "CODE" 35179; "FOR" X:= .5, 1.0, 1.5, 2.0, 2.5 "DO" "BEGIN" NONEXP BESS I(X, 2, I); OUTPUT(61, "("/,4BZ.D,3(B-.12D"-D)")",X, I[0], I[1], I[2]) "END" "END" 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 13 PRINTS THE FOLLOWING RESULTS: .5 .152410938577" 1 .273100970821" 1 .124481482186" 2 1.0 .114446307981" 1 .163615348626" 1 .441677005233" 1 1.5 .958210053295" 0 .124316587355" 1 .261576455136" 1 2.0 .841568215071" 0 .103347684707" 0 .187504506214" 1 2.5 .759548690328" 0 .900174423908" 0 .147968822945" 1 SUBSECTION: NONEXP BESS K01. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" NONEXP BESS K01(X, K0, K1); "VALUE" X; "REAL" X, K0, K1; "CODE" 35178; THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTIONS; X > 0; K0: ; EXIT: K0 HAS THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER 0 WITH ARGUMENT X MULTIPLIED BY EXP(X); K1: ; EXIT: K1 HAS THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER 1 MULTIPLIED BY EXP(X). PROCEDURES USED: BESS K01 = CP35173. REQUIRED CENTRAL MEMORY: NO ARRAYS ARE USED. RUNNING TIME: DEPENDS ON THE VALUE OF X; BECAUSE OF THE STRONG INTERDEPENDENCE OF THE BESS K01 ( = CP35173) AND NONEXP BESS K01 THE READER IS REFERRED TO THE TABLE OF RUNNING TIMES BELONGING TO BESS K01. 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 14 METHOD AND PERFORMANCE: FOR THE COMPUTATION OF K0 AND K1 THREE DIFFERENT METHODS ARE USED DEPENDING ON THE VALUE OF X: FOR 0 < X <= 1.5 K0 AND K1 ARE COMPUTED BY MEANS OF MULTIPLICATION OF THE MODIFIED BESSEL FUNCTIONS OF ORDER ZERO AND ONE (SEE DESCRIPTION OF K0) BY EXP(X); FOR 1.5 < X <= 5 K0 AND K1 ARE COMPUTED BY THE EVALUATION OF THEIR INTEGRAL REPRESENTATIONS (SEE [1], P. 376, FORMULA 9.6.23) BY MEANS OF THE TRAPEZOIDAL RULE (SEE [2]); FOR X > 5 K0 AND K1 ARE COMPUTED BY MEANS OF A FINITE CHEBYSHEV SERIES EXPANSION (SEE [3], P. 339 AND [4]). EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X, K1; "PROCEDURE" NONEXP BESS K01(X, K0, K1); "CODE" 35178; "FOR" X:= .5, 1.0, 1.5, 2.0, 2.5 "DO" "BEGIN" NON EXP BESS K01(X, K0, K1); OUTPUT(61,"("/,4BZ.D,2(5B-.14D"-ZD)")", X, K0, K1) "END" "END" PRINTS THE FOLLOWING RESULTS: .5 .15241093857739" 1 .27310097082118" 1 1.0 .11444630798069" 1 .16361534862633" 1 1.5 .95821005329496" 0 .12431658735525" 1 2.0 .84156821507078" 0 .10334768470687" 1 2.5 .75954869032810" 0 .90017442390788" 0 SUBSECTION: NONEXP BESS K. CALLING SEQUENCE: THE HEADING OF THE PROCEDURE READS: "PROCEDURE" NONEXP BESS K(X, N,K); "VALUE" X, N; "INTEGER" N; "REAL" X; "ARRAY" K; "CODE" 35179; 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 15 THE MEANING OF THE FORMAL PARAMETERS IS: X: ; THE ARGUMENT OF THE BESSEL FUNCTIONS; X > 0; N: ; THE UPPER BOUND OF THE INDICES OF THE ARRAY K; N >= 0; K: ; "ARRAY" K[0:N]; EXIT: K[I] POSSESSES THE VALUE OF THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF ORDER I (I = 0, ..., N) MULTIPLIED BY EXP(X). PROCEDURES USED: NONEXP BESS K01 = CP 35178. REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE USED. METHOD AND PERFORMANCE: K[0] AND K[1] ARE COMPUTED BY USING NONEXP BESS K01 (CP 35178), WHILE K[2], ..., K[N] ARE COMPUTED ACCORDING TO THE RECURRENCE RELATION K[I+1]=K[I]+(2*I/X)*K[I], I>=2 (SEE [1], P. 376, FORMULA 9.6.26). EXAMPLE OF USE: THE PROGRAM "BEGIN" "REAL" X; "ARRAY" K[0:1]; "PROCEDURE" NONEXP BESS K(X, N, K); "CODE" 35179; "FOR" X:= .5, 1.0, 1.5, 2.0 "DO" "BEGIN" NONEXP BESS K(X, 2, K); OUTPUT(61, "("/,Z.D,3(5B.14D"D)")",X,K) "END" "END" PRINTS THE FOLLOWING RESULTS: .5 .15241093857739"1 .27310097082118"1 .12448148218621"2 1.0 .11444630798069"1 .16361534862633"1 .44167700523334"1 1.5 .95821005329496"0 .12431658735525"1 .26157645513649"1 2.0 .84156821507078"0 .10334768470687"U .18750450621395"1 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 16 SOURCE TEXT(S): "CODE" 35170; "REAL" "PROCEDURE" BESS I0(X); "VALUE" X; "REAL" X; "IF" X= 0 "THEN" BESS I0:=1 "ELSE" "IF" ABS(X) < = 15.0 "THEN" "BEGIN" "REAL" Z, DENOMINATOR, NUMERATOR; Z:= X*X; NUMERATOR:= (Z*(Z*(Z*(Z*(Z*(Z*(Z* (Z*(Z*(Z*(Z*(Z*(Z*(Z* .21058 07228 90567 "-22 +.38071 52423 45326 "-19) +.47944 02575 48300 "-16) +.43512 59712 62668 "-13) +.30093 11271 12960 "-10) +.16022 46793 95361 "-07) +.65485 83700 96785 "-05) +.20259 10841 43397 "-02) +.46307 62847 21000 "+00) +.75433 73289 48189 "+02) +.83079 25418 09429 "+04) +.57166 11305 63785 "+06) +.21641 55723 61227 "+08) +.35664 44822 44025 "+09) +.14404 82982 27235 "+10); DENOMINATOR:= (Z*(Z* (Z-.30764 69126 82801 "04) +.34762 63324 05882 "07) -.14404 82982 27235 "10); BESS I0:= -NUMERATOR/DENOMINATOR; "END" "ELSE" "BEGIN" "REAL" "PROCEDURE" NONEXP BESS I0(X); "CODE" 35175; BESS I0:= EXP(ABS(X)) * NONEXP BESS I0(X) "END" "EOP" 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 17 "CODE" 35171; "REAL" "PROCEDURE" BESS I1(X); "VALUE" X; "REAL" X; "IF" X=0 "THEN" BESS I1:=0 "ELSE" "IF" ABS(X) <= 15.0 "THEN" "BEGIN" "REAL" Z, DENOMINATOR, NUMERATOR; Z:= X*X; DENOMINATOR:= Z*(Z-.22258 36740 00860 "4) +.13629 35930 52499 "7; NUMERATOR:= (Z*(Z*(Z*(Z*(Z*(Z*(Z* (Z*(Z*(Z*(Z*(Z*(Z*(Z* .20717 57672 32792 "-26 +.25709 19055 84414 "-23) +.30627 92836 56135 "-20) +.26137 27721 58124 "-17) +.17846 93614 10091 "-14) +.96362 88915 18450 "-12) +.41006 89068 47159 "-09) +.13545 52288 41096 "-06) +.33947 28903 08516 "-04) +.62472 61951 27003 "-02) +.80614 48788 21295 "-00) +.68210 05679 80207 "+02) +.34106 97522 84422 "+04) +.84070 57728 77836 "+05) +.68146 79652 62502 "+06); BESS I1:= X*(NUMERATOR/DENOMINATOR) "END" "ELSE" "BEGIN" "REAL" "PROCEDURE" NONEXP BESS I1(X); "CODE" 35176; BESS I1:= EXP(ABS(X))*NONEXP BESS I1(X) "END"; "EOP" "CODE" 35172; "PROCEDURE" BESS I(X, N, I); "VALUE" X, N; "INTEGER" N; "REAL" X; "ARRAY" I; "IF" X = 0 "THEN" "BEGIN" I[0]:= 1; "FOR" N:= N "STEP" - 1 "UNTIL" 1 "DO" I[N]:= 0; "END" "ELSE" "BEGIN" "REAL" EXPX; "PROCEDURE" NONEXP BESS I(X, N, NONEXPI); "CODE" 35177; EXPX:= EXP(ABS(X)); NONEXP BESS I(X, N, I); "FOR" N:= N "STEP" - 1 "UNTIL" 0 "DO" I[N]:= I[N] * EXPX "END" BESS I; "EOP" 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 18 "CODE" 35173; "PROCEDURE" BESS K01(X, K0, K1); "VALUE" X; "REAL" X, K0, K1; "IF" X <= 1.5 "THEN" "BEGIN" "INTEGER" K; "REAL" C, D, R, S, SUM0, SUM1, T, TERM, T0, T1; SUM0:= D:= LN(2/X) -.5772156649015328606; SUM1:= C:= -1 -2 * D; R:= TERM:= 1; T:= X * X/4; "FOR" K:= 1,K+1 "WHILE" ABS(T0/SUM0) + ABS(T1/SUM1) > "-15 "DO" "BEGIN" TERM:= T * TERM * R * R; D:= D + R; C:= C - R; R:= 1/(K+1); C:= C - R; T0:= TERM * D; T1:= TERM * C * R; SUM0:= SUM0 + T0; SUM1:= SUM1 + T1 "END"; K0:= SUM0; K1:= (1 + T * SUM1) / X "END" "ELSE" "BEGIN" "REAL" EXPX; "PROCEDURE" NONEXP BESS K01(X, K0, K1); "CODE" 35178; EXPX:= EXP(- X); NONEXP BESS K01(X, K0, K1); K1:= EXPX * K1; K0:= K0 * EXPX "END" BESS K01; "EOP" "CODE" 35174; "PROCEDURE" BESS K(X, N, K); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" K; "BEGIN" "INTEGER" I; "REAL" K0, K1, K2; "PROCEDURE" BESS K01(X, K0, K1); "CODE" 35173; BESS K01(X, K0, K1); K[0]:= K0; "IF" N > 0 "THEN" K[1]:= K1; X:= 2 / X; "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO" "BEGIN" K[I]:= K2:= K0 + X * (I-1)* K1; K0:= K1; K1:= K2 "END" "END" BESS K; "EOP" "CODE" 35175; "REAL" "PROCEDURE" NONEXP BESS I0(X); "VALUE" X; "REAL" X; "IF" X= 0 "THEN" NONEXP BESS I0:=1 "ELSE" "IF" ABS(X) <= 15.0 "THEN" "BEGIN" "REAL" "PROCEDURE" BESS I0(X); "CODE" 35170; NONEXP BESS I0:= EXP(-ABS(X))*BESS I0(X) "END" "ELSE" "BEGIN" "REAL" SQRTX, AR, BR, BR1, BR2, Z, Z2, NUMERATOR, DENOMINATOR; X:=ABS(X); SQRTX:= SQRT(X); "COMMENT" 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 19 ; BR1:= BR2:= 0; Z:= 30/X-1; Z2:= Z+Z; "FOR" AR:= .24392 60769 778, -.11559 19781 04435 "3, +.78403 42490 05088 "4, -.14346 46313 13583 "6 "DO" "BEGIN" BR:= Z2*BR1-BR2+AR; BR2:= BR1; BR1:= BR "END"; NUMERATOR:= Z*BR1-BR2+.34651 98333 57379 "6; BR1:= BR2:= 0; "FOR" AR:= 1, -.32519 73333 69824 "3, +.20312 84361 00794 "5, -.36184 77792 19653 "6 "DO" "BEGIN" BR:= Z2*BR1 - BR2 + AR; BR2:= BR1; BR1:= BR "END"; DENOMINATOR:= Z*BR1 - BR2 +.86566 52748 32055 "6; NONEXP BESS I0:= (NUMERATOR/DENOMINATOR)/SQRTX; "END"; "EOP" "CODE" 35176; "REAL" "PROCEDURE" NONEXP BESS I1(X); "VALUE" X; "REAL" X; "IF" X=0 "THEN" NONEXP BESS I1:= 0 "ELSE" "IF" ABS(X)> 15.0 "THEN" "BEGIN" "INTEGER" SIGNX ; "REAL" AR, BR, BR1, BR2, Z, Z2, SQRTX, DENOMINATOR, NUMERATOR; SIGNX:= SIGN(X); X:= ABS(X); SQRTX:= SQRT(X); Z:= 30/X - 1; Z2 := Z + Z; BR1:= BR2:= 0; "FOR" AR:= +.14940 52814 740 "+1, -.36202 64202 42263 "+3, +.22054 97222 60336 "+5, -.40892 80849 44275 "+6 "DO" "BEGIN" BR:= Z2 * BR1 - BR2 + AR; BR2:= BR1; BR1:= BR "END"; NUMERATOR:= Z * BR1 -BR2 +.10277 66923 71524 "7; BR1:= BR2:= 0; "FOR" AR:= 1, -.63100 32005 51590 "3, +.49681 19495 33398 "5, -.10042 54281 33695 "7 "DO" "BEGIN" BR:= Z2 * BR1 - BR2 + AR; BR2:= BR1; BR1:=BR "END"; DENOMINATOR:= Z * BR1 - BR2 +.26028 87678 9105 "7; NONEXP BESS I1:= ((NUMERATOR/DENOMINATOR)/SQRTX) * SIGN X "END" "ELSE" "BEGIN" "REAL" "PROCEDURE" BESS I1(X); "CODE" 35171; NONEXP BESS I1:= EXP(-ABS(X))*BESS I1(X) "END"; "EOP" 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 20 "CODE" 35177; "PROCEDURE" NONEXP BESS I(X, N, I); "VALUE" X, N; "INTEGER" N; "REAL" X; "ARRAY" I; "IF" X = 0 "THEN" "BEGIN" I[0]:= 1; "FOR" N:= N "STEP" - 1 "UNTIL" 1 "DO" I[N]:= 0 "END" "ELSE" "BEGIN" "INTEGER" K; "REAL" X2, R, S; "BOOLEAN" NEGATIVE; "INTEGER" "PROCEDURE" START(X,N,T); "CODE" 35185; NEGATIVE:= X < 0; X:= ABS(X); R:= S:= 0; X2:= 2/X; K:= START(X,N,1); "FOR" K:= K "STEP" - 1 "UNTIL" 1 "DO" "BEGIN" R:= 1 / (R + X2 * K); S:= R * (2 + S); "IF" K <= N "THEN" I[K]:= R "END"; I[0]:= R:= 1 / (1 + S); "IF" NEGATIVE "THEN" "BEGIN" "FOR" K:= 1 "STEP" 1 "UNTIL" N "DO" I[K]:= R:= - R * I[K] "END" "ELSE" "FOR" K:=1 "STEP" 1 "UNTIL" N "DO" I[K]:= R:= R * I[K]; "END" NONEXP BESS I; "EOP" "CODE" 35178; "PROCEDURE" NONEXP BESS K01(X, K0, K1);"VALUE" X;"REAL" X, K0, K1; "IF" X <= 1.5 "THEN" "BEGIN" "REAL" EXPX; "PROCEDURE" BESS K01(X, K0, K1); "CODE" 35173; EXPX:= EXP(X); BESS K01(X, K0, K1); K0:= K0 * EXPX; K1:= EXPX * K1 "END" "ELSE" "IF" X <= 5 "THEN" "BEGIN" "INTEGER" R; "REAL" T2, FAC, S1, S2, TERM1, TERM2, SQRTEXPR, EXPH2, X2; S1:= .5; S2:=0; R:= 0; X2:= X + X; EXPH2:= 1 / SQRT(5 * X); "FOR" FAC:= .90483741803596, .67032004603564, .40656965974060, .20189651799466, .82084998623899"-1, .27323722447293"-1, .74465830709243"-2, .16615572731739"-2, .30353913807887"-3, .45399929762485"-4, .55595132416500"-5, .55739036926944"-6, .45753387694459"-7, .307487987958650"-8, .16918979226151"-9, .76218651945127"-11, .28111852987891"-12, .84890440338729"-14, .2098791048793"-15, .42483542552916"-17 "DO" "BEGIN" R:= R + 1; T2:= R * R / 10; SQRTEXPR:= SQRT(T2 / X2 + 1); TERM1:= FAC / SQRTEXPR; TERM2:= FAC * SQRTEXPR * T2; S1:= S1 + TERM1; S2:= S2 + TERM2 "END"; "COMMENT" 1SECTION : 6.9.2 (DECEMBER 1978) PAGE 21 ; K0:= EXPH2 * S1; K1:= EXPH2 * S2 * 2 "END" "ELSE" "BEGIN" "INTEGER" R; "REAL" BR, BR1, BR2, CR, CR1, CR2, DR, ERMIN1, ERPLUS1, ER, F0, F1, EXPX, Y, Y2; Y:= 10 / X - 1; Y2:= Y + Y; R:= 30; BR1:= BR2:= CR1:= CR2:= ERPLUS1:= ER:= 0; "FOR" DR:= .27545" - 15, -.172697" - 14, .1136042 " - 13, -.7883236 " -13, .58081063 " -12, -.457993622 " -11, .3904375576 " -10, -.36454717921 " - 9, .379299645568 " - 8, -.450473376411 " - 7, .63257510850049 " - 6, -.11106685196665" - 4, .26953261276272 " - 3, -.11310504646928" - 1 "DO" "BEGIN" R:= R - 2; BR:= Y2 * BR1 - BR2 + DR; CR:= CR1 * Y2 - CR2 + ER; ERMIN1:= R * DR + ERPLUS1; ERPLUS1:= ER; ER:= ERMIN1; BR2:= BR1; BR1:= BR; CR2:= CR1; CR1:= CR "END"; F0:= Y * BR1 - BR2 + .9884081742308258; F1:= Y * CR1 - CR2 + ER / 2; EXPX:= SQRT(1.5707963267949 / X); K0:= F0:= F0 * EXPX; K1:= (1 + .5 / X) * F0 + (10 / X / X) * EXPX * F1 "END" K0; "EOP" "CODE" 35179; "PROCEDURE" NONEXP BESS K(X, N, K); "VALUE" X, N; "REAL" X; "INTEGER" N; "ARRAY" K; "BEGIN" "INTEGER" I; "REAL" K0, K1, K2; "PROCEDURE" NONEXP BESS K01(X, K0, K1); "CODE" 35178; NONEXP BESS K01(X, K0, K1); K[0]:= K0; "IF" N> 0 "THEN" K[1]:= K1; X:= 2 / X; "FOR" I:= 2 "STEP" 1 "UNTIL" N "DO" "BEGIN" K[I]:= K2:= K0 + X * (I-1)* K1; K0:= K1; K1:= K2 "END" "END" NONEXP BESS K; "EOP" 1SECTION : 4.2.2 (DECEMBER 1979) PAGE 1 AUTHOR : P.W. HEMKER. CONTRIBUTOR : F.GROEN. INSTITUTE : MATHEMATICAL CENTRE. RECEIVED : 740620. BRIEF DESCRIPTION : TRICUB COMPUTES THE DEFINITE INTEGRAL OF A FUNCTION OF TWO VARIABLES OVER A TRIANGULAR DOMAIN. KEYWORDS : INTEGRATION, QUADRATURE, MOREDIMENSIONAL QUADRATURE, CUBATURE, DEFINITE INTEGRAL. CALLING SEQUENCE : THE HEADING OF THE PROCEDURE READS : "REAL" "PROCEDURE" TRICUB ( XI, YI, XJ, YJ, XK, YK, F, RE, AE ); "VALUE" XI, YI, XJ, YJ, XK, YK, RE, AE; "REAL" XI, YI, XJ, YJ, XK, YK, RE, AE; "REAL" "PROCEDURE" F; "CODE" 32075; TRICUB := THE COMPUTED VALUE OF THE DEFINITE INTEGRAL OF THE FUNCTION F ( X, Y ) OVER THE TRIANGULAR DOMAIN T WITH VERTICES ( XI, YI ), ( XJ, YJ ) AND ( XK, YK ). 1SECTION : 4.2.2 (DECEMBER 1979) PAGE 2 THE MEANING OF THE FORMAL PARAMETERS IS : XI, YI: ; ENTRY : THE COORDINATES OF THE FIRST VERTEX OF THE TRIANGULAR DOMAIN OF INTEGRATION; XJ, YJ: ; ENTRY : THE COORDINATES OF THE SECOND VERTEX OF THE TRIANGULAR DOMAIN OF INTEGRATION; XK, YK: ; ENTRY : THE COORDINATES OF THE THIRD VERTEX OF THE TRIANGULAR DOMAIN OF INTEGRATION; REMARK: THE ALGORITHM IS SYMMETRICAL IN THE VERTICES; THIS IMPLIES THAT THE RESULT OF THE PROCEDURE (ON ALL COUNTS) IS INVARIANT FOR ANY PERMUTATION OF THE VERTICES. F : ; THE HEADING OF THIS PROCEDURE READS: "REAL" "PROCEDURE" F ( X, Y ); "REAL" X, Y; THIS PROCEDURE DEFINES THE INTEGRAND; AE, RE: ; ENTRY: THE REQUIRED ABSOLUTE AND RELATIVE ERROR RESPECTIVELY. ONE SHOULD TAKE FOR "AE" AND "RE" VALUES WHICH ARE GREATER THAN THE ABSOLUTE AND RELATIVE ERROR IN THE COMPUTATION OF THE INTEGRAND F. PROCEDURES USED : NONE. REQUIRED CENTRAL MEMORY : THE PROCESS IS PROGRAMMED RECURSIVELY. AT EACH RECURSION LEVEL 43 REAL NUMBERS ARE USED. HOWEVER, FOR ANY PROPERLY CHOSEN VALUES OF RE AND AE THE RECURSION DEPTH WILL NOT EXCEED THE NUMBER OF BITS IN A REAL'S MANTISSA. RUNNING TIME : DEPENDS STRONGLY ON THE INTEGRAL TO COMPUTE. METHOD AND PERFORMANCE : A NESTED SEQUENCE OF CUBATURE RULES OF ORDER 2, 3, 4 AND 5 IS APPLIED. IF THE DIFFERENCE BETWEEN THE RESULT WITH THE 4-TH DEGREE RULE AND THE RESULT WITH THE 5-TH DEGREE RULE IS TOO LARGE, THEN THE TRIANGLE IS DIVIDED INTO FOUR CONGRUENT TRIANGLES. THIS PROCESS IS APPLIED RECURSIVELY IN ORDER TO OBTAIN AN ADAPTIVE CUBATURE ALGORITHM. 1SECTION : 4.2.2 (OCTOBER 1975) PAGE 3 REFERENCES : [1].P.W.HEMKER, A SEQUENCE OF NESTED CUBATURE RULES, MATH.CENTRE, AMSTERDAM, REPORT NW 3/73. EXAMPLE OF USE: THE FOLLOWING PROGRAM EVALUATES THE INTEGRAL OF F ( X, Y ) = COS ( X ) * COS ( Y ) OVER THE TRIANGLE T WITH VERTICES ( 0, 0 ), ( O, PI/2 ) AND ( PI/2, PI/2 ). ON EACH LINE ARE LISTED : A: THE REQUIRED RELATIVE AND ABSOLUTE PRECISION; B: THE COMPUTED VALUE OF THE INTEGRAL; C: THE NUMBER OF CALLS OF THE FUNCTION F, "BEGIN" "REAL" "PROCEDURE" TRICUB(A,B,C,D,E,F,G,H,I); "CODE" 32075; "INTEGER" N,C,I,K; "REAL" PI,ACC,R,S; "REAL" "PROCEDURE" E(X,Y); "REAL" X,Y; "BEGIN" C:= C+1; "IF" C> 20000 "THEN" "GOTO" CC; E:= COS(X) * COS(Y); "END" E; PI:= 3.14159265359; "FOR" ACC:= "-1,"-2,"-3,"-4,"-5,"-6,"-7,"-8,"-9,"-10,"-11 "DO" "BEGIN" C:=0; OUTPUT(61,"("+.D"+ZD,2B,+.14D"+2ZD,2B,10ZD,/")", ACC, TRICUB(0,0,0,PI/2,PI/2,PI/2,E,ACC,ACC) ,C); "END"; CC: OUTPUT(61,"("*")"); "END" RESULTS: +.1" +0 +.50063973801970" +0 7 +.1" -1 +.50063973801970" +0 7 +.1" -2 +.50063973801970" +0 7 +.1" -3 +.49999110261504" +0 10 +.1" -4 +.49999848959226" +0 13 +.1" -5 +.49999848959226" +0 13 +.1" -6 +.49999997378240" +0 43 +.1" -7 +.49999999209792" +0 133 +.1" -8 +.49999999893172" +0 313 +.1" -9 +.49999999985571" +0 733 +.1"-10 +.49999999998692" +0 1723 1SECTION : 4.2.2 (DECEMBER 1979) PAGE 4 SOURCE TEXT(S): 0"CODE" 32075; "REAL" "PROCEDURE" TRICUB(XI,YI,XJ,YJ,XK,YK,G,RE,AE); "VALUE" XI,YI,XJ,YJ,XK,YK,RE,AE; "REAL" XI,YI,XJ,YJ,XK,YK,RE,AE; "REAL" "PROCEDURE" G; "BEGIN" "REAL" SURF,SURFMIN,XZ,YZ,GI,GJ,GK; "REAL" "PROCEDURE" INT(AX1,AY1,AF1,AX2,AY2,AF2,AX3,AY3,AF3, BX1,BY1,BF1,BX2,BY2,BF2,BX3,BY3,BF3, PX,PY,PF); "VALUE" BX1,BY1,BF1,BX2,BY2,BF2,BX3,BY3,BF3,PX,PY,PF; "REAL" BX1,BY1,BF1,BX2,BY2,BF2,BX3,BY3,BF3,PX,PY,PF, AX1,AY1,AF1,AX2,AY2,AF2,AX3,AY3,AF3; "BEGIN" "REAL" E,I3,I4,I5,A,B,C,SX1,SY1,SX2,SY2,SX3,SY3, CX1,CY1,CF1,CX2,CY2,CF2,CX3,CY3,CF3, DX1,DY1,DF1,DX2,DY2,DF2,DX3,DY3,DF3; A:= AF1 + AF2 + AF3; B:= BF1 + BF2 + BF3; I3:= 3 * A + 27 * PF + 8 * B; E:= ABS(I3) * RE + AE; "IF" SURF < SURFMIN "OR" ABS(5 * A + 45 * PF - I3) < E "THEN" INT:= I3 * SURF "ELSE" "BEGIN" CX1:= AX1 + PX; CY1:= AY1 + PY; CF1:= G(CX1,CY1); CX2:= AX2 + PX; CY2:= AY2 + PY; CF2:= G(CX2,CY2); CX3:= AX3 + PX; CY3:= AY3 + PY; CF3:= G(CX3,CY3); C:= CF1 + CF2 + CF3; I4:= A + 9 * PF + 4 * B + 12 * C; "IF" ABS(I3 - I4) < E "THEN" INT:= I4 * SURF "ELSE" "BEGIN" SX1:= .5 * BX1; SY1:= .5 * BY1; DX1:= AX1 + SX1; DY1:= AY1 + SY1; DF1:= G(DX1,DY1); SX2:= .5 * BX2; SY2:= .5 * BY2; DX2:= AX2 + SX2; DY2:= AY2 + SY2; DF2:= G(DX2,DY2); SX3:= .5 * BX3; SY3:= .5 * BY3; DX3:= AX3 + SX3; DY3:= AY3 + SY3; DF3:= G(DX3,DY3); I5:= (51 * A + 2187 * PF + 276 * B + 972 * C - 768 * (DF1 + DF2 + DF3))/63; "COMMENT" 1SECTION : 4.2.2 (OCTOBER 1975) PAGE 5 ; "IF" ABS(I4 - I5) < E "THEN" INT:= I5 * SURF "ELSE" "BEGIN" SURF:= .25 * SURF; INT:= INT(SX1,SY1,BF1,SX2,SY2,BF2,SX3,SY3,BF3, DX1,DY1,DF1,DX2,DY2,DF2,DX3,DY3,DF3, PX,PY,PF) + INT(AX1,AY1,AF1,SX3,SY3,BF3,SX2,SY2,BF2,DX1,DY1,DF1, AX1 + SX2,AY1 + SY2,G(AX1 + SX2,AY1 + SY2), AX1 + SX3,AY1 + SY3,G(AX1 + SX3,AY1 + SY3), .5 * CX1,.5 * CY1,CF1) + INT(AX2,AY2,AF2,SX3,SY3,BF3,SX1,SY1,BF1,DX2,DY2,DF2, AX2 + SX1,AY2 + SY1,G(AX2 + SX1,AY2 + SY1), AX2 + SX3,AY2 + SY3,G(AX2 + SX3,AY2 + SY3), .5 * CX2,.5 * CY2,CF2) + INT(AX3,AY3,AF3,SX1,SY1,BF1,SX2,SY2,BF2,DX3,DY3,DF3, AX3 + SX2,AY3 + SY2,G(AX3 + SX2,AY3 + SY2), AX3 + SX1,AY3 + SY1,G(AX3 + SX1,AY3 + SY1), .5 * CX3,.5 * CY3,CF3); SURF:= 4 * SURF "END" "END" "END" "END" INT; SURF:= 0.5 * ABS(XJ * YK - XK * YJ + XI * YJ - XJ * YI + XK * YI - XI * YK); SURFMIN:= SURF*RE; RE:= 30*RE; AE:= 30*AE/SURF; XZ:= (XI + XJ + XK)/3; YZ:= (YI + YJ + YK)/3; GI:= G(XI,YI); GJ:= G(XJ,YJ); GK:= G(XK,YK); XI:= XI*.5; YI:= YI*.5; XJ:= XJ*.5; YJ:= YJ*.5; XK:= XK*.5; YK:= YK*.5; TRICUB:= INT(XI,YI,GI,XJ,YJ,GJ,XK,YK,GK, XJ+XK,YJ+YK,G(XJ+XK,YJ+YK), XK+XI,YK+YI,G(XK+XI,YK+YI), XI+XJ,YI+YJ,G(XI+XJ,YI+YJ), .5 * XZ,.5 * YZ,G(XZ,YZ))/60 "END" TRICUB; "EOP" 1SECTION : 5.2.1.3.1 (FEBRUARY 1979) PAGE 1 AUTHOR : B. VAN DOMSELAAR. INSTITUTE: MATHEMATICAL CENTRE, AMSTERDAM. RECEIVED: 750601. BRIEF DESCRIPTION: PEIDE ESTIMATES UNKNOWN VARIABLES IN A SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS; THE UNKNOWN VARIABLES MAY APPEAR NONLINEAR BOTH IN THE DIFFERENTIAL EQUATIONS AND ITS INITIAL VALUES; A SET OF OBSERVED VALUES OF SOME COMPONENTS OF THE SOLUTION OF THE DIFFERENTIAL EQUATIONS MUST BE GIVEN; KEYWORDS: PARAMETER ESTIMATION, DIFFERENTIAL EQUATIONS, INITIAL VALUE PROBLEM, DATA FITTING. CALLING SEQUENCE: THE HEADING OF THIS PROCEDURE IS: "PROCEDURE" PEIDE(N, M, NOBS, NBP, PAR, RV, BP, JTJINV, IN, OUT, DERIV, JAC DFDY, JACDFDP, CALL YSTART, DATA, MONITOR); "VALUE" N,M,NOBS; "INTEGER" N,M,NOBS,NBP; "ARRAY" PAR,RV,JTJINV,IN,OUT; "INTEGER" "ARRAY" BP; "PROCEDURE" CALL YSTART,DATA,MONITOR; "BOOLEAN" "PROCEDURE" DERIV,JAC DFDY,JAC DFDP; "CODE" 34444; THE MEANING OF THE FORMAL PARAMETERS IS: N: ; THE NUMBER OF DIFFERENTIAL EQUATIONS; M: ; THE NUMBER OF UNKNOWN VARIABLES; NOBS: ; THE NUMBER OF OBSERVATIONS; NOBS SHOULD SATISFY NOBS>=M; NBP: ; ENTRY: THE NUMBER OF BREAK-POINTS; IF NO BREAK-POINTS ARE USED THEN SET NBP=0; EXIT: WITH NORMAL TERMINATION OF THE PROCESS NBP=0; OTHERWISE, IF THE PROCESS HAS BEEN BROKEN OFF (SEE OUT[1]), THE VALUE OF NBP IS THE NUMBER OF BREAK- POINTS USED BEFORE THE PROCESS BROKE OFF; 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 2 PAR: ; "ARRAY" PAR[1 : M+NBP]; ENTRY: PAR[1:M] SHOULD CONTAIN AN INITIAL APPROXIMATION TO THE REQUIRED PARAMETER VECTOR; EXIT: PAR[1:M] CONTAINS THE CALCULATED PARAMETER VECTOR; IF OUT[1]>0 AND NBP>0 THEN PAR[M+1:M+NBP] CONTAINS THE VALUES OF THE NEWLY INTRODUCED PARAMETERS BEFORE THE PROCESS BROKE OFF; RV: ; "ARRAY" RV[1 : NOBS+NBP]; EXIT: RV[1:NOBS] CONTAINS THE RESIDUAL VECTOR AT THE CALCULATED MINIMUM; IF OUT[1]>0 AND NBP>0 THEN RV[NOBS+1:NOBS+NBP] CONTAINS THE ADDITIONAL CONTINUITY REQUIREMENTS AT THE BREAK-POINTS BEFORE THE PROCESS BROKE OFF; BP: ; "INTEGER" "ARRAY" BP[0 : NBP]; ENTRY: BP[I], I=1,...,NBP, SHOULD CORRESPOND TO THE INDEX OF THAT TIME OF OBSERVATION WHICH WILL BE USED AS A BREAK-POINT (1<=BP[I]<=NOBS); THE BREAK-POINTS HAVE TO BE ORDERED SUCH THAT BP[I]<=BP[J] IF I<=J; EXIT: WITH NORMAL TERMINATION OF THE PROCESS BP[1:NBP] CONTAINS NO INFORMATION; OTHERWISE, IF OUT[1]>0 AND NBP>0 THEN BP[I], I=1,...,NBP, CONTAINS THE INDEX OF THAT TIME OF OBSERVATION WHICH WAS USED AS A BREAK-POINT BEFORE THE PROCESS BROKE OFF; JTJINV: ; "ARRAY" JTJINV[1 : M, 1 : M]; EXIT: THE INVERSE OF THE MATRIX J' * J WHERE J DENOTES THE MATRIX OF PARTIAL DERIVATIVES DRV[I] / DPAR[K] (I=1,...,NOBS ; K=1,...,M) AND J' DENOTES THE TRANSPOSE OF J; THIS MATRIX CAN BE USED IF ADDITIONAL INFORMATION ABOUT THE RESULT IS REQUIRED; E.G. STATISTICAL DATA SUCH AS THE COVARIANCE MATRIX, CORRELATION MATRIX AND CONFIDENCE INTERVALS CAN EASILY BE CALCULATED FROM JTJINV AND OUT[2]; IN: ; "ARRAY" IN[0 : 6]; ENTRY: IN THIS ARRAY THE USER SHOULD GIVE SOME DATA TO CONTROL THE PROCESS; IN[0]: THE MACHINE PRECISION; FOR THE CYBER 73 A SUITABLE VALUE IS "-14; IN[1]: THE RATIO: THE MINIMAL STEPLENGTH FOR THE INTEGRATION OF THE DIFFERENTIAL EQUATIONS DIVIDED BY THE DISTANCE BETWEEN TWO NEIGHBOURING OBSERVATIONS; MOSTLY, A SUITABLE VALUE IS "-4; IN[2]: THE RELATIVE LOCAL ERROR BOUND FOR THE INTEGRATION PROCESS; THIS VALUE SHOULD SATISFY IN[2]<=IN[3]; THIS PARAMETER CONTROLS THE ACCURACY OF THE NUMERICAL INTEGRATION; MOSTLY, A SUITABLE VALUE IS IN[3]/100; 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 3 IN[3], IN[4]: THE RELATIVE AND THE ABSOLUTE TOLERANCE FOR THE DIFFERENCE BETWEEN THE EUCLIDEAN NORM OF THE ULTIMATE AND PENULTIMATE RESIDUAL VECTOR RESPECTIVELY; THE PROCESS IS TERMINATED IF THE IMPROVEMENT OF THE SUM OF SQUARES IS LESS THAN IN[3] * (SUM OF SQUARES) + IN[4] * IN[4]; THESE TOLERANCES SHOULD BE CHOSEN IN ACCORDANCE WITH THE RELATIVE, RESP. ABSOLUTE ERRORS IN THE OBSERVATIONS; NOTE THAT THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR IS DEFINED AS THE SQUARE ROOT OF THE SUM OF SQUARES; IN[5]: THE MAXIMUM NUMBER OF TIMES THAT THE INTEGRATION OF THE DIFFERENTIAL EQUATIONS IS PERFORMED; IN[6]: A STARTING VALUE USED FOR THE RELATION BETWEEN THE GRADIENT AND THE GAUSS-NEWTON DIRECTION (SEE [1]); IF THE PROBLEM IS WELL CONDITIONED THEN A SUITABLE VALUE FOR IN[6] WILL BE 0.01; IF THE PROBLEM IS ILL CONDITIONED THEN IN[6] SHOULD BE GREATER, BUT THE VALUE OF IN[6] SHOULD SATISFY: IN[0] < IN[6] <= 1/IN[0]; OUT: ; "ARRAY" OUT[1 : 7]; EXIT : IN ARRAY OUT SOME BY-PRODUCTS ARE DELIVERED; OUT[1]: THIS VALUE GIVES INFORMATION ABOUT THE TERMINATION OF THE PROCESS; OUT[1]=0: NORMAL TERMINATION; IF OUT[1]>0 THEN THE PROCESS HAS BEEN BROKEN OFF AND THIS MAY OCCUR BECAUSE OF THE FOLLOWING REASONS: OUT[1]=1: THE NUMBER OF INTEGRATIONS PERFORMED EXCEEDED THE NUMBER GIVEN IN IN[5]; OUT[1]=2: THE DIFFERENTIAL EQUATIONS ARE VERY NONLINEAR; DURING AN INTEGRATION THE VALUE OF IN[1] WAS DECREASED BY A FACTOR 10000 AND IT IS ADVISED TO DECREASE IN[1], ALTHOUGH THIS WILL INCREASE COMPUTING TIME; OUT[1]=3: A CALL OF DERIV DELIVERED THE VALUE FALSE; OUT[1]=4: A CALL OF JAC DFDY DELIVERED THE VALUE FALSE; OUT[1]=5: A CALL OF JAC DFDP DELIVERED THE VALUE FALSE; OUT[1]=6: THE PRECISION ASKED FOR CAN NOT BE ATTAINED; THIS PRECISION IS POSSIBLY CHOSEN TOO HIGH, RELATIVE TO THE PRECISION IN WHICH THE RESIDUAL VECTOR IS CALCULATED (SEE IN[3]); OUT[2]: THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR CALCULATED WITH VALUES OF THE UNKNOWNS DELIVERED; 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 4 OUT[3]: THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR CALCULATED WITH THE INITIAL VALUES OF THE UNKNOWN VARIABLES; OUT[4]: THE NUMBER OF INTEGRATIONS PERFORMED, NEEDED TO OBTAIN THE CALCULATED RESULT; IF OUT[4]=1 AND OUT[1]>0 THEN THE MATRIX JTJINV CAN NOT BE USED; OUT[5]: THE MAXIMUM NUMBER OF TIMES THAT THE REQUESTED LOCAL ERROR BOUND WAS EXCEEDED IN ONE INTEGRATION; IF IT IS A LARGE NUMBER, IT MAY BE BETTER TO DECREASE THE VALUE OF IN[1]; OUT[6]: THE IMPROVEMENT OF THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR IN THE LAST ITERATION STEP OF THE PROCESS OF MARQUARDT; OUT[7]: THE CONDITION NUMBER OF J' * J , I.E. THE RATIO OF ITS LARGEST TO SMALLEST EIGENVALUES; DERIV: ; THIS PROCEDURE DEFINES THE RIGHT HAND SIDE OF THE DIFFERENTIAL EQUATIONS; THE HEADING OF THIS PROCEDURE SHOULD BE: "BOOLEAN" "PROCEDURE" DERIV(PAR, Y, T, DF); "VALUE" T; "REAL" T; "ARRAY" PAR,Y,DF; ENTRY: PAR,Y,T; PAR[1:M] CONTAINS THE CURRENT VALUES OF THE UNKNOWNS AND SHOULD NOT BE ALTERED; Y[1:N] CONTAINS THE SOLUTIONS OF THE DIFFERENTIAL EQUATIONS AT TIME T AND SHOULD NOT BE ALTERED; EXIT: "ARRAY" DF[1 : N]; AN ARRAY ELEMENT DF[I] SHOULD CONTAIN THE RIGHT HAND SIDE OF THE I-TH DIFFERENTIAL EQUATION; AFTER A SUCCESSFUL CALL OF DERIV, THE BOOLEAN PROCEDURE SHOULD DELIVER THE VALUE TRUE; HOWEVER, IF DERIV DELIVERS THE VALUE FALSE, THEN THE PROCESS IS TERMINATED (SEE OUT[1]); HENCE, PROPER PROGRAMMING OF DERIV MAKES IT POSSIBLE TO AVOID CALCULATION OF THE RIGHT HAND SIDE WITH VALUES OF THE UNKNOWN VARIABLES WHICH CAUSE OVERFLOW IN THE COMPUTATION; JAC DFDY: ; THE HEADING OF THIS PROCEDURE SHOULD BE: "BOOLEAN" "PROCEDURE" JAC DFDY(PAR, Y, T, FY); "VALUE" T; "REAL" T; "ARRAY" PAR,Y,FY; ENTRY: PAR,Y,T; SEE DERIV; EXIT: "ARRAY" FY[1 : N,1 : N]; AN ARRAY ELEMENT FY[I,J] SHOULD CONTAIN THE PARTIAL DERIVATIVE OF THE RIGHT HAND SIDE OF THE I-TH DIFFERENTIAL EQUATION WITH RESPECT TO Y[J], I.E. DF[I]/DY[J]; THE BOOLEAN VALUE SHOULD BE ASSIGNED TO THIS PROCEDURE IN THE SAME WAY AS IT IS DONE FOR THE VALUE OF DERIV; JAC DFDP: ; THE HEADING OF THIS PROCEDURE SHOULD BE: "BOOLEAN" "PROCEDURE" JAC DFDP(PAR, Y, T, FP); "VALUE" T; "REAL" T; "ARRAY" PAR,Y,FP; 1SECTION : 5.2.1.3.1 (FEBRUARY 1979) PAGE 5 ENTRY: PAR,Y,T; SEE DERIV; EXIT: "ARRAY" FP[1 : N,1 : M]; AN ARRAY ELEMENT FP[I,J] SHOULD CONTAIN THE PARTIAL DERIVATIVE OF THE RIGHT HAND SIDE OF THE I-TH DIFFERENTIAL EQUATION WITH RESPECT TO PAR[J], I.E. DF[I]/DPAR[J]; THE BOOLEAN VALUE SHOULD BE ASSIGNED TO THIS PROCEDURE IN THE SAME WAY AS IT IS DONE FOR THE VALUE OF DERIV; CALL YSTART: ; THIS PROCEDURE DEFINES THE INITIAL VALUES OF THE INITIAL VALUE PROBLEM; THE HEADING OF THIS PROCEDURE SHOULD BE: "BOOLEAN" "PROCEDURE" CALL YSTART(PAR, Y, YMAX); "ARRAY" PAR,Y,YMAX; ENTRY: PAR; PAR[1:M] CONTAINS THE CURRENT VALUES OF THE UNKNOWN VARIABLES AND SHOULD NOT BE ALTERED; EXIT: Y,YMAX; Y[1:N] SHOULD CONTAIN THE INITIAL VALUES OF THE CORRESPONDING DIFFERENTIAL EQUATIONS; THE INITIAL VALUES MAY BE FUNCTIONS OF THE UNKNOWN VARIABLES PAR; IN THAT CASE, THE INITIAL VALUES OF DY/DPAR ALSO HAVE TO BE SUPPLIED; NOTE THAT DY[I]/DPAR[J] CORRESPONDS WITH Y[5*N+J*N+I] (I=1,...,N , J=1,...,M); YMAX[I], I=1,...,N, SHOULD CONTAIN A ROUGH ESTIMATE TO THE MAXIMAL ABSOLUTE VALUE OF Y[I] OVER THE INTEGRATION INTERVAL; DATA: ; THIS PROCEDURE TAKES THE DATA TO FIT INTO THE PROCEDURE PEIDE; THE HEADING OF THIS PROCEDURE SHOULD BE: "PROCEDURE" DATA(NOBS, TOBS, OBS, COBS); "VALUE" NOBS; "INTEGER" NOBS; "ARRAY" TOBS,OBS; "INTEGER" "ARRAY" COBS; ENTRY: NOBS; NOBS HAS THE SAME MEANING AS IN PEIDE; EXIT: "ARRAY" TOBS[0 : NOBS]; THE ARRAY ELEMENT TOBS[0] SHOULD CONTAIN THE TIME, CORRESPONDING TO THE INITIAL VALUES OF Y GIVEN IN THE PROCEDURE CALL YSTART; AN ARRAY ELEMENT TOBS[I], 1<=I<=NOBS, SHOULD CONTAIN THE I-TH TIME OF OBSERVATION; THE OBSERVATIONS HAVE TO BE ORDERED SUCH THAT TOBS[I]<=TOBS[J] IF I<=J; "INTEGER" "ARRAY" COBS[1:NOBS]; AN ARRAY ELEMENT COBS[I] SHOULD CONTAIN THE COMPONENT OF Y OBSERVED AT TIME TOBS[I]; NOTE THAT 1<=COBS[I]<=N; "ARRAY" OBS[1:NOBS]; AN ARRAY ELEMENT OBS[I] SHOULD CONTAIN THE OBSERVED VALUE OF THE COMPONENT COBS[I] OF Y AT THE TIME TOBS[I]; 1SECTION : 5.2.1.3.1 (FEBRUARY 1979) PAGE 6 MONITOR: ; THIS PROCEDURE CAN BE USED TO OBTAIN INFORMATION ABOUT THE COURSE OF THE ITERATION PROCESS; IF NO INTERMEDIATE RESULTS ARE DESIRED, A DUMMY PROCEDURE SATISFIES; THE HEADING OF THIS PROCEDURE SHOULD BE: "PROCEDURE" MONITOR(POST,NCOL,NROW,PAR,RV,WEIGHT,NIS); "VALUE" POST,NCOL,NROW,WEIGHT,NIS; "INTEGER" POST,NCOL,NROW,WEIGHT,NIS; "ARRAY" PAR,RV; INSIDE PEIDE, THE PROCEDURE MONITOR IS CALLED AT TWO DIFFERENT PLACES AND THIS IS DENOTED BY THE VALUE OF POST: POST=1: MONITOR IS CALLED AFTER AN INTEGRATION OF THE DIFFERENTIAL EQUATIONS; AT THIS PLACE ARE AVAILABLE: THE CURRENT VALUES OF THE UNKNOWN VARIABLES PAR[1:NCOL], WHERE NCOL=M+NBP, THE CALCULATED RESIDUAL VECTOR RV[1:NROW], WHERE NROW=NOBS+NBP, AND THE VALUE OF NIS, WHICH IS THE NUMBER OF INTEGRATION STEPS PERFORMED DURING THE SOLUTION OF THE LAST INITIAL VALUE PROBLEM; POST=2: MONITOR IS CALLED BEFORE A MINIMIZATION OF THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR WITH THE PROCEDURE MARQUARDT (SEE SECTION 5.1.3.1.3) IS STARTED; AVAILABLE ARE THE CURRENT VALUES OF PAR[1:NCOL] AND THE VALUE OF THE WEIGHT, WITH WHICH THE CONTINUITY REQUIREMENTS AT THE BREAK- POINTS ARE ADDED TO THE ORIGINAL LEAST SQUARES PROBLEM. DATA AND RESULTS: SEE REF[1]. PROCEDURES USED: INIVEC = CP31010, INIMAT = CP31011, MULVEC = CP31020, MULROW = CP31021, DUPVEC = CP31030, DUPMAT = CP31035, VECVEC = CP34010, MATVEC = CP34011, ELMVEC = CP34020, SOL = CP34051, DEC = CP34300, MARQUARDT = CP34440. 1SECTION : 5.2.1.3.1 (FEBRUARY 1979) PAGE 7 REQUIRED CENTRAL MEMORY : IN THE BODY OF PEIDE (3 + NBP) * NOBS + N * (13 + N + 7 * M + 7 * NBP) ARRAY ELEMENTS ARE DECLARED. METHOD AND PERFORMANCE: PEIDE ESTIMATES UNKNOWN VARIABLES IN THE SYSTEM OF DIFFERENTIAL EQUATIONS DY/DT (T, PAR) = F (T, Y, PAR), BY USING A SET OF OBSERVED VALUES OF Y; THE UNKNOWN VARIABLES PAR ARE OBTAINED IN THE LEAST SQUARES SENSE; AN ELEMENT OF THE RESIDUAL VECTOR IS DEFINED BY THE CALCULATED VALUE OF Y MINUS ITS OBSERVED VALUE; THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR IS MINIMIZED BY THE ITERATION PROCESS OF MARQUARDT; THE DIFFERENTIAL EQUATIONS ARE SOLVED BY THE INTEGRATION PROCESS OF GEAR; A MULTIPLE SHOOTING TECHNIQUE HAS BEEN IMPLEMENTED TO IMPROVE BAD STARTING VALUES OF THE UNKNOWNS; IF THIS TECHNIQUE IS USED, ONE HAS TO GIVE SOME BREAK-POINTS, I.E. TIMES OF OBSERVATIONS WHERE A NEW INITIAL VALUE PROBLEM SHOULD BE STARTED; THE NEW INITIAL VALUES OF Y BECOME EXTRA UNKNOWN VARIABLES AND THE CONTINUITY REQUIREMENTS AT THE BREAK-POINTS ARE ADDED WITH SOME WEIGHTING FACTOR TO THE LEAST SQUARES PROBLEM; FOR DETAILS SEE REF[1]. REFERENCES: [1]: B. VAN DOMSELAAR, NONLINEAR PARAMETER ESTIMATION IN INITIAL VALUE PROBLEMS, MATH. CENTRE, AMSTERDAM (TO APPEAR). EXAMPLE OF USE: THE PARAMETERS PAR[1:3] IN THE DIFFERENTIAL EQUATIONS DY[1]/DT = - (1 - Y[2]) * Y[1] + EXP(PAR[2]) * Y[2], DY[2]/DT = EXP(PAR[1]) * ((1 - Y[2]) * Y[1] - (EXP(PAR[2])+ +EXP(PAR[3])) * Y[2]), WITH 23 OBSERVATIONS OF Y[2], MAY BE OBTAINED BY THE FOLLOWING PROGRAM, THAT CONSISTS OF 1: A CODE PROCEDURE WHICH TAKES CARE OF THE OUTPUT OF THE EXAMPLE PROGRAM. IT ALSO INTERPRETS THE NUMERICAL DATA THAT CAN BE USED TO OBTAIN STATISTICAL RESULTS; 2: THE USERS PROGRAM IN WHICH THE PROBLEM EXAMPLE IS DEFINED. 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 8 "CODE" 34445; "PROCEDURE" COMMUNICATION(POST,FA,N,M,NOBS,NBP,PAR,RES,BP,JTJINV, IN,OUT,WEIGHT,NIS); "VALUE" POST,FA,N,M,NOBS,NBP,WEIGHT,NIS; "INTEGER" POST,N,M,NOBS,NBP,WEIGHT,NIS; "REAL" FA; "ARRAY" PAR,RES,BP,JTJINV,IN,OUT; "BEGIN" "INTEGER" I,J; "REAL" C; "ARRAY" CONF[1:M]; "REAL" "PROCEDURE" VECVEC(L,U,S,A,B); "CODE" 34010; "IF" POST=5 "THEN" "BEGIN" OUTPUT(61,"("*,/,10B,"("THE FIRST RESIDUAL VECTOR")",//,16B, "("I")",4B,"("RES[I]")",/")"); "FOR" I:=1 "STEP" 1 "UNTIL" NOBS "DO" OUTPUT(61,"("15B,ZD,2B,+.4D"+ZD,/")",I,RES[I]); "END" "ELSE" "IF" POST=3 "THEN" "BEGIN" OUTPUT(61,"("*,/, "("THE EUCLIDEAN NORM OF THE RESIDUAL VECTOR:")", .7D"+ZD,2/,5B,"("CALCULATED PARAMETERS")",/")", SQRT(VECVEC(1,NOBS,0,RES,RES))); "FOR" I:=1 "STEP" 1 "UNTIL" M "DO" OUTPUT(61,"("9B,+.7D"+ZD,/")",PAR[I]); OUTPUT(61,"("/, "("NUMBER OF INTEGRATION STEPS PERFORMED: ")",ZZD,//")",NIS); "END" "ELSE" "IF" POST=4 "THEN" "BEGIN" "IF" NBP=0 "THEN" OUTPUT(61,"("*,//,5B, "("THE MINIMIZATION IS STARTED WITHOUT BREAK-POINTS")"")") "ELSE" "BEGIN" OUTPUT(61,"("*,5/,20B, "("THE MINIMIZATION IS STARTED WITH W E I G H T =")",ZD, 3/")",WEIGHT); OUTPUT(61,"("/,5B, "("THE EXTRA PARAMETERS ARE THE OBSERVATIONS:")"")"); "FOR" I:=1 "STEP" 1 "UNTIL" NBP "DO" OUTPUT(61,"("8B,ZD,2B")",BP[I]); "END"; OUTPUT(61,"("6/,10B, "("STARTING VALUES OF THE PARAMETERS")",/")"); "FOR" I:=1 "STEP" 1 "UNTIL" M "DO" OUTPUT(61,"("20B,+.7D"+ZD,/")",PAR[I]); OUTPUT(61,"("//, "("REL. TOLERANCE FOR THE EUCL. NORM OF THE RES. VECTOR:")" ,B,.7D"+ZD,/, "("ABS. TOLERANCE FOR THE EUCL. NORM OF THE RES. VECTOR:")" ,B,.7D"+ZD,/,"("RELATIVE STARTING VALUE OF LAMBDA")",19B, "(":")",B,.7D"+ZD")",IN[3],IN[4],IN[6]) "END" "ELSE" "IF" POST=1 "THEN" "BEGIN" OUTPUT(61,"("10B,"("STARTING VALUES OF THE PARAMETERS")",/")"); "FOR" I:=1 "STEP" 1 "UNTIL" M "DO" OUTPUT(61,"("20B,+.7D"+ZD,/")",PAR[I]); "COMMENT" 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 9 ; OUTPUT(61,"("2/,"("NUMBER OF EQUATIONS")",3B,"(":")",ZD,/, "("NUMBER OF OBSERVATIONS:")",ZD,2/, "("MACHINE PRECISION")",30B,"(":")",+.D"+ZD,/, "("RELATIVE LOCAL ERROR BOUND FOR INTEGRATION")",5B,"(":")",+.D"+ZD,/, "("RELATIVE TOLERANCE FOR RESIDUE")",17B,"(":")",+.2D"+ZD,/, "("ABSOLUTE TOLERANCE FOR RESIDUE")",17B,"(":")",+.2D"+ZD,/, "("MAXIMUM NUMBER OF INTEGRATIONS TO PERFORM")",6B,"(":")",ZZD,/, "("RELATIVE STARTING VALUE OF LAMBDA")",14B,"(":")",+.2D"+ZD,/, "("RELATIVE MINIMAL STEPLENGTH")",20B,"(":")",+.2D"+ZD,/")", N,NOBS,IN[0],IN[2],IN[3],IN[4],IN[5],IN[6],IN[1]); "IF" NBP=0 "THEN" OUTPUT(61,"("//, "("THERE ARE NO BREAK-POINTS")"")") "ELSE" "BEGIN" OUTPUT(61,"("//, "("BREAK-POINTS ARE THE OBSERVATIONS :")"")"); "FOR" I:=1 "STEP" 1 "UNTIL" NBP "DO" OUTPUT(61,"("ZZD,B")",BP[I]) "END"; OUTPUT(61,"("//, "("THE ALPHA-POINT OF THE F-DISTIBUTION :")", ZD.DD")",FA); "END" "ELSE" "IF" POST=2 "THEN" "BEGIN" OUTPUT(61,"("*")"); "IF" OUT[1]=0 "THEN" OUTPUT(61,"("2/, "("NORMAL TERMINATION OF THE PROCESS")"")") "ELSE" "IF" OUT[1]=1 "THEN" OUTPUT(61,"("2/, "("NUMBER OF INTEGRATIONS ALLOWED WAS EXCEEDED")"")") "ELSE" "IF" OUT[1]=2 "THEN" OUTPUT(61,"("2/, "("MINIMAL STEPLENGTH WAS DECREASED FOUR TIMES")"")") "ELSE" "IF" OUT[1]=3 "THEN" OUTPUT(61,"("2/, "("A CALL OF DERIV DELIVERED FALSE")"")") "ELSE" "IF" OUT[1]=4 "THEN" OUTPUT(61,"("2/, "("A CALL OF JAC DFDY DELIVERED FALSE ")"")") "ELSE" "IF" OUT[1]=5 "THEN" OUTPUT(61,"("2/, "("A CALL OF JAC DFDP DELIVERED FALSE ")"")") "ELSE" "IF" OUT[1]=6 "THEN" OUTPUT(61,"("2/, "("PRECISION ASKED FOR MAY NOT BE ATTAINED")"")"); "IF" NBP=0 "THEN" OUTPUT(61,"("2/, "("LAST INTEGRATION WAS PERFORMED WITHOUT BREAK-POINTS")"")") "ELSE" "BEGIN" OUTPUT(61,"("2/, "("THE PROCESS STOPPED WITH BREAK-POINTS: ")"")"); "FOR" I:=1 "STEP" 1 "UNTIL" NBP "DO" OUTPUT(61,"("ZZD,B")",BP[I]) "END"; OUTPUT(61,"("4/, "("EUCL. NORM OF THE LAST RESIDUAL VECTOR :")",.7D"+ZD,/, "("EUCL. NORM OF THE FIRST RESIDUAL VECTOR:")",.7D"+ZD,/, "("NUMBER OF INTEGRATIONS PERFORMED")",7B,"(":")",ZZD,/, "("LAST IMPROVEMENT OF THE EUCLIDEAN NORM :")",.7D"+ZD,/, "("CONDITON NUMBER OF J'*J")",15B,"(":")",.7D"+ZD,/, "("LOCAL ERROR BOUND WAS EXCEEDED (MAXIM.):")",ZZD,7/")", OUT[2],OUT[3],OUT[4],OUT[6],OUT[7],OUT[5]); "COMMENT" 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 10 ; "COMMENT" STATISTICS FOR THE PARAMETERS; OUTPUT(61,"("//,B,"("PARAMETERS")",12B,"("CONFIDENCE INTERVAL")", /")"); "FOR" I:=1 "STEP" 1 "UNTIL" M "DO" "BEGIN" CONF[I]:=SQRT(M*FA*JTJINV[I,I]/(NOBS-M))*OUT[2]; OUTPUT(61,"("+.7D"+ZD,12B,+.7D"+ZD,/")",PAR[I],CONF[I]); "END"; C:="IF" NOBS=M "THEN" 0 "ELSE" OUT[2]*OUT[2]/(NOBS-M); OUTPUT(61,"("5/,"("CORRELATION MATRIX")",11B,"("COVARIANCE MATRIX")", /")"); "FOR" I:=1 "STEP" 1 "UNTIL" M "DO" "BEGIN" "FOR" J:=1 "STEP" 1 "UNTIL" M "DO" "BEGIN" "IF" I=J "THEN" OUTPUT(61,"("29B")"); "IF" I>J "THEN" OUTPUT(61,"("+.7D"+ZD,B")", JTJINV[I,J]/SQRT(JTJINV[I,I]*JTJINV[J,J])) "ELSE" OUTPUT(61,"("+.7D"+ZD,B")",JTJINV[I,J]*C) "END"; OUTPUT(61,"("/")"); "END"; OUTPUT(61,"("*")"); OUTPUT(61,"("3/,10B,"("THE LAST RESIDUAL VECTOR")",//,15B, "("I")",4B,"("RES[I]")",/")"); "FOR" I:=1 "STEP" 1 "UNTIL" NOBS "DO" OUTPUT(61,"("14B,ZD,2B,+.4D"+ZD,/")",I,RES[I]) "END" "END" COMMUNICATION; "EOP" THE USER PROGRAM READS: "BEGIN" "INTEGER" I,M,N,NOBS,NBP; "REAL" TIME,FA; "ARRAY" PAR[1:6],RES[1:26],JTJINV[1:3,1:3],IN[0:6],OUT[1:7]; "INTEGER" "ARRAY" BP[0:3]; "PROCEDURE" PEIDE(N,M,NO,NB,P,R,BP,J,I,O,D,JDY,JDP,CY,DA,MO); "CODE" 34444; "PROCEDURE" COMMUNICATION(P,F,M,N,NO,NP,PA,R,BP,J,I,O,W,NI); "CODE" 34445; "BOOLEAN" "PROCEDURE" JAC DFDP(PAR,Y,X,FP); "REAL" X; "ARRAY" PAR,Y,FP; "BEGIN" "REAL" Y2; Y2:=Y[2]; FP[1,1]:=FP[1,3]:=0; FP[1,2]:=Y2*EXP(PAR[2]); FP[2,1]:=EXP(PAR[1])*(Y[1]*(1-Y2)-(EXP(PAR[2])+EXP(PAR[3]))*Y2); FP[2,2]:=-EXP(PAR[1]+PAR[2])*Y2; FP[2,3]:=-EXP(PAR[1]+PAR[3])*Y2; JAC DFDP:="TRUE" "END" JAC DFDP 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 11 ; "PROCEDURE" DATA(NOBS,TOBS,OBS,COBS); "VALUE" NOBS; "INTEGER" NOBS; "ARRAY" TOBS,OBS,COBS; "BEGIN" "INTEGER" I; TOBS[0]:=0; OUTPUT(61,"("*,4/,4B,"("THE OBSERVATIONS WERE:")", //,B,"("I")",3B,"("TOBS[I]")",3B,"("COBS[I]")",3B, "("OBS[I]")",/")"); "FOR" I:=1 "STEP" 1 "UNTIL" NOBS "DO" "BEGIN" INPUT(60,"("3(N)")",TOBS[I],COBS[I],OBS[I]); OUTPUT(61,"("ZD,3B,ZD.4D,6B,D,6B,.4D,/")",I,TOBS[I],COBS[I], OBS[I]) "END" "END" DATA; "PROCEDURE" CALL YSTART(PAR,Y,YMAX); "ARRAY" PAR,Y,YMAX; "BEGIN" Y[1]:=YMAX[1]:=YMAX[2]:=1; Y[2]:=0 "END" CALL YSTART; "BOOLEAN" "PROCEDURE" DERIV(PAR,Y,X,DF); "REAL" X; "ARRAY" PAR,Y,DF; "BEGIN" "REAL" Y2; Y2:=Y[2]; DF[1]:=-(1-Y2)*Y[1]+EXP(PAR[2])*Y2; DF[2]:=EXP(PAR[1])*((1-Y2)*Y[1]-(EXP(PAR[2])+EXP(PAR[3]))*Y2); DERIV:="TRUE" "END" DERIV; "BOOLEAN" "PROCEDURE" JAC DFDY(PAR,Y,X,FY); "REAL" X; "ARRAY" PAR,Y,FY; "BEGIN" FY[1,1]:=-1+Y[2]; FY[1,2]:=EXP(PAR[2])+Y[1]; FY[2,1]:=EXP(PAR[1])*(1-Y[2]); FY[2,2]:=-EXP(PAR[1])*(EXP(PAR[2])+EXP(PAR[3])+Y[1]); JAC DFDY:="TRUE" "END" JAC DFDY; "PROCEDURE" MONITOR(POST,NCOL,NROW,PAR,RES,WEIGHT,NIS); "VALUE" POST,NCOL,NROW,WEIGHT,NIS; "INTEGER" POST,NCOL,NROW,WEIGHT,NIS; "ARRAY" PAR,RES;; OUTPUT(61,"("2/,30B,"("E S C E P - PROBLEM")",3/")"); M:= 3; N:=2; NOBS:=23; NBP:=3; PAR[1]:=LN(1600); PAR[2]:=LN(.8); PAR[3]:=LN(1.2); IN[0]:="-14; IN[3]:="-4; IN[4]:="-4; IN[5]:=50; IN[6]:="-2; IN[1]:="-4; IN[2]:="-5; BP[1]:=17; BP[2]:=19; BP[3]:=21; "COMMENT" 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 12 ; FA:=4.94; "COMMENT" FA DENOTES THE ALPHA-POINT OF THE FISHER-DISTRIBUTION; COMMUNICATION(1,FA,N,M,NOBS,NBP,PAR,RES,BP,JTJINV,IN,OUT,0,0); TIME:=CLOCK; PEIDE(N,M,NOBS,NBP,PAR,RES,BP,JTJINV,IN,OUT,DERIV,JAC DFDY,JAC DFDP, CALL YSTART,DATA,MONITOR); TIME:=CLOCK-TIME; COMMUNICATION(2,FA,N,M,NOBS,NBP,PAR,RES,BP,JTJINV,IN,OUT,0,0); OUTPUT(61,"("3/,5B, "("THE CALCULATION IN PEIDE CONSUMED")",B,ZZD.DD,2B, "("SECONDS")",*")",TIME) "END" "EOP" THIS PROGRAM DELIVERS: E S C E P - PROBLEM STARTING VALUES OF THE PARAMETERS +.7377759" +1 -.2231436" +0 +.1823216" +0 NUMBER OF EQUATIONS : 2 NUMBER OF OBSERVATIONS:23 MACHINE PRECISION :+.1"-13 RELATIVE LOCAL ERROR BOUND FOR INTEGRATION :+.1" -4 RELATIVE TOLERANCE FOR RESIDUE :+.10" -3 ABSOLUTE TOLERANCE FOR RESIDUE :+.10" -3 MAXIMUM NUMBER OF INTEGRATIONS TO PERFORM : 50 RELATIVE STARTING VALUE OF LAMBDA :+.10" -1 RELATIVE MINIMAL STEPLENGTH :+.10" -3 BREAK-POINTS ARE THE OBSERVATIONS : 17 19 21 THE ALPHA-POINT OF THE F-DISTIBUTION : 4.94 THE OBSERVATIONS WERE: 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 13 I TOBS[I] COBS[I] OBS[I] 1 0.0002 2 .1648 2 0.0004 2 .2753 3 0.0006 2 .3493 4 0.0008 2 .3990 5 0.0010 2 .4322 6 0.0012 2 .4545 7 0.0014 2 .4695 8 0.0016 2 .4795 9 0.0018 2 .4862 10 0.0020 2 .4907 11 0.0200 2 .4999 12 0.0400 2 .4998 13 0.0600 2 .4998 14 0.0800 2 .4998 15 0.1000 2 .4998 16 1.0000 2 .4986 17 2.0000 2 .4973 18 5.0000 2 .4936 19 10.0000 2 .4872 20 15.0000 2 .4808 21 20.0000 2 .4743 22 25.0000 2 .4677 23 30.0000 2 .4610 NORMAL TERMINATION OF THE PROCESS LAST INTEGRATION WAS PERFORMED WITHOUT BREAK-POINTS EUCL. NORM OF THE LAST RESIDUAL VECTOR :.1430776" -3 EUCL. NORM OF THE FIRST RESIDUAL VECTOR:.1331071" +1 NUMBER OF INTEGRATIONS PERFORMED : 12 LAST IMPROVEMENT OF THE EUCLIDEAN NORM :.2223694" -4 CONDITON NUMBER OF J'*J :.2582882" +3 LOCAL ERROR BOUND WAS EXCEEDED (MAXIM.): 37 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 14 PARAMETERS CONFIDENCE INTERVAL +.6907670" +1 +.3209313" -3 -.1003941" -1 +.1687774" -3 -.4605292" +1 +.1942501" -2 CORRELATION MATRIX COVARIANCE MATRIX +.6949857" -8 +.1407628" -8 -.9129848" -8 +.3851320" +0 +.1922119" -8 -.1414245" -7 -.2170393" +0 -.6392889" +0 +.2546094" -6 THE LAST RESIDUAL VECTOR I RES[I] 1 +.1748" -5 2 -.2905" -4 3 +.2814" -4 4 -.3879" -4 5 +.3069" -4 6 +.3101" -4 7 -.2019" -4 8 -.3887" -5 9 +.1052" -4 10 +.1391" -4 11 -.5109" -4 12 +.2384" -4 13 -.1156" -5 14 -.2616" -4 15 -.5116" -4 16 +.2244" -4 17 +.6794" -4 18 -.1418" -4 19 +.2087" -4 20 -.1980" -4 21 -.3476" -4 22 -.2245" -4 23 +.1886" -4 THE CALCULATION IN PEIDE CONSUMED 108.57 SECONDS 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 15 SOURCE TEXT(S): 0"CODE" 34444; "PROCEDURE" PEIDE(N,M,NOBS,NBP,PAR,RES,BP,JTJINV,IN,OUT,DERIV,JAC DFDY, JAC DFDP, CALL YSTART,DATA,MONITOR); "VALUE" N,M,NOBS; "INTEGER" N,M,NOBS,NBP; "ARRAY" PAR,RES,JTJINV,IN,OUT; "INTEGER" "ARRAY" BP; "PROCEDURE" CALL YSTART,DATA,MONITOR; "BOOLEAN" "PROCEDURE" DERIV,JAC DFDY,JACDFDP; "BEGIN" "INTEGER" I,J,EXTRA,WEIGHT,NCOL,NROW,AWAY,NPAR,II,JJ,MAX, NFE,NIS; "REAL" EPS,EPS1,XEND,C,X,T,HMIN,HMAX,RES1,IN3,IN4,FAC3,FAC4; "ARRAY" AUX[1:3],OBS[1:NOBS],SAVE[-38:6*N],TOBS[0:NOBS], YP[1:NBP+NOBS,1:NBP+M],YMAX[1:N],Y[1:6*N*(NBP+M+1)],FY[1:N,1:N], FP[1:N,1:M+NBP]; "INTEGER" "ARRAY" COBS[1:NOBS]; "BOOLEAN" FIRST,SEC,CLEAN; "PROCEDURE" INIVEC(L,U,A,X); "CODE" 31010; "PROCEDURE" INIMAT(L1,U1,L2,U2,A,X); "CODE" 31011; "PROCEDURE" MULVEC(L,U,S,A,B,X); "CODE" 31020; "PROCEDURE" MULROW(L,U,I,J,A,B,X); "CODE" 31021; "PROCEDURE" DUPVEC(L,U,S,A,B); "CODE" 31030; "PROCEDURE" DUPMAT(L1,U1,L2,U2,A,B); "CODE" 31035; "REAL" "PROCEDURE" VECVEC(L,U,S,A,B); "CODE" 34010; "REAL" "PROCEDURE" MATVEC(L,U,I,A,B); "CODE"34011; "PROCEDURE" ELMVEC(L,U,S,A,B,X); "CODE"34020; "PROCEDURE" SOL(A,N,P,B); "CODE" 34051; "PROCEDURE" DEC(A,N,AUX,P); "CODE" 34300; "PROCEDURE" MARQUARDT(M,N,P,R,C,F,J,I,O); "CODE" 34440; "REAL" "PROCEDURE" INTERPOL(STARTINDEX,JUMP,K,TOBSDIF); "VALUE" STARTINDEX,JUMP,K,TOBSDIF; "INTEGER" STARTINDEX,JUMP,K; "REAL" TOBSDIF; "BEGIN" "INTEGER" I; "REAL" S,R; S:=Y[STARTINDEX]; R:=TOBSDIF; "FOR" I:=1 "STEP" 1 "UNTIL" K "DO" "BEGIN" STARTINDEX:=STARTINDEX+JUMP; S:=S+Y[STARTINDEX]*R; R:=R*TOBSDIF "END"; INTERPOL:=S "END" INTERPOL; "PROCEDURE" JAC DYDP(NROW,NCOL,PAR,RES,JAC,LOCFUNCT); "VALUE" NROW,NCOL; "INTEGER" NROW,NCOL; "ARRAY" PAR,RES,JAC; "PROCEDURE" LOCFUNCT; "BEGIN" DUPMAT(1,NROW,1,NCOL,JAC,YP) "END" JACOBIAN; 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 16 ; "BOOLEAN" "PROCEDURE" FUNCT(NROW,NCOL,PAR,RES); "VALUE" NROW,NCOL; "INTEGER" NROW,NCOL; "ARRAY" PAR,RES; "BEGIN" "INTEGER" L,K,KNEW,FAILS,SAME,KPOLD,N6,NNPAR,J5N, COBSII; "REAL" XOLD,HOLD,A0,TOLUP,TOL,TOLDWN,TOLCONV,H,CH,CHNEW, ERROR,DFI,TOBSDIF; "BOOLEAN" EVALUATE,EVALUATED,DECOMPOSE,CONV; "ARRAY" A[0:5],DELTA,LAST DELTA,DF,Y0[1:N],JACOB[1:N,1:N]; "INTEGER" "ARRAY" P[1:N]; "REAL" "PROCEDURE" NORM2(AI); "REAL" AI; "BEGIN" "REAL" S,A; S:= "-100; "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" A:= AI/YMAX[I]; S:= S + A * A "END"; NORM2:= S "END" NORM2; "PROCEDURE" RESET; "BEGIN" "IF" CH < HMIN/HOLD "THEN" CH:= HMIN/HOLD "ELSE" "IF" CH > HMAX/HOLD "THEN" CH:= HMAX/HOLD; X:= XOLD; H:= HOLD * CH; C:= 1; "FOR" J:= 0 "STEP" N "UNTIL" K*N "DO" "BEGIN" "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" Y[J+I]:= SAVE[J+I] * C; C:= C * CH "END"; DECOMPOSE:="TRUE" "END" RESET; "PROCEDURE" ORDER; "BEGIN" C:= EPS * EPS; J:= (K-1) * (K + 8)/2 - 38; "FOR" I:= 0 "STEP" 1 "UNTIL" K "DO" A[I]:= SAVE[I+J]; J:= J + K + 1; TOLUP := C * SAVE[J]; TOL := C * SAVE[J + 1]; TOLDWN := C * SAVE[J + 2]; TOLCONV:= EPS/(2 * N * (K + 2)); A0:= A[0]; DECOMPOSE:= "TRUE"; "END" ORDER; "PROCEDURE" EVALUATE JACOBIAN; "BEGIN" EVALUATE:= "FALSE"; DECOMPOSE:= EVALUATED:= "TRUE"; "IF" "NOT" JAC DFDY(PAR,Y,X,FY) "THEN" "BEGIN" SAVE[-3]:=4; "GOTO" RETURN "END"; "END" EVALUATE JACOBIAN; 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 17 ; "PROCEDURE" DECOMPOSE JACOBIAN; "BEGIN" DECOMPOSE:= "FALSE"; C:= -A0 * H; "FOR" J:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" JACOB[I,J]:= FY[I,J] * C; JACOB[J,J]:= JACOB[J,J] + 1 "END"; DEC(JACOB,N,AUX,P) "END" DECOMPOSE JACOBIAN; "PROCEDURE" CALCULATE STEP AND ORDER; "BEGIN" "REAL" A1,A2,A3; A1:= "IF" K <= 1 "THEN" 0 "ELSE" 0.75 * (TOLDWN/NORM2(Y[K*N+I])) ** (0.5/K); A2:= 0.80 * (TOL/ERROR) ** (0.5/(K + 1)); A3:= "IF" K >= 5 "OR" FAILS ^= 0 "THEN" 0 "ELSE" 0.70 * (TOLUP/NORM2(DELTA[I] - LAST DELTA[I]))** (0.5/(K+2)); "IF" A1 > A2 "AND" A1 > A3 "THEN" "BEGIN" KNEW:= K-1; CHNEW:= A1 "END" "ELSE" "IF" A2 > A3 "THEN" "BEGIN" KNEW:= K ; CHNEW:= A2 "END" "ELSE" "BEGIN" KNEW:= K+1; CHNEW:= A3 "END" "END" CALCULATE STEP AND ORDER; "IF" SEC "THEN" "BEGIN" SEC:="FALSE"; "GOTO" RETURN "END"; NPAR:=M; EXTRA:=NIS:=0; II:=1; JJ:="IF" NBP=0 "THEN" 0 "ELSE" 1; N6:=N*6; INIVEC(-3,-1,SAVE,0); INIVEC(N6+1,(6+M)*N,Y,0); INIMAT(1,NOBS+NBP,1,M+NBP,YP,0); T:=TOBS[1]; X:=TOBS[0]; CALL YSTART(PAR,Y,YMAX); HMAX:=TOBS[1]-TOBS[0]; HMIN:=HMAX*IN[1]; EVALUATE JACOBIAN; NNPAR:=N*NPAR; NEW START: K:= 1; KPOLD:=0; SAME:= 2; ORDER; "IF" "NOT" DERIV(PAR,Y,X,DF) "THEN" "BEGIN" SAVE[-3]:=3; "GOTO" RETURN "END"; H:=SQRT(2 * EPS/SQRT(NORM2 (MATVEC(1,N,I,FY,DF)))); "IF" H > HMAX "THEN" H:= HMAX "ELSE" "IF" H < HMIN "THEN" H:= HMIN; XOLD:= X; HOLD:= H; CH:= 1; "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" SAVE[I]:=Y[I]; SAVE[N+I]:=Y[N+I]:=DF[I]*H "END"; FAILS:= 0; "COMMENT" 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 18 ; "FOR" L:= 0 "WHILE" X < XEND "DO" "BEGIN" "IF" X + H <= XEND "THEN" X:= X + H "ELSE" "BEGIN" H:= XEND-X; X:= XEND; CH:= H/HOLD; C:= 1; "FOR" J:= N "STEP" N "UNTIL" K*N "DO" "BEGIN" C:= C* CH; "FOR" I:= J+1 "STEP" 1 "UNTIL" J+N "DO" Y[I]:= Y[I] * C "END"; SAME:= "IF" SAME<3 "THEN" 3 "ELSE" SAME+1; "END"; "COMMENT" PREDICTION; "FOR" L:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" "FOR" I:= L "STEP" N "UNTIL" (K-1)*N+L "DO" "FOR" J:= (K-1)*N+L "STEP" -N "UNTIL" I "DO" Y[J]:= Y[J] + Y[J+N]; DELTA[L]:= 0 "END"; EVALUATED:= "FALSE"; "COMMENT" CORRECTION AND ESTIMATION LOCAL ERROR; "FOR" L:= 1,2,3 "DO" "BEGIN" "IF" "NOT" DERIV(PAR,Y,X,DF) "THEN" "BEGIN" SAVE[-3]:=3; "GOTO" RETURN "END"; "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" DF[I]:= DF[I] * H - Y[N+I]; "IF" EVALUATE "THEN" EVALUATE JACOBIAN; "IF" DECOMPOSE "THEN" DECOMPOSE JACOBIAN; SOL(JACOB,N,P,DF); CONV:= "TRUE"; "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" DFI:= DF[I]; Y[ I]:= Y[ I] + A0 * DFI; Y[N+I]:= Y[N+I] + DFI; DELTA[I]:= DELTA[I] + DFI; CONV:= CONV "AND" ABS(DFI) < TOLCONV * YMAX[I] "END"; "IF" CONV "THEN" "BEGIN" ERROR:= NORM2(DELTA[I]); "GOTO" CONVERGENCE "END" "END"; "COMMENT" ACCEPTANCE OR REJECTION; "IF" "NOT" CONV "THEN" "BEGIN" "IF" "NOT" EVALUATED "THEN" EVALUATE:= "TRUE" "ELSE" "BEGIN" CH:=CH/4; "IF" H<4*HMIN "THEN" "BEGIN" SAVE[-1]:= SAVE[-1]+10; HMIN:=HMIN/10; "COMMENT" 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 19 ; "IF" SAVE[-1]>40 "THEN" "GOTO" RETURN "END" "END"; RESET "END" "ELSE" CONVERGENCE: "IF" ERROR > TOL "THEN" "BEGIN" FAILS:= FAILS + 1; "IF" H > 1.1 * HMIN "THEN" "BEGIN" "IF" FAILS > 2 "THEN" "BEGIN" RESET; "GOTO" NEW START "END" "ELSE" "BEGIN" CALCULATE STEP AND ORDER; "IF" KNEW ^= K "THEN" "BEGIN" K:= KNEW; ORDER "END"; CH:= CH * CHNEW; RESET "END" "END" "ELSE" "BEGIN" "IF" K = 1 "THEN" "BEGIN" "COMMENT" VIOLATE EPS CRITERION; SAVE[-2]:= SAVE[-2] + 1; SAME:= 4; "GOTO" ERROR TEST OK "END"; K:=1; RESET; ORDER; SAME:= 2 "END" "END" "ELSE" ERROR TEST OK: "BEGIN" FAILS:= 0; "FOR" I:= 1 "STEP" 1 "UNTIL" N "DO" "BEGIN" C:= DELTA[I]; "FOR" L:= 2 "STEP" 1 "UNTIL" K "DO" Y[L*N+I]:= Y[L*N+I] + A[L] * C; "IF" ABS(Y[I]) > YMAX[I] "THEN" YMAX[I]:= ABS(Y[I]) "END"; SAME:= SAME-1; "IF" SAME= 1 "THEN" DUPVEC(1,N,0,LAST DELTA,DELTA) "ELSE" "IF" SAME= 0 "THEN" "BEGIN" CALCULATE STEP AND ORDER; "IF" CHNEW > 1.1 "THEN" "BEGIN" "IF" K ^= KNEW "THEN" "BEGIN" "IF" KNEW > K "THEN" MULVEC(KNEW*N+1,KNEW*N+N,-KNEW*N,Y,DELTA, A[K]/KNEW); K:= KNEW; ORDER "END"; SAME:= K+1; "IF" CHNEW * H > HMAX "THEN" CHNEW:= HMAX/H; "COMMENT" 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 20 ; H:= H * CHNEW; C:= 1; "FOR" J:= N "STEP" N "UNTIL" K*N "DO" "BEGIN" C:= C * CHNEW; MULVEC(J+1,J+N,0,Y,Y,C) "END"; DECOMPOSE:="TRUE" "END" "ELSE" SAME:= 10 "END" OF A SINGLE INTEGRATION STEP OF Y; NIS:=NIS+1; "COMMENT" START OF A INTEGRATION STEP OF YP; "IF" CLEAN "THEN" "BEGIN" HOLD:=H; XOLD:=X; KPOLD:=K; CH:=1; DUPVEC(1,K*N+N,0,SAVE,Y) "END" "ELSE" "BEGIN" "IF" H^=HOLD "THEN" "BEGIN" CH:=H/HOLD; C:=1; "FOR" J:=N6+NNPAR "STEP" NNPAR "UNTIL" KPOLD*NNPAR+N6 "DO" "BEGIN" C:=C*CH; "FOR" I:=J+1 "STEP" 1 "UNTIL" J+NNPAR "DO" Y[I]:=Y[I]*C "END"; HOLD:=H "END"; "IF" K>KPOLD "THEN" INIVEC(N6+K*NNPAR+1,N6+K*NNPAR+NNPAR,Y,0); XOLD:= X; KPOLD:= K; CH:= 1; DUPVEC(1,K*N+N,0,SAVE,Y); EVALUATE JACOBIAN; DECOMPOSE JACOBIAN; "IF" "NOT" JAC DFDP(PAR,Y,X,FP) "THEN" "BEGIN" SAVE[-3]:=5; "GOTO" RETURN "END"; "IF" NPAR>M "THEN" INIMAT(1,N,M+1,NPAR,FP,0); "COMMENT" PREDICTION; "FOR" L:=0 "STEP" 1 "UNTIL" K-1 "DO" "FOR" J:=K-1 "STEP" -1 "UNTIL" L "DO" ELMVEC(J*NNPAR+N6+1,J*NNPAR+N6+NNPAR,NNPAR,Y,Y,1); "COMMENT" CORRECTION; "FOR" J:=1 "STEP" 1 "UNTIL" NPAR "DO" "BEGIN" J5N:=(J+5)*N; DUPVEC(1,N,J5N,Y0,Y); "FOR" I:=1 "STEP" 1 "UNTIL" N "DO" DF[I]:= H*(FP[I,J]+MATVEC(1,N,I,FY,Y0)) -Y[NNPAR+J5N+I]; SOL(JACOB,N,P,DF); "FOR" L:=0 "STEP" 1 "UNTIL" K "DO" "BEGIN" I:=L*NNPAR+J5N; ELMVEC(I+1,I+N,-I,Y,DF,A[L]) "END" "END" "END"; "COMMENT" 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 21 ; "FOR" L:=0 "WHILE" X>=T "DO" "BEGIN" "COMMENT" CALCULATION OF A ROW OF THE JACOBIAN MATRIX AND AN ELEMENT OF THE RESIDUAL VECTOR; TOBSDIF:=(TOBS[II]-X)/H; COBSII:=COBS[II]; RES[II]:=INTERPOL(COBSII,N,K,TOBSDIF)-OBS[II]; "IF" "NOT" CLEAN "THEN" "BEGIN" "FOR" I:=1 "STEP" 1 "UNTIL" NPAR "DO" YP[II,I]:=INTERPOL(COBSII+(I+5)*N,NNPAR,K, TOBSDIF); "COMMENT" INTRODUCING OF BREAK-POINTS; "IF" BP[JJ]^=II "THEN" "ELSE" "IF" FIRST "AND" ABS(RES[II])0 "THEN" "BEGIN" "FOR" I:=1 "STEP" 1 "UNTIL" N "DO" "BEGIN" Y[I]:=INTERPOL(I,N,K,TOBSDIF); "FOR" J:=1 "STEP" 1 "UNTIL" NPAR "DO" Y[I+(J+5)*N]:=INTERPOL(I+(J+5)*N,NNPAR,K, TOBSDIF) "END"; "FOR" L:=1 "STEP" 1 "UNTIL" EXTRA "DO" "BEGIN" COBSII:=COBS[BP[NPAR-M+L]]; Y[COBSII]:=PAR[NPAR+L]; "FOR" I:=1 "STEP" 1 "UNTIL" NPAR+EXTRA "DO" Y[COBSII+(5+I)*N]:=0; "COMMENT" 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 22 ; INIVEC(1+NNPAR+(L+5)*N,NNPAR+(L+6)*N,Y,0); Y[COBSII+(5+NPAR+L)*N]:=1 "END"; NPAR:=NPAR+EXTRA; EXTRA:=0; X:=TOBS[II-1]; EVALUATE JACOBIAN; NNPAR:=N*NPAR; "GOTO" NEW START "END" "END" "END" STEP; RETURN: "IF" SAVE[-2]>MAX "THEN" MAX:=SAVE[-2]; FUNCT:=SAVE[-1]<=40 "AND" SAVE[-3]=0; "IF" "NOT" FIRST "THEN" MONITOR(1,NCOL,NROW,PAR,RES,WEIGHT,NIS) "END" FUNCT; I:= -39; "FOR" C:= 1,1,9,4,0,2/3,1,1/3,36,20.25,1,6/11, 1,6/11,1/11,84.028,53.778,0.25,.48,1,.7,.2,.02, 156.25, 108.51, .027778, 120/274, 1, 225/274, 85/274, 15/274, 1/274, 0, 187.69, .0047361 "DO" "BEGIN" I:= I + 1; SAVE[I]:= C "END"; DATA(NOBS,TOBS,OBS,COBS); WEIGHT:=1; FIRST:=SEC:="FALSE"; CLEAN:=NBP>0; AUX[2]:="-12; EPS:=IN[2]; EPS1:="10; XEND:=TOBS[NOBS]; OUT[1]:=0; BP[0]:=MAX:=0; "COMMENT" SMOOTH INTEGRATION WITHOUT BREAK-POINTS; "IF" "NOT" FUNCT(NOBS,M,PAR,RES) "THEN" "GOTO" ESCAPE; RES1:=SQRT(VECVEC(1,NOBS,0,RES,RES)); NFE:=1; "IF" IN[5]=1 "THEN" "BEGIN" OUT[1]:=1; "GOTO" ESCAPE "END"; "IF" CLEAN "THEN" "BEGIN" FIRST:="TRUE"; CLEAN:="FALSE"; FAC3:=SQRT(SQRT(IN[3]/RES1)); FAC4:=SQRT(SQRT(IN[4]/RES1)); EPS1:=RES1*FAC4; "IF" "NOT" FUNCT(NOBS,M,PAR,RES) "THEN" "GOTO" ESCAPE; FIRST:="FALSE" "END" "ELSE" NFE:=0; NCOL:=M+NBP; NROW:=NOBS+NBP; SEC:="TRUE"; IN3:=IN[3]; IN4:=IN[4]; IN[3]:=RES1; "BEGIN" "REAL" W; "ARRAY" AID[1:NCOL,1:NCOL]; WEIGHT:=AWAY:=0; OUT[4]:=OUT[5]:=W:=0; "COMMENT" 1SECTION : 5.2.1.3.1 (OCTOBER 1975) PAGE 23 ; "FOR" WEIGHT:=(SQRT(WEIGHT)+1)**2 "WHILE" WEIGHT^=16 "AND" NBP>0 "DO" "BEGIN" "IF" AWAY=0 "AND" W^=0 "THEN" "BEGIN" "COMMENT" IF NO BREAK-POINTS WERE OMITTED THEN ONE FUNCTION EVALUATION IS SAVED; W:=WEIGHT/W; "FOR" I:=NOBS+1 "STEP" 1 "UNTIL" NROW "DO" "BEGIN" "FOR" J:=1 "STEP" 1 "UNTIL" NCOL "DO" YP[I,J]:=W*YP[I,J]; RES[I]:=W*RES[I] "END"; SEC:="TRUE"; NFE:=NFE-1 "END"; IN[3]:=IN[3]*FAC3*WEIGHT; IN[4]:=EPS1; MONITOR(2,NCOL,NROW,PAR,RES,WEIGHT,NIS); MARQUARDT(NROW,NCOL,PAR,RES,AID,FUNCT,JAC DYDP,IN,OUT); "IF" OUT[1]>0 "THEN" "GOTO" ESCAPE; "COMMENT" THE RELATIVE STARTING VALUE OF LAMBDA IS ADJUSTED TO THE LAST VALUE OF LAMBDA USED; AWAY:=OUT[4]-OUT[5]-1; IN[6]:=IN[6] * 5**AWAY * 2**(AWAY-OUT[5]); NFE:=NFE+OUT[4]; W:=WEIGHT; EPS1:=(SQRT(WEIGHT)+1)**2*IN[4]*FAC4; AWAY:=0; "COMMENT" USELESS BREAK-POINTS ARE OMITTED; "FOR" J:=1 "STEP" 1 "UNTIL" NBP "DO" "BEGIN" "IF" ABS(OBS[BP[J]]+RES[BP[J]]-PAR[J+M])