SymPy subs().evalf() vs evalf(subs=dict(...))

[s5s]

Hi,

First, apologies if this has been discussed before. I'm browsing through the tutorial and encountered what appears to be an inaccuracy in the documentation. In particular, I am reading Basic Operations section which states that

To numerically evaluate an expression with a Symbol at a point, we might use subs followed by evalf, but it is more efficient and numerically stable to pass the substitution to evalf using the subs flag, which takes a dictionary of Symbol: point pairs.

 However, when I test this in jupyter lab using "%timeit" it appears the reverse is true. I'll avoid posting images, but I got the following results:

expr = cos(x) - 2*sin(exp(x))

%timeit expr.subs(x, 0).evalf()

92.9 µs ± 628 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit expr.evalf(subs=dict(x=0))

264 µs ± 7.21 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

expr = cos(2*x)

%timeit expr.evalf(subs={x: 2.4})

73.8 µs ± 3.28 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit expr.subs(x, 2.4).evalf()

35.1 µs ± 434 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Perhaps the documentation needs to elaborate more what it means by "stable" and "efficient" and under what conditions using one is preferable over the other because from the tests I performed, it appears using subs().evalf() is the way to go. The docs are either incorrect or out of date.

Best regards

[Oscar Benjamin]

Hi,

The evalf function is intended to compute the result to any given
specified precision and return a result such that it is correct to
almost every digit. However evalf needs access to the full expression
to be able to do that. If you first substitute a Float in because the
expression will partially or fully evaluate before evalf gets a
chance:

In [5]: expr = cos(2*x)

In [6]: expr.evalf(50, subs={x:2.4})
Out[6]: 0.087498983439446392365833650247410931894111629665389

In [7]: expr.subs(x, 2.4).evalf(50)
Out[7]: 0.087498983439446398335803678492084145545959472656250

Note that eval isn't doing anything in the second example because the
subs already evaluated to a Float (ignoring our request for 50
digits):

In [10]: expr.subs(x, 2.4)
Out[10]: 0.0874989834394464

The first example gives an accurate result for cos(2x) based on the
true value of the float 2.4 (which is not exactly equal to 12/5). If
we really want to calculate cos(2x) for x=12/5 with minimal rounding
error we should use:

In [9]: expr.subs(x, Rational('2.4')).evalf(50)
Out[9]: 0.087498983439446569320215257649487633957449890596100

--
Oscar

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/03831493-431b-45e7-809f-6cc691cab191n%40googlegroups.com.

[Davide Sandona']

Hello Oscar,

why is the float approximation of Rational('2.4') more precise of Float('2.4')?

Davide.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxTTpyuVr%3DxhVOXyxPPUXyb9%2B1cCE2n9s2_Yqeig6u-aWw%40mail.gmail.com.

[Aaron Meurer]

Probably the word "efficient" shouldn't be there. The key thing is
that it can potentially be more accurate, because it is able to be
better at avoiding loss of significance due to numerical cancellation.
For example

>>> (x + y + z).subs({x: 1, y: 1e-100, z: -1})
0
>>> (x + y + z).evalf(subs={x: 1, y: 1e-100, z: -1})
1.00000000000000e-100

Aaron Meurer

[Aaron Meurer]

Rational('2.4') creates the exact rational number 12/5. It can then be
evaluated to as many digits as needed. Float('2.4') creates the
floating point version of this rational with 15 digits of precision
(the default). This corresponds to this rational number

>>> Rational(Float('2.4'))
5404319552844595/2251799813685248

If you later wanted to evaluate this to more digits, you would be
evaluating that rational number to more digits, not 12/5.

>>> Rational('2.4').evalf(50)
2.4000000000000000000000000000000000000000000000000
>>> Float('2.4').evalf(50)
2.3999999999999999111821580299874767661094665527344

Aaron Meurer

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAO%3D1Z08%2B_28O_y2zumiA9AB4844Hb9%3DFbXws-qLbW3%2BNgymwDw%40mail.gmail.com.