# Problem with numerical integral

Is it normal that the second integral can't be computed ?

show(numerical_integral((x)^(-2/3),0,1)); show(numerical_integral((1-x)^(-2/3),0,1));

(2.99999987765,2.57881950224×10−06)
(+infinity,NaN)


Thanks

edit retag close merge delete

Sort by » oldest newest most voted

The following works for me:

sage: numerical_integral( lambda x: x^(-2/3), (0,1) )
(2.99999987764845, 2.578819502242948e-06)
sage: numerical_integral( lambda x: (1-x)^(-2/3), (0,1) )
0.0 cannot be raised to a negative power
0.0 cannot be raised to a negative power
0.0 cannot be raised to a negative power
0.0 cannot be raised to a negative power
0.0 cannot be raised to a negative power
0.0 cannot be raised to a negative power
(2.999988120409371, 3.768990092235661e-05)


Here, numerical_integral gets a (true pythonic) function, a "callable object" as its first argument. There are seldom problems with such an input. Expressions instead of callables are sometimes a problem for the numerical_integral.

Note: Related numerical integral, that goes to maxima and uses expressions:

sage: (x^(-2/3)).nintegral(x,0,1)
(3.000000000000006, 3.01980662698043e-14, 231, 0)
sage: ((1-x)^(-2/3)).nintegral(x,0,1)
(3.000000000002707, 8.78142003557514e-12, 231, 0)

more

Well... Let's do it symbolically :

sage: ((1-x)^(-2/3)).integrate(x)
-3*(-x + 1)^(1/3)


This is indeed a correct primitive :

sage: ((1-x)^(-2/3)).integrate(x).diff(x)
(-x + 1)^(-2/3)


Therefore, on might compute :

sage: ((1-x)^(-2/3)).integrate(x).limit(x=0)
-3
sage: ((1-x)^(-2/3)).integrate(x).limit(x=1)
0


And, by the way :

sage: ((1-x)^(-2/3)).integrate(x,0,1)
3


But, indeed :

sage: numerical_integral((1-x)^(-2/3),0,1)
(inf, nan)


I can but conclude that this is a bug of numerical_integral...

Note that, misery loving company, Sage is not alone in having problems with this integration :

sage: maxima.quad_qagp((1-x)^(-2/3),x,0,1,[1])
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-51-039b3b49b5dd> in <module>()

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __call__(self, *args, **kwds)
627
628     def __call__(self, *args, **kwds):
--> 629         return self._parent.function_call(self._name, list(args), kwds)
630
631     def _instancedoc_(self):

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in function_call(self, function, args, kwds)
574                                        [s.name() for s in args],
575                                        ['%s=%s'%(key,value.name()) for key, value in kwds.items()])
--> 576         return self.new(s)
577
578     def _function_call_string(self, function, args, kwds):

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in new(self, code)
343
344     def new(self, code):
--> 345         return self(code)
346
347     ###################################################################

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __call__(self, x, name)
278
279         if isinstance(x, string_types):
--> 280             return cls(self, x, name=name)
281         try:
282             return self._coerce_from_special_method(x)

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/maxima.pyc in __init__(self, parent, value, is_name, name)
1159             True
1160         """
-> 1161         ExpectElement.__init__(self, parent, value, is_name=False, name=None)
1162
1163     def display2d(self, onscreen=True):

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/expect.pyc in __init__(self, parent, value, is_name, name)
1375         else:
1376             try:
-> 1377                 self._name = parent._create(value, name=name)
1378             # Convert ValueError and RuntimeError to TypeError for
1379             # coercion to work properly.

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in _create(self, value, name)
474     def _create(self, value, name=None):
475         name = self._next_var_name() if name is None else name
--> 476         self.set(name, value)
477         return name
478

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/maxima.pyc in set(self, var, value)
1005         else:
-> 1006             self._eval_line(cmd)
1007             #self._sendline(cmd)
1008             #self._expect_expr()

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/maxima.pyc in _eval_line(self, line, allow_use_file, wait_for_prompt, reformat, error_check, restart_if_needed)
796         out = self._before()        # input echo + output prompt + output
797         if error_check:
--> 798             self._error_check(line, out)
799         if not reformat:
800             return out

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/maxima.pyc in _error_check(self, cmd, out)
936         m = r.search(out)
937         if not m is None:
--> 938             self._error_msg(cmd, out)
939
940     def _error_msg(self, cmd, out):

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/maxima.pyc in _error_msg(self, cmd, out)
953                 Principal Value
954         """
--> 955         raise TypeError("Error executing code in Maxima\nCODE:\n\t%s\nMaxima ERROR:\n\t%s"%(cmd, out.replace('-- an error.  To debug this try debugmode(true);','')))
956
957     ###########################################

TypeError: Error executing code in Maxima
CODE:
Maxima ERROR:

quad_qagp: Cannot numerically evaluate 1/(1-_SAGE_VAR_x)^(2/3) at 1.0
-- an error. To debug this try: debugmode(true);

sage: mathematica.NIntegrate((1-x)^(-2/2),[x,0,1])
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-52-af13341319ac> in <module>()
----> 1 mathematica.NIntegrate((Integer(1)-x)**(-Integer(2)/Integer(2)),[x,Integer(0),Integer(1)])

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __call__(self, *args, **kwds)
627
628     def __call__(self, *args, **kwds):
--> 629         return self._parent.function_call(self._name, list(args), kwds)
630
631     def _instancedoc_(self):

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in function_call(self, function, args, kwds)
574                                        [s.name() for s in args],
575                                        ['%s=%s'%(key,value.name()) for key, value in kwds.items()])
--> 576         return self.new(s)
577
578     def _function_call_string(self, function, args, kwds):

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in new(self, code)
343
344     def new(self, code):
--> 345         return self(code)
346
347     ###################################################################

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __call__(self, x, name)
278
279         if isinstance(x, string_types):
--> 280             return cls(self, x, name=name)
281         try:
282             return self._coerce_from_special_method(x)

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/expect.pyc in __init__(self, parent, value, is_name, name)
1375         else:
1376             try:
-> 1377                 self._name = parent._create(value, name=name)
1378             # Convert ValueError and RuntimeError to TypeError for
1379             # coercion to work properly.

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in _create(self, value, name)
474     def _create(self, value, name=None):
475         name = self._next_var_name() if name is None else name
--> 476         self.set(name, value)
477         return name
478

/usr/local/sage-8/local/lib/python2.7/site-packages/sage/interfaces/mathematica.pyc in set(self, var, value)
533         out = self._eval_line(cmd, allow_use_file=True)
534         if len(out) > 8:
--> 535             raise TypeError("Error executing code in Mathematica\nCODE:\n\t%s\nMathematica ERROR:\n\t%s"%(cmd, out))
536
537     def get(self, var, ascii_art=False):

TypeError: Error executing code in Mathematica
CODE:
sage11=NIntegrate[sage7,sage1];
Mathematica ERROR:
NIntegrate::slwcon:
Numerical integration converging too slowly; suspect one of the following:
singularity, value of the integration is 0, highly oscillatory integrand,
or WorkingPrecision too small.

NIntegrate::ncvb:
NIntegrate failed to converge to prescribed accuracy after 9
recursive bisections in x near {x} =
{0.99999999999999999999999999998369273566250875422125492455831169218}
. NIntegrate obtained 149.573 and 23.2473
for the integral and error estimates.


But :

sage: import mpmath
mpf('2.9999987601927596')


It would be interesting to test other numerical integratin libraries available in Sage (and, more generally, Python...). And to report this bug...

more