Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to solve a system of recurrence relations as follows.

RSolve[{a[p] == 1 + (2/n) b[p - 1] + ((p - 1)/n) a[p - 1], 
        b[p] == 1 + (p/n) b[p - 1], 
        a[1] == 1 + 2/n, 
        b[0] == 1}, 
     {a[p], b[p]}, p]

Having done this, I would like to set $p=n-1$ and then plot the resulting functions $a[n-1]$ and $b[n-1]$ for $n=1, ..., 1000$.

I am failing at the first hurdle as Mathematica version 8 just returns the RSolve line when I press shift-enter.

Update. Following a request for an example using Maple 16. I do the following.

rsolve({a(1) = 1+2/n, a(p) = 1+(p-1)*a(p-1)/n+2*b(p-1)/n, b(0) = 1, b(p) = 1+p*b(p-1)/n}, {a(p), b(p)});
subs(p = n-1, %);

I then copy the result for a(n-1) and plot it

plot([seq([n, (1/n)^(n-1)*n*GAMMA(n-1)*(Sum((n+2*(1/n)^p1*exp(n)*GAMMA(p1+1, n))/(n^2*(1/n)^(p1+1)*GAMMA(p1+1)), p1 = 1 .. n-2)+1+2/n)], n = 2 .. 200)], style = point);

This gives

enter image description here

share|improve this question
1  
First issue is the syntax : RSolve with capital S. –  b.gatessucks Jan 1 '13 at 19:30
    
Strongly related mathematica.stackexchange.com/q/16919/193 –  belisarius Jan 1 '13 at 19:41
    
Sorry... copy and paste error. I have the same problem with RSolve with a capital S. Is there a way to restart so that there are no existing function/variable settings or is there some other reason it doesn't appear to do anything? –  lip1 Jan 1 '13 at 19:48
1  
solving for b[p] seperately produces a solution for b : solb = RSolve[{b[p] == 1 + (p/n)*b[p - 1], b[0] == 1}, b, p][[1, 1]]. Substituting this in the other equation seems to be more complicated though sola = RSolve[{a[p] == 1 + (2/n)*b[p - 1] + ((p - 1)/n)*a[p - 1], a[1] == 1 + 2/n} /. solb, a, p] –  Thies Heidecke Jan 1 '13 at 19:53

1 Answer 1

In two steps, since they can be decoupled:

h1 = RSolve[{
   b[p] == 1 + (p/n) b[p - 1],
   b[0] == 1}, {b[p]}, p]
h2 = RSolve[{
   a[p] == Evaluate[b[p] /. h1[[1]]] + ((p - 1)/n) a[p - 1],
   a[1] == 1 + 2/n}, {a[p]}, p]

DiscretePlot[{b[i] /. Evaluate[h1 /. p -> n - 1 /. n -> i + 1], Sqrt[i]}, {i, 300}]
DiscretePlot[ a[i] /. Evaluate[h2 /. p -> n - 1 /. n -> i + 1], {i, 300}]

Comparing b[i]with Sqrt[i]

Mathematica graphics

And this is a[i]

Mathematica graphics

Edit

BTW, answering your comment

Thanks. How can I now do the substitution and plotting? If there isn't a bug then it should look like $sqrt(n)$

here is a reasonable approximation to b[i] using Sqrt[]:

y[i_] := Sqrt@(1.78081659706051 + 1.57076044099885*i) - 0.332147990949542; 
DiscretePlot[Evaluate[b[p] /. h1[[1]] /. p -> n - 1 /. n -> i + 1] - y[i], {i,10000}, 
             PlotRange -> All]

Mathematica graphics

Or, if you want an exact result:

b[n] -> E^n n ExpIntegralE[1 - n, n]
share|improve this answer
    
Thanks. How can I now do the substitution and plotting? If there isn't a bug then it should look like $\sqrt{n}$. –  lip1 Jan 1 '13 at 19:56
    
Thanks but isn't it possible to use the answer to RSolve[{ a[p] == Evaluate[b[p] /. h[[1]]] + ((p - 1)/n) a[p - 1], a[1] == 1 + 2/n}, {a[p]}, p] without having to retype the answer for a[n]? –  lip1 Jan 1 '13 at 20:04
    
@lip1 Look at it now –  belisarius Jan 1 '13 at 20:13
    
Thanks very much. –  lip1 Jan 1 '13 at 22:30
    
@belisarius Don't you think Mathematica should be able to solve it the way you did ? –  b.gatessucks Jan 2 '13 at 7:59

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.