ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Fri, 31 Aug 2018 21:38:39 -0500Conflicting Sage vs Wolfram evaluation of a limit?http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/<s> >Why are the following computed limits different (1 by Sage, 0 by Wolfram), and which (if either) is correct?
</s>
**EDIT**: Increasing the numerical precision in Wolfram produces a plot that *strongly suggests* that the limit is indeed $0$, which it had already computed. Presumably, Sage is computing the wrong limit simply because of inadequate numerical precision, so the question is now ...
>How can I increase the numerical precision in Sage, so that `limit()` and `plot()` will produce the correct results (i.e., the limit should be $0$ and the plot should show a stable approach to $0$)?
**Sage**: (you can cut/paste/execute this code [here](http://sagecell.sagemath.org/))
#in()=
f(x) = exp(-x^2/2)/sqrt(2*pi)
F(x) = (1 + erf(x/sqrt(2)))/2
num1(a,w) = (a+w)*f(a+w) - a*f(a)
num2(a,w) = f(a+w) - f(a)
den(a,w) = F(a+w) - F(a)
V(a,w) = 1 - num1(a,w)/den(a,w) - (num2(a,w)/den(a,w))^2
assume(w>0); print(limit(V(a,w), a=oo))
plot(V(a,1),a,0,8)
#out()=
1 #<--------- computed limit = 1
[![enter image description here][1]][1]
**Wolfram**: (you can execute this code [here](https://sandbox.open.wolframcloud.com/app/objects/0e2860d3-6c86-4d61-a9cf-e97fcf88c3b5#sidebar=compute))
#in()=
f[x_]:=Exp[-x^2/2]/Sqrt[2*Pi]
F[x_]:=(1 + Erf[x/Sqrt[2]])/2
num1[a_,w_] := (a+w)*f[a+w] - a*f[a]
num2[a_,w_] := f[a+w] - f[a]
den[a_,w_] := F[a+w] - F[a]
V[a_,w_] := 1 - num1[a,w]/den[a,w] - (num2[a,w]/den[a,w])^2
Assuming[w>0, Limit[V[a,w], a -> Infinity]]
Plot[V[a, 10], {a, 0, 100}, WorkingPrecision -> 128]
#out()=
0 (* <--------- computed limit = 0 *)
[![enter image description here][2]][2]
(This is supposed to compute the limit, as a -> oo, of the [variance of a standard normal distribution when truncated to the interval (a,a+w)](https://en.wikipedia.org/wiki/Truncated_normal_distribution#Moments).)
[1]: https://i.stack.imgur.com/eRAjI.png
[2]: https://i.stack.imgur.com/WTSf9.pngWed, 29 Aug 2018 09:17:24 -0500http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/Comment by John Palmieri for <p><s> >Why are the following computed limits different (1 by Sage, 0 by Wolfram), and which (if either) is correct?
</s></p>
<p><strong>EDIT</strong>: Increasing the numerical precision in Wolfram produces a plot that <em>strongly suggests</em> that the limit is indeed $0$, which it had already computed. Presumably, Sage is computing the wrong limit simply because of inadequate numerical precision, so the question is now ...</p>
<blockquote>
<p>How can I increase the numerical precision in Sage, so that <code>limit()</code> and <code>plot()</code> will produce the correct results (i.e., the limit should be $0$ and the plot should show a stable approach to $0$)?</p>
</blockquote>
<p><strong>Sage</strong>: (you can cut/paste/execute this code <a href="http://sagecell.sagemath.org/">here</a>) </p>
<pre><code>#in()=
f(x) = exp(-x^2/2)/sqrt(2*pi)
F(x) = (1 + erf(x/sqrt(2)))/2
num1(a,w) = (a+w)*f(a+w) - a*f(a)
num2(a,w) = f(a+w) - f(a)
den(a,w) = F(a+w) - F(a)
V(a,w) = 1 - num1(a,w)/den(a,w) - (num2(a,w)/den(a,w))^2
assume(w>0); print(limit(V(a,w), a=oo))
plot(V(a,1),a,0,8)
#out()=
1 #<--------- computed limit = 1
</code></pre>
<p><a href="https://i.stack.imgur.com/eRAjI.png"><img alt="enter image description here" src="https://i.stack.imgur.com/eRAjI.png"></a></p>
<p><strong>Wolfram</strong>: (you can execute this code <a href="https://sandbox.open.wolframcloud.com/app/objects/0e2860d3-6c86-4d61-a9cf-e97fcf88c3b5#sidebar=compute">here</a>)</p>
<pre><code>#in()=
f[x_]:=Exp[-x^2/2]/Sqrt[2*Pi]
F[x_]:=(1 + Erf[x/Sqrt[2]])/2
num1[a_,w_] := (a+w)*f[a+w] - a*f[a]
num2[a_,w_] := f[a+w] - f[a]
den[a_,w_] := F[a+w] - F[a]
V[a_,w_] := 1 - num1[a,w]/den[a,w] - (num2[a,w]/den[a,w])^2
Assuming[w>0, Limit[V[a,w], a -> Infinity]]
Plot[V[a, 10], {a, 0, 100}, WorkingPrecision -> 128]
#out()=
0 (* <--------- computed limit = 0 *)
</code></pre>
<p><a href="https://i.stack.imgur.com/WTSf9.png"><img alt="enter image description here" src="https://i.stack.imgur.com/WTSf9.png"></a></p>
<p>(This is supposed to compute the limit, as a -> oo, of the <a href="https://en.wikipedia.org/wiki/Truncated_normal_distribution#Moments">variance of a standard normal distribution when truncated to the interval (a,a+w)</a>.)</p>
http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/?comment=43541#post-id-43541I don't know how to get the plot to use extra precision. You can evaluate `V` with as high precision as you like: `V(20, 1).n(100)` will give 100 bits of precision, while `V(20,1).n(digits=100)` will give 100 digits of precision. For plotting, you could play with the `adaptive_tolerance` and `adaptive_recursion` arguments, although I haven't been able to use those to get anything good in this case. (It actually looks from the code like plotting uses `float(...)` on its arguments, without an option to specify precision. But I am not an expert in the plotting code in Sage.)Fri, 31 Aug 2018 13:10:01 -0500http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/?comment=43541#post-id-43541Answer by John Palmieri for <p><s> >Why are the following computed limits different (1 by Sage, 0 by Wolfram), and which (if either) is correct?
</s></p>
<p><strong>EDIT</strong>: Increasing the numerical precision in Wolfram produces a plot that <em>strongly suggests</em> that the limit is indeed $0$, which it had already computed. Presumably, Sage is computing the wrong limit simply because of inadequate numerical precision, so the question is now ...</p>
<blockquote>
<p>How can I increase the numerical precision in Sage, so that <code>limit()</code> and <code>plot()</code> will produce the correct results (i.e., the limit should be $0$ and the plot should show a stable approach to $0$)?</p>
</blockquote>
<p><strong>Sage</strong>: (you can cut/paste/execute this code <a href="http://sagecell.sagemath.org/">here</a>) </p>
<pre><code>#in()=
f(x) = exp(-x^2/2)/sqrt(2*pi)
F(x) = (1 + erf(x/sqrt(2)))/2
num1(a,w) = (a+w)*f(a+w) - a*f(a)
num2(a,w) = f(a+w) - f(a)
den(a,w) = F(a+w) - F(a)
V(a,w) = 1 - num1(a,w)/den(a,w) - (num2(a,w)/den(a,w))^2
assume(w>0); print(limit(V(a,w), a=oo))
plot(V(a,1),a,0,8)
#out()=
1 #<--------- computed limit = 1
</code></pre>
<p><a href="https://i.stack.imgur.com/eRAjI.png"><img alt="enter image description here" src="https://i.stack.imgur.com/eRAjI.png"></a></p>
<p><strong>Wolfram</strong>: (you can execute this code <a href="https://sandbox.open.wolframcloud.com/app/objects/0e2860d3-6c86-4d61-a9cf-e97fcf88c3b5#sidebar=compute">here</a>)</p>
<pre><code>#in()=
f[x_]:=Exp[-x^2/2]/Sqrt[2*Pi]
F[x_]:=(1 + Erf[x/Sqrt[2]])/2
num1[a_,w_] := (a+w)*f[a+w] - a*f[a]
num2[a_,w_] := f[a+w] - f[a]
den[a_,w_] := F[a+w] - F[a]
V[a_,w_] := 1 - num1[a,w]/den[a,w] - (num2[a,w]/den[a,w])^2
Assuming[w>0, Limit[V[a,w], a -> Infinity]]
Plot[V[a, 10], {a, 0, 100}, WorkingPrecision -> 128]
#out()=
0 (* <--------- computed limit = 0 *)
</code></pre>
<p><a href="https://i.stack.imgur.com/WTSf9.png"><img alt="enter image description here" src="https://i.stack.imgur.com/WTSf9.png"></a></p>
<p>(This is supposed to compute the limit, as a -> oo, of the <a href="https://en.wikipedia.org/wiki/Truncated_normal_distribution#Moments">variance of a standard normal distribution when truncated to the interval (a,a+w)</a>.)</p>
http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/?answer=43542#post-id-43542I asked about plotting with high precision on the Google group [sage-support](https://groups.google.com/d/msg/sage-support/mz7DaVyhC0c/McWpIHU8AwAJ) and got the a few suggestions. One that works for me is, rather than plot the function, instead generate pairs `(a, V(a,1))` and draw the plot connecting those points. Combining this with an appropriate call to `numerical_approx` (also known as just `n`) works well:
L = [(x, V(x,1).n(1000)) for x in sxrange(0, 30, 0.01)]
plot(line(L))
generates
![image description](/upfiles/15357465267144159.png)Fri, 31 Aug 2018 15:16:39 -0500http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/?answer=43542#post-id-43542Comment by res0001 for <p>I asked about plotting with high precision on the Google group <a href="https://groups.google.com/d/msg/sage-support/mz7DaVyhC0c/McWpIHU8AwAJ">sage-support</a> and got the a few suggestions. One that works for me is, rather than plot the function, instead generate pairs <code>(a, V(a,1))</code> and draw the plot connecting those points. Combining this with an appropriate call to <code>numerical_approx</code> (also known as just <code>n</code>) works well:</p>
<pre><code>L = [(x, V(x,1).n(1000)) for x in sxrange(0, 30, 0.01)]
plot(line(L))
</code></pre>
<p>generates</p>
<p><img alt="image description" src="/upfiles/15357465267144159.png"></p>
http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/?comment=43546#post-id-43546Thanks -- that is good to know! Now I just need to know how to get Sage to compute the correct limit (i.e. `Limit[V[a,w], a -> Infinity]` should output $0$ instead of its current output of $1$).Fri, 31 Aug 2018 21:38:39 -0500http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/?comment=43546#post-id-43546Answer by kcrisman for <p><s> >Why are the following computed limits different (1 by Sage, 0 by Wolfram), and which (if either) is correct?
</s></p>
<p><strong>EDIT</strong>: Increasing the numerical precision in Wolfram produces a plot that <em>strongly suggests</em> that the limit is indeed $0$, which it had already computed. Presumably, Sage is computing the wrong limit simply because of inadequate numerical precision, so the question is now ...</p>
<blockquote>
<p>How can I increase the numerical precision in Sage, so that <code>limit()</code> and <code>plot()</code> will produce the correct results (i.e., the limit should be $0$ and the plot should show a stable approach to $0$)?</p>
</blockquote>
<p><strong>Sage</strong>: (you can cut/paste/execute this code <a href="http://sagecell.sagemath.org/">here</a>) </p>
<pre><code>#in()=
f(x) = exp(-x^2/2)/sqrt(2*pi)
F(x) = (1 + erf(x/sqrt(2)))/2
num1(a,w) = (a+w)*f(a+w) - a*f(a)
num2(a,w) = f(a+w) - f(a)
den(a,w) = F(a+w) - F(a)
V(a,w) = 1 - num1(a,w)/den(a,w) - (num2(a,w)/den(a,w))^2
assume(w>0); print(limit(V(a,w), a=oo))
plot(V(a,1),a,0,8)
#out()=
1 #<--------- computed limit = 1
</code></pre>
<p><a href="https://i.stack.imgur.com/eRAjI.png"><img alt="enter image description here" src="https://i.stack.imgur.com/eRAjI.png"></a></p>
<p><strong>Wolfram</strong>: (you can execute this code <a href="https://sandbox.open.wolframcloud.com/app/objects/0e2860d3-6c86-4d61-a9cf-e97fcf88c3b5#sidebar=compute">here</a>)</p>
<pre><code>#in()=
f[x_]:=Exp[-x^2/2]/Sqrt[2*Pi]
F[x_]:=(1 + Erf[x/Sqrt[2]])/2
num1[a_,w_] := (a+w)*f[a+w] - a*f[a]
num2[a_,w_] := f[a+w] - f[a]
den[a_,w_] := F[a+w] - F[a]
V[a_,w_] := 1 - num1[a,w]/den[a,w] - (num2[a,w]/den[a,w])^2
Assuming[w>0, Limit[V[a,w], a -> Infinity]]
Plot[V[a, 10], {a, 0, 100}, WorkingPrecision -> 128]
#out()=
0 (* <--------- computed limit = 0 *)
</code></pre>
<p><a href="https://i.stack.imgur.com/WTSf9.png"><img alt="enter image description here" src="https://i.stack.imgur.com/WTSf9.png"></a></p>
<p>(This is supposed to compute the limit, as a -> oo, of the <a href="https://en.wikipedia.org/wiki/Truncated_normal_distribution#Moments">variance of a standard normal distribution when truncated to the interval (a,a+w)</a>.)</p>
http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/?answer=43532#post-id-43532There are two separate problems here.
- Sage appears to calculate the wrong symbolic limit.
- Plotting expressions of this kind has inherent numerical instability, which Sage could make easier to control.
You seem to have almost purposely chosen a formula which would have high numerical instability! Any time you have denominators and numerators both going to zero and you start plugging in for tiny changes $\Delta x$ you have the chance at running into trouble.
Unfortunately, `limit(V(a,w), a=+oo, algorithm='sympy')` doesn't seem to produce a result (perhaps Sympy doesn't handle these sorts of limits). If there is independent confirmation of this error, a ticket should be opened.
However, the issue is symbolics, not numerical precision. You may find some luck with [adapting some plot refinement parameters](http://doc.sagemath.org/html/en/reference/plotting/sage/plot/plot.html#sage.plot.plot.adaptive_refinement), or possibly inputting numbers of higher precision - something like `M = [(x,V(x,1)) for x in [RealField(1000)(y) for y in [0,0.001..8]]]; points(M)` but I have to admit I was unable to get that to avoid the instability either.Thu, 30 Aug 2018 08:53:40 -0500http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/?answer=43532#post-id-43532Answer by tmonteil for <p><s> >Why are the following computed limits different (1 by Sage, 0 by Wolfram), and which (if either) is correct?
</s></p>
<p><strong>EDIT</strong>: Increasing the numerical precision in Wolfram produces a plot that <em>strongly suggests</em> that the limit is indeed $0$, which it had already computed. Presumably, Sage is computing the wrong limit simply because of inadequate numerical precision, so the question is now ...</p>
<blockquote>
<p>How can I increase the numerical precision in Sage, so that <code>limit()</code> and <code>plot()</code> will produce the correct results (i.e., the limit should be $0$ and the plot should show a stable approach to $0$)?</p>
</blockquote>
<p><strong>Sage</strong>: (you can cut/paste/execute this code <a href="http://sagecell.sagemath.org/">here</a>) </p>
<pre><code>#in()=
f(x) = exp(-x^2/2)/sqrt(2*pi)
F(x) = (1 + erf(x/sqrt(2)))/2
num1(a,w) = (a+w)*f(a+w) - a*f(a)
num2(a,w) = f(a+w) - f(a)
den(a,w) = F(a+w) - F(a)
V(a,w) = 1 - num1(a,w)/den(a,w) - (num2(a,w)/den(a,w))^2
assume(w>0); print(limit(V(a,w), a=oo))
plot(V(a,1),a,0,8)
#out()=
1 #<--------- computed limit = 1
</code></pre>
<p><a href="https://i.stack.imgur.com/eRAjI.png"><img alt="enter image description here" src="https://i.stack.imgur.com/eRAjI.png"></a></p>
<p><strong>Wolfram</strong>: (you can execute this code <a href="https://sandbox.open.wolframcloud.com/app/objects/0e2860d3-6c86-4d61-a9cf-e97fcf88c3b5#sidebar=compute">here</a>)</p>
<pre><code>#in()=
f[x_]:=Exp[-x^2/2]/Sqrt[2*Pi]
F[x_]:=(1 + Erf[x/Sqrt[2]])/2
num1[a_,w_] := (a+w)*f[a+w] - a*f[a]
num2[a_,w_] := f[a+w] - f[a]
den[a_,w_] := F[a+w] - F[a]
V[a_,w_] := 1 - num1[a,w]/den[a,w] - (num2[a,w]/den[a,w])^2
Assuming[w>0, Limit[V[a,w], a -> Infinity]]
Plot[V[a, 10], {a, 0, 100}, WorkingPrecision -> 128]
#out()=
0 (* <--------- computed limit = 0 *)
</code></pre>
<p><a href="https://i.stack.imgur.com/WTSf9.png"><img alt="enter image description here" src="https://i.stack.imgur.com/WTSf9.png"></a></p>
<p>(This is supposed to compute the limit, as a -> oo, of the <a href="https://en.wikipedia.org/wiki/Truncated_normal_distribution#Moments">variance of a standard normal distribution when truncated to the interval (a,a+w)</a>.)</p>
http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/?answer=43518#post-id-43518I do not see any difference, except that wolfram cut part of the graph where your function takes higer values.
You can change the range of y-values to be between -0.05 and 0.15: just replace the last line of the Sage code by:
plot(V(a,1),a,0,8,ymax=0.15,ymin=-0.05)
You will get:
![this image](/upfiles/15355596853838028.png)Wed, 29 Aug 2018 11:22:00 -0500http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/?answer=43518#post-id-43518Comment by res0001 for <p>I do not see any difference, except that wolfram cut part of the graph where your function takes higer values.</p>
<p>You can change the range of y-values to be between -0.05 and 0.15: just replace the last line of the Sage code by:</p>
<pre><code>plot(V(a,1),a,0,8,ymax=0.15,ymin=-0.05)
</code></pre>
<p>You will get:</p>
<p><img alt="this image" src="/upfiles/15355596853838028.png"></p>
http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/?comment=43519#post-id-43519The point is that the limit is computed to be 1 by Sage and 0 by Wolfram, in spite of the unstable behavior indicated by the plots. (I've edited to clarify.)Wed, 29 Aug 2018 12:10:33 -0500http://ask.sagemath.org/question/43517/conflicting-sage-vs-wolfram-evaluation-of-a-limit/?comment=43519#post-id-43519