(Credit to Ralf Fritzsch)
[6.60] Faster, more accurate exponential integrals
A post on the comp.sys.hp48 newsgroup announced that the HP-49G did not accurately find a
numerical value for this integral
¶
−3.724
∞
e
−x
x
dx
Except for a sign change, this is identical to the special function called the exponential integral:
Ei
(
x
)
=−
¶
−x
∞
e
−x
x
dx
The TI-89 / TI-92 Plus also cannot find accurate values for this integral, and it is no wonder, since at
the integrand singularity at x = 0, the limit is +infinity from the right and -infinity from the left. One
expedient solution is a routine specifically designed to find the exponential integral:
ei(x)
Func
©(x) exponential integral Ei
©13jul02/dburkett@infinet.com
local eps,euler,fpmin,maxit,k,fact,prev,sum1,term
6⁻12→eps © Desired relative accuracy
.57721566490153→euler © Euler's constant
100→maxit © Maximum iterations allowed for convergence
1⁻900→fpmin © Number near floating-point minimum
if x≤0:return "ei arg error" © Argument must be > 0
if x<fpmin then © Handle very small arguments
return euler+ln(x)
elseif x≤⁻ln(eps) then © Use power series for x < 25.840
0→sum1
1→fact
for k,1,maxit
fact*x/k→fact
fact/k→term
sum1+term→sum1
if term<eps*sum1
return sum1+ln(x)+euler
endfor
return "ei failed series" © Return error string if convergence fails
else © Use asymptotic expansion for large arguments
0→sum1
1→term
for k,1,maxit
term→prev
term*k/x→term
if term<eps:goto l1
if term<prev then
sum1+term→sum1
else
sum1-prev→sum1
goto l1
endif
endfor
lbl l1
return ℯ^(x)*(1+sum1)/x
6 - 114
To Next Page
Loading...
Need help?
General
