Approximation using Laurent polynomial I would like to obtain an approximation of a function $1/(a x + b y)$ around $x=1$ in terms of a Laurent polynomial. I do not need the complete series; in fact perhaps the 1st order terms in $x$ and $x^{-1}$ may be sufficient for my purposes, leading to something of the form $\alpha_{-1} x^{-1} + \alpha_0 + \alpha_1 x$. However a Taylor expansion is not what I need. What is the correct way of obtaining such an expansion?
The next step is to expand in $y$ around 1 as well, obtaining something of the form $\alpha_{-1} x^{-1} + \beta_{-1} y^{-1} + \alpha_0 + \alpha_1 x + \beta_1 y$. Is this reasonable?
 A: Expanding around $x=1$ and $y=1$ seem to be the wrong approaches to take here. The standard approach is to factor the larger (in magnitude) of $ax$ and $by$ out from the denominator, and use the geometric series summation formula in reverse (equivalently, the binomial expansion formula). So:
$$\begin{align} \frac{1}{ax + by} & = \frac{1}{ax} \sum_{n=0}^\infty \left(-\frac{by}{ax}\right)^n \ \forall\, |ax| > |by|,\ \mathrm{and} \\
& = \frac{1}{by} \sum_{n=0}^\infty \left(-\frac{ay}{by}\right)^n \ \forall\, |ax| < |by| \end{align}.$$
If you want to expand about a different point, just shift your variables and apply the same formulae. For example:
$$\begin{align} x' &= x - 1,\ \mathrm{and} \\
y' & = y + \frac{a}{b}. \end{align}$$
To more literally apply repeated expansions, as requested in the original post, you'll apply a binomial expansion to each term in one of the above formulae after applying another shift/change of variables. The formula for more general binomial expansions is:
$$\begin{align}\frac{1}{(ax + b)^n} &= \frac{1}{(ax)^n} \sum_{m=0}^\infty \frac{\Gamma(-n+1)}{m! \Gamma(-n-m+1)} \left(\frac{b}{ax}\right)^m\ \forall\, |ax| > |b|,\ \mathrm{and} \\
& = \frac{1}{b^n} \sum_{m=0}^\infty \frac{\Gamma(-n+1)}{m! \Gamma(-n-m+1)} \left(\frac{ax}{b}\right)^m\ \forall\, |ax| < |b|. \end{align}$$
A: I didn't meet the Laurent polynomials of two variables yet, but it can be used the next considerations.
It is possible to build "the standard" Laurent polynomial only for complex analytic function of one variable $z=x+iy$. That gives
$$F(x+iy)=u(x,y)+iv(x,y),$$
so $u(x,y)$ and $v(x,y)$ must satisfy to Cauchy–Riemann equations
$$\dfrac{\partial{u}}{\partial{x}}=\dfrac{\partial{v}}{\partial{y}},\quad \dfrac{\partial{u}}{\partial{y}}=-\dfrac{\partial{v}}{\partial{x}}.$$
Assuming
$$u(x,y) = \mathcal RF(x+iy)=\dfrac1{ax+by},$$
we have 
$$\dfrac{\partial{v}}{\partial{y}} = -\dfrac{a}{(ax+by)^2},\quad \dfrac{\partial{v}}{\partial{x}} = \dfrac{b}{(ax+by)^2},\quad a=b=0.$$
That means that there are not Laurent polynomials for $f(x,y).$
Moreover, it can not be represented as real (or imaginary) part of this series.
At the same time, for the complex function
$$F(z)=1/z$$
$$F(x+iy)=\dfrac{x}{x^2+y^2}-i\dfrac{y}{x^2+y^2},$$
so the function
$$u(x,y)=\dfrac{x}{x^2+y^2}$$
can be presented as the real part of Laurent polynomial.
Within the real analysis, Taylor series of two variables in form
$$f(x,y)=f(x_0,y_0) + f'_x\biggr|_{x=x_0,y=y_0}(x-x_0) + f'_y\biggr|_{x=x_0,y=y_0}(x-x_0)(y-y_0)$$
$$+ \dfrac12f''_{xx}\biggr|_{x=x_0,y=y_0}(x-x_0)^2 + f''_{xy}\biggr|_{x=x_0,y=y_0}(x-x_0)(y-y_0) + \dfrac12f''_{yy}\biggr|_{x=x_0,y=y_0}(y-y_0)^2$$
$$+ \dfrac16f'''_{xxx}\biggr|_{x=x_0,y=y_0}(x-x_0)^3 + \dfrac12f'''_{xxy}\biggr|_{x=x_0,y=y_0}(x-x_0)^2(y-y_0) + \dfrac12f'''_{xyy}\biggr|_{x=x_0,y=y_0}(x-x_0)(y-y_0)^2+ \dfrac16f'''_{yyy}\biggr|_{x=x_0,y=y_0}(y-y_0)^3+\dots.$$
It's easy to see that Taylor series of $n$th order use the full set of partial derivatives to $n$th order.
To obtain the series in negative powers of $(x-x_0)$ and $(y-y_0),$ we can present $f(x,y)$ as $g\left(\dfrac1{x-x_0},\dfrac1{y-y_0}\right),$ where
$$g(r,s) = \dfrac{rs}{as+br+ax_0+by_0},$$ 
and use the above Taylor series in form of Maclaurin series.
It is easy to see that
$$g(0,0)=g'_r(0,0)=g'_s(0,0)=g''_{rr}(0,0)=g''_{ss}(0,0)=g''''_{rrr}(0,0)=g'''_{sss}(0,0).$$
Calculating the second
$$g''_{rs}\biggr|_{r=0,s=0} = \dfrac1{ax_0+by_0}$$
and the third
$$g'''_{rrs}\biggr|_{r=0,s=0} = -\dfrac{2b_0^2y_0}{(ax_0+by_0)^2}, \quad g'''_{rrs}\biggr|_{r=0,s=0} = -\dfrac{2a_0^2x_0}{(ax_0+by_0)^2}$$
derivatives, we can get "the Laurent" polynomial in form
$$\dfrac1{ax+by} = \dfrac1{(ax_0+by_0)^2(x-x_0)(y-y_0)}\left(ax_0+by_0 -\dfrac{2b_0^2y_0}{x-x_0}- \dfrac{2a_0^2x_0}{y-y_0}+\dots\right)$$
To choice the reqired form of series, it follows to consider the actual area. Maybe, the Maclaurin series on $((x-x_0)^2+(y-y_0)^2-r^2)$ or similar can satisfy to requirements.  
