This question already has an answer here:

I would like to compute a[1000] of the following recursive sequence

a[1] = 1;
a[n_] := a[n - 1] + 1/a[n - 1]^2

but it takes a long time for computing even small numbers like a[30], however, if I used a last computed term like a[20] as the value of the first term, then Mathematica compute the a[30] faster than before. is there any method for computing this kind of values more faster?


marked as duplicate by Lukas Lang, AccidentalFourierTransform, Roman, march, Carl Lange May 23 at 9:02

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.


Even with memoization, this paricular problem becomes very slow very quickly, since mathematica struggles with saving a huge numer of numerator and denominator digits for rational numbers that grow fairly slowly in magnitude overall.

One answer to this problem may be to do the recursion completely numerically:

aN[1] = N[1, 10000];
aN[n_] := aN[n] = aN[n - 1] + 1/aN[n - 1]^2

The overall error $\Delta a$ grows with rounding error $\Delta r=10^{-10000}$ as $\Delta a=(1-2/a_{n-1}^3)\Delta r$, so considering that $a_n>2$ already at $n>2$, the above numerical approach should have constant calculation time for each step and have no appreciable errors for $n$ up to $n<9970$ (for 30 digits of precision remaining).

Indeed the following plot is generated in a second or two, and looks good:

ListLinePlot[Table[aN[i], {i, 1, 10000}], PlotRange -> All]

enter image description here

If you want to go even higher, increase the number of digits of precision accordingly. This will make the calculation a bit slower but give a bit more digits of precision.


Also, note that the resulting function of $n$ becomes increasingly flat as $n$ grows. Once the increase in a[n] versus a[n-1] becomes negligible, you can approximate the relation a[n]-a[n-1]=1/a[n-1]^2 as a'[n]=1/a[n]^2 and solve

a[n] /. DSolve[a'[n] == 1/a[n]^2 && a[1] == 1, a[n], n][[1]]

(-2 + 3 n)^(1/3)

This is the function which your sequence asymptotes to for n->Infinity:

    ListLinePlot[Table[aN[i]/(-2 + 3 i)^(1/3), {i, 200, 10000}], PlotStyle -> Green, PlotRange -> All], 
    ListLinePlot[Table[1, {i, 200, 10000}], PlotStyle -> Red, PlotRange -> All]
, PlotRange -> All]

enter image description here


Since you only want a numerical solution, Nest or NestList might be good alternatives:

ListPlot[NestList[# + 1/#^2 &, 1., 1000]]

enter image description here

Or to find the value of a[1000]:

Nest[# + 1/#^2 &, 1., 1000]

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