# Limit as n -> ∞ for RSolve with Multivariate Recurrence Relation

I'm trying to solve the following when n -> ∞:

RSolve[{
a[n + 1] == a[n] + x (1 - a[n] - b[n]),
b[n + 1] == b[n] + y (1 - a[n] - b[n]),
a[0] == 0, b[0] == 0},
{a[n], b[n]}, n]


After a FullSimplify, I get the following for a[n] (and similarly for b[n]):

It's clear that as n gets larger that this expression converges so long as 0 < x < 1 and 0 < y < 1, but simply plugging in n = ∞ I get "Sum does not converge" as an error. For instance when I set n = 20 and x = 1/3 and y = 1/4 I get the following:

547679985297149068793/
958439998111868780544


For n = 50:

130006259285744996749134383652034872751260675268341857/
227510953750053744333189631883564158123178815081414656


In floating point notation these both are ~0.571429, which is approximately 4/7 (thanks @hftf), suggesting the answer is (1/y)(1/x + 1/y). How does one get RSolve to return this result?

-
You mean, 4/7? –  hftf Apr 15 at 21:23
Did you do the simplification with assumptions on x,y, or (especially) n? If not, that could help. –  Daniel Lichtblau Apr 15 at 21:24
@hftf ...yes, haha. So I guess the answer is just (1/y)/(1/y + 1/x) for a[n]? But how to find that from RSolve... –  Guillochon Apr 15 at 21:31
@DanielLichtblau I tried, but the two Sums remain, regardless of the assumptions. The best I was able to do is get rid of the UnitStep[n-1] from the numerators when assuming n >= 1. –  Guillochon Apr 15 at 21:34

Look at the unsimplified UnitStep nonsense in Mathematica’s result. For example, the expression
(1 - x - y)^K[1] UnitStep[-1 + K[1]] + UnitStep[-K[1]] really just ought to be (1 - x - y)^K[1] instead. So let’s write a unitStepSimplify function to simplify expressions containing that, and apply it:

unitStepSimplify = # /. {UnitStep[-i_] + f_ UnitStep[i_ - 1] -> f} &;

closedForms = RSolve[{
a[n + 1] == a[n] + x (1 - a[n] - b[n]),
b[n + 1] == b[n] + y (1 - a[n] - b[n]),
a[0] == 0, b[0] == 0}, {a[n], b[n]}, n];

simplifiedClosedForms =
FullSimplify[closedForms // unitStepSimplify,
Assumptions -> n >= 1]

{{a[n] -> (x - x (1 - x - y)^n)/(x + y),
b[n] -> (y - (1 - x - y)^n y)/(x + y)}}


Nice. To find the n → ∞ limit:

% /. (f_ -> r_ ) :> (f -> Limit[r, n -> ∞, Assumptions :> {0 < x < 1, 0 < y < 1}]

{{a[n] -> x/(x + y), b[n] -> y/(x + y)}}


There you go:

$$a(n)=\frac x{x+y}\qquad\qquad b(n)=\frac y{x+y}$$

-