Right, this error comes up when using floating point numbers. Things will be fixed in that regard in Mateusz's branch. I am concerned that that integral is almost definitely wrong, though. This is the final integral:
In [18]: Integral(curve_expr, (t, 0, oo)).subs(vals) Out[18]: ∞ ⌠ ⎮ ⎧ 0 for t - 478.515625⋅π < 0 ⎮ ⎪ ⎮ ⎪ 13.2075145209219⋅π ⎮ ⎨──────────────────────── for t - 478.515625⋅π ≥ 0 dt ⎮ ⎪ 2 ⎮ ⎪(0.000871222⋅t + 0.995) ⎮ ⎩ ⌡ 0 In [21]: print Integral(curve_expr, (t, 0, oo)).subs(vals) Integral(Piecewise((0, t - 478.515625*pi < 0), (13.2075145209219*pi/(0.000871222*t + 0.995)**2, t - 478.515625*pi >= 0)), (t, 0, oo)) It is clearly positive for t >= 478*pi. Aaron Meurer On Fri, May 31, 2013 at 9:10 AM, <[email protected]> wrote: > I don't. 0 sounds unlikely but may be correct given the maths. > > Just realised that again the error reported wasn't produced by the posted > code. I get the error shown in my previous post when I try to evaluate the > integral using a value for n other than 1 or 2 e.g. 1.1. Given that my > ultimate aim is to explore the influence of this n parameter on the shown > integral, do you think that the reported error indicates that I should be > looking to rework my formulations or is it simply indicative of a software > bug? > > Regards, > > Will > > > On Wednesday, May 29, 2013 4:10:12 AM UTC+1, Aaron Meurer wrote: >> >> Do you know what the right answer should be? It gives 0 in Mateusz's >> new-polys branch. >> >> Aaron Meurer >> >> On Tue, May 28, 2013 at 11:30 AM, Will Furnass >> <[email protected]> wrote: >> > Thanks Aaron. My code now looks as follows (T is now a lambda function >> > and >> > the integral is evaluated over t=0..oo): >> > >> > ... >> > T = lambda t: Piecewise((0, t < 0), (k * D * pi * dx * expr8 / Q, t >= >> > 0)) >> > >> > curve_expr = integrate(T(t - (x/V)), (x, 0, L)) >> > vals = { >> > P : 0.00022, >> > k : 0.5, >> > tau_s : 0.01, >> > tau_a : 2.0, >> > Q : 0.01, >> > D : 0.125, >> > L : 1225, >> > n : 2, >> > dx : 1} >> > Integral(curve_expr, (t, 0, oo)).subs(vals).doit() >> > >> > I now get the following error, regardless of whether I use the git >> > master or >> > mattpap's new-polys branch: >> > >> > AttributeError Traceback (most recent call >> > last) >> > <ipython-input-7-55341b44f4e2> in <module>() >> > 1 #plot(turb_curve_expr.subs(vals), (t, 0, 3600)) >> > ----> 2 Integral(turb_curve_expr, (t, 0, oo)).subs(vals).doit() >> > >> > /usr/local/lib/python2.7/dist-packages/sympy/integrals/integrals.pyc in >> > doit(self, **hints) >> > 886 antideriv = None >> > 887 else: >> > --> 888 antideriv = self._eval_integral(function, >> > xab[0], >> > meijerg=meijerg1, risch=risch, conds=conds) >> > 889 if antideriv is None and meijerg1 is True: >> > 890 ret = try_meijerg(function, xab) >> > >> > /usr/local/lib/python2.7/dist-packages/sympy/integrals/integrals.pyc in >> > _eval_integral(self, f, x, meijerg, risch, conds) >> > 1113 # Piecewise antiderivatives need to call special >> > integrate. >> > 1114 if f.func is Piecewise: >> > -> 1115 return f._eval_integral(x) >> > 1116 >> > 1117 # let's cut it short if `f` does not depend on `x` >> > >> > >> > /usr/local/lib/python2.7/dist-packages/sympy/functions/elementary/piecewise.pyc >> > in _eval_integral(self, x) >> > 197 def _eval_integral(self, x): >> > 198 from sympy.integrals import integrate >> > --> 199 return Piecewise(*[(integrate(e, x), c) for e, c in >> > self.args]) >> > 200 >> > 201 def _eval_interval(self, sym, a, b): >> > >> > /usr/local/lib/python2.7/dist-packages/sympy/utilities/decorator.pyc in >> > threaded_func(expr, *args, **kwargs) >> > 31 func(expr.rhs, *args, >> > **kwargs)) >> > 32 else: >> > ---> 33 return func(expr, *args, **kwargs) >> > 34 >> > 35 return threaded_func >> > >> > /usr/local/lib/python2.7/dist-packages/sympy/integrals/integrals.pyc in >> > integrate(*args, **kwargs) >> > 1587 >> > 1588 if isinstance(integral, Integral): >> > -> 1589 return integral.doit(deep=False, meijerg=meijerg, >> > conds=conds, risch=risch) >> > 1590 else: >> > 1591 return integral >> > >> > /usr/local/lib/python2.7/dist-packages/sympy/integrals/integrals.pyc in >> > doit(self, **hints) >> > 886 antideriv = None >> > 887 else: >> > --> 888 antideriv = self._eval_integral(function, >> > xab[0], >> > meijerg=meijerg1, risch=risch, conds=conds) >> > 889 if antideriv is None and meijerg1 is True: >> > 890 ret = try_meijerg(function, xab) >> > >> > /usr/local/lib/python2.7/dist-packages/sympy/integrals/integrals.pyc in >> > _eval_integral(self, f, x, meijerg, risch, conds) >> > 1235 try: >> > 1236 if conds == 'piecewise': >> > -> 1237 h = heurisch_wrapper(g, x, hints=[]) >> > 1238 else: >> > 1239 h = heurisch(g, x, hints=[]) >> > >> > /usr/local/lib/python2.7/dist-packages/sympy/integrals/heurisch.pyc in >> > heurisch_wrapper(f, x, rewrite, hints, mappings, retries, degree_offset, >> > unnecessary_permutations) >> > 124 >> > 125 res = heurisch(f, x, rewrite, hints, mappings, retries, >> > degree_offset, >> > --> 126 unnecessary_permutations) >> > 127 if not isinstance(res, Basic): >> > 128 return res >> > >> > /usr/local/lib/python2.7/dist-packages/sympy/integrals/heurisch.pyc in >> > heurisch(f, x, rewrite, hints, mappings, retries, degree_offset, >> > unnecessary_permutations) >> > 422 >> > 423 u_split = _splitter(denom) >> > --> 424 v_split = _splitter(Q) >> > 425 >> > 426 polys = list(v_split) + [ u_split[0] ] + special.keys() >> > >> > /usr/local/lib/python2.7/dist-packages/sympy/integrals/heurisch.pyc in >> > _splitter(p) >> > 399 return (c_split[0], q * c_split[1]) >> > 400 >> > --> 401 q_split = _splitter(cancel(q / s)) >> > 402 >> > 403 return (c_split[0]*q_split[0]*s, >> > c_split[1]*q_split[1]) >> > >> > /usr/local/lib/python2.7/dist-packages/sympy/integrals/heurisch.pyc in >> > _splitter(p) >> > 387 >> > 388 if _derivation(y) is not S.Zero: >> > --> 389 c, q = p.as_poly(y).primitive() >> > 390 >> > 391 q = q.as_expr() >> > >> > AttributeError: 'NoneType' object has no attribute 'primitive' >> > >> > >> > Regards, >> > >> > Will >> > >> > >> > On Tuesday, 28 May 2013 15:43:44 UTC+1, Aaron Meurer wrote: >> >> >> >> Are you sure you're using the git head? The T(t - (x/V)) in the >> >> integral is not allowed (we stopped making SymPy objects arbitrarily >> >> callable). I'm not sure if you meant for that to be a substitution or >> >> a multiplication, so I didn't try it further. Chances are it's fixed >> >> in https://github.com/sympy/sympy/pull/2126 though. >> >> >> >> Aaron Meurer >> >> >> >> On Tue, May 28, 2013 at 9:23 AM, Will Furnass >> >> <[email protected]> wrote: >> >> > Hi, >> >> > >> >> > I'm fairly new to sympy. I'm trying to evalute an integral >> >> > analytically >> >> > over the range 0..oo but when I run the last line in the script below >> >> > I >> >> > get >> >> > the following error (using the sympy git head). >> >> > >> >> > CoercionFailed: can't convert DMP([[1], []], ZZ, ZZ[_a1,_b1]) of type >> >> > <class >> >> > 'sympy.polys.polyclasses.DMP'> from ZZ[_a1,_b1] to RR >> >> > >> >> > Anyone got any ideas as to what this means? >> >> > >> >> > Regards, >> >> > >> >> > Will Furnass >> >> > >> >> > >> >> > ###### >> >> > >> >> > D, dx, k, L, n, P, Q, t, tau_a, tau_s, x= symbols('D dx k, L, n, P, >> >> > Q, >> >> > t, >> >> > \\tau_a, \\tau_s x') >> >> > >> >> > V = Q * 4 / (pi * (D**2)) >> >> > >> >> > expr1 = (tau_a - tau_s)**n * k >> >> > expr2 = (t * P * (tau_a - tau_s)**n * n) - (t * P * (tau_a - >> >> > tau_s)**n) >> >> > + (k >> >> > * tau_a) - (k * tau_s) >> >> > expr3 = 1 / (n-1) >> >> > expr4 = P * (tau_a - tau_s)**n >> >> > expr5 = (t * P * (tau_a - tau_s)**n * n) - (t * P * (tau_a - >> >> > tau_s)**n) >> >> > + (k >> >> > * tau_a) - (k * tau_s) >> >> > expr6 = ((expr1 / expr2) ** expr3) * expr4 / expr5 >> >> > expr7 = - P * exp(-t * P / k ) * (-tau_a + tau_s) / k >> >> > expr8 = Piecewise((expr6, Ne(n,1)), (expr7, Eq(n,1))) >> >> > >> >> > T = Piecewise((0, t < 0), (k * D * pi * dx * expr8 / Q, t >= 0)) >> >> > >> >> > curve_expr = integrate(T(t - (x/V)), (x, 0, L)) >> >> > >> >> > vals = { >> >> > P : 0.00022, >> >> > k : 0.5, >> >> > tau_s : 0.01, >> >> > tau_a : 2.0, >> >> > Q : 0.01, >> >> > D : 0.125, >> >> > L : 1225, >> >> > n : 3.0, >> >> > dx : 1} >> >> > >> >> > integrate(curve_expr.subs(vals)).eval() >> >> > >> >> > -- >> >> > You received this message because you are subscribed to the Google >> >> > Groups >> >> > "sympy" group. >> >> > To unsubscribe from this group and stop receiving emails from it, >> >> > send >> >> > an >> >> > email to [email protected]. >> >> > To post to this group, send email to [email protected]. >> >> > Visit this group at http://groups.google.com/group/sympy?hl=en-US. >> >> > For more options, visit https://groups.google.com/groups/opt_out. >> >> > >> >> > >> > >> > -- >> > You received this message because you are subscribed to the Google >> > Groups >> > "sympy" group. >> > To unsubscribe from this group and stop receiving emails from it, send >> > an >> > email to [email protected]. >> > To post to this group, send email to [email protected]. >> > Visit this group at http://groups.google.com/group/sympy?hl=en-US. >> > For more options, visit https://groups.google.com/groups/opt_out. >> > >> > > > -- > You received this message because you are subscribed to the Google Groups > "sympy" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to [email protected]. > To post to this group, send email to [email protected]. > Visit this group at http://groups.google.com/group/sympy?hl=en-US. > For more options, visit https://groups.google.com/groups/opt_out. > > -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To post to this group, send email to [email protected]. Visit this group at http://groups.google.com/group/sympy?hl=en-US. For more options, visit https://groups.google.com/groups/opt_out.
