python - SymPy 1.0 中的集成错误

Question

我正在尝试用 python 语言编写我的 mathcad 模型，但我遇到了一些错误。集成函数应如下所示:

在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

所以，当我尝试执行它时，我在

    return sp.integrate

功能:

    $ 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)

Answer 1

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

输出:

-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)))

您可以替换此表达式中的常量以获得如下数值结果(此处，我使用 t=4 的示例值):

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

0.0154626407404632

另一种选择是使用 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