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

Fri Oct 12 16:11:10 EDT 2012

> -----Original Message-----
> From: "Claas H. Köhler" [mailto:claas.koehler at dlr.de]
> Sent: Friday, October 12, 2012 12:39 PM
> To: scipy-user at scipy.org
> Subject: Re: [SciPy-User] SciPy-User Digest, Vol 110, Issue 21
> 
> 
> 
> On 12/10/12 12:48, The Helmbolds wrote:
> >> 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 er
> ror function returned.
> > I do not think SciPy's erf is supposed to be called with a complex
> number.
> According to the docs it is. Otherwise I would expect to see a domain
> error, similar to erfc.
> 
> Regards
> Claas
> 
> 
> >
> > For example:
> >>>> special.erf(1j)
> > 1.6504257587975431j                  # Wrong answer!

This is the right answer for erf(1j).  And the behavior you detail below is exactly how ufuncs work.  You get a  scalar back if you provide one, and get an array back if you provide a sequence. 

Eric.

> >>>> 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
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User at scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user
> >
> 
> --
> Deutsches Zentrum für Luft- und Raumfahrt e.V. (DLR)
> Institut für Methodik der Fernerkundung | Experimentelle Verfahren |
> Münchner Str | 82234 Weßling
> 
> Claas H. Köhler
> Telefon 08153 28-1274 | Telefax 08153 28-1337 | claas.koehler at dlr.de
> 
> www.DLR.de/EOC
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user