Hi Raymond san, all,
I am looking at the expintegral.lisp line 534 where the continued fraction
or power series
is chosen.
According to DLMF Section 6.18 (http://dlmf.nist.gov/6.18 which is an
updated publication of A&S),
"For small or moderate values of x and z, the expansion in power series
(?6.6) or in series of spherical Bessel functions (?6.10(ii)) can be used.
....
Lastly, the continued fraction (6.9.1) can be used if z is
bounded away from the origin. Convergence becomes slow when z is near the
negative real axis, however."
So, I would like to propose:
if z is close to origin or z is close to the negative real axis, we expand
in power series,
otherwise, we expand in continued fraction. The actual patch comes below.
It would be great if someone could verify (or modify if necessary) this
patch.
Thanks and best regards,
Yasuaki Honda, Chiba Japan
[iMacG5-2F-2:maxima-git-development/maxima/src] yasube% git diff
expintegral.lisp
diff --git a/src/expintegral.lisp b/src/expintegral.lisp
index 6bf7d82..e27339e 100644
--- a/src/expintegral.lisp
+++ b/src/expintegral.lisp
@@ -531,7 +531,7 @@
(format t "~& : z = ~A~%" z))
(cond
- ((and (>= (realpart z) 0) (> (abs z) 1.0))
+ ((not (or (< (abs z) 5.0) (and (< (realpart z) 0.0) (< (abs
(imagpart z)) 2.0))))
;; We expand in continued fractions.
(when *debug-expintegral*
(format t "~&We expand in continued fractions.~%"))
2012?4?3?16:24 ???? <yasuaki.honda at gmail.com>:
> Hi Raymond san,
>
> Thanks for looking into this problem. When you create a fix and want to
> test it,
> drop me the code and I will test it.
>
> Thanks and best regards,
> Yasuaki Honda, Chiba, Japan
>
> 2012?4?3?0:52 Raymond Toy <toy.raymond at gmail.com>:
>
>
>>
>> 2012/4/2 ???? <yasuaki.honda at gmail.com>
>>
>>> Dear all,
>>>
>>> It seems that the Maxima function expintegral_ei() behaves unexpectedly
>>> for
>>> some complex arguments. (There should be some bugs in the
>>> implementation).
>>>
>>> Specifically, considering a line log(20.0)/2+%i*t on the complex plane
>>> and t gets
>>> increasing, such as t=30, 40, 50.
>>> the values can
>>> (%i22) expintegral_ei(log(20.0)/2+30*%i);
>>> (%o22) 3.116181583331404*%i-.1466880831183789
>>>
>>> (%i23) expintegral_ei(log(20.0)/2+40*%i);
>>> (%o23) 3.205076248461513*%i+.1490949857975915
>>>
>>> (%i24) expintegral_ei(log(20.0)/2+50*%i);
>>> (%o24) 1766.649087960532-1098.205534849491*%i
>>>
>>> So, (%o24) indicates that somewhere betwee log(20.0)/2+40*%i and
>>> log(20.0)/2+50*%i,
>>> there is a point where expintegral_ei() starts unstable ( not converge).
>>> This function,
>>> however, should converge for these values. (See below for the values
>>> obtained using
>>> mpmath in python).
>>>
>>> Thanks for reporting this. It looks like it's really an issue with
>> expintegral_e. Maxima compute expintegral_ei(z) using expintegral_e(1,-z),
>> and the expIntegral_e appears to have problems because it's using a
>> series. I think if the continued fraction is used instead, the problem
>> goes away.
>>
>> I'll look into it.
>>
>> Ray
>>
>>
>>
>