Ask Your Question
1

Problem with numerical integral

asked 2018-04-10 10:13:28 +0100

user111 gravatar image

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 flag offensive close merge delete

2 Answers

Sort by » oldest newest most voted
4

answered 2018-04-10 12:49:57 +0100

dan_fulea gravatar image

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)
edit flag offensive delete link more
4

answered 2018-04-10 10:59:39 +0100

Emmanuel Charpentier gravatar image

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

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2018-04-10 10:13:28 +0100

Seen: 1,165 times

Last updated: Apr 10 '18