Integrationsfehler in SymPy 1.0

Question

Ich versuche, mein Mathcad-Modell in Python-Sprache zu schreiben, aber ich bekomme einen Fehler. Die Integrationsfunktion sollte wie folgt aussehen:

In Python ich einen solchen Code

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 

Also schrieb, als ich es auszuführen bin versucht, hatte ich solche Fehler in

Funktion:

$ 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 scheint Schwierigkeiten zu haben, das Integral mit Gleitkommakoeffizienten zu bewerten (in diesem Fall). Es kann jedoch das Integral in geschlossener Form finden, wenn die Konstanten des Integrandenausdrucks symbolisch sind.

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 

Ausgang:

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

Sie können die Konstanten in diesem Ausdruck ersetzen numerische Ergebnisse zu erhalten, wie folgt (hier verwende ich ein Beispiel Wert von t=4):

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

0.0154626407404632 

Eine weitere Option ist quad von scipy zu verwenden (wieder mit 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