Refining the brackets

Once bracketed, there are several classical ways of refining a root of a non-linear or non-analytical function. Among them, the most robust is the bissection (Press et al. (1992)). It always gives the correct answer if the root is correctly bracketed and if the function is continuous. However, it is not the quickest way in most of the situations. Like Herrmann (1994), we implemented an algorithm that mixes the bissection method and a Lagrange polynomial fit. The Lagrange polynomial is best constructed using the iterative Neville's algorithm (Press et al. (1992)). Corrections have been brought to Herrmann's algorithm to achieve better performances.

Figure 3.6: Method for refining roots. (a) to (f) are successive steps of refinement. The thick plain curve is the unknown theoretical curve. The thin plain line is the polynomial fit from already calculated samples (black dots). The grey dot is the new root computed from polynomial fit. The grey rectangles show the zoom area of the next step.

Practically, from the two initial bracketing values ($V_1$ and $V_2$ ), a third point $V_3$ is calculated by bissection. $R_1$ to $R_3$ are the values of $R_{1212}(z_0)$ for $V_1$ to $V_3$ , respectively. According to the sign of $R_3$ , either $V_1$ or $V_2$ is replaced by the value $V_3$ , swapping them eventually afterwards to keep $V_1<V_2$ . It is the state described by figure 3.6(a), where $V_1$ and $V_2$ are represented by black dots. A Lagrange polynomial is constructed on those two samples (a line in this case), shown by the thin plain line. From the intersection of the polynomial with axis $y=0$ , a new $V_3$ is deduced, it is shifted by a tenth of $V_2-V_1$ towards the limit with the highest $R_{1212}$ value, and $R_3$ is re-calculated. If $R_3$ is located between $R_1$ and $R_2$ , the function is bijective inside $[V_1, V_2]$ and hence invertible between $V_1$ and $V_2$ . This is not the case in figure 3.6(b). Consequently, the algorithm returns to bissection to generate a new sample $V_3$ from the current $V_1$ and $V_2$ . As in the first step, either $V_1$ or $V_2$ is replaced by the value $V_3$ (new brackets are shown in figure 3.6(c) by black dots). New samples (grey dots) are generated from the polynomial fit and integrated into it. In figures 3.6(d) to 3.6(f), the degree increases at each step as new samples are added.

To efficiently calculate the root of a Lagrange polynomial $P(V)$ , the axis X and Y are swapped during its construction. The coordinates of the samples are swapped so that $P'(R)$ fits $(R_1,V_1)$ , $(R_2,V_2)$ , .... The current estimate for the root is $V=P'(0)$ . If $R_1$ and $R_2$ differ from a factor 10 or more, the Neville's algorithm may fail and it is better to return to the bissection until reducing the ratio. To avoid a quick return to bissection, the newly generated point has to be on the side of the true root where either $R_1$ or $R_2$ is maximum. For doing so, the calculated value is considered as the true root, and it is shifted by a tenth of the current bracket interval towards the boundary having the highest value for the function. In most cases, when the true function is locally weakly non-linear, the limit with the highest function value is replaced by the new shifted estimate of the root at the next iteration, keeping the order of magnitude of $R_1$ and $R_2$ in the same range. This is the minor modification we brought to Herrmann's algorithm (1994), but it has a major influence over the global performances, as the full power of the polynomial fit is used. The maximum degree of the polynomial is set to 19, because its coefficients are stored in a static vector of limited length for efficiency. In a contradictory way, finding the perfect root ^3.3 is a problem because the inferior and superior brackets are lost as no sign can indicate on which side of the true root a new sample is. It may ruin the root search for the next modes. This case is checked, and the computation of the function is redone at a slightly distinct value (minus one tenth of the current bracketed interval).

With this method, only 4 to 6 iterations are generally necessary to obtain a $10^{-7}$ relative precision. One iteration corresponds to one evaluation of the numerical function $R_{1212}(z_0)$ . In the original code written by Herrmann (1994), the degree of the polynomial never increases over 2 or 3, quickly returning to bissection. More than 10 iterations are necessary to achieve the same precision. Our code rarely returns to bissection, increasing the degree of the polynomial at each iteration. Together with the removal of all file Input/Output, it has been possible to drop the time consumption by a factor 5 to 6.