# Conflicting Sage vs Wolfram evaluation of a limit?

>Why are the following computed limits different (1 by Sage, 0 by Wolfram), and which (if either) is correct?

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)

#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


Wolfram: (you can execute this code here)

#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 *)


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

edit retag close merge delete

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

( 2018-08-31 13:10:01 -0600 )edit

Sort by » oldest newest most voted

I asked about plotting with high precision on the Google group sage-support 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

more

Thanks -- 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$).

( 2018-08-31 21:38:39 -0600 )edit

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

more

I 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:

more

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

( 2018-08-29 12:10:33 -0600 )edit