Efficient Techniques for Locating Roots of Complex Functions
Written on
Chapter 1: Introduction to Root Finding
In scientific programming, there are occasions when one must deal with a smooth, real-valued function f(t) that is costly to compute. Specifically, it may be necessary to determine the points where this function equals zero. I faced this challenge while developing an open-source astronomy library, aiming to calculate moonrise and moonset times. This process involves determining the Earth's rotation and the Moon's orbit to establish when the moon appears to rise or set, based on an observer's geographical position.
In this scenario, f(t) represents the Moon's apparent angular altitude above or below the horizon at time t. To find the moments of moonrise or moonset, it is essential to identify the value of t such that f(t)=0. Mathematically, this means locating a root of the function f.
Calculating f(t) involves determining the Moon's position, which requires extensive floating-point calculations. Thus, I aimed to minimize the frequency of evaluating f(t) while refining the value of t to achieve an accuracy of a few seconds.
Section 1.1: Initial Approach to Root Finding
One effective method for tackling this problem is to select a time t₀ before the root and another time t₁ after it. This approach assumes that the interval is small enough to contain only a single root. For instance, I know that there can’t be more than one moonrise within an hour.
If it is established that f(t₀)<0 and f(t₁)>0, this indicates that there is at least one root within the interval marked by the arrow. Given that f is continuous, at least one value of t exists between t₀ and t₁ for which f(t)=0. A continuous function cannot transition from a negative to a positive value without passing through zero.
To make further progress, we evaluate f at a midpoint t₂ within the interval: If f(t₂)<0, we have gained new information, narrowing the root’s location to the interval t₂<t<t₁, effectively halving the uncertainty. Conversely, if f(t₂)>0, the search can be constrained to the left half-interval t₀<t<t₂.
Continuing this process of halving the interval is known as binary search. This method can be repeated until the time interval is sufficiently small.
Subsection 1.1.1: The Role of Quadratic Interpolation
While binary search is a viable method, under certain conditions, quadratic interpolation can be significantly faster for locating a function's root. These conditions include:
- The function is smooth and continuous.
- The function can be closely approximated by a parabolic curve over a small interval t₀<t<t₁.
- There exists exactly one value of t between t₀ and t₁ where the function value equals zero.
This method is similar to binary search, but instead of evaluating the function at two points, we assess it at three equally spaced values. The unique quadratic function that fits these three points approximates the expensive function f over the defined interval.
For visualization, refer to the following video on quadratic interpolation:
Once the parabolic formula is derived, the quadratic formula can be applied to find the value of t where the parabola intersects the t-axis.
Section 1.2: Deriving the Quadratic Function
To derive the parabola, we first substitute the time variable t in the range t₀≤t≤t₁ with an abstract parameter x, which simplifies to the range −1≤x≤+1.
Verify that: - If t=t₀, then x=−1. - If t=t₂, then x=0. - If t=t₁, then x=+1.
Next, we define a generic parabolic function p(x) that must pass through the evaluated points for f. We establish three linear equations reflecting this requirement.
To solve for the coefficients Q, R, and S that ensure the parabola p(x) intersects the three evaluated points for f, we derive a system of linear equations.
Before proceeding, it is crucial to check for special cases. If Q=0, this indicates a linear function instead of a parabola, represented as p(x) = Rx + S. The root can be found as x = −S/R. If Q=0, we also need to check for R=0 to prevent division by zero—this scenario indicates a horizontal line without a zero-crossing.
For cases where Q≠0, we can utilize the quadratic formula to find the roots of the curve p(x).
To ensure accurate results, check if the radicand R²−4QS is negative. If R²−4QS < 0, no real roots exist for p(x), meaning the parabola does not cross zero, and the solver should report a failure.
You will need to decide how to handle cases where R²−4QS = 0. In my astronomy application, it was reasonable to disregard this as a valid moonrise or moonset.
Chapter 2: Python Implementation
By evaluating the expensive function only three times, assuming it behaves sufficiently parabolically within the interval, we can accurately approximate t such that f(t)=0.
However, the accuracy of t may vary. Experimentation with specific functions and interval sizes is essential.
A hybrid approach may also be beneficial: initially applying binary search to narrow the interval, followed by quadratic interpolation. My algorithm, found in the Python implementation of the Astronomy Engine, uses a hybrid method that intelligently guesses small intervals for quadratic interpolation and falls back to binary search when necessary. This method has proven to be as accurate as a pure binary search while significantly reducing computation time.
For practical implementation, here’s a sample of Python code that showcases using quadratic interpolation to approximate roots. You may use the function QuadraticInterpolate as a basis for your projects.