# ODE Event detection for calculating multiple roots of continuous sinusoidal equation

I found a paper [1] that has a method for computing rise and set times of a satellite given a closed form solution. It is a complicated sinusoidal function and the paper has a method to calculate a single rise set over a single period using newton rhapson method. I want to find a way to calculate all zeros over a given interval of time.

Here is the equation: $$E=\alpha-\cos^{-1}\Bigg[\frac{G+\mathbf{P}\cdot\mathbf{\bar{Z}}ae}{a\big[(\mathbf{P}\cdot\mathbf{\bar{Z}})^2+(\mathbf{Q}\cdot\mathbf{\bar{Z}})^2(1-e^2)\big]^\frac{1}{2}}\Bigg]$$

Here is a picture of one of the initial intervals of the plot

I know Newton's approach will work given a good initial guess. My first thought was to create a loop and when it finds a zero to move forward to another set, but I find it hard to believe there is not a method already developed.

I have come across ODE event detection in MATLAB, or a lot of other solutions that seem to be for polynomials but wanted to see if there were better approaches.

1. RISE AND SET TIME OF A SATELLITE ABOUT AN OBLATE PLANET P. R. ESCOBAL AIAA Journal 1963 1:10, 2306-2310 DOI: 10.2514/3.2057
• Do you know the approximate interval of your roots? If yes, I would give you a simple bisection method that worked really well for me in finding the roots of Confluent Hypergeometric function. Commented Oct 4, 2019 at 20:20
• @AloneProgrammer : Proposing bisection to find roots is similar to proposing bubblesort as a simple and reliable method to sort long lists. Use better bracketing methods, like regula falsi (illinois, which is also very simple to program and reasonably fast), Dekker's fzeroin, Muller's inverse quadratic interpolation or, combining most of the features of the previous methods, Brent's method. Commented Oct 6, 2019 at 6:27
• @LutzL Method you mentioned have issues for convergence, accuracy, etc. But, bisection method is the only one that guarantees that if there is only one root in [a,b] range it will find it for 100% sure. I’ll be really grateful if you could bring an example that there is one root for a function in [a,b] range but bisection can’t find it. For finding multiple roots, you need to know the intervals of the roots to isolate each root in a given range and then use bisection to find it and for finding the all of them you need to just traverse the range and make sure there is one root in each interval. Commented Oct 6, 2019 at 13:18
• @AloneProgrammer : The mentioned methods are all bracketing methods, which means they all keep an ever shrinking interval containing a sign change (and thus root) of the function. This is the same guarantee that the bisection method gives. Under the same starting conditions, all these methods are faster than bisection at finding a root, except in very constructed cases where the speed reduces to about the bisection case. Commented Oct 6, 2019 at 14:18
• None of these methods seem to answer the key question of automating the process over continuous intervals with multiple roots. The process of root finding for single root over a defined interval is well studied, as you all have pointed out with many robust algorithms Commented Oct 7, 2019 at 16:49

Assuming it is not too expensive to evaluate the function, I would recommend using the chebfun toolbox (matlab), as shown here: https://www.chebfun.org/docs/guide/guide03.html

It builds an approximation of the function using Chebyshev polynomials and then finds the polynomial’s roots using a highly efficient root finder. It will find all roots in one fell swoop. Methods like Newton’s are less than ideal because it is hard to avoid converging to a root you have already found.

• the computation of the equation isn't expensive for small amounts, but the trade off is for zero finding i need tons of computations. I have been exploring these more and am very interested. Commented Oct 9, 2019 at 1:00

I don't know about a generalized approach, but satellite around Earth can not physically orbit around it in less than about 87 minutes; at an altitude of 100 km the period is about 5190 seconds, but of course so low it will burn up quickly!

That means that no matter where you are on Earth, those peaks can only reach maxima that are greater than zero (i.e. above the horizon) 16 times a day.

Therefore if you start your Newton root finder on an evenly spaced grid of time points spaced every 20 minutes or so, they will (probably) find every rise and set time.

There will of course be a high rate of redundancies, but this can be easily vectorized, then sorted, so I think it's going to be fine.