[SciPy-User] SciPy-User Digest, Vol 110, Issue 21
Moore, Eric (NIH/NIDDK) [F]
eric.moore2@nih....
Fri Oct 12 15:11:10 CDT 2012
> -----Original Message-----
> From: "Claas H. Köhler" [mailto:claas.koehler@dlr.de]
> Sent: Friday, October 12, 2012 12:39 PM
> To: scipy-user@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@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@dlr.de
>
> www.DLR.de/EOC
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
More information about the SciPy-User
mailing list