PolynomialDivisionFailed in integrate #18386
Comments
|
It works with In [5]: x = sympy.Symbol("x")
...: a = sympy.Symbol("a")
...: f = x**2 + a/(x + Rational(11, 10))
...: g = abs(x)
...: sympy.integrate((g - f)**2, (x, -1, 0))
Out[5]:
2
100⋅a 6⋅a 11⋅a⋅log(11) 1
────── - ─── + ──────────── + ──
11 5 50 30 |
|
i ran into the same problem, i've attached an interactive python notebook with which to reproduce my case: reproduction.ipynb.txt (GitHub doesn't allow uploading my workaround (i didn't find this ticket until after) was to use sympy variables instead of floats, do the integration and only then substitute in the real values. this makes it of course more cumbersome but at least it works. this was with SymPy 1.10.1 and python 3.8.13 (if you prefer to have an own ticket for my case i can also raise a separate one; but i presume based on the description that it's the same as it's the same error and also caused by floats). |
I think it's good to add the example here. It would be better just to paste a short demonstration as a code block like in the OP though (rather than attaching a notebook). |
|
ok, there you go: python codeimport sympy as sym
import scipy.constants as const
from IPython.display import display,Math
a = 2.4e-2 # m
Q1 = 10e-9 # C
Q2 = -10e-9 # C
x = sym.symbols('x', real=True, positive=True)
P1 = 2*a
P2 = 1*a
r1 = sym.sqrt(a**2 + x**2)
E1 = 1/(4*const.pi*const.epsilon_0) * Q1 / r1**2 * 1
r2 = sym.sqrt(a**2 + (-3*a + x)**2)
E2 = 1/(4*const.pi*const.epsilon_0) * Q2 / r2**2 * -1 # direction is inversed
E = E1 + E2
display(Math('E = ' + sym.latex(E)))
U = E.expand().integrate((x,P1,P2))
display(Math(r'U \approx %.2f~\mathrm V' %U))
qe = -const.elementary_charge
dW = qe * U
display(Math(r'\Delta W \approx %.2f \cdot 10^{-16}~\mathrm J' %(dW * 1e16)))error message & stack trace |
|
Thanks. Really it's better if you just strip out all the unrelated stuff though in which case your example is just this: from sympy import *
x = symbols('x')
p = ((179.751035845223*x**2 - 12.9420745808561*x + 0.569451281557668)
/(1.0*x**4 - 0.144*x**3 + 0.006336*x**2 - 8.2944e-5*x + 3.31776e-6))
integrate(p, x)A simple workaround is to convert the floats to rational: In [6]: nsimplify(p).integrate().n()
Out[6]:
⎛ 2 9 ⎞ ⎛ 2 18⋅x 18 ⎞
- 8.9031339031339e-14⋅log⎜x + ─────⎟ + 8.9031339031339e-14⋅log⎜x - ──── + ────⎟ + 3744.8132467754
⎝ 15625⎠ ⎝ 125 3125⎠
⎛125⋅x⎞ ⎛125⋅x ⎞
9⋅atan⎜─────⎟ + 3744.81324677547⋅atan⎜───── - 3⎟
⎝ 3 ⎠ ⎝ 3 ⎠ |
|
The question is really if integrate should convert to rational internally to handle cases like this. |
|
The long division of polynomials fails when the leading coefficient of a remainder is nonzero because of a rounding error. If f and g have leading coefficients a and b, then it is expected that sympy/sympy/polys/densearith.py Lines 1450 to 1455 in ab414d1 (There are other places in densearith.py where this could also be done.) |
They all seem to follow def dup_sub_mul_strip(f, g, c, i, K):
"""Compute r = f - g*c*x**i where the leading term of f should cancel"""
G = dup_mul_term(g, c, i, K)
r = dup_sub(f, G, K)
if not K.is_exact and dup_degree(r) == dup_degree(f):
# remove leading term left over by rounding error
r = dup_strip(r[1:])
return rWhat happens if the new leading term should also cancel though? |
That is probably not easy to detect. Maybe the dup/dmp methods should convert the coefficients for the computation. |
|
Maybe def dup_sub_mul_strip(f, g, K):
"""Compute r = f - g where the leading terms should cancel"""
r = dup_sub(f, g, K)
if not K.is_exact and dup_degree(r) == dup_degree(f):
# remove leading term left over by rounding error
r = dup_strip(r[1:])
return r
I think conversion is better done at a higher level if it needs to be done. Probably stripping the leading term is the best that can be done if conversion is not used. |
Depends on the
1.1constant. Some values error, some don't.I am using
sympy 1.5.1on arch linux.The text was updated successfully, but these errors were encountered: