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

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.

• 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? – NickD Mar 5 '18 at 17:30

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

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