# Calculating N-th step of recurrent sequence

I am figuring out how to calculate the $n$-th step of a square root approximation sequence: $$x_{n+1} = \frac12\left(x_{n} + \frac{a}{x_n}\right),$$ where $a$ is a number the root of which is to be found and $x_1$ is any positive real.

Trying to follow Mathematica's language lean towards functional paradigm, I have encountered an evaluation issue: namely, how to hold evaluation up until the function is called and particular $n$ value is known and more exactly how to evaluate nested Holds upon a call. Here is my code:

f[a_, 1] = 1/2 (1 + a);
f[a_, n_] = 1/2 Hold[(f[a, n - 1] + a/f[a, n - 1])];
mysqrt[a_] = With[{N = 3}, f[a, N]] // ReleaseHold // N;
mysqrt


This is the result I'm getting.

• More methods on recursion can be found here Dec 21 '15 at 5:47

Here's one direct way to define the recursion:

Clear[f];
f[n_, a_] := f[n, a] = 1/2 (f[n - 1, a] + a/f[n - 1, a]);
f[0, a_] := 1;


For example, f[10,1] gives 1 and

f[10, 2] // N

1.41421


which is probably pretty close to Sqrt.

• Aha, so my mistake was using "=" instead of ":=". Dec 20 '15 at 17:38
• And consequently that none of the Hold and ReleaseHold stuff is needed. Dec 20 '15 at 18:36

another way:

NestList[N[1/2 (# + a/#) &], 1, 5] /. a -> 16
{1, 8.5, 5.19118, 4.13666, 4.00226, 4.}


You can use RSolve whith the boundary condition x0

x0 = 1;
f[N_, A_, x0_] := Re[RSolve[{x[n + 1] == 1/2 (x[n] + a/x[n]), x == x0},
x[n],  n][[1, 1, 2]] /. {n -> N, a -> A}]
f[10, 2., x0]


You can also use FixedPoint or FixedPointList

fpl[n_Integer, a_, x0_: 1] :=
FixedPointList[(# + a/#)/2 &, x0, n]

fpl[10, 2.]

(*  {1, 1.5, 1.4166666666666665,
1.4142156862745097,
1.4142135623746899,
1.414213562373095,
1.414213562373095}  *)


In this case it converged within machine precision in fewer than 10 iterations.

The problem with the solution returned by RSolve[] (as used in another answer) is that it uses transcendental functions when the solution can be expressed entirely algebraically. Here's what RSolve[] returns:

RSolve[{x[k + 1] == (x[k] + a/x[k])/2, x == x0}, x, k]
{{x -> Function[{k}, Sqrt[a]*Coth[2^k ArcCoth[x0/Sqrt[a]]]]}}


For comparison, here is the closed form for the $k$-th iterate of the Newton-Raphson iteration for the square root, due to Kalantari and Kalantari:

nrsqrt[x_, x0_, k_Integer?NonNegative] := With[{r = Sqrt[x], e = 2^k},
r ((x0 + r)^e + (x0 - r)^e)/((x0 + r)^e - (x0 - r)^e)]


Compare:

Table[nrsqrt[2, 1, k] // Simplify, {k, 0, 5}]
{1, 3/2, 17/12, 577/408, 665857/470832, 886731088897/627013566048}

With[{x = 2, x0 = 1, n = 5}, NestList[(# + x/#)/2 &, x0, n]]
{1, 3/2, 17/12, 577/408, 665857/470832, 886731088897/627013566048}


The caveat, of course, is that nrsqrt[] is slower than the conventional method of using Nest[]/NestList[], but the closed form can still be useful in some applications (e.g. the study of fractional iterates).