[SciPy-User] SciPy-User Digest, Vol 110, Issue 21

The Helmbolds helmrp@yahoo....
Fri Oct 12 05:48:25 CDT 2012


> On 09/10/12 19:12, Pauli Virtanen wrote:
>>  09.10.2012 19:28, "Claas H. K?hler" kirjoitti:
>>>  I have a question regarding the error function scipy.special.erf:
>>> 
>>>  Is it intended, that the erf of an imaginary argument yields a 
> non-vanishing real-part?
>>> 
>>>  I get e.g.
>>>  erf(1j)= 1.6504257587975431j
>>>  erf(5j)= (1+8298273879.8992386j)
>>> 
>>>  The first result is what I would expect in accordance with Wolfram 
> alpha. The second result, however,
>>>  has a real part of unity. As far as I know, the real part of erf should 
> always vanish for purely
>>>  imaginary numbers.
>>> 
>>>  Any support would be appreciated.
>> 
>>  The reason here is that the ye olde complex erf Fortran implementation
>>  that Scipy has uses the asymptotic expansion (Abramowitz & Stegun
>>  7.1.23) to compute large-argument values. The asymptotic series is for
>>  erfc, and one always gets Re erf = 1 along the imaginary axis.
>> 
>>  Of course, this is somewhat naive. While it does produce reasonable
>>  relative accuracy as a complex number, the accuracy of the real and
>>  imaginary parts separately is not necessarily OK near the imaginary axis.
>> 
>>  The issue with Scipy here is twofold -- first, there are no better
>>  existing special function libraries we could use, or at least I'm not
>>  aware of them. Second, writing these from scratch takes time and
>>  expertise and nobody has so far volunteered to do any work in this
>>  direction.
>> 
> Thanks for the quick response!
> 
> The bottom line is that erf is actually not (correctly) implemented for complex 
> arguments, if I
> understand you correctly.
> 
> I suspect there are good reasons to provide a function which is known to yield 
> incorrect results, so
> that throwing a type error is not an option? (This is what erfc does on my 
> machine)
> 
> However, adding a warning when called with complex arguments could be helpful to 
> prevent naiive use
> as in my case. Adding this important piece of information to the docs would not 
> harm either, from my
> point of view.
> 
> In any case, thanks for the quick support.
> 
> Regards
> Claas

On my system, I get the correct answers if I'm careful about the call to erf. 
If I call erf with a single real value, I get the ordinary (not the complex) error function value. 
If I call erf with a NumPy array or a Python sequence, I get the complex error function returned.
I do not think SciPy's erf is supposed to be called with a complex number. 

For example:
>>> special.erf(1j)
1.6504257587975431j                  # Wrong answer!
>>> special.erf((0,1))
array([ 0.        ,  0.84270079])        # Right answer.

Two more examples:
>>> for y in range(-10, 11):
 temp = special.erf((0,y))
 print y, temp                                 # Calling with a sequence, returns a NumPy array

-10 [ 0.         -1.]
-9 [ 0.           -1.]
-8 [ 0.           -1.]
-7 [ 0.           -1.]
-6 [ 0.           -1.]
-5 [ 0.           -1.]
-4 [ 0.         -0.99999998]
-3 [ 0.         -0.99997791]
-2 [ 0.         -0.99532227]
-1 [ 0.         -0.84270079]
0 [ 0.           0.]
1 [ 0.          0.84270079]
2 [ 0.          0.99532227]
3 [ 0.          0.99997791]
4 [ 0.          0.99999998]
5 [ 0.          1.]
6 [ 0.          1.]
7 [ 0.          1.]
8 [ 0.          1.]
9 [ 0.          1.]
10 [ 0.         1.]


OTOH--------------------------------------------------------------------------------------------
>>> for y in range(-10, 11):
 temp = special.erf(y)
 print y, temp                            # Calling with a (scalar)  real value returns a (scalar) real value.

-10    -1.0
-9     -1.0
-8     -1.0
-7     -1.0
-6     -1.0
-5     -0.999999999998
-4     -0.999999984583
-3     -0.999977909503
-2     -0.995322265019
-1     -0.84270079295
0      0.0
1     0.84270079295
2     0.995322265019
3     0.999977909503
4     0.999999984583
5     0.999999999998
6     1.0
7     1.0
8     1.0
9     1.0
10   1.0

Bob and Paula H        


More information about the SciPy-User mailing list