Finding Limits in several variables

Question

Is there a way to find a limit of a multivariable function, like

$$\lim_{(x,y)\to (0,0)} f(x,y)$$ with Mathematica?

When $f$ is continuous, we can use $$\lim_{(x,y)\to (0,0)} f(x,y)=\lim_{(x,0)\to (0,0)} f(x,0)$$ or something similar, which could be implemented with code like

Limit[Limit[f[x,y],y->0],x->0]

When $f$ is not continuous, like for example $$f(x,y)=\begin{cases} \frac{xy}{x^2+y^2} & (x,y)\neq (0,0)\\ 0 & \text{else} \end{cases}$$ we could make at least a small control with something like

Table[Limit[f[x,y] /. y -> k * x, x -> 0], {k, 0, 5}]

which gives us a list of the limits $$\lim_{x\to 0} f(x,k\cdot x)$$ Taking $f$ to be $$ f(x,y)=\begin{cases} 1 & \exists \alpha \in \mathbb{Q} \text{ such that } x=\alpha \cdot y\\ 0 & \text{else}\\ \end{cases}$$ gives us a function where this won't work.

So, is there a way to compute those limits (or getting the result that the limits don't exist) in Mathematica?

Answer 1

Taking a limit depends on the path used to approach that limit.

Consider the function in the question:

f[x_, y_] := Piecewise[{{x y / (x^2 + y^2), x != 0 && y != 0}}, 0];
base = Plot3D[f[x, y], {x, -1, 1}, {y, -1, 1}, MeshStyle->Opacity[0.2], PlotStyle->Opacity[0.5]]

(A plot of its graph, saved here as base, appears in subsequent figures.)

For instance, we may approach the origin along any ray given by a nonzero vector $(u,v)$; such rays can be parameterized by $t \to t(u,v)$:

c[t_, {u_, v_}] := t {u, v}

The value $t=0$ corresponds to the origin. To find the limit, we need to "lift" the approach path up to the "elevations" determined by $f$. When the original path is parameterized as $(x(t), y(t))$, then the lift is parameterized by $\left(x(t), y(t), f(x(t), y(t))\right)$:

lift[t_, f_, c_, opts___] := With[{x = c[t, opts]}, Append[x, f @@ x]]

(In this definition I have provided a mechanism to pass the ray vector $(u,v)$ as an option to lift.) Let's graph some of these lifts for various rays, using (of course) ParametricPlot3D:

trailRays = ParametricPlot3D[Evaluate@lift[t, f, c, #], {t, -1, 1}, 
    PlotStyle -> Thickness[0.01], ColorFunction -> Function[{x, y, z, u}, Hue[u]]] & /@
 {{1, 0}, {1,1}, {1, -2}};
Show[base, trailRays, Boxed -> False]

Especially when you can manipulate this plot in Mathematica, it is evident how the lifts of the various curves approach different limiting elevations at the origin.

Here is a more interesting way to approach the origin: spiral in. This time the value $t=\infty$ corresponds to the limit at the origin:

b[t_] := {Cos[t], Sin[t]} /Sqrt[t]

Let's plot (part of) its lift:

trail = ParametricPlot3D[Evaluate@lift[t, f, b], {t, 1, 30}, 
  PlotStyle -> Thickness[0.01], ColorFunction -> Function[{x, y, z, u}, Hue[u]]];
Show[base, trail, Boxed -> False]

As we spiral in toward the origin, the lift swoops up and down as the curve passes back and forth past all the incoming rays infinitely many times, never approaching a definite limit. This is a partial plot of the elevation as a function of the parameter $t$ along the curve; the graduated hues match those of the preceding plot:

Plot[f @@ b[t], {t, 0, 30}, PlotStyle -> Thick, AxesLabel -> {t, f}, ColorFunction -> (Hue[#1] &)]

(Analogous plots along the rays would be uninteresting, because along any ray through the origin, the value of $f$ does not vary at all!)

We may confirm visual impressions by applying Limit. The whole point is that the limit is taken along a curve, so it involves only a single (real) variable. Using the work we have already done, this is easy. Thus:

Limit[f @@ c[t, {u, v}], t -> 0]

$\begin{array}{ll} \{ & \begin{array}{ll} \frac{u v}{u^2+v^2} & (\text{Im}[u]\neq 0\|\text{Re}[u]\neq 0)\&\&(\text{Im}[v]\neq 0\|\text{Re}[v]\neq 0) \\ 0 & \text{True} \end{array} \end{array}$

There is a definite limit along each ray whose value is given here in terms of the ray vector's coordinates $(u,v)$. How about approaching the origin along the spiraling curve $b$?

Limit[f @@ b[t], t -> Infinity] // FullSimplify

Mathematica evaluates and simplifies f @@ b[t], but otherwise it--correctly--cannot obtain any limit and so just spits out another Limit expression.

---

Note that to study limits in more than one dimension, it does not suffice to study limits along rays (or lines) only. One can construct "nastier" functions f which have no definite limit when the origin is approached along a line, but do have definite limits when approached along particular spirals (or other curves of your choice). For instance, take our spiraling path $b$. At every point $(x,y)$ not at the origin we may locate two "arms" of the spiral, a nearest one and a next-nearest one, at distances $d_0$ and $d_1$, respectively. Let $f(x,y)$ be $d_0^2/d_1^3$. Because these distance functions are continuous and $d_1$ is never zero, $f$ is continuous. At all points equidistant between two arms (and there are many of these spiraling into the origin), $d_0=d_1$ and so $f = 1/d_1$, which clearly grows unbounded as the origin is approached. Any ray into the origin will hit infinitely many such points. But if we stick along one of the original arms of the spiral, $d_0$ is constantly $0$ and so, therefore, is $f$, whence its limiting value along that arm is $0$.

Illustrating this are plots of $r^2 g(r, \theta)$ and $g(r, \theta)$ for the function $g = 4t(1-t)$ (defined in polar coordinates $(r,\theta)$) where $t =\mod(\frac{1}{r} - \frac{\theta}{2\pi})$ has a fairly innocent definition that nevertheless reproduces the qualitative structure of the preceding description, including the lack of any limits at the origin along rays. The left plot lowers the values of $g$ near the origin so we can see the structure; the right one shows the spiraling "fences" created by the graph of $g$. Underneath each plot is show (in black) the locus of points where $g=0$: the limiting value of $g$ at $(0,0)$ along this curve clearly is zero.

Answer 2

Edit: Updated to handle an arbitrary number of variables.

Here's a fairly naive way, trusting in Mathematica's power. Find the maximum and minimum of f in a square of "radius" d about the point being approached. If the maximum and minimum approach the same value as d tends to 0, then by the "squeeze" theorem, the value is the limit. One caveat: Maximize and Minimize are not guaranteed to find the global extrema (according to the documentation). If they miss, then the following might say a limit exists, when it doesn't.

lim[f_, vars_List -> point_List] := Module[{lim0, d, max, min},
  max = Maximize[{f, 0 < d < 1 && And @@ MapThread[-d <= #1 - #2 <= d &, {vars, point}]},
    vars, Reals];
  If[Head[max] === Maximize,(*Maximize failed to find a max*)
   Return[$Failed],
       max = Simplify[First@max, d > 0];
       min = Minimize[{f, 0 < d < 1 && And @@ MapThread[-d <= #1 - #2 <= d &, {vars, point}]},
         vars, Reals];
       If[Head[min] === Minimize,(*Minimize failed to find a min*)
        Return[$Failed],
    min = Simplify[First@min, d > 0]
    ]];
  lim0 = Limit[max, d -> 0];
  If[lim0 == Limit[min, d -> 0], lim0, "Does not exist."]
  ]

Examples:

lim[(x y)/(x^2 + y^2), {x, y} -> {0, 0}]
(* "Does not exist." *)

lim[1 + (x^2 y^2)/(x^2 + y^4), {x, y} -> {0, 0}]
(* 1 *)

lim[(x y^2)/(x^2 + x^3 + y^4), {x, y} -> {0, 0}]
(* "Does not exist." *)

lim[(x y z)/(x^2 + y^2 + z^2), {x, y, z} -> {0, 0, 0}]
(* 0 *)

It returns $Failed on Sin[x] + y. :)

Answer 3

If a limit exists that is independent of path when restricting to real variables, then, on an especially good day, Wolfram|Alpha might tell you what that limiting value might be. Here are some possibly interesting examples from a selection of one day's worth of multivariate limit queries.

WolframAlpha["Limit[(x^2+y^2)^(x^2*y^2),{x -> 0,y -> 0}]"]
WolframAlpha["Limit[(x*y^2)/(x^2 + y^4), {x -> 0, y -> 0}]"]
WolframAlpha["Limit[(y^x - 1)/(x*(y - 1)), {x -> 1, y -> 1}]]}]"]
WolframAlpha["Limit[y*Log[x^2 + y^2], {y -> 0, x -> 0}]"]
WolframAlpha["Limit[Sin[x*y]/(Sin[x]*Sin[y]), {x -> 0, y -> 0}]"]
WolframAlpha["Limit[(Sqrt[y] - Sqrt[Tan[x*y]])/
    (Sqrt[y] + Sqrt[Tan[x*y]]), {x -> 2, y -> 0}]"]
WolframAlpha["Limit[Log[1 + 2*x - y]/(y - 2*x), {x -> 0, y -> 0}]"]

This in no way contradicts what @whuber notes. For most purposes limits are path-dependent creatures. Also the methods used in W|A are a bit, shall I say, on the heuristic side. (Not likely to give an incorrect limit, but could well state none exists when in fact one does.)

--- edit ---

I thought I should say a bit more about limits in multivariate real space. If only to convince myself that I actually might know anything about them.

Here is some code for the case of bivariate limits at the origin of functions that are amenable to solving for critical points. Rational functions of modest degree tend to work in finite time.

critpts[f_, vars_, r_] := Module[
  {lam, polys, gb},
  polys = 
   Join[Numerator[
     Together[
      Grad[f, vars] - lam*Grad[vars.vars, vars]]], {vars.vars - r}];
  gb = GroebnerBasis[polys, Join[vars, {r}], {lam}, 
    MonomialOrder -> EliminationOrder];
  gb]

bivariateLimit[rat_] := Catch[Module[
   {vars = Variables[rat], cpts, r, curves, realcurves, vals, lims},
   If[Length[vars] != 2, Throw[$Failed]];
   cpts = critpts[rat, vars, r];
   curves = Solve[cpts == 0, vars, Cubics -> False, Quartics -> False];
   realcurves = Select[curves, FreeQ[Chop[# /. r -> .111], Complex] &];
   vals = rat /. realcurves;
   lims = Limit[vals, r -> 0];
   Union[lims]
   ]]

It returns values that will include the extremal limiting values along any approach path to $(0, 0)$ in $R^2$. It is not hard to alter the code to allow for approaching different points.

Here is an example.

bivariateLimit[(x^4 - y^2 + 3*x^2*y - x^2)/(x^2 + y^2)]

(* Out[18]= {-1} *)

As there is only one value, -1, it has a limiting value thereof.

The example, and the idea behind the method although not the particulars, are from

C. Cadavid, S. Molina, and J. D. Vélez. Limits of quotients of bivariate real analytic functions. Journal of Symbolic Computation 50:197-207. 2013.

This is, as best I can tell, an area that has seen but little attention in the literature. (Perhaps because, as has been pointed out, these limits don't really exist. Except when they do.)

--- end edit ---

--- edit #2 ---

I should add some remarks about this method.

(1) In theory it is not really limited to the bivaraite case. But at three variables it will flat-out hang.

(2) As @Michael E2 notes, it will balk on inputs that are ostensibly bivariate but evaluate to constant or univariate functions. This is of course readily redressed (but not here and not tonight).

(3) (Also from a comment that I decided I should get in here).

When it is a function of one variable, in disguise, then it will fail. Example that messes up:

bivariateLimit[(x - y)^2/((x - y)^4 + (x - y)^2)]

During evaluation of In[21]:= Power::infy: Infinite expression 1/0 encountered. >>

During evaluation of In[21]:= Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >>

During evaluation of In[21]:= Power::infy: Infinite expression 1/0 encountered. >>

During evaluation of In[21]:= Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >>

(* Out[21]= {1, Indeterminate} *)
 bivariateLimit[(x - y)^2/((x - y)^4 + (x - y)^2)]

A production code version would need to trap for these..

--- end edit #2 ---

Answer 4

As I said in my comment, the limit in the multidimensional space doesn't exist if the linear limit depends on the way you approach the point of interest (which can be chosen as the origin unless it's infinity). But if a limit exists, then of course you can find it by approaching from any direction.

A counter-example is

$$\frac{x}{x+y}$$

for $x, y\to 0$, which generically "should" be $0$. But there is an exception for the line $y=0$ along which a linear limit gives $1$. And you can get other values by following other paths to the origin, see below.

Here is a limiting procedure that identifies generic values automatically if they exist, under the condition that we approach the limit radially. There are many other ways of approach, but this choice will always work if a true limit exists. The function radialLimit works in arbitrary number of dimensions, specified by the list of variables supplied in the second argument:

radialLimit[f_, vars_List] := Module[{r, angles, localVars, n},
  n = Length[vars];
  angles = Array[C, n];
  localVars = 
   Reverse@Table[
     r Times @@ (Sin /@ angles[[1 ;; i - 1]]) Cos[angles[[i]]], {i, 1,
       n - 1}];
  localVars = 
   Insert[localVars, 
    First[localVars] /. Cos[angles[[n - 1]]] -> Sin[angles[[n - 1]]], 
    2] ; {Limit[f /. Thread[vars -> localVars], r -> 0, 
     Assumptions -> angles \[Element] Reals] // FullSimplify, 
   Thread[vars -> localVars] /. r -> "r"}
  ]

It uses hyperspherical coordinates (the generalization of polar coordinates), and returns the limit that was found, together with the mapping between the original function variables and the hyperspherical angles which are named C[i]. The limit is taken by doing $r \to 0$ where $r$ is the radial distance from the origin. If a unique limit does not exist, you'll get the dependence of the result on these angles. For example, here is the function used in Daniel's answer:

radialLimit[(x^4 - y^2 + 3*x^2*y - x^2)/(x^2 + y^2), {x, y}]

The limit is unique because no dependence on angles remains after simplification. The previous non-unique example yields

radialLimit[x/(x + y), {x, y}]

The example in whuber's post would be treated as follows:

f[x_, y_] := Piecewise[{{x y/(x^2 + y^2), x != 0 && y != 0}}, 0];

radialLimit[f[x, y], {x, y}]

With this function, you can now also calculate limits in arbitrary dimensions, e.g., 3D:

radialLimit[Sin[x^2 + y^2 + z^2]/(x^2 + y^2 + z^2), {x, y, z}]

Another example from Michael's answer:

radialLimit[(x y^2)/(x^2 + x^3 + y^4), {x, y}]

Michael's function states that the limit doesn't exist. Here the dependence on the path again shows up. What this shows is that getting an angle-independent result in the radial approach doesn't guarantee that a unique limit exists. Only the converse is true: you'll find an existing limit easily this way.