Devise a Newton iteration formula for computing $\sqrt[3]{R}$ where $R>0$. Perform a graphical analysis of your function $f(x)$ to determine

Question

Devise a Newton iteration formula for computing $\sqrt[3]{R}$ where $R>0$. Perform a graphical analysis of your function $f(x)$ to determine the starting values for which the iteration will converge.

I know that $\sqrt[3]{R}$ is the root of $f(x)=x^3-R$, with which, the iteration formula we need is $x_{n+1}=x_n-(\frac{x_n^3-R}{3x_n^2})$, here I have the graphs for when $R=1,2,3$. But I'm not sure where to start the method so that convergence is assured, it will be in the $[0,1]$ interval since according to the graph the roots go through there? Thank you very much.

Answer 1

You can prove by induction that the Newton iteration sequence is decreasing and always greater than $\sqrt[3]{R}$ if you start at $x_1 > \sqrt[3]{R}$:

• If $x_n > \sqrt[3]{R}$ then $x_{n+1}-x_n=-(\frac{x_n^3-R}{3x_n^2}) < 0$, hence $x_{n+1} < x_n$.
• As $f$ is convex for $x > \sqrt[3]{R}$, you get that $x_{n+1} > \sqrt[3]{R}$.

As a decreasing lower bounded sequence converges, if you start with $x_1 > \sqrt[3]{R}$, the Newton iteration sequence converges. As $f$ is continuous, it can in that case only converge to $l= \sqrt[3]{R}$.

Answer 2

The typical Newton iteration requires a division at every step which is time-consuming. When finding the square root, the normal way is to iterate for $D^{-1/2}$ and in a later stage to effectively multiply by $D$. In this answer we will solve $$y=D^2-x^{-3}=0$$ Then we have $$y^{\prime}=3x^{-4}$$ So that $$x_{n+1}=\frac{x_n}3\left(4-x_n(Dx_n)^2\right)$$ In the end we will find $$Dx_{\infty}=D\times D^{-2/3}=\sqrt[3]{D}$$ That division by $3$ in the iterative step will be computed as a floating-point multiplication by $1.0/3$. The algorithm will converge if $0<x_0<\left(\frac4D\right)^{2/3}$ because the graph or $y(x)$ is increasing and concave down, so if the first iteration doesn't cause $x_1$ to be negative it will be in the convergence zone, $0<x_n<D^{-2/3}$.

There is a bit of a trick to finding a good starting value $x_0$. Since $$D^{-2/3}=2^{-\frac23\log_2(D)}$$ If we thought about the $\text{IEEE-754}$ double precision floating point representation of $D$ as an integer, then multiplying it by $-\frac23$ would have the effect of raising $D$ to the $-\frac23$ power. Of course the exponent has a bias so some magic number will have to be added before the integer, considered again as an $\text{IEEE-754}$ double precision floating point number is close enough to $D^{-2/3}$ to be a good starting value $x_0$. The multiplication by $\frac23$ will be carried out in fixed point arithmetic via an integer scale factor. This will get close enough that convergence is assured within $4$ iterations, so no testing for convergence will be carried out. Here is the program:

module cbrtmod
   use ISO_FORTRAN_ENV
   implicit none
   private
   public cbrt
   contains
      function cbrt(x) bind(C,name='cbrt') result(Dx)
         real(REAL64) Dx
         real(REAL64), value :: x
         integer(INT64), parameter :: &
            magic = int(Z'000000009fd60dd7',INT64), &
            scale = int(Z'00000000aaaaaaaa',INT64)
         real(REAL64), parameter :: third = 1.0_REAL64/3
         real(REAL64) D, xx, corr, Dx3, x3
         integer j
         xx = transfer((magic-shiftr(transfer(abs(x),magic),32))*scale,xx)
         Dx = x*xx
         corr = 4-(Dx)**2*xx
         Dx3 = Dx*third
         do j = 1, 3
            x3 = xx*third
            Dx = Dx3*corr
            xx = x3*corr
            corr = 4-(Dx)**2*xx
            Dx3 = Dx*third
         end do
         Dx = sign(Dx3*corr,x)
      end function cbrt
end module cbrtmod

program cbrt_test
   use ISO_FORTRAN_ENV
   use cbrtmod
   implicit none
   real(REAL64) x, r, maxerr
   integer i
   maxerr = 0
   do i = 100, 800
      x = i*0.01_REAL64
      r = sign(abs(x)**(1.0_REAL64/3),x)
      maxerr = max(maxerr,abs(cbrt(x)-r))
      x = -x
      r = sign(abs(x)**(1.0_REAL64/3),x)
      maxerr = max(maxerr,abs(cbrt(x)-r))
   end do
   write(*,*) 'maxerr = ',maxerr
end program cbrt_test

Maximum error was about $8.881784197001252E-016$.