The roots have multiplicity
The situation for the given function is that the roots are at the same time maxima of the function, that is, they have multiplicity $2$, as
$$
f(x)=\left(1-2\sin^2\frac{33\pi}{2x}\right)\left(1-2\sin^2\frac{\pi x}{2}\right)-1
$$
so after expanding $-f(x)$ is a sum of squares minus the product of these terms. Methods that are developed for finding single roots will either slow down or fail to converge at roots of higher multiplicity. Newton's and Halley's method slow down.
There are many local extrema
Another problem with applying Newton is that this function has many local maxima and minima at small $x$ due to the first factor. There the derivative is zero, so that the Newton step, considered as function of $x$, has as many poles. Any improved method based on Newton's method will have as many or more poles, even if locally around the roots of $f$ the convergence is better.
Note that at a double root, where locally $f(x)=c(x-r)^2$, the Newton step maps $x$ to $\frac{x+r}2$ and the Halley step to $\frac{x+2r}3$. In the plots, this is somewhat visible around the roots $x=3$ and $x=11$.
Modifying Newton's method
Knowing this and the possibility of a double root, one can change the Newton step to alternate steps of single and double step size. Then at simple roots the single step will reduce the distance to the root quadratically while the following double step will overshoot the root, however with a smaller step size. At a double root the single step will reduce the distance by half, while the following double step will restore quadratic convergence. In each case, the "wrong" step does not make the situation worse, while the "right" step proceeds with the expected quadratic convergence.
Finding roots inside intervals
If an interval is small enough for a given function, then it either has no root inside the interval or it is contained in the basin of attraction of the root inside. Finding a subdivision of a given interval that is fine enough is again a heuristic task.
- When performing the iteration, if it leaves the given small interval then it has failed with a high probability of no root inside.
- Another failure condition is that the iteration enters a cycle. There might be a root inside the span of the cycle, but for simplicity let the iteration fail if after a small number of iterations the step size is not small relative to the interval length. With a good probability this means that the subdivision is not fine enough
- The convergence should now be at least linear, reducing the step size be one half each step. To guard against strange floating point effects, stop based on the iteration count after a number of iterations that theoretically should be sufficient to reach the desired accuracy.
- Of course, also stop if the desired accuracy is reached.
As a python code, this could look like
def find_roots(method,a,b,segments=10):
seg = np.linspace(a,b,segments+1);
for k in range(segments):
ak, bk = seg[k:k+2];
#print "searching for roots in",[ak,bk]
x = (ak+bk)/2;
count = 0;
while ak<=x<=bk and count < 50:
count += 1;
xold, x = x, method(x);
#print x
if count==2 and abs(x-xold)>1e-1*(bk-ak): break;
if abs(x-xold)<1e-8:
y,_,_ = f(x)
print "found root x=%.15f with f(x)=%.8e in %d iterations"%(x,y,count);
break;
Called as find_roots(method,2,12,segments=14) this returns the results
find roots with Newton step
found root x=3.000000007315551 with f(x)=-3.77475828e-15 in 23 iterations
found root x=10.999999991701889 with f(x)=-3.33066907e-16 in 23 iterations
find roots with Halley step
found root x=3.000000004913715 with f(x)=-1.66533454e-15 in 15 iterations
found root x=10.999999999234854 with f(x)=0.00000000e+00 in 16 iterations
find roots with Newton plus double Newton step
found root x=2.999999999980970 with f(x)=0.00000000e+00 in 4 iterations
found root x=10.999999999997232 with f(x)=0.00000000e+00 in 3 iterations
Note that in the last method, each iteration contains two Newton steps. If one counts the effort in function evaluations, then Newton gets a factor of $2$, Halley a factor of $3$, and the double step method a factor of $4$, giving the first two methods a similar complexity.
Appendix: More code
The method steps are standard implementations
def Newton_f(x): vf, df, _ = f(x); return x-vf/df
def Halley_f(x): vf, df, ddf = f(x); return x-(vf*df)/(df**2-0.5*vf*ddf)
def TwoStep_f(x):
vf,df,_ = f(x);
x = x - vf/df;
vf,df,_ = f(x);
return x - 2*vf/df;
The function implementation provides also the first and second derivative à la algorithmic differentiation (AD) in forward mode
def f(x):
v1 = 33*np.pi/x; dv1 = -v1/x; ddv1 = -2*dv1/x;
v2 = np.cos(v1); v3 = np.sin(v1);
dv2 = -v3*dv1; dv3 = v2*dv1;
ddv2 = -dv3*dv1-v3*ddv1; ddv3 = dv2*dv1+v2*ddv1;
v4 = np.pi*x; dv4 = np.pi; ddv4 = 0;
v5 = np.cos(v4); v6 = np.sin(v4);
dv5 = -v6*dv4; dv6 = v5*dv4;
ddv5 = -dv6*dv4-v6*ddv4; ddv6 = dv5*dv4+v5*ddv4;
return v2*v5-1, dv2*v5+v2*dv5, ddv2*v5+2*dv2*dv5+v2*ddv5;
The call of the root finder procedure is
names = ["Newton step", "Halley step", "Newton plus double Newton step"]
for k, method in enumerate([Newton_f, Halley_f, TwoStep_f]):
print "find roots with %s"%names[k];
find_roots(method,2,12,segments=14)
Best Answer
It is practical to recall that $z\mapsto\frac{1-z}{1+z}$ is an involution. In particular $$ \frac{1+\alpha}{1-\alpha}+\frac{1+\beta}{1-\beta}+\frac{1+\gamma}{1-\gamma} $$ is the sum of the reciprocal of the roots of $p\left(\frac{1-x}{1+x}\right)=-\frac{x^3-x^2+7x+1}{(1+x)^3}$, which is also the sum of the reciprocal of the roots of $x^3-x^2+7x+1$. By Vieta's formulas, this is $\color{red}{-7}$ (option 3.).