Calculating N-th step of recurrent sequence

Question

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[10]

This is the result I'm getting.

More methods on recursion can be found here – garej Dec 21 '15 at 05:47 — garej, Dec 21 '15 at 05:47

score 5 · Answer 1 · answered Dec 20 '15 at 15:44

5

another way:

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

answered Dec 20 '15 at 15:44

score 4 · Accepted Answer · answered Dec 20 '15 at 15:40

4

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[2].

answered Dec 20 '15 at 15:40

bill s

68,936
4
101
191

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

score 3 · Answer 3 · edited Feb 18 '19 at 00:59

3

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[0] == x0},
                            x[n],  n][[1, 1, 2]] /. {n -> N, a -> A}]
f[10, 2., x0]

edited Feb 18 '19 at 00:59

J. M.'s missing motivation

124,525
11
401
574

answered Dec 20 '15 at 15:51

Stratus

2,942
12
24

score 2 · Answer 4 · answered Dec 20 '15 at 17:10

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.

score 1 · Answer 5 · answered Feb 21 '16 at 16:27

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[0] == 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).

Calculating N-th step of recurrent sequence

5 Answers5