[scilab-Users] intg: results differ substantially from those from Wolfram Alpha, which are correct?
michael.baudin at scilab.org
michael.baudin at scilab.org
Wed Jun 1 09:07:35 CEST 2011
Hi,
Serge is absolutely right: any other splitting of the integral
leads to the result.
In my experiments, Wolfram Alpha does not return
-0.111. Notice that Wolfram Alpha can also be wrong. The result depends
on the way to ask the question :
*
http://www.wolframalpha.com/input/?i=integrate+%281%2Fsqrt%282*pi%29%29*log%28abs%281%2B0.1*x%29%29*exp%28-%28x^2%29%2F2%29+from+-100+to+-100
[1] [2]
returns 0.
*
http://www.wolframalpha.com/input/?i=integrate&a=*C.integrate-_*Calculator.dflt-&f2=%281%2Fsqrt%282*pi%29%29*log%28abs%281%2B0.1*x%29%29*exp%28-%28x^2%29%2F2%29&f=Integral.integrand_%281%2Fsqrt%282*pi%29%29*log%28abs%281%2B0.1*x%29%29*exp%28-%28x^2%29%2F2%29&f3=-100&f=Integral.rangestart_-100&f4=100&f=Integral.rangeend_100&a=*FVarOpt.1-_**-.***Integral.variable---.**Integral.rangestart-.*Integral.rangeend
[3]--- [4]
returns -0.00507764.
Bruno Pincon warned me that the
integrand is singular at x=-a/sigma. By chance, this singularity is not
seen by the integrator, which does not sample at this node. Moreover,
the neighbors of this singularity have a y-value which is too small to
be influential in the result.
-->a = 1;
-->sigma = 0.1;
-->s =
-a/sigma;
-->ieee(2);
-->y=f(s, a,sigma)
y =
- Inf
-->xnear =
nearfloat("succ", s);
-->y=f(xnear, a,sigma)
y =
- 2.827D-21
-->xnear = nearfloat("pred", s);
-->y=f(xnear, a,sigma)
y =
-
2.773D-21
Splitting along this discontinuity returns the correct
result :
-->out1=intg(-100,s,list(f,1,.1)); // Wolfram =
-3.95999x10-23
-->out2=intg(s,100,list(f,1,.1)); // Wolfram =
-0.00507764
-->I = out1 + out2 // Wolfram = -0.00507764
I =
-
0.0050776
But splitting along another point also returns the correct
result.
-->out1=intg(-100,10,list(f,1,.1));
-->out2=intg(10,100,list(f,1,.1));
-->I = out1 + out2
I =
-
0.0050776
I guess that splitting along any point for which the
function is non-zero allows to retrieve the result...
We must be
warned that it is easy to fool any deterministic integration algorithm.
Indeed, the integration algorithm behind the intg function is based on
an adaptive algorithm, based on deterministic integration nodes. The
algorithm is there :
http://www.netlib.org/quadpack/dqagse.f [5]
The
first step is to use a 21 point rule :
http://www.netlib.org/quadpack/dqk21.f [6]
If this rule identifies a
zero integral, the algorithm returns. Hence, there is an easy way to
fool this algorithm, as suggested by Kahan in a 1980 paper ("Handheld
calculator evaluates integrals"). In the first step, we call the
integrator with a function which returns y=f(x)=0. The integrator
measures zero function values, measures a zero integral and returns.
Meanwhile, we have recorded the integration nodes. Then, we integrate a
special function, which is equal to zero at the integration nodes and is
non-zero, positive, elsewhere. The integrator samples at the same
points, and, as before, returns a zero integral: but this time it is
wrong. This idea is explored in the following script.
a = -10;
b =
10;
function y=foolme(x)
global _storex_
y=0
mprintf("x=%fn",x)
_storex_($+1)=x
endfunction
global _storex_
_storex_ = [];
intg(a,b,foolme)
function y=reallyfoolme(x)
global _storex_
r =
_storex_
y = prod(x-r)^2
endfunction
intg(a,b,reallyfoolme)
x =
linspace(a,b,1000);
for i = 1 : 1000
y(i)=reallyfoolme(x(i));
end
plot(x,y,"b-")
Best regards,
Michael Baudin
Le 24/05/2011 16:51,
Serge Steer a écrit : Le 24/05/2011 09:34, Ginters Bušs a écrit : -----
Message d'origine -----
De : Ginters Bušs
Date : 20/05/2011 13:09:
Dear all,
Let's integrate:
function y=f(x, a,
sigma),y=(1/sqrt(2*%pi))*log(abs(a+sigma*x))*exp(-(x^2)/2),endfunction
out=intg(-1e+2,1e+2,list(f,1,.1))
out=8.605D-49
but Wolfram Alpha
gives out= -0.111
which is a totally different answer.
I've noticed
that intg and integrate incline to give values close to zero when
boundaries tend to infinity. So, I trust Wolfram Alpha more. How to get
around the apparent mistakes in intg, integrate (particularly, I'm
interested in indefinite integrals)?
Gin.
On Sat, May 21, 2011 at
1:20 PM, Samuel GOUGEON wrote:
Hello,
The following thread on
Bugzilla may also feed the present discussion:
http://bugzilla.scilab.org/show_bug.cgi?id=5728 [8]
Regards
Samuel
Thanks. I like your insistance in that discussion, Sam.
Although a
"known-limitation-of-current-algorithm" might be a correct description
of the current algorithm, it nevertheless does not help much when one
wants the algorithm to give the correct solution and observes that a guy
like Wolfram has apparently has done something about this issue: for y
given function, Scilab's intg/integrate gives ok value for bounds
(-50,50), but screws up for (-100,100) and larger bounds, which I
consider still pretty narrow. On the contrary, Wolfram Alpha
http://www.wolframalpha.com/input/?i=integration&a=*C.integration-_*Calculator.dflt-&f2=%281%2Fsqrt%282*pi%29%29*log%28abs%281%2B.1*x%29%29*exp%28-%28x
[9]^2%29%2F2%29&x=0&y=0&f=Integral.integrand_%281%2Fsqrt%282*pi%29%29*log%28abs%281%2B.1*x%29%29*exp%28-%28x^2%29%2F2%29&f3=-infinity&f=Integral.rangestart_-infinity&f4=%2Binfinity&f=Integral.rangeend_%2Binfinity&a=*FVarOpt.1-_**-.***Integral.variable---.**Integral.rangestart-.*Integral.rangeend---
gives ok values for bounds (-50,50), (-100,100), (-1000,1000) and it
allows to directly use infinite bounds and get ok result (interim bounds
like (-1e9,1e9) are not ok). I consider this a superior solution. If
Wolfram can, why shouldn't Scilab try? The status quo might hurt
Scilab's reputation in the domain in long run.
Gin.
WolframAlpha is
based on a symbolic computation engine that furnishes much more tools
for integration. In particular it is possible for it doing serie
expansion near the infinity.
If one remarks that for x>35 and
xintg(-35,35,myfun)
ans =
- 0.0050776
but without symbolic
computation facilities it is difficult to automatically determine that
such integration interval reduction can be done...
Serge Steer
INRIA
--
Michaël Baudin
Ingénieur de
développement
michael.baudin at scilab.org
[10]
-------------------------
Consortium Scilab - Digiteo
Domaine de
Voluceau - Rocquencourt
B.P. 105 - 78153 Le Chesnay Cedex
Tel. : 01 39
63 56 87 - Fax : 01 39 63 55 94
Links:
------
[1]
http://www.wolframalpha.com/input/?i=integrate+%281%2Fsqrt%282*pi%29%29*log%28abs%281%2B0.1*x%29%29*exp%28-%28x%5E2%29%2F2%29+from+-100+to+-100
[2]
http://www.wolframalpha.com/input/?i=integrate+%281%2Fsqrt%282*pi%29%29*log%28abs%281%2B0.1*x%29%29*exp%28-%28x%5E2%29%2F2%29+from+-100+to+-100
[3]
http://www.wolframalpha.com/input/?i=integrate&a=*C.integrate-_*Calculator.dflt-&f2=%281%2Fsqrt%282*pi%29%29*log%28abs%281%2B0.1*x%29%29*exp%28-%28x%5E2%29%2F2%29&f=Integral.integrand_%281%2Fsqrt%282*pi%29%29*log%28abs%281%2B0.1*x%29%29*exp%28-%28x%5E2%29%2F2%29&f3=-100&f=Integral.rangestart_-100&f4=100&f=Integral.rangeend_100&a=*FVarOpt.1-_**-.***Integral.variable---.**Integral.rangestart-.*Integral.rangeend
[4]
http://www.wolframalpha.com/input/?i=integrate&a=*C.integrate-_*Calculator.dflt-&f2=%281%2Fsqrt%282*pi%29%29*log%28abs%281%2B0.1*x%29%29*exp%28-%28x%5E2%29%2F2%29&f=Integral.integrand_%281%2Fsqrt%282*pi%29%29*log%28abs%281%2B0.1*x%29%29*exp%28-%28x%5E2%29%2F2%29&f3=-100&f=Integral.rangestart_-100&f4=100&f=Integral.rangeend_100&a=*FVarOpt.1-_**-.***Integral.variable---.**Integral.rangestart-.*Integral.rangeend---
[5]
http://www.netlib.org/quadpack/dqagse.f
[6]
http://www.netlib.org/quadpack/dqk21.f
[7]
mailto:Samuel.Gougeon at univ-lemans.fr
[8]
http://bugzilla.scilab.org/show_bug.cgi?id=5728
[9]
http://www.wolframalpha.com/input/?i=integration&a=*C.integration-_*Calculator.dflt-&f2=%281%2Fsqrt%282*pi%29%29*log%28abs%281%2B.1*x%29%29*exp%28-%28x
[10]
mailto:michael.baudin at scilab.org
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20110601/ce4152a8/attachment.htm>
More information about the users
mailing list