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)
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>()
----> 1 maxima.quad_qagp((Integer(1)-x)**(-Integer(2)/Integer(3)),x,Integer(0),Integer(1),[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/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)
1004 self._batch(cmd, batchload=True)
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:
sage40 : quad_qagp(sage34,sage35,sage36,sage37,sage39)$
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
sage: mpmath.quad(lambda x:(1-x)^(-2/3),[0,1])
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...
Asked: 2018-04-10 10:13:28 +0100
Seen: 1,200 times
Last updated: Apr 10 '18
Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.