6
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ More methods on recursion can be found here $\endgroup$
    – garej
    Dec 21 '15 at 5:47
4
$\begingroup$

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

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

another way:

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

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]
$\endgroup$
2
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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