NUMAL Section 6.5.1
BEGIN SECTION : 6.5.1 (December, 1979)
AUTHOR(S) : H.FIOLET, N.TEMME.
INSTITUTE : MATHEMATICAL CENTRE.
RECEIVED: 740628.
BRIEF DESCRIPTION :
THIS SECTION CONTAINS FOUR PROCEDURES :
A.
EI CALCULATES THE EXPONENTIAL INTEGRAL DEFINED AS FOLLOWS (SEE
ALSO REF[1] , EQ. (5.1.1)) : EI(X) = INTEGRAL (EXP(T)/T DT) FROM
T=-INFINITY TO T=X , WHERE THE INTEGRAL IS TO BE INTERPRETED AS THE
CAUCHY PRINCIPAL VALUE. ALSO THE RELATED FUNCTION E1(X), DEFINED BY
THE INTEGRAL (EXP(-T)/T DT) FROM T= X TO T= INFINITY, FOR POSITIVE
X (REF[1], EQ.(5.1.2)) CAN EASILY BE OBTAINED BY THE RELATION
E1(X) = - EI(-X). FOR X=0 THE INTEGRAL IS UNDEFINED AND THE
PROCEDURE WILL CAUSE OVERFLOW.
B.
EIALPHA CALCULATES A SEQUENCE OF INTEGRALS OF THE FORM
INTEGRAL( EXP(-X*T)*T**I DT )
FROM T=1 TO T= INFINITY,
WHERE X IS POSITIVE AND I = 0,...,N.
(SEE ALSO REF[1], EQ. (5.1.5)).
C.
ENX COMPUTES A SEQUENCE OF INTEGRALS E(N,X),
N=N1, N1+1,...,N2, WHERE X>0 AND N1, N2 ARE POSITIVE INTEGERS WITH
N2>=N1; E(N,X) IS DEFINED AS FOLLOWS:
E(N,X)= THE INTEGRAL FROM 1 TO INFINITY OF EXP(-X * T)/ T**N DT;
(SEE ALSO REF[1], EQ.(5.1.4));
D.
NONEXPENX COMPUTES A SEQUENCE OF INTEGRALS
EXP(X)*E(N,X), N=N1, N1+1,...,N2, WHERE X>0 AND N1, N2 ARE POSITIVE
INTEGERS WITH N2>=N1; E(N,X) IS DEFINED UNDER C).
KEYWORDS :
EXPONENTIAL INTEGRAL,
SPECIAL FUNCTIONS.
SUBSECTION : EI.
CALLING SEQUENCE :
THE HEADING OF THE PROCEDURE READS:
"REAL" "PROCEDURE" EI(X);
"VALUE" X;"REAL" X;
"CODE" 35080;
EI: DELIVERS THE VALUE OF THE EXPONENTIAL INTEGRAL;
THE MEANING OF THE FORMAL PARAMETER IS :
X: <ARITHMETIC EXPRESSION>;
THE ARGUMENT OF THE INTEGRAL.
PROCEDURES USED :
CHEPOLSUM = CP31046 ,
POL = CP31040 ,
JFRAC = CP35083 .
LANGUAGE : ALGOL 60.
METHOD AND PERFORMANCE :
THE INTEGRAL IS CALCULATED BY MEANS OF THE RATIONAL CHEBYSHEV
APPROXIMATIONS GIVEN IN REFERENCES [1] AND [2]. ONLY RATIOS OF
POLYNOMIALS WITH EQUAL DEGREE L ARE CONSIDERED. BELOW,THE DIFFERENT
INTERVALS ARE LISTED, TOGETHER WITH THE CORRESPONDING DEGREE L AND
THE NUMBER OF CORRECT DIGITS OF THE APPROXIMATIONS :
[-INFINITY,-4] 6 15.1
[-4,-1] 7 16.9
[-1, 0] 5 18.5
[ 0, 6] 7 15.2
[ 6,12] 7 15.1
[12,24] 7 15.0
[24,+INFINITY] 7 15.9 .
VARIOUS TESTS SHOWED A RELATIVE ACCURACY OF AT LEAST "-13, EXEPT
IN THE NEIGHBOURHOOD OF X=.37250 , THE ZERO OF THE INTEGRAL, WHERE
ONLY AN ABSOLUTE ACCURACY OF .3"-13 IS REACHED . IN SOME OF THE
INTERVALS , THE RATIONAL FUNCTIONS ARE EXPRESSED EITHER AS RATIOS
OF FINITE SUMS OF CHEBYSHEV POLYNOMIALS OR AS J-FRACTIONS, SINCE
THE ORIGINAL FORMS ARE POORLY CONDITIONED.
REFERENCES : SEE REFERENCES [1], [2] AND [3] OF THE PROCEDURE
NONEXPENX (THIS SECTION).
EXAMPLE OF USE :
"BEGIN"
"COMMENT" THE COMPUTATION OF E1(.5);
OUTPUT(61,"("N")",-EI(-.5))
"END"
DELIVERS :
+5.5977359477616"-001 .
SUBSECTION : EIALPHA.
CALLING SEQUENCE :
THE HEADING OF THE PROCEDURE READS :
"PROCEDURE" EIALPHA(X,N,ALPHA);
"VALUE" N,X;"INTEGER" N;"REAL" X;"ARRAY" ALPHA;
"CODE" 35081;
THE MEANING OF THE FORMAL PARAMETERS IS :
X: <ARITHMETIC EXPRESSION>;
THE REAL X OCCURING IN THE INTEGRAND.
N: <ARITHMETIC EXPRESSION>;
UPPER BOUND FOR THE INTEGER I OCCURING IN THE INTEGRAND;
ALPHA: <ARRAY IDENTIFIER>;
"ARRAY" ALPHA[0:N];
THE VALUE OF THE INTEGRAL(EXP(-X*T)*T**I DT) FROM T=1 TO
T=INFINITY IS STORED IN ALPHA[I].
PROCEDURES USED : NONE.
RUNNING TIME : CIRCA ( 6 + N * .8 ) * "-4 SEC.
LANGUAGE : ALGOL 60.
METHOD AND PERFORMANCE :
THE INTEGRAL IS CALCULATED BY MEANS OF THE RECURSION FORMULA
A[N]:=A[0] + N * A[N-1] / X, WITH A[0]:= EXP(-X)/X. FOR X CLOSE TO
ZERO, EIALPHA MIGHT CAUSE OVERFLOW, SINCE THE VALUE OF THE INTEGRAL
IS INFINITE FOR X=0. THE PROCEDURE IS NOT PROTECTED AGAINST THIS
TYPE OF OVERFLOW. THE MINIMAL VALUE FOR THE ARGUMENT X DEPENDS ON
THE PARAMETER N :
N=20 X CIRCA "-14
N=15 X CIRCA "-18
N=10 X CIRCA "-28
N= 5 X CIRCA "-53
THE RECURSION FORMULA IS STABLE AND VARIOUS TESTS EXECUTED ON THE
CD CYBER 7228 SHOWED A RELATIVE ACCURACY OF AT LEAST .2"-12.
EXAMPLE OF USE :
"BEGIN"
"INTEGER" K;"REAL" "ARRAY" A[0:5];
EIALPHA(.25,5,A);
"FOR" K:=0 "STEP" 1 "UNTIL" 5 "DO"
OUTPUT(61,"("DBBB,N,/")",K,A[K]);
"END"
DELIVERS :
0 +3.1152031322856"+000
1 +1.5576015661428"+001
2 +1.2772332842371"+002
3 +1.5357951442168"+003
4 +2.4575837510601"+004
5 +4.9151986541516"+005 .
REFERENCES: SEE REFERENCE [1] OF THE PROCEDURE NONEXPENX(THIS SECTION).
SUBSECTION: ENX.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS :
"PROCEDURE" ENX(X, N1, N2, A);
"VALUE" X, N1, N2; "REAL" X; "INTEGER" N1, N2; "ARRAY" A;
"CODE" 35086;
THE MEANING OF THE FORMAL PARAMETERS IS :
X : <ARITHMETIC EXPRESSION>;
ENTRY: THE (REAL) POSITIVE X OCCURING IN THE INTEGRAND;
N1, N2: <ARITHMETIC EXPRESSION>;
ENTRY: LOWER AND UPPER BOUND, RESPECTIVELY, OF THE INTEGER
N OCCURING IN THE INTEGRAND;
A: <ARRAY IDENTIFIER>;
"ARRAY" A[N1:N2];
EXIT: THE VALUE OF THE INTEGRAL(EXP(-X * T)/T**I DT) FROM
T=1 TO T= INFINITY IS STORED IN A[I].
PROCEDURES USED:
EI = CP35080,
NONEXPENX = CP35087.
RUNNING TIME:
DEPENDS STRONGLY ON THE VALUES OF X, N1, AND N2, WITH A MAXIMUM
OF ROUGHLY ( 5 + .1 * NUMBER OF NECESSARY ITERATIONS ) MSEC.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
SEE METHOD AND PERFORMANCE OF THE PROCEDURE NONEXPENX(THIS SECTION)
SUBSECTION: NONEXPENX.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS :
"PROCEDURE" NONEXPENX(X, N1, N2, A);
"VALUE" X, N1, N2; "REAL" X; "INTEGER" N1, N2; "ARRAY" A;
"CODE" 35087;
THE MEANING OF THE FORMAL PARAMETERS IS :
X: <ARITHMETIC EXPRESSION>;
ENTRY: THE (REAL) POSITIVE X OCCURING IN THE INTEGRAND;
N1, N2: <ARITHMETIC EXPRESSION>;
ENTRY: LOWER AND UPPER BOUND, RESPECTIVELY, OF THE INTEGER
I OCCURING IN THE INTEGRAND;
A: <ARRAY IDENTIFIER>;
"ARRAY" A[N1:N2];
EXIT: THE VALUE OF EXP(X) * INTEGRAL(EXP(-X*T)/T**I DT)
FROM T=1 TO T=INFINITY IS STORED IN A[I].
PROCEDURES USED:
ENX = CP35086.
RUNNING TIME:
DEPENDS STRONGLY ON THE VALUES OF X, N1, AND N2, WITH A MAXIMUM
OF ROUGHLY ( 5 + .1 * NUMBER OF NECESSARY ITERATIONS) MSEC.
LANGUAGE: ALGOL 60.
METHOD AND PERFORMANCE:
THE SEQUENCE OF INTEGRALS IS GENERATED BY MEANS OF THE RECURRENCE
RELATION:
E(N+1,X) = (EXP(-X) - X * E(N,X))/N.
FOR REASONS OF STABILITY THE RECURSION STARTS WITH E(N0,X), WHERE
N0=ENTIER(X+.5), (SEE ALSO REF[5]). THE INTEGRALS ARE THEN COMPUTED
BY BACKWARD RECURRENCE IF N<N0 AND BY FORWARD RECURRENCE IF N>N0.
TO OBTAIN THE STARTING VALUES E(N0,X) OF THE RECURSION THE
FOLLOWING CASES ARE DISTINGUISHED:
A) N0 = 1: THE PROCEDURE EI IS USED (THIS SECTION);
B) N0<=10: A TAYLOR EXPANSION ABOUT X=N0 IS USED, WHICH MADE IT
NECESSARY TO STORE THE VALUES OF E(N,N) IN THE PROCEDURE
FOR N= 2, 3,...,10;
C) N0 >10: THE FOLLOWING CONTINUED FRACTION IS USED:
EXP(X)*E(N,X) = 1/(X+N/(1+1/(X+(N+1)/(1+...)))),
(SEE ALSO REF[4], EQ.(2.3));
THE CASES A) AND B) ARE TREATED IN ENX, WHILE NONEXPENX EVALUATES
THE CONTINUED FRACTION IN CASE C).
ENX CALLS FOR NONEXPENX IN CASE C).
NONEXPENX CALLS FOR ENX IN THE CASES A) AND B).
VARIOUS TESTS SHOWED A RELATIVE ACCURACY OF AT LEAST 5"-14.
REFERENCES:
[1].M.ABRAMOWITZ AND I.A.STEGUN.
HANDBOOK OF MATHEMATICAL FUNCTIONS.
DOVER PUBLICATIONS, INC. NEW YORK (1965).
[2] W.J.CODY AND H.C.THACHER, JR.
RATIONAL CHEBYSHEV APPROXIMATIONS FOR THE EXPONENTIAL INTEGRAL
E1(X).
MATH. COMP. 22 (JULY 1968), 641-649.
[3] W.J.CODY AND H.C.THACHER, JR.
CHEBYSHEV APPROXIMATIONS FOR THE EXPONENTIAL INTEGRAL EI(X).
MATH. COMP. 23 (APRIL 1969), 289-303.
[4].W.GAUTSCHI.
EXPONENTIAL INTEGRALS.
CACM, DECEMBER 1973, P.761-763.
[5].W.GAUTSCHI.
RECURSIVE COMPUTATION OF CERTAIN INTEGRALS.
JACM, VOL.8, 1961, P.21-40.
EXAMPLE OF USE:
IN THE FOLLOWING PROGRAM WE COMPUTE THE VALUES OF
E(40,1.1), E(41,1.1), E(42,1.1) AND EXP(X)*E(1,50.1).
"BEGIN"
"INTEGER" I;
"REAL" "ARRAY" A[40:42], B[1:1];
ENX(1.1, 40, 42, A);
"FOR" I:= 40, 41, 42 "DO"
OUTPUT(61,"("4B,"("E(")",DD,"(",1.1)= ")",N/")",I,A[I]);
NONEXPENX(50.1, 1, 1, B);
OUTPUT(61,"("/,4B,"("EXP(50.1)*E(1,50.1)= ")",N")",B[1]);
"END"
THIS PROGRAM DELIVERS:
E(40,1.1)= +8.2952134128634"-003
E(41,1.1)= +8.0936587235982"-003
E(42,1.1)= +7.9016599781006"-003
EXP(50.1)*E(1,50.1)= +1.9576696324723"-002
SOURCE TEXT(S):
"CODE" 35080;
"CODE" 35081;
"CODE" 35086;
"CODE" 35087;