SymPy 1.0中的集成错误

I'm trying to write my mathcad model in python language, but I get some mistake. The integration function should look like this:

我正在尝试用python语言编写我的mathcad模型，但是我犯了一些错误。集成功能应如下所示：

In python I wrote such code

在python中我写了这样的代码

    from __future__ import division
    import sympy as sp
    import numpy as np
    import math
    from pylab import *

    print(sp.__version__)

    s  =  sp.Symbol('s')
    x = sp.symbols('x')

    t_start = 11
    t_info = 1
    t_transf = 2
    t_stat_analyze = 3
    t_repeat = 3.2
    P = 0.1

    def M1(s):
        return P/(t_info*t_start*t_stat_analyze*t_transf*(1 - (-P + 1)/(t_repeat*t_transf*(s + 1/t_repeat)*(s + 1/t_transf)))*(s + 1/t_info)*(s + 1/t_start)*(s + 1/t_stat_analyze)*(s + 1/t_transf)**2) + P/(        t_info*t_start*t_stat_analyze*t_transf*(1 - (-P + 1)/(t_repeat*t_transf*(s + 1/t_repeat)*(s + 1/t_transf)))*(s + 1/t_info)*(s + 1/t_start)*(s + 1/t_stat_analyze)**2*(s + 1/t_transf)) + P/(t_info*t_start*t_stat_analyze*t_transf*(1 - (-P +         1)/(t_repeat*t_transf*(s + 1/t_repeat)*(s + 1/t_transf)))*(s + 1/t_info)*(s + 1/t_start)**2*(s + 1/t_stat_analyze)*(s + 1/t_transf)) + P/(t_info*t_start*t_stat_analyze*t_transf*(1 - (-P + 1)/(t_repeat*t_transf*(s + 1/t_repeat)*(s + 1/        t_transf)))*(s + 1/t_info)**2*(s + 1/t_start)*(s + 1/t_stat_analyze)*(s + 1/t_transf)) - P*(-(-P + 1)/(t_repeat*t_transf*(s + 1/t_repeat)*(s + 1/t_transf)**2) - (-P + 1)/(t_repeat*t_transf*(s + 1/t_repeat)**2*(s + 1/t_transf)))/(        t_info*t_start*t_stat_analyze*t_transf*(1 - (-P + 1)/(t_repeat*t_transf*(s + 1/t_repeat)*(s + 1/t_transf)))**2*(s + 1/t_info)*(s + 1/t_start)*(s + 1/t_stat_analyze)*(s + 1/t_transf))

    def M2(s):
        return 2*P*((s + 1/t_transf)**(-2) + 1/((s + 1/t_stat_analyze)*(s + 1/t_transf)) + (s + 1/t_stat_analyze)**(-2) + 1/((s + 1/t_start)*(s + 1/t_transf)) + 1/((s + 1/t_start)*(s + 1/t_stat_analyze)) + (s + 1/t_start)**(-2) + 1/((s + 1/        t_info)*(s + 1/t_transf)) + 1/((s + 1/t_info)*(s + 1/t_stat_analyze)) + 1/((s + 1/t_info)*(s + 1/t_start)) + (s + 1/t_info)**(-2) - (P - 1)*((s + 1/t_transf)**(-2) + 1/((s + 1/t_repeat)*(s + 1/t_transf)) + (s + 1/t_repeat)**(-2))/(        t_repeat*t_transf*(1 + (P - 1)/(t_repeat*t_transf*(s + 1/t_repeat)*(s + 1/t_transf)))*(s + 1/t_repeat)*(s + 1/t_transf)) - (P - 1)*(1/(s + 1/t_transf) + 1/(s + 1/t_repeat))/(t_repeat*t_transf*(1 + (P - 1)/(t_repeat*t_transf*(s + 1/        t_repeat)*(s + 1/t_transf)))*(s + 1/t_repeat)*(s + 1/t_transf)**2) - (P - 1)*(1/(s + 1/t_transf) + 1/(s + 1/t_repeat))/(t_repeat*t_transf*(1 + (P - 1)/(t_repeat*t_transf*(s + 1/t_repeat)*(s + 1/t_transf)))*(s + 1/t_repeat)*(s + 1/        t_stat_analyze)*(s + 1/t_transf)) - (P - 1)*(1/(s + 1/t_transf) + 1/(s + 1/t_repeat))/(t_repeat*t_transf*(1 + (P - 1)/(t_repeat*t_transf*(s + 1/t_repeat)*(s + 1/t_transf)))*(s + 1/t_repeat)*(s + 1/t_start)*(s + 1/t_transf)) - (P - 1)*(1/(        s + 1/t_transf) + 1/(s + 1/t_repeat))/(t_repeat*t_transf*(1 + (P - 1)/(t_repeat*t_transf*(s + 1/t_repeat)*(s + 1/t_transf)))*(s + 1/t_info)*(s + 1/t_repeat)*(s + 1/t_transf)) + (P - 1)**2*(1/(s + 1/t_transf) + 1/(s + 1/t_repeat))**2/(        t_repeat**2*t_transf**2*(1 + (P - 1)/(t_repeat*t_transf*(s + 1/t_repeat)*(s + 1/t_transf)))**2*(s + 1/t_repeat)**2*(s + 1/t_transf)**2))/(t_info*t_start*t_stat_analyze*t_transf*(1 + (P - 1)/(t_repeat*t_transf*(s + 1/t_repeat)*(s + 1/        t_transf)))*(s + 1/t_info)*(s + 1/t_start)*(s + 1/t_stat_analyze)*(s + 1/t_transf))

    T_realyze = M1(0)
    D = M2(0)-M1(0)**2

    alpha = T_realyze**2/D
    myu = T_realyze/D

    def F(t):
        if t<0:
            return 0
        else:
            return sp.integrate((myu**alpha)/(sp.gamma(alpha)*(x**(alpha-1))*sp.exp(myu*x)), (x, 0, t))

    t=arange(0, 200, 1)
    for i in t:
        print(F(i))
        i = i+1

So, when I'm trying to execute it, I had such error in

所以，当我试图执行它时，我有这样的错误

    return sp.integrate

function:

功能：

    $ python2.7 nta.py
    1.0
    ('T_realyze  =  ', 63.800000000000026)
    ('D  =  ', 2696.760000000001)
    ('alpha  =  ', 1.5093816283243602)
    ('myu  =  ', 0.02365801925273291)
    0
    ('myu*x  =  ', 0.0236580192527329*x)
    ('sp.exp(myu*x)', exp(0.0236580192527329*x))
    0
    1
    ('myu*x  =  ', 0.0236580192527329*x)
    ('sp.exp(myu*x)', exp(0.0236580192527329*x))
    Traceback (most recent call last):
      File "nta.py", line 48, in <module>
        print(F(i))
      File "nta.py", line 43, in F
        return sp.integrate((myu**alpha)/(sp.gamma(alpha)*(x**(alpha-1))*sp.exp(myu*x)), (x, 0, t))
      File "/root/anaconda2/lib/python2.7/site-packages/sympy/integrals/integrals.py", line 1280, in integrate
        risch=risch, manual=manual)
      File "/root/anaconda2/lib/python2.7/site-packages/sympy/integrals/integrals.py", line 486, in doit
        conds=conds)
      File "/root/anaconda2/lib/python2.7/site-packages/sympy/integrals/integrals.py", line 887, in _eval_integral
        h = heurisch_wrapper(g, x, hints=[])
      File "/root/anaconda2/lib/python2.7/site-packages/sympy/integrals/heurisch.py", line 130, in heurisch_wrapper
        unnecessary_permutations)
      File "/root/anaconda2/lib/python2.7/site-packages/sympy/integrals/heurisch.py", line 657, in heurisch
        solution = _integrate('Q')
      File "/root/anaconda2/lib/python2.7/site-packages/sympy/integrals/heurisch.py", line 646, in _integrate
        numer = ring.from_expr(raw_numer)
      File "/root/anaconda2/lib/python2.7/site-packages/sympy/polys/rings.py", line 371, in from_expr
        raise ValueError("expected an expression convertible to a polynomial in %s, got %s" % (self, expr))
    ValueError: expected an expression convertible to a polynomial in Polynomial ring in _x0, _x1, _x2, _x3 over RR[_A0,_A1,_A2,_A3,_A4,_A5,_A6,_A7,_A8,_A9,_A10,_A11,_A12,_A13,_A14,_A15,_A16,_A17,_A18,_A19,_A20,_A21,_A22,_A23,_A24,_A25,_A26,_A27,_A28,_A29,_A30,_A31,_A32,_A33,_A34] with lex order, got 0.50938162832436*_x3**2.96316463805253*(_A0 + _A10*_x0*_x1 + 2*_A11*_x1*_x3 + _x0**2*_A12 + _A14*_x0*_x2 + _A2*_x0 + 2*_A20*_x0*_x3 + _A24*_x1*_x2 + _x2**2*_A27 + 2*_A28*_x3 + _x1**2*_A30 + 3*_x3**2*_A31 + 2*_A6*_x2*_x3 + _A8*_x2 + _A9*_x1) + 1.50938162832436*_x3**4.92632927610506*(_A10*_x1*_x3 + 2*_A12*_x0*_x3 + _A13*_x1*_x2 + _A14*_x2*_x3 + 2*_A15*_x0 + _A16*_x2 + _x2**2*_A18 + _A2*_x3 + _x3**2*_A20 + _A21 + _x1**2*_A3 + 2*_A33*_x0*_x2 + _A34*_x1 + 3*_x0**2*_A5 + 2*_A7*_x0*_x1) - _A10*_x0*_x3 - _x3**2*_A11 - _A13*_x0*_x2 - _x2**2*_A17 - 2*_A19*_x1*_x2 - _A22 - _A24*_x2*_x3 - 2*_A25*_x1 - 3*_x1**2*_A29 - 2*_A3*_x0*_x1 - 2*_A30*_x1*_x3 - _A34*_x0 - _A4*_x2 - _x0**2*_A7 - _A9*_x3 + _x2*_x3 + 0.0236580192527329*_x2*(_A13*_x0*_x1 + _A14*_x0*_x3 + _A16*_x0 + 2*_A17*_x1*_x2 + 2*_A18*_x0*_x2 + _x1**2*_A19 + 2*_A23*_x2 + _A24*_x1*_x3 + 3*_x2**2*_A26 + 2*_A27*_x2*_x3 + _A32 + _x0**2*_A33 + _A4*_x1 + _x3**2*_A6 + _A8*_x3)

1 个解决方案

#1

Sympy appears to have difficulties evaluating the integral with floating point coefficients (in this case). However, it can find the integral in closed form when the constants of the integrand expression are symbolic.

Sympy似乎很难评估浮点系数的积分（在这种情况下）。但是，当被积函数表达式的常数是符号时，它可以以闭合形式找到积分。

a, b, c, t = sp.symbols('a,b,c,t', positive = True)
f = sp.Integral(a * sp.exp(-c*x)/(x**b),(x,0,t)).doit()
print f

Output:

输出：

-a*(-b*c**b*gamma(-b + 1)*lowergamma(-b + 1, 0)/(c*gamma(-b + 2)) + c**b*gamma(-b + 1)*lowergamma(-b + 1, 0)/(c*gamma(-b + 2))) + a*(-b*c**b*gamma(-b + 1)*lowergamma(-b + 1, c*t)/(c*gamma(-b + 2)) + c**b*gamma(-b + 1)*lowergamma(-b + 1, c*t)/(c*gamma(-b + 2)))

You can substitute the constants in this expression to get numerical results as follows (here, I use an example value of t=4):

您可以替换此表达式中的常量来获得如下的数值结果（这里，我使用t = 4的示例值）：

f.subs({a:(myu**alpha)/sp.gamma(alpha), b:(alpha-1), c:myu, t:4}).n()

0.0154626407404632

Another option is to use quad from scipy (again using t=4):

另一个选择是使用scipy中的quad（再次使用t = 4）：

from scipy.integrate import quad

quad(lambda x: (myu**alpha)/(sp.gamma(alpha)*(x**(alpha-1))*sp.exp(myu*x)), 0 ,4)[0]

0.015462640740458165

#1