2
$\begingroup$

I have a trigonometric function; for instance $$f(x)=\left(\cos{\frac{33}{x}\pi}\right) (\cos{x \pi})-1$$

I wanted to know the zeroes of this particular function, so I thought I could look into some root-finding algorithms (Newton's, Halley's, Secant...). However, they don't seem to work as $f'(x)=0$ at the roots of $f(x)$, so all those methods are not guaranteed to converge.

So, I was thinking, is there some type of root-finding algorithm for this particular trigonometric equation? Or at least transform this equation into one that the roots would go through the x-axis rather than "bouncing" off it, so Newton's method would apply.

Also, I am focused on roots $>1$ and $<33$.

Note: Although the given example can be solved with trigonometric techniques, I am specifically looking for numerical methods. The example was chosen to make it easy to check the roots. I can generalize it to say for any $$f(x)=\left(\cos{\frac{n}{x}\pi}\right) (\cos{x \pi})-1$$ and an interval $$[a,b]$$ where there is only one root in that interval, is there a way to use numerical methods that are guaranteed to converge at the root to find that root?

$\endgroup$
7
  • $\begingroup$ $x=33$ is a root $\endgroup$
    – Vasili
    Jul 19, 2019 at 16:54
  • $\begingroup$ Ok yes, I know that, but how would you find the roots $3$ and $11$ numerically? $\endgroup$
    – DUO Labs
    Jul 19, 2019 at 16:55
  • $\begingroup$ Why do you need to find them numerically? You just need the argument of each cosine to be an integer multiple of $\pi$. Hence roots are divisors of $33$ such that $33/x\in\mathbb{Z}$. $\endgroup$ Jul 19, 2019 at 16:56
  • $\begingroup$ Yeah, but I picked this example so when found a method, I can check that it is correct, as the roots are easy to find. $\endgroup$
    – DUO Labs
    Jul 19, 2019 at 16:57
  • 1
    $\begingroup$ Could you give the code and iteration results for a method you tried? What does "converge to a wrong value" mean, what is wrong in that value? $\endgroup$ Jul 19, 2019 at 17:27

5 Answers 5

5
$\begingroup$

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.

enter image description here

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)
$\endgroup$
19
  • $\begingroup$ This is a great answer! +1 from me. Just one question, what are the segments for? Also, I just want one root, the first one. Any way to modify to get it? Also, what is method? But overall great answer! $\endgroup$
    – DUO Labs
    Jul 20, 2019 at 13:58
  • $\begingroup$ As said in the motivation, you need that the "working" interval are small enough so that when one contains a root, it is contained in the basin of attraction of that root (for the method used). Otherwise the first step from the selected midpoint might lead outside the interval, or a cycle of the method might be contained inside the interval. Usually it is not possible to find the first root, you will need to determine all roots and then in post-processing determine the smallest. Or you need an enclosing interval where you know that only one root is contained. $\endgroup$ Jul 20, 2019 at 14:05
  • $\begingroup$ Yeah I have a root where only one root is contained. Also, which python is this? $\endgroup$
    – DUO Labs
    Jul 20, 2019 at 14:07
  • $\begingroup$ This is python 2.7, for python 3 you need to adapt the print statements. You will need implementations of x_next=method(x) which are the single steps of the methods used. $\endgroup$ Jul 20, 2019 at 14:19
  • $\begingroup$ I'm sorry if I don't understand, what implementation of method did you use? $\endgroup$
    – DUO Labs
    Jul 20, 2019 at 14:22
2
$\begingroup$

We have $$\frac{33\pi}{x}=2\pi k,$$ where $k\in\mathbb Z$ and $$x\pi=2\pi n,$$ where $n\in\mathbb Z$.

We obtain: $$33=4kn,$$ which is impossible.

Also, there is a case $$\cos\frac{33\pi}{x}=\cos{\pi x}=-1.$$ Here we obtain: $$33=(1+2k)(1+2n).$$ Can you end it now?

$\endgroup$
0
1
$\begingroup$

Your task is equivalent to solving $\cos{\frac{33}{x}\pi}=\cos{x \pi}=1$ or $\cos{\frac{33}{x}\pi}=\cos{x \pi}=-1$.

The first equation results in $\frac{33}{x}\pi=2\pi n$, $x=\frac{33}{2n}$ -not a solution because $\cos \frac{33}{2n}\pi \ne 1, n \in Z, n \ne 0$.

The second equation results in $\frac{33}{x}\pi=\pi(1+2n)$, $x=\frac{33}{2n+1}$. Now we have $\cos \frac{33}{2n+1}\pi =-1$, or $\frac{33}{2n+1}\pi=(2k+1)\pi$, $k \in Z$. We can rewrite the last equation as $(2n+1)(2k+1)=33$ which gives us solutions $(3,11)$, $(-3,-11)$, $(1,33)$, $(-1,-33)$.

As for using numerical methods, there may be difficulty with using Newton method because the function and it's derivative have a lot of points of discontinuity and derivative may have a point of discontinuity where the value of function is zero.

$\endgroup$
1
  • $\begingroup$ Yeah I know that, but are there other methods I can use? $\endgroup$
    – DUO Labs
    Jul 19, 2019 at 21:30
1
$\begingroup$

For $\cos(x)\cos(y)$ to be equal to $1$, either both $\cos(x)$ and $\cos(y)$ must be equal to $1$ or both equal to $-1$. This is because the range of $\cos(x)$ is $[-1, 1]$. This means we want to solve $$\cos(x\pi) = 1, \cos\left(\frac{33}{x}\pi\right) = 1$$

and $$\cos(x\pi) = -1, \cos\left(\frac{33}{x}\pi\right) = -1$$

Tackling the first case first, for $\cos(t)$ to equal $1$, $t$ must be $2\pi k$, with $k$ an integer. This means $x = 2k_1$ is an integer and $x = \frac{33}{2k_2}$ is an integer. This cannot happen as $33$ has no even divisors.

For the second case, for $\cos(t)$ to equal $-1$, $t$ must equal $\pi + 2\pi k$. This means $x = 1+2k_1$ and $\frac{33}{x} = 1 + 2k_2$. For $33/x$ to be an integer, $x$ must be equal to $\pm 1, \pm 3, \pm 11, \pm 33$. $x$ and $33/x$ for all of these $x$ is odd.

Therefore, the solutions are $x = \pm 1, \pm 3, \pm 11, \pm 33$.

$\endgroup$
0
$\begingroup$

In a general manner, if you want to find the zero of $f(x)=0$ knowing that the solution is such that $a < x <b$, a good algorithm is used in subroutine $\color{red}{\text{rtsafe}}$ from Numerical Recipes (have a look here for the source code in C).

Basically, what it does is to combine bisection steps (whenever Newton method would make the iterate to be out of the given bounds - these are permanently updated) and Newton steps.

$\endgroup$
2
  • $\begingroup$ Bisection won't work because f(x)<=0. $\endgroup$
    – miracle173
    Jul 20, 2019 at 6:17
  • $\begingroup$ @miracle173. I was addressing the general problem of $f(x)=0$ not this specific one. $\endgroup$ Jul 20, 2019 at 7:53

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.