1
$\begingroup$

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.

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ You have to realize what Newton's method does geometrically: it takes the current guess, and draws a tangent to the graph at that point; then it calculates where that tangent intersects the x axis and takes that intersection as the next iteration point. Based on this, can you see where you might get into trouble? Hint: what happens if the tangent is horizontal? Where is it horizontal? Also what happens when you start with $x_0 < -1$ say? $\endgroup$ – NickD Mar 5 '18 at 17:30
1
$\begingroup$

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}$.

$\endgroup$
1
$\begingroup$

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

$\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.