# Limiting form of a polynomial expression

When simplifying an expression by hand, a trick that is often used is to remove terms that are lower powers of the independent variable, for instance, as $x \rightarrow \infty$,

$x^2 + x$

becomes simply

$x^2$

I am trying to replicate this functionality in Mathematica. Limit seems like an obvious choice, but

Limit[x^2 + x, x -> ∞]


simply returns $∞$ (not terribly helpful). Another option is Series, but it seems that one has to pass the highest power in the expression to get the desired result, with

Normal@Series[x^2 + x, {x, ∞, -2}]


returning the expected $x^2$. However, this does not work in general for more complicated expressions, when the highest power may not be so obvious. Consider the following example:

$\frac{x^2+x}{x+1}$

We expect that the limiting form of this expression as $x \rightarrow \infty$ should be $x$ (the top becomes $x^2$, the bottom becomes $x$, and this then simplifies to $x$). However, if you just use the highest power available, Series returns 0 (you have to actually use Normal@Series[x^2 + x, {x, ∞, -1}]).

My question is, can anyone produce a generalized piece of code that can produce a limiting expression? It would be good if the procedure handled the other extreme as well, when the variable in question approaches 0 (in which case $x^2 + x$ would simplify to $x$).

Edit: It seems that the main issue for the solutions presented so far is that the form of the expression is not consistent. These operations are relatively easy to perform by hand, but sometimes tedious, which is why it would be nice to have Mathematica do it for me! Here's a pathological example that I can reduce by hand and for which the accepted solution should be able to return the correct result:

$\frac{\sqrt{x^3 + x}}{x^5 + x + 1} \frac{x^2}{x^{10} + x} + \frac{x^5 + x}{x}$

(The right answer is $x^4$)

-
That solution seems to fail for the example I provided: $\frac{x^2+x}{x+1}$, it returns $O[x^2]$ instead of the expected $O[x]$. – Guillochon Aug 28 '12 at 17:05
Something like 3 x^2 + x // Extract[#, Exponent[#, x, Max]] &? – J. M. Aug 28 '12 at 17:11
@J.M. Unfortunately, that too does not work for my fraction example. – Guillochon Aug 28 '12 at 17:17
In that case, a combination of Cancel[] and Apart[] might do you good... – J. M. Aug 28 '12 at 17:19
@Guillochon That is easily fixed by a Simplify. I get $O[x]^4$ for the second (as do both of Daniel's approaches). Could you show your working as to how you get $x^3$? – R. M. Aug 28 '12 at 17:52

The solution to this requires just a minor modification of my answer here:

ClearAll[limitingTerm1]
SetAttributes[limitingTerm1, HoldAll]
limitingTerm1[expr_, var_Symbol, func_] := var^Exponent[#, var, func] &@Simplify[expr]


This gives the right answer for both the examples provided by the OP:

limitingTerm1[(x^2 + x)/(x + 1), x, Max]
(* x *)

limitingTerm1[Sqrt[x^3 + x]/(x^5 + x + 1) x^2/(x^10 + x) + (x^5 + x)/x, x, Max]
(* x^4 *)


One can now generalize this to polynomials with arbitrary coefficients. It works as a two-pass approach, where the second pass is applied to the result of the first if its denominator still contains an x term. This seems to work very well in examples I've tested, but it might be possible to find edge cases where it doesn't work. Nevertheless, some modifications and extensions along these lines should help solve it.

ClearAll[limitingTerm2]
SetAttributes[limitingTerm2, HoldAll]
limitingTerm2[expr_, var_Symbol, func_] := Module[{term},
term := With[{pow = Exponent[#, var, func]}, Coefficient[#, var, pow] var^pow] &;
If[FreeQ[Denominator@#, x], #, term@Apart@#] &@term@Simplify[expr]
]

limitingTerm2[(a x^2 + bx)/(d + c x), x, Max]
(* (a x)/c *)


It also gives the correct answer for the other examples above.

-
For the second example, if I pass Min, it returns $O[x]^1$. Shouldn't it be $O[x]^0$? The first term reduces to $x^{3/2}$, the second to 1, so the limiting order is 0th, no? – Guillochon Aug 28 '12 at 18:09
@Guillochon Yes, I had removed the constant term for the previous usage, which was left behind in this. The above definition is far simpler and should give you 0 if you pass Min for the second example. – R. M. Aug 28 '12 at 18:17
Great...I guess the last bit of functionality I would like is to include any coefficients, if they exist. So for instance $\frac{ax^2 + bx}{cx + d}$ should return $\frac{a}{c}x$. – Guillochon Aug 28 '12 at 18:21
This is a great answer, as it works the same way one would calculate the asymptotic expansion by hand. In the case of rational polynomials it's very straight forward, but for more complicated functions, one has to calculate the divergent part first, then feed the convergent part to Series. – Pragabhava Aug 28 '12 at 18:23
@Guillochon See my edit – R. M. Aug 28 '12 at 18:42

Method I

(1) Take series at infinity to get a (Laurent) polynomial.

(2) Find largest exponent.

(3) Find corresponding coefficient.

mainTerm[expr_, x_] := Module[
{approxpoly, prec = 1, j = 0, expon, coef},
approxpoly = Normal[Series[expr, {x, Infinity, prec}]];
While[approxpoly === 0 && j <= 5,
j++;
prec = 2*prec + j;
approxpoly = Normal[Series[expr, {x, Infinity, prec}]];
];
approxpoly = PowerExpand@approxpoly;
expon = Exponent[approxpoly, x];
coef = Coefficient[approxpoly, x, expon];
coef*x^expon
]


Example:

In[83]:= mainTerm[(x^3 + 1)/(x + 1), x]

(* Out[83]= x^2 *)


Method II

(1) Take the series at infinity to get a Laurent polynomial.

(2) Homogenize it wrt a new variable.

(3) Set that new variable to zero.

homogenize[poly_, vars_, new_] :=
Module[{degfunc, totdeg, j},
degfunc = Apply[Plus, Map[Function[x, Exponent[#, x] &], vars]];
degfunc = Distribute[degfunc, Function, Plus];
totdeg = Max[Map[degfunc, Apply[List, poly]]];
Apply[Plus,
Table[poly[[j]]*new^(totdeg - degfunc[poly[[j]]]), {j,
Length[poly]}]]]

mainTerm2[expr_, x_] := Module[
{approxpoly, prec = 1, j = 0, hpoly, y},
approxpoly = Normal[Series[expr, {x, Infinity, prec}]];
While[approxpoly === 0 && j <= 5,
j++;
prec = 2*prec + j;
approxpoly = Normal[Series[expr, {x, Infinity, prec}]];
];
hpoly = homogenize[approxpoly, {x}, y];
hpoly /. y -> 0
]


Example:

In[91]:= mainTerm2[(x^3 + 1)/(x + 1), x]

(* Out[91]= x^2 *)


--- edit ---

For getting the asymptotic behavior at the origin you can do as follows.

(1) Replace x by 1/x. (2) Get the asymptotic term at infinity. (3) Transform back with x->1/x.

mainTermAtOrigin[expr_, x_] :=
mainTerm[expr /. x -> 1/x, x] /. x -> 1/x


Example:

mainTermAtOrigin[(x^2 + 1)/(x + 1) - 1/x, x]

(* Out[106]= -(1/x) *)


--- end edit ---

--- edit 2 ---

After adding a PowerExpand to deal appropriately with fractional power issues, we can wotk with Puisieux expansions. Here is an example.

ee = Sqrt[x^3 + x]/(x^5 + x + 1)*x^2/(x^10 + x^7) + (x^5 + x)/x;

In[170]:= mt = mainTermAtOrigin[ee, x]

(* Out[170]= (1/x)^(9/2) *)


We can test that this is correct by checking whether the expression divided by it approaches unity as x goes to zero.

In[171]:= Limit[ee/mt, x -> 0]

(* Out[171]= 1 *)

-
Verified that this works for my pathological example. How about as $x \rightarrow 0$? – Guillochon Aug 28 '12 at 17:52
Also, any particular reason for the j <= 5? – Guillochon Aug 28 '12 at 17:57
The limitation on j is just to prevent it from going forever on pathological examples. It will of course break in such cases. I did not put in error checking but that isn't too hard. – Daniel Lichtblau Aug 28 '12 at 18:04
See edit for handling at origin – Daniel Lichtblau Aug 28 '12 at 18:11
@Guillochon You can't expect it to work straight forward. When $x \rightarrow 0$ you are looking for the smallest power, while $x \rightarrow \infty$, you are looking for the greatest. If you are only interested on the order, a simple solution for rational functions would be something like (Exponent[Numerator[#],x]-Exponent[Denominator[#],x]) & @ rational[x] for the case $x \rightarrow \infty$. In the case $x \rightarrow 0$, (Exponent[Numerator[#],x,Min]-Exponent[Denominator[#],x,Min]) & @ rational[x] – Pragabhava Aug 28 '12 at 18:38

This works for me:

Normal@Series[((x^2 + x ) / (x + 1) ) , {x, Infinity, 1}] -> x


or this..

Normal@Series[((x^2 + x ) / (x + 1) /. x -> 1/y ) , {y, 0, 1}] /.  y -> 1/x


Note this example is a case where there is not a Taylor series, but Series[] is being smart for you and returning a negative power solution, which should give a clue why a fully general method is going to be problematic.

Edit another approach:

n = 0;
While[ (lim = Limit[(1/x)^n (x^3 + x ) / (x + 1)  , x -> Infinity]) ==  Infinity , n++];
lim x^(n)


Obvoulsly you'd want to limit n to prevent runaway.

I'm suspecting your question isn't mathematically well posed. Suppse the "limit" is a Log[x], do you expect your 'general' approach to pick that up?

-
Unfortunately, this solution is no better than the series solution I prescribed, as it only works for that specific fraction. Change the $x^2$ to $x^3$ and it no longer works. I'm looking for a solution that can simplify arbitrary polynomial expressions. – Guillochon Aug 28 '12 at 17:30
No, only polynomial expressions. I have edited the question title. – Guillochon Aug 28 '12 at 18:01