# Finding Limits in several variables

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?

-
With Limit, you're always restricted to a line in the larger space, and you can't make statements about the existence of the limit in the sense of the higher-dimensional space. For that you have to show the independence of the result on the direction of the line. If you intentionally set up a function to have different limits along different lines, I dont see what else you can do with Mathematica. –  Jens Mar 17 at 20:56
@Jens "Intentionally set up" makes it sound like we're discussing a mere pathology. However, if I recall correctly (after 30 years of not thinking about it), the Edge of the Wedge Theorem of multivariate complex analysis concerns an important natural class of functions (i.e., holomorphic in some domain) which potentially have precisely these bad limiting properties at their boundaries, unless the limiting curves are carefully controlled: and therein lies the need to restrict to wedges. –  whuber Mar 18 at 17:10
@Jens, After re-reading the original version, I still have the impression that you and the OP may be thinking about this question differently. In the example, although a unique limit does not exist, there is nevertheless a set of possible limits that can be achieved; it is also the set of accumulation points of $f$ at $(0,0)$. The question asks--in a clear manner, in my view--how to compute that set or at least to establish that a unique limit does not exist. (My answer does not fully address this--it was offered more as an extended comment. Michael E2 probably comes closest.) –  whuber Mar 18 at 17:42
Thank you @Daniel. As I confessed earlier, I'm dredging up memories a third of a century old :-). My recollection is that the proof of the theorem relies on analysis of how limits are approached at the boundary within carefully restricted curves. (There may be multiple proofs by now; my source for this is Steven Krantz's SCV book.) Regardless of the application, the point was to indicate that these issues about limits are not mere matters of "pathology" but can be important in Analysis generally. –  whuber Mar 23 at 23:58
@whuber I suspect you are on target about that proof (my memories are a mere 20+ years dated, but I have even less of an excuse since it was my area at one time). I certainly agree that limiting behaviors in the multivariate setting are of importance, in a way that is not always evocative of the term "pathological". –  Daniel Lichtblau Mar 24 at 0:05

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.

-
This is one very nice and illustrative answer. +1 –  halirutan Mar 17 at 22:25

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\[Rule]0,y\[Rule]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[
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. Vele'z. 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 --- - What is an especially good day? Full moon? Payday? – Rolf Mertig Mar 17 at 23:18 @Rolf Mertig Any day that Alpha limit code behaves is a good day.. – Daniel Lichtblau Mar 17 at 23:23 +1 Nice. Interesting: fails on bivariateLimit[(x^2 + y^2)/(x^2 + y^2)] but not on bivariateLimit[(2 x^2 + 2 y^2)/(x^2 + y^2)]. The use of GroebnerBasis suggests it will only work on rational functions. – Michael E2 Mar 18 at 2:00 @Michael E2 See edit #2 (and I'm deleteing my comment response since that covers it). – Daniel Lichtblau Mar 18 at 2:36 @DanielLichtblau Your example limit is trivial to calculate: FullSimplify[(x^4-y^2+3*x^2*y-x^2)/(x^2+y^2)/.{x->r Cos[\[Phi]],y->r Sin[\[Phi]]}]/.r->0 yields -1. This doesn't seem to be very deep to me. What am I missing? – Jens Mar 18 at 3:38 show 2 more comments 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. :)

-

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];



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}]


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