<p>Hi,<br /> <br /> Serge is absolutely right: any other splitting of the integral leads to the result.<br /> <br /> 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 :<br /> <br /> * <a class="moz-txt-link-freetext" href="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">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</a> <a class="moz-txt-link-rfc2396E" href="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"></a> <br /> <br /> returns 0. <br /> <br /> * <a class="moz-txt-link-freetext" href="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">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</a>--- <a class="moz-txt-link-rfc2396E" href="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---"></a> <br /> <br /> returns -0.00507764. <br /> <br /> 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.<br /> <br /> -->a = 1;<br /> -->sigma = 0.1;<br /> -->s = -a/sigma;<br /> -->ieee(2);<br /> -->y=f(s, a,sigma)<br /> y =<br /> - Inf <br /> -->xnear = nearfloat("succ", s);<br /> -->y=f(xnear, a,sigma)<br /> y =<br /> - 2.827D-21 <br /> -->xnear = nearfloat("pred", s);<br /> -->y=f(xnear, a,sigma)<br /> y =<br /> - 2.773D-21 <br /> <br /> Splitting along this discontinuity returns the correct result :<br /> <br /> -->out1=intg(-100,s,list(f,1,.1)); // Wolfram = -3.95999x10-23<br /> -->out2=intg(s,100,list(f,1,.1)); // Wolfram = -0.00507764<br /> -->I = out1 + out2 // Wolfram = -0.00507764<br /> I =<br /> - 0.0050776 <br /> <br /> But splitting along another point also returns the correct result. <br /> <br /> -->out1=intg(-100,10,list(f,1,.1)); <br /> -->out2=intg(10,100,list(f,1,.1)); <br /> -->I = out1 + out2 <br /> I =<br /> - 0.0050776 <br /> <br /> I guess that splitting along any point for which the function is non-zero allows to retrieve the result...<br /> <br /> 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 :<br /> <br /> <a href="http://www.netlib.org/quadpack/dqagse.f">http://www.netlib.org/quadpack/dqagse.f</a><br /> <br /> The first step is to use a 21 point rule :<br /> <br /> <a href="http://www.netlib.org/quadpack/dqk21.f">http://www.netlib.org/quadpack/dqk21.f</a><br /> <br /> 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.<br /> <br /> a = -10;<br /> b = 10;<br /> function y=foolme(x)<br /> global _storex_<br /> y=0<br /> mprintf("x=%f\n",x)<br /> _storex_($+1)=x<br /> endfunction<br /> <br /> global _storex_<br /> _storex_ = [];<br /> intg(a,b,foolme)<br /> <br /> function y=reallyfoolme(x)<br /> global _storex_<br /> r = _storex_<br /> y = prod(x-r)^2<br /> endfunction<br /> intg(a,b,reallyfoolme)<br /> <br /> x = linspace(a,b,1000);<br /> for i = 1 : 1000<br /> y(i)=reallyfoolme(x(i));<br /> end<br /> plot(x,y,"b-")<br /> <br /> <br /> <br /> <br /> Best regards,<br /> <br /> Michael Baudin<br /> <br /> Le 24/05/2011 16:51, Serge Steer a écrit :</p>
<blockquote cite="mid:4DDBC5FA.5000605@inria.fr">Le 24/05/2011 09:34, Ginters Bušs a écrit :
<blockquote cite="mid:BANLkTincfJTS0mzzmfZ+gHJpqpZy8KkOrQ@mail.gmail.com">----- Message d'origine -----<br /> De : Ginters Bušs<br /> Date : 20/05/2011 13:09:<br />
<blockquote class="gmail_quote" style="border-left: 1px solid #cccccc; margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">Dear all,<br /> <br /> Let's integrate:<br /> <br /> function y=f(x, a, sigma),y=(1/sqrt(2*%pi))*log(abs(a+sigma*x))*exp(-(x^2)/2),endfunction<br /> <br /> out=intg(-1e+2,1e+2,list(f,1,.1))<br /> <br /> out=8.605D-49<br /> <br /> but Wolfram Alpha gives out= -0.111<br /> <br /> which is a totally different answer.<br /> <br /> 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)?<br /> <br /> Gin.<br /> <br /> <br /></blockquote>
<br /> <br />
<div class="gmail_quote">On Sat, May 21, 2011 at 1:20 PM, Samuel GOUGEON <span><<a href="mailto:Samuel.Gougeon@univ-lemans.fr">Samuel.Gougeon@univ-lemans.fr</a>></span> wrote:<br />
<blockquote class="gmail_quote" style="border-left: 1px solid #cccccc; margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;"> Hello,<br /> The following thread on Bugzilla may also feed the present discussion:<br /> <a href="http://bugzilla.scilab.org/show_bug.cgi?id=5728" target="_blank">http://bugzilla.scilab.org/show_bug.cgi?id=5728</a><br /> Regards<br /> Samuel<br /> <br />
<div>
<div class="h5"></div>
</div>
</blockquote>
</div>
Thanks. I like your insistance in that discussion, Sam. <br /> <br /> 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 <br /> <br /> <a href="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">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</a>^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---<br /> <br /> 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.<br /> <br /> Gin.</blockquote>
<br /> 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.<br /> <br /> If one remarks that for x>35 and x<-35 the function evaluates to zero it it possible to reduce the integration interval and obtain a good result<br /> -->intg(-35,35,myfun)<br /> ans =<br /> <br /> - 0.0050776 <br /> <br /> but without symbolic computation facilities it is difficult to automatically determine that such integration interval reduction can be done...<br /> <br /> Serge Steer<br /> INRIA<br /> <br /> <br /> <br /> <br /></blockquote>
<p><br /> <br /></p>
<pre class="moz-signature">--
Michaël Baudin
Ingénieur de développement
<a class="moz-txt-link-abbreviated" href="mailto:michael.baudin@scilab.org">michael.baudin@scilab.org</a>
-------------------------
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
</pre>