# An apparently “simple” limit?

Let $c$ and $h$ be real values. I was using Mathematica to compute the limit $(h \rightarrow \infty)$ of the following expression:

$$\frac{\left(h^2 +c^2 h^2 + \sqrt{4 h^2+\left(h^2+c^2 h^2\right)^2}\right)^2}{4\left(h^2 + \frac{1}{4}\left(h^2 +c^2 h^2 + \sqrt{4 h^2+\left(h^2+c^2 h^2\right)^2}\right)^2\right)}\text{//FullSimplify}$$ which is equal to $$=\frac{2}{4+(1+c^2)^2 h^2 - (1+c^2)\sqrt{4h^2+(1+c^2)^2 h^4}}$$

Just as the rough estimate anyone can make by counting powers of $h$ in the numerator and denominator predicts, the Mathematica function Limit[f[h], h -> ∞] gives $1$ in the nonsimplified case and gives $0$ in the simplified case.

Now, I am a little confused about whether the correct limit is $1$ or $0$ ? And is the function Limit supposed to work that way as to not being able to distinguish the correct result from the wrong one ?

Here is the Mathematica code

a = (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2/(
4 (h^2 +
1/4 (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2)) //
FullSimplify;

Limit[a, h -> Infinity]


as you know, you can add or remove the //FullSimplify at the end of the definition for $a$ to switch between the two cases.

UPDATE:

As mentioned in the comments below, I have reported this as a bug to Wolfram support. They wrote back, were very polite and thankful, mentioned that this issue has been forwarded to their developer team and recommended to use Assumptions in order to be sure that everything goes well and additionally even speed up the computation.

EDIT:

Fun fact - 5 years later, if I try to evaluate the above limit in Mathematica

$Version  "11.2.0 for Microsoft Windows (64-bit) (September 11, 2017)" I get exactly the same discrepancy in the two limits as described above. Which means, Wolfram has ignored this bug in the Limit routine in the meantime. • @Artes: OK, I will report this as a bug. – Kagaratsch Jan 22 '13 at 14:08 • Just to make sure - I am right in believing that 0 is the right answer, aren't I? – Kagaratsch Jan 22 '13 at 14:09 • @Kagaratsch: 1 is the correct result.... – Fabian Jan 22 '13 at 14:25 • @Kagaratsch Since the WRI people say they are aware of it, they'll say it is not a bug. Nevertheless one might expect a warning that some inconsitencies would appear. Since there is no warning one can encounter serious troubles. However Simplify and FullSimplify work well if one uses assumptions therein or in Limit. – Artes Jan 22 '13 at 17:20 • Actually Limit has been considerably improved in the time frame under question. This particular bug has not been fixed, but inductive reasoning from a sample size of 1 is probably not a statistically sound notion. – Daniel Lichtblau Jun 20 '18 at 14:00 ## 4 Answers There's a bit more to the story. Mathematica treats variables as complex by default, and I for one have had trouble figuring out how Limit figures out how to treat variables such as c in this case. Some analysis First, let's examine a0 (= a in OP) with the assumption thatc is real: a0 = (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2 / (4 (h^2 + 1/4 (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2)); aR = FullSimplify[a0, h > 0 && c \[Element] Reals] (* -> (h + c^2 h + Sqrt[4 + (1 + c^2)^2 h^2])/(2 Sqrt[ 4 + (1 + c^2)^2 h^2]) *)  $$\frac{\sqrt{\left(c^2+1\right)^2 h^2+4}+c^2 h+h}{2 \sqrt{\left(c^2+1\right)^2 h^2+4}}$$ It simplifies a little if we split the fraction up: 1/2 + Factor[aR - 1/2] (* -> 1/2 + ((1 + c^2) h)/(2 Sqrt[4 + (1 + c^2)^2 h^2]) *)  $$\frac{\left(c^2+1\right) h}{2 \sqrt{\left(c^2+1\right)^2 h^2+4}}+\frac{1}{2}$$ And let's look at it if c is an imaginary number k I: aI = FullSimplify[a0 /. c -> k I, h > 0 && k \[Element] Reals] (* -> (h - h k^2 + Sqrt[4 + h^2 (-1 + k^2)^2])/(2 Sqrt[4 + h^2 (-1 + k^2)^2]) *)  $$\frac{\sqrt{h^2 \left(k^2-1\right)^2+4}-h k^2+h}{2 \sqrt{h^2 \left(k^2-1\right)^2+4}}$$ Again it simplifies if we split the fraction up: 1/2 + Factor[aI - 1/2] (* -> 1/2 - (h (-1 + k) (1 + k))/(2 Sqrt[4 + h^2 (-1 + k^2)^2]) *)  $$\frac{1}{2}-\frac{h (k-1) (k+1)}{2 \sqrt{h^2 \left(k^2-1\right)^2+4}}$$ We see if c is +/-I (if k is +/-1), the fraction is 0, so the limit is 1/2. If k > 1 or k < -1, then the limit of the fraction is -1/2, so the limit of a0 is 0. And if -1 < k < 1, then the limit of the fraction is +1/2, so the limit of a0 is 1. If c is real, then up above one can see the limit of the fraction is 1/2, so the limit of a0 is 1. Here are some examples: Table[Print["c = ", c, " --> ", HoldForm[Limit[a0, h -> Infinity]], " = ", Limit[a0, h -> Infinity]], {c, {1, I/2, I, 2 I}}];  c = 1 --> Limit[a0, h -> Infinity] = 1 c = I/2 --> Limit[a0, h -> Infinity] = 1 c = I --> Limit[a0, h -> Infinity] = 1/2 c = 2 I --> Limit[a0, h -> Infinity] = 0  ## How to find the limit sought One should probably use assumptions: a1 = Simplify[a]; a2 = FullSimplify[a]; Block[{$Assumptions = c \[Element] Reals}, Limit[a0, h -> Infinity]]
Block[{$Assumptions = c \[Element] Reals}, Limit[a1, h -> Infinity]] Block[{$Assumptions = c \[Element] Reals}, Limit[a2, h -> Infinity]]
(* -> 1, 1, 1 *)


Prompted to investigate further by various comments, I have something to add, including what the actual limit for a generic complex $c$ is (which one might expect Mathematica could return as the correct answer). This goes beyond the current version of the OP's question, which stipulates that $c$ be real, but I hope it sheds some light on the issues Limit has dealing with this function.

What should the limit be?

Working it out by hand, I get the limit for a complex $c$ to be

$$\begin{cases} 1 & {\rm Im}(c)^2-{\rm Re}(c)^2<1 \\ {1}/{2} & {\rm Im}(c)^2-{\rm Re}(c)^2=1 \\ 0 & {\rm Im}(c)^2-{\rm Re}(c)^2>1 \end{cases}$$

Here is a plot of the real part of a0 vs. c when h = 100. The imaginary parts are nearly zero, and the real part is nearly equal to the limit. You can see the discontinuity forming along the red cliff. The limits at six test values are shown by the green spheres. We will use these test values for c again later.

testvalues = {1, I/2, I, 2 I, 1 + 2 I, 2 + 2 I};

Show[ParametricPlot3D[
Evaluate[{Re[c], Im[c], Re[a0]} /. {c -> r E^(I q), h -> 100.}],
{r, 0, 3}, {q, 0, 2 π}, PlotPoints -> {50, 150}, MaxRecursion -> 5,
Exclusions -> None, MeshFunctions -> {#2^2 - #1^2 &, #4 &, #5 &},
MeshStyle -> {Red, GrayLevel[0.5], GrayLevel[0.5]},
Mesh -> {{-1, 0, 1, 2}; {0.97, 1.03}, 11, {0, π/2, π, (3 π)/2}},
MeshShading -> {{{Yellow, Red, Lighter[Purple, 0.5]}}},
BoundaryStyle -> GrayLevel[0.5], Lighting -> "Neutral", PlotRange -> Full,
AxesLabel -> {Re[c], Im[c], ""}],
Graphics3D[{Darker@Green,
Sphere[{Re[#], Im[#], Limit[a0 /. {c -> #}, h -> Infinity]}, 0.07] & /@
testvalues}]
] There are spikes (poles) where the denominator of a0 is zero:

Solve[((4 (h^2 + 1/4 (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2)) /. h -> 100) == 0, c] // N
(* -> {{c -> -0.0099995 + 1.00005 I}, {c -> 0.0099995 - 1.00005 I}, {c -> -0.0099995 - 1.00005 I}, {c -> 0.0099995 + 1.00005 I}} *)


How could Mathematica get different answers?

First, it is not a problem with Simplify (a1 above) or FullSimplify (a2 above). All give the same correct limits with definite values for c:

Limit[a0 /. {c -> #}, h -> Infinity] & /@ testvalues
Limit[a1 /. {c -> #}, h -> Infinity] & /@ testvalues
Limit[a2 /. {c -> #}, h -> Infinity] & /@ testvalues

(* -> {1, 1, 1/2, 0, 0, 1}
{1, 1, 1/2, 0, 0, 1}
{1, 1, 1/2, 0, 0, 1} *)


I came up with ways to arrive algebraicallly at the different answers for the limits of a0 and a2 for a generic c. I doubt this is how Mathematica arrives at its answers, but it shows how the mistake might arise.

To find a limit at infinity, we can replace h by 1/k and let k approach 0. In calculus we learn to simplify the expression until we can plug in k = 0.

b0 = Simplify[a0 /. h -> 1/k, k > 0]
b2 = Simplify[a2 /. h -> 1/k, k > 0]

(* -> (1 + c^2 + Sqrt[1 + 2 c^2 + c^4 + 4 k^2])^2/(2 (1 + c^4 + 4 k^2 + Sqrt[1 + 2 c^2 + c^4 + 4 k^2] + c^2 (2 + Sqrt[1 + 2 c^2 + c^4 + 4 k^2]))) *)
(* -> (2 k^2)/(1 + c^4 + 4 k^2 - Sqrt[1 + 2 c^2 + c^4 + 4 k^2] - c^2 (-2 + Sqrt[1 + 2 c^2 + c^4 + 4 k^2])) *)

b0 // TeXForm // Print
b2 // TeXForm // Print


$$\frac{\left(c^2+\sqrt{c^4+2 c^2+4 k^2+1}+1\right)^2}{2 \left(c^4+c^2 \left(\sqrt{c^4+2 c^2+4 k^2+1}+2\right)+\sqrt{c^4+2 c^2+4 k^2+1}+4 k^2+1\right)}$$

$$\frac{2 k^2}{c^4-c^2 \left(\sqrt{c^4+2 c^2+4 k^2+1}-2\right)-\sqrt{c^4+2 c^2+4 k^2+1}+4 k^2+1}$$

So far, so good:

b0 == b2 // Simplify
(* -> True *)


For b0, setting k to zero yields 1 the same as Limit[a0, h -> Infinity]

b0 /. k -> 0 // TeXForm


$$\frac{\left(c^2+\sqrt{c^4+2 c^2+1}+1\right)^2}{2 \left(c^4+\left(\sqrt{c^4+2 c^2+1}+2\right) c^2+\sqrt{c^4+2 c^2+1}+1\right)}$$

which we can see by multiplying out the numerator and simplifying:

b0 /. k -> 0 // Expand // Together
(* -> 1 *)


Surprisingly, Simplify does something else and returns

b0 /. k -> 0 // Simplify // TeXForm


$$\frac{c^2+\sqrt{\left(c^2+1\right)^2}+1}{2 c^2+2}$$

whose value depends on c:

% /. {c -> #} & /@ testvalues
(* -> {1, 1, Indeterminate, 0, 0, 1} *)


Turning to b2, setting k to zero yields 0, the same as Limit[a2, h -> Infinity]

b2 /. k -> 0
(* -> 0 *)


As with a0 and a2, the limits of b0 and b2 at definite values of c are correct:

Limit[b0 /. {c -> #}, k -> 0] & /@ testvalues
Limit[b2 /. {c -> #}, k -> 0] & /@ testvalues
(* -> {1, 1, 1/2, 0, 0, 1}
{1, 1, 1/2, 0, 0, 1} *)


## A related example

Change 1+c^2 to 1-c^2:

aIm = a0 /. c -> I c
(* -> (h^2 - c^2 h^2 + Sqrt[4 h^2 + (h^2 - c^2 h^2)^2])^2/(4 (h^2 + 1/4 (h^2 - c^2 h^2 + Sqrt[4 h^2 + (h^2 - c^2 h^2)^2])^2)) *)


$$\frac{\left(-c^2 h^2+\sqrt{\left(h^2-c^2 h^2\right)^2+4 h^2}+h^2\right)^2}{4 \left(\frac{1}{4} \left(-c^2 h^2+\sqrt{\left(h^2-c^2 h^2\right)^2+4 h^2}+h^2\right)^2+h^2\right)}$$

Now I'm primarily interest in real c and the limit at h -> Infinity. It is similar to previous one above:

$$\begin{cases} 1 & \left| c\right| <1 \\ {1}/{2} & \left| c\right| =1 \\ 0 & \left| c\right| >1 \end{cases}$$

Mathematica gives one of the answers by default

Limit[aIm, h -> Infinity]
(* -> 1 *)


At definite values of c, we get the right answers

Limit[aIm /. {c -> #}, h -> Infinity] & /@ {0, 1, 2}
(* -> {1, 1/2, 0} *)


With various assumptions, one cannot get an answer that depends on c.

Limit[aIm, h -> Infinity, Assumptions -> c \[Element] Reals]
(* -> 1 *)


The correct answer appears in some cases:

Limit[aIm, h -> Infinity, Assumptions -> c > 1 || c < -1]
Limit[aIm, h -> Infinity, Assumptions -> c == 1]
Limit[aIm, h -> Infinity, Assumptions -> c == -1]
Limit[aIm, h -> Infinity, Assumptions -> -1 < c < 1]
(* -> 0, 1/2, 1/2, 1 *)


These are wrong; the limit is 1/2 in both cases:

Limit[aIm, h -> Infinity, Assumptions -> c \[Element] Reals && Abs[c] == 1]
Limit[aIm, h -> Infinity, Assumptions -> c == -1 || c == 1]
(* -> 0, 0 *)


## Conclusion

Whatever definite value I set c to, I always got the correct limit, but Limit does not handle c as a symbol correctly. With some careless algebraic operations, I was able to derive the answers Limit gave for original function and the FullSimplified version. Even with Assumptions that yield a unique limit, Limit does not always give the right answer. I would be surprised if it were mathematically impossible to develop an algorithm for finding a large class of limits including the OP's, given what can be done for integrals and by standard calculus techniques. I've never found a lot of use for Limit, but perhaps it is because of its limitations.

• I was composing my response and so missed yours until after I posted. What I had to say is pretty much what you said (use assumptions). – Daniel Lichtblau Jan 22 '13 at 15:35
• @DanielLichtblau I'm not surprised that assumptions affect the limits but I'd call it a bug that we have different values of limits with Simplify and FullSimplify. – Artes Jan 22 '13 at 16:07
• @Artes No, these discrepancies all follow from the same underlying weaknesses in handling of radicals. The discrepancies themselves are not bugs because the Simplify/FullSimplify rewritings are not incorrect. – Daniel Lichtblau Jan 22 '13 at 16:12
• @DanielLichtblau Nonetheless I'm not pleased with this issue. I would expect a warning that one should make assumptions on the variables rather than simply different values from limits. Moreover you can see that Limit[Simplify[...],h->Infinity] yields different results in M8 than in M7 and M9. – Artes Jan 22 '13 at 16:18
• @Kagaratsch They're pretty much equivalent. With Block you can temporarily reset global variables. $Assumptions is used by many functions, including Limit. All or most such functions provide an Assumptions option to do the same thing. Using the option seems easier here, but I've got used to using $Assumptions inside a Block when testing things. You can cut and paste things into it, and my assumptions are up front and not buried in the options. You can use Block to wrap several commands that all have the same assumptions, for instance. – Michael E2 Jan 22 '13 at 19:24

It is a bug in Series. Note that

a := (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2/( 4 (h^2 + 1/4 (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2))
b = FullSimplify[a]
Series[a,{h,Infinity,0}]
(* Out: 1 + O[1/h]^2 *)
Series[b,{h,Infinity,0}]
(* Out: O[1/h]^2 *)


The fact is that for $h\to\infty$ there are two terms cancelling each other in the denominator of $b$, one of which is a square root.

In fact further assumptions make it work

Assuming[c >= 0, Simplify@Series[b, {h, Infinity, 0}]]
(* Out: 1 + O[1/h]^2 *)
Assuming[c <= 0, Simplify@Series[b, {h, Infinity, 0}]]
(* Out: 1 + O[1/h]^2 *)


This is a known limitation in Series and Limit. Series does not handle roots in a flawless manner. For example, here is an expansion at the branch point (zero) that is only "half" correct.

In:= Series[Sqrt[x^2], {x,0,2}]

3
Out= x + O[x]


This is fine for re(x)>0, but not so good for re(x)<0 in a neighborhood of the origin.

Limit has the ability to handle this situation if provided sufficiently strong assumptions. For this example one could do as follows.

a = (h^2 + c^2 h^2 +
Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2/(4 *(h^2 +
1/4 (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2));

a2 = FullSimplify[a];

Limit[a2, h -> Infinity, Assumptions -> Element[c, Reals]]

(* Out= 1 *)


This next one is actually not so good.

Limit[a, h -> Infinity]

(* Out= 1 *)


That is to say, while the original result was what the poster wanted, it is only by accident of haw some symbolic radicals were handled. Here is why it is not fully correct: for "many" values of c the limit is not 1.

Limit[{a, a2} /. c -> 2 I, h -> Infinity]

(* Out= {0, 0} *)


Indeed, there is an interplay between real and imaginary parts that determines the limiting behavior.

Limit[{a, a2} /. c -> 1 - 2 I, h -> Infinity]

(* Out= {0, 0} *)

Limit[{a, a2} /. c -> 3 - 2 I, h -> Infinity]

(* Out= {1, 1} *)


As for bug reports, this one is in the data base and has been for quite some time. It is unlikely to change unless and until there is some redesign both of Series and Limit.

• Yes, I figured we must have been writing simultaneously. As you said, our solution is basically the same: use assumptions. – Michael E2 Jan 22 '13 at 15:39

Let's define :

a1 = (h^2 + c^2 h^2 +
Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2/(4 (h^2 +
1/4 (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2));

a2 = a1 // Simplify

a3 = a1 // FullSimplify Mathematica 7 and 9

Limit[ #, h -> Infinity]& /@ {a1, a2, a3}

{1, 1, 0}


while assuming that c is a real number :

Limit[#, h -> Infinity, Assumptions -> c ∈  Reals] & /@ {a1, a2, a3}

{1, 1, 1}


Mathematica 8

Limit[#, h -> Infinity]& /@ {a1, a2, a3} assuming that c is a real number :

Limit[#, h -> Infinity, Assumptions -> c ∈  Reals] & /@ {a1, a2, a3}

{1, 1, 1}


Conclusions

All versions (7, 8, 9) return the same expressions respectively for a1, a2 and a3.

There is a serious problem in Limit at least since the version 7 ( seems like WRI emploees are aware of it). The way we simplify the expressions shouldn't affect its limits. So there should be a warning on possible inconsistencies.

• Thank you, I will report it. I accepted the answer of Fabian since he has fewer points than you and therefore would benifit more from it. Hope that is alright. – Kagaratsch Jan 22 '13 at 14:36
• @Kagaratsch it is recommended that you select the answer that actually helps you the most, or when you cannot decide the oldest answer. – Mr.Wizard Jan 22 '13 at 14:46
• @Mr.Wizard: In this case both answers correctly identified the problem as bug, so I was undecided. However, the answer of Fabian came 3 minutes earlier, so I am glad to know that apparently I did the right thing. – Kagaratsch Jan 22 '13 at 14:50