Интегрирование кусочно с иррациональным показателем дает ошибку - PullRequest
0 голосов
/ 26 мая 2018

При попытке интегрировать кусочную функцию с иррациональным показателем я наткнулся на следующую проблему:

import sympy as sym

x, y = sym.symbols(['x', 'y'])
pdf = sym.Piecewise((0, x < 1),
                    (x ** -2.5, x < 2),
                    (0, True))
cdf = sym.integrate(pdf, (x, -sym.oo, y))

Это дает следующую ошибку:

ValueError                                Traceback (most recent call last)
<ipython-input-3-c1f460aac84e> in <module>()
      5                     (x ** -2.5, x < 2),
      6                     (0, True))
----> 7 cdf = sym.integrate(pdf, (x, -sym.oo, y))
      8 print(cdf)

~/.local/lib/python3.6/site-packages/sympy/integrals/integrals.py in integrate(*args, **kwargs)
   1293     if isinstance(integral, Integral):
   1294         return integral.doit(deep=False, meijerg=meijerg, conds=conds,
-> 1295                              risch=risch, manual=manual)
   1296     else:
   1297         return integral

~/.local/lib/python3.6/site-packages/sympy/integrals/integrals.py in doit(self, **hints)
    551                                        for f in integrals])
    552                         try:
--> 553                             evalued = Add(*others)._eval_interval(x, a, b)
    554                             function = uneval + evalued
    555                         except NotImplementedError:

~/.local/lib/python3.6/site-packages/sympy/functions/elementary/piecewise.py in _eval_interval(self, sym, a, b)
    229                         rep = b
    230                         val = e._eval_interval(sym, mid, b)
--> 231                         val += self._eval_interval(sym, a, mid)
    232                     elif (a > upper) == True:
    233                         mid = upper

~/.local/lib/python3.6/site-packages/sympy/functions/elementary/piecewise.py in _eval_interval(self, sym, a, b)
    282                     super(Piecewise, expr)._eval_interval(sym, Max(a, int_a), Min(b, int_b)))
    283             else:
--> 284                 ret_fun += expr._eval_interval(sym, Max(a, int_a), Min(b, int_b))
    285         return mul * ret_fun
    286 

~/.local/lib/python3.6/site-packages/sympy/core/expr.py in _eval_interval(self, x, a, b)
    832             else:
    833                 domain = Interval(b, a)
--> 834             singularities = list(solveset(self.cancel().as_numer_denom()[1], x, domain = domain))
    835             for s in singularities:
    836                 if a < s < b:

~/.local/lib/python3.6/site-packages/sympy/solvers/solveset.py in solveset(f, symbol, domain)
    919         return result
    920 
--> 921     return _solveset(f, symbol, domain, _check=True)
    922 
    923 

~/.local/lib/python3.6/site-packages/sympy/solvers/solveset.py in _solveset(f, symbol, domain, _check)
    709             result += solns
    710     else:
--> 711         lhs, rhs_s = inverter(f, 0, symbol)
    712         if lhs == symbol:
    713             # do some very minimal simplification since

~/.local/lib/python3.6/site-packages/sympy/solvers/solveset.py in <lambda>(f, rhs, symbol)
    677     else:
    678         inverter_func = invert_complex
--> 679     inverter = lambda f, rhs, symbol: inverter_func(f, rhs, symbol, domain)
    680 
    681     result = EmptySet()

~/.local/lib/python3.6/site-packages/sympy/solvers/solveset.py in invert_real(f_x, y, x, domain)
    116     the domain to ``S.Reals`` before inverting.
    117     """
--> 118     return _invert(f_x, y, x, domain)
    119 
    120 

~/.local/lib/python3.6/site-packages/sympy/solvers/solveset.py in _invert(f_x, y, x, domain)
     98 
     99     if domain.is_subset(S.Reals):
--> 100         x1, s = _invert_real(f_x, FiniteSet(y), x)
    101     else:
    102         x1, s = _invert_complex(f_x, FiniteSet(y), x)

~/.local/lib/python3.6/site-packages/sympy/solvers/solveset.py in _invert_real(f, g_ys, symbol)
    184             else:
    185                 if not base.is_positive:
--> 186                     raise ValueError("x**w where w is irrational is not "
    187                                      "defined for negative x")
    188                 return _invert_real(base, res, symbol)

ValueError: x**w where w is irrational is not defined for negative x

Но если я возьму толькочасть с показателем степени и интегрировать ее, она отлично работает:

pdf = x ** -2.5
cdf = sym.integrate(pdf, (x, -sym.oo, y))
print(cdf)

Это дает:

-0.666666666666667*_y**(-1.5)

Вопрос:
Как я могу интегрироватьпервая кусочная функция?

То, что я пытался:
Указание y неотрицательным:

y = sym.Symbol('y', nonnegative=True)

Так я не получу первую ошибку, но результатстановится довольно странно:

0.666666666666667*Min(1, y)**(-1.5) - 0.666666666666667*Min(2, y)**(-1.5)

, который я не могу использовать позже:

eq = sym.Eq(x, cdf)
inverse_solutions = sym.solve(eq, y, rational=False)  

дает:

NotImplementedError: multiple generators [Min(1, y)**1.5, Min(2, y)**1.5]
No algorithms are implemented to solve equation x - 0.666666666666667*Min(1, y)**(-1.5) + 0.666666666666667*Min(2, y)**(-1.5)

1 Ответ

0 голосов
/ 26 мая 2018

Похоже, я нашел обходной путь.
Мне просто пришлось упростить pdf на sympy.nsimplify, и ошибка исчезла:

import sympy as sym

x, y = sym.symbols(['x', 'y'])
pdf = sym.Piecewise((0, x < 1),
                    (x ** -2.5, x < 2),
                    (0, True))
pdf = sym.nsimplify(pdf)
print(pdf)
cdf = sym.integrate(pdf, (x, -sym.oo, y))
print(cdf)

дает:

Piecewise((0, x < 1), (x**(-5/2), x < 2), (0, True))
Piecewise((0, y < 1), (2/3 - 2/(3*y**(3/2)), y < 2), (-sqrt(2)/6 + 2/3, True))
...