Ошибка интеграции в SymPy 1.0

Question

Я пытаюсь написать свою модель 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

У Симпы, кажется, трудности с вычислением интеграла с коэффициентами с плавающей запятой (в данном случае). Однако он может найти интеграл в замкнутой форме, когда константы выражения подынтегральной функции являются символическими.

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 

Другой вариант заключается в использовании quad из SciPy (опять-таки с помощью 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