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

Question

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?

Answer 1

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