An apparently "simple" limit?

Question

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.

Answer 1

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

Then the answers are correct.

EDIT - Addendum

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.

Answer 2

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

Answer 3

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[4]:= Series[Sqrt[x^2], {x,0,2}]                                              

                3
Out[4]= 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[80]= 1 *)

This next one is actually not so good.

Limit[a, h -> Infinity]

(* Out[63]= 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[81]= {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[85]= {0, 0} *)

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

(* Out[86]= {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.

Answer 4

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.