Algorithm 1048: A C++ Class for Robust Linear Barycentric Rational Interpolation

Barycentric rational interpolation (BRI) is an efficient and numerically robust interpolation method that performs well even in the case of equidistant nodes where polynomial interpolation can fail badly. Because of these favourable properties, BRI has recently become popular not only in mathematics, where it is used for the definition of quadrature rules [6, 22] and for solving partial differential equations [26, 27], but also in more applied contexts. For example, the method proposed by Floater and Hormann [13] is used in statistics [2], engineering [32], as well as in the medical [12], aerospace [21, 25], and chemical sciences [23, 24].

Given a set of \(n+1\) interpolation nodes \(X_{n}=(x_{0},x_{1},\dots,x_{n})\) with \(x_{0} \lt x_{1} \lt \dots \lt x_{n}\) and some associated data \(Y_{n}=(y_{0},y_{1},\dots,y_{n})\), any rational interpolant \(r\colon\mathbb{R}\to\mathbb{R}\) that interpolates \(y_{i}\) at \(x_{i}\) can be expressed in its second barycentric form asfor some non-zero weights \(W_{n}=(w_{0},w_{1},\dots,w_{n})\) [7]. We focus on linear BRI, meaning that the denominator does not depend on the data \(Y_{n}\), and we use weights \(w_{i}\) that guarantee the absence of poles in the domain \([x_{0},x_{n}]\). The first to introduce a method that is linear and prevents poles was Berrut [3], who defined the weights as

Afterwards, Floater and Hormann [13] proposed a new linear barycentric interpolant by considering a parameter \(d\in\{0,1,\dots,n\}\) and the barycentric weights

Actually, this method is a generalization of Berrut's case since, for \(d=0\), the weights in Equation (3) are exactly the ones in Equation (2). Furthermore, the polynomial interpolant [8] can be expressed as in Equation (1) with weightswhich can be also obtained from Equation (3) for \(d=n\).

The barycentric formula (1) is widely used to evaluate the interpolant \(r\) as it can be implemented with an efficient \(O(n)\) algorithm (see Algorithm 1) and the \(w_{i}\) can be rescaled by a common factor to prevent overflow and underflow errors [8]. Moreover, it has been proven [13] that this family of rational interpolants guarantees a fast convergence rate on the order of \(O(h^{d+1})\), where \(h\) denotes the maximum distance between consecutive nodes. However, there exists another mathematically equivalent formula to evaluate \(r\), namely the first barycentric formwhere

While this formula is slightly inferior in terms of efficiency because it requires \(O(nd)\) operations to be computed straightforwardly (see Algorithm 2), there exists a more efficient way of computing its denominator in \(O(n)\) operations. The resulting algorithm (see Algorithm 3) has the same computational cost as Equation (1), but it is slightly less numerically precise than the standard implementation [15]. Moreover, it may happen that Equation (5) is more numerically stable than Equation (1) for some configurations of nodes. In fact, it is known that the forward stability of both Equations (1) and (5) depends on the condition number \(\kappa(x;X_{n},Y_{n})\) related to the problem and on a second function that is different for both cases [15]. On the one hand, for the second barycentric form, the bound on the relative forward error is affected by the Lebesgue function \(\Lambda_{n}(x;X_{n})\), which is known to grow logarithmically in \(n\) and exponentially in \(d\) for equidistant nodes [10, 11]. On the other hand, the numerical stability of the first form depends on a function \(\Gamma_{d}(x;X_{n})\), which is independent of \(n\) (see Section 6), but the only proven upper bound so far depends on \(\mu^{d}\), where \(\mu\) is the mesh ratio of the nodes \(X_{n}\) [15]. However, numerical experiments indicate that \(\Gamma_{d}\) is usually considerably smaller than \(\Lambda_{n}\), suggesting that the first barycentric form may provide a more stable solution when the Lebesgue function related to \(X_{n}\) is large. Finally, the choice of the integer \(d\) is very important because, contrary to the approximation order, the numerical stability of the BRI may be compromised when \(d\) grows. For equidistant, quasi-equidistant [18], or conformally shifted nodes [5], this class of rational interpolants exhibits excellent approximation properties due to the favorable behavior of the Lebesgue function. In particular, for equidistant nodes, Güttel and Klein [16] provide insights on how to choose \(d\) with respect to \(n\) in order to preserve both the convergence rate and the stability. However, for cases like Chebyshev nodes, where polynomial interpolation is effective, these interpolants may become unstable under certain circumstances [4], so it is then better to use the polynomial interpolant in barycentric rational form by choosing \(d=n\).

The goal of this work is to present our class, named BRI, which contains all the necessary functions for handling a barycentric rational interpolant with all its features through a robust and efficient implementation. One main feature of the BRI class is that it provides the user with an “intelligent” function for evaluating \(r\), which computes the best result in terms of numerical stability, efficiency, or both by autonomously choosing one of the three evaluation options (Algorithms 1–3). The other main feature is the possibility to choose a “guarded” version of the evaluation algorithms, which avoids potential overflow and underflow errors by adaptively scaling the numerator and the denominator of \(r\). The functionality of our class exceeds those of currently available libraries [1, 9, 28, 30], which simply implement the second barycentric form in Equation (1) without any further considerations.

After a brief overview of all the functions contained in the BRI class (Section 2), we present some preliminary results needed to handle overflow and underflow errors in the guarded version (Section 3). Then, we discuss the implementations of the barycentric weights \(w_{i}\) in Equation (3) (Section 4) and of the interpolant \(r\) (Section 5). Finally, we describe the algorithms related to the computation of the functions \(\kappa\), \(\Lambda_{n}\), and \(\Gamma_{d}\) (Section 6).

The BRI class arises from the idea of providing the user with a class that holds both variables and functions related to BRI with a simple and intuitive interface and in an efficient programming environment. This idea results in a C++ class template, which comes with the advantage of having all input variables and function outputs defined as generic types. In addition to the basic data types available in C++, such as float or double, the BRI class supports arbitrary precision thanks to the compatibility with the Multiple Precision Floating-Point Reliable (MPFR) library [14] and the MPFR C++ interface [17]. The BRI class does not have any other dependencies other than the C++ standard libraries and, if arbitrary precision is used, the MPFR library.

To initialize a new instance of this class, the user has to specify the integer \(d\), the nodes \(X_{n}\), and the data \(Y_{n}\), and optionally their type and the type of the output results, which are both considered to be double by default. The vectors \(X_{n}\) and \(Y_{n}\) are stored internally as private vectors, but the class provides the user with some functions to access or modify them. Moreover, since they can be passed in different ways, the class includes different constructors. In particular, the user can pass both nodes and data by reference or they can be read from external files by specifying the filenames. Another constructor variant allows the user to pass one of the three keywords UNIFORM, CHEBYSHEV, or EXTENDED_CHEBYSHEV, together with the values \(a\), \(b\), and \(n\). In this case, the class automatically generates \(X_{n}\) as a vector of uniform nodesor Chebyshev nodesor extended Chebyshev nodesover the interval \([a,b]\). Of course, it would be possible to extend the class to allow for other predefined sets of interpolation nodes. As for the data \(Y_{n}\), the user can also call the constructor with a pointer to a procedure f, which implements some function \(f\colon\mathbb{R}\to\mathbb{R}\), so that the data is automatically generated as \(f(X_{n})\) by the class. After having defined these initial parameters, the class allows the user to modify them at any moment. Furthermore, the constructors create the vector \(W_{n}\) of barycentric weights (3), which get updated automatically whenever the nodes \(X_{n}\) are changed. Like nodes and data, the weights are private and can be accessed using suitable functions.

The core of the BRI class is the evaluation of the barycentric rational interpolant \(r\), defined by the input \(X_{n}\), \(Y_{n}\), and \(d\), through the eval function. The latter takes \(x\) as input, which can be either a single evaluation point or a vector of points, and outputs \(r(x)\), which in turn is a point or a vector, respectively. If the user wants to specify explicitly which algorithm to use for the evaluation, the eval function has an optional second parameter, which is one of the following keywords: SECOND to evaluate \(r\) with the second barycentric form (Algorithm 1), FIRST_DEF to use the standard implementation of the first barycentric form (Algorithm 2), or FIRST_EFF for its more efficient variant (Algorithm 3). Without this second parameter, the code uses by default the algorithm that best balances efficiency and numerical stability. Likewise, the user can call the functions numerator and denominator to evaluate the numerator and denominator of \(r\), respectively. In case of the denominator, one of the three aforementioned keywords can be added to explicitly ask for the denominator to be computed by the corresponding algorithms.

The BRI class also provides three flags: the stability flag for controlling the numerical stability of the result, the efficiency flag for indicating a preference for the fastest evaluation routine whenever the first barycentric form is used, and the guard flag for activating additional checks that prevent overflow and underflow. The user can turn these flags on or off at any moment. In the following sections, we present in detail how the code decides which algorithm guarantees the numerically most stable or most efficient result if the respective flag is activated and how the guarding mechanism avoids possible overflow or underflow errors without affecting the numerical accuracy of the result.

Finally, the class also allows to evaluate all the functions that are related to the numerical stability of the barycentric rational interpolant, namely the condition number \(\kappa\), the Lebesgue function \(\Lambda_{n}\), and the function \(\Gamma_{d}\). They can be called either with an input \(x\), which again can be an evaluation point or a vector of evaluation points, or without any arguments, upon which the algorithm searches for the maximum of the function using Newton's method.

We consider the binary floating-point number system [31], which is characterized by a precision \(t\in\mathbb{N}\) and two exponents \(E_{\smash{\min}},E_{\max}\in\mathbb{Z}\). For example, the standard float type in C++ represents a single precision floating-point number with \(t=23\), \(E_{\smash{\min}}=-125\), \(E_{\max}=128\), while the standard double type represents a double precision floating-point number with \(t=52\), \(E_{\smash{\min}}=-1021\), \(E_{\max}=1024\). Note that the IEEE standard 754 [20] specifies these numbers slightly differently in terms of the number of digits in the significand \(p=t+1\), the maximum exponent \(\mathit{emax}=E_{\max}-1\), and the minimum exponent \(\mathit{emin}=1-\mathit{emax}=E_{\smash{\min}}-1\). In such a system, every element of the set \(\mathbb{F}\) of floating point numbers is \(0\) or has the form \(\pm\mu\times 2^{E}\), with a mantissa (or significand) \(\mu\in[\mu_{\smash{\min}},\mu_{\max}]=[2^{-1},1-2^{-t}]\) and an exponent \(E\in\{E_{\smash{\min}},E_{\smash{\min}}+1,\dots,E_{\max}\}\), and the smallest and largest positive normal floating-point numbers are

There also exist so-called subnormal floating-point numbers, but they have reduced precision and are not relevant in our context as all the rescale operations are performed with respect to the interval \([F_{\smash{\min}},F_{\max}]\).

In the following proposition, we present some basic facts about two arbitrary positive floating-point numbers that will be used frequently later.

Let \(M,J,K\in\mathbb{N}\) and consider, for any \(i=0,1,\dots,M\), the numbers \(a_{i,j}=\alpha_{i,j}\times 2^{A_{i,j}}\in\mathbb{F}\) for \(j=1,2,\dots,J\) and \(b_{i,k}=\beta_{i,k}\times 2^{B_{i,k}}\in\mathbb{F}\) for \(k=1,2,\dots,K\). Denoting the products of these numbers bywhich can be expressed in floating point encoding asandsuppose we want to compute the quantity

Even if we know that \(f\in[F_{\smash{\min}},F_{\max}]\), it may happen that the algorithm that implements \(f\) runs into overflow or underflow errors in some of its intermediate steps. However, we can try to avoid this problem by appropriately rescaling the intermediate floating-point values by some constant \(C\), so that they are kept far from the overflow and underflow regions. Moreover, if we use a rescale factor of the type \(2^{E}\), for some \(E\in\mathbb{Z}\), then this operation modifies only the exponent of each floating-point number without changing the mantissa, meaning that we are not introducing any additional rounding error. Therefore, the goal now is to properly define the exponent \(E\) so that we are sure that \(\hat{f}=2^{E}f\) can be safely computed without any overflow or underflow error.

As for the barycentric weights in Equation (3), the simplest method to implement them is to follow their definition, thus obtaining an algorithm that computes all \(w_{i}\) in \(O(nd^{2})\) operations. However, Hormann and Schaefer [19] proposed a pyramid algorithm that requires only \(O(nd)\) operations and has the same precision [15]. This algorithm starts with the valuesand iteratively definesfor \(l=d-1,d-2,\dots,0\), tacitly assuming that \(v^{l}_{i}=0\) for \(i \lt 0\) and \(i \gt n-l\). The barycentric weights are then given as

To save memory, the BRI class uses the vector of weights also to store the intermediate values \(v^{l}_{i}\) in \(w_{i}\), which requires the assignment in Equation (21) to be performed in reverse order (see Algorithm 5).

After having initialized a new instance of the BRI class with the variables \(X_{n}\), \(Y_{n}\), and \(d\), the vector \(W_{n}\) is automatically computed as described above, which is efficient and robust. However, there are cases, albeit extreme, in which the Weights function runs into overflow or underflow errors.

It is precisely for situations like these that we provide the guard flag, which, if activated, calls code that tries to prevent overflow and underflow errors. The basic idea is to rescale intermediate floating-point values, so that they are kept far from the overflow and underflow regions. To do this, we multiply the values \(v_{i}^{l}\) in Equation (21) for \(i=0,1,\dots,n-l\) in each step of the pyramid algorithm with a suitable common constant \(2^{C_{l}}\), with \(C_{l}\in\mathbb{Z}\). Of course, we need to define these constants appropriately to make sure that these rescaling operations are safe and do not in turn cause any overflow or underflow errors. To this end, we shall first work out intervals \(I_{l}=\bigl{[}\mu_{\smash{\min}}2^{L_{l}},\mu_{\max}2^{U_{l}}\bigr{]}\) that are guaranteed to contain all \(v^{l}_{i}\). Before delving into these details, it is worth recalling that such rescaling operations do not change the second barycentric formula in Equation (1).

Considering the floating-point representation of each value \(v_{i}^{l+1}=\nu_{i}\times 2^{V_{i}}\) and the differences \(x_{i+l+1}-x_{i}=\xi_{i}\times 2^{X_{i}}\), for some \(l\in\{0,1,\dots,d-1\}\) and \(i=0,1,\dots,n-l\), we know from Proposition 4.4 that \(L_{l+1}\leq V_{i}\leq U_{l+1}\) and, from Equations (25) and (27), that \(H_{\smash{\min}}+E_{l}\leq X_{i}\leq H_{\max}+E_{l}+1\). Moreover, we can observe that the expression for the values \(v_{i}^{l}\) in Equation (21) fits into the more general formula Equation (14) for \(M=J=K=1\). This means that \(L_{l}\) and \(U_{l}\) in Equation (23) are exactly \(L\) and \(U\) in Equation (18), where \(R_{\smash{\min}}=L_{l+1}\), \(R_{\max}=U_{l+1}\), \(S_{\smash{\min}}=H_{\smash{\min}}+E_{l}\), and \(S_{\max}=H_{\max}+E_{l}+1\). Consequently, by Corollary 3.3, if \(U_{l}-L_{l} \lt E_{\max}-E_{\smash{\min}}\), then we can defineand compute \(2^{C_{l}}v_{i}^{l}\) for \(i=0,1,\dots,n-l\) without any overflow or underflow error using Algorithm 4.

Therefore, the “guarded” version of the pyramid algorithm (see Algorithm 6) rescales the values \(v_{i}^{l}\) by the constant \(2^{C_{l}}\) in each step \(l=d-1,d-2,\dots,0\). Therefore, if we first compute \(L_{l}\) and \(U_{l}\) as in Equation (23) and find that \(U_{l}-L_{l} \lt E_{\max}-E_{\smash{\min}}\), then we can be sure that computing the values \(2^{C_{l}}v^{l}_{i}\) with Algorithm 4 will not cause any overflow or underflow. If this latter condition is not satisfied for some \(\hat{l}\in\{0,1,\dots,d-1\}\), it means that we cannot be sure anymore that all \(v_{i}^{l}\) are normal floating-point numbers. Therefore, we try to repeat the same procedure by first re-defining

If it fails again, then the algorithm proceeds normally by updating the weights as in lines 3–7 of Algorithm 5. If the calculations are successful for every \(i=0,1,\dots,n-\hat{l}\) without overflow nor underflow errors, then we define \(L_{\hat{l}}\) and \(U_{\hat{l}}\) as in Equation (28), and we continue with the remaining iterations. Otherwise, the code stops and reports an error message to the user. At the end, the weights are rescaled with respect to the constant \(2^{C_{w}}\), where

For example, considering the setting of Example 4.1, we can solve the overflow problem with Algorithm 6, which returns the weights rescaled by the constant \(C_{w}=2^{-127}\).

As already mentioned in Section 2, the evaluation of the barycentric rational interpolant \(r\) can be realized by using the second rational barycentric form with Algorithm 1 or the first form with either Algorithm 2 or 3. The BRI class allows the user to choose one of the three algorithms by calling the eval function with two input parameters: the first is the evaluation point or vector \(x\) and the second is one of the three keywords FIRST_DEF, FIRST_EFF, or SECOND. However, the second parameter can be omitted, obtaining Algorithm 7 that autonomously chooses the best algorithm to use. Let us now focus on the latter case.

By default, the class is always called Algorithm 1 because it is generally the one that best combines numerical stability and efficiency of the method. On the other hand, if the user asks for the most stable solution by turning on the stability flag, then the code first evaluates \(\operatorname{\kappa}(x)\) and \(\Lambda_{n}(x)\) and decides, based on these values, which is the best algorithm to use. In particular, it is known that the condition number affects the upper bound on the relative forward error of both the first and the second barycentric form [15], so if \(\operatorname{\kappa}(x) \gt 10^{3}\), then the class warns the user about this and continues without interruption. It is further known that the stability of Equation (1) depends on the function \(\Lambda_{n}\), which grows with \(n\), whereas the stability of Equation (5) depends on the function \(\Gamma_{d}\), whose upper bound is independent of \(n\). However, numerical experiments show that the estimated upper bound on \(\Gamma_{d}\) seems to be always much larger than \(\Gamma_{d}\) itself. For this reason, if \(\Lambda_{n}(x) \gt 10^{2}\), then the result is computed with Algorithm 2. If it happens that also \(\Gamma_{d}(x) \gt 10^{2}\), then the user simply receives a warning without the code stopping. Finally, if both the stability and the efficiency flags are active and it turns out that the most stable solution is given by the first barycentric form, then it is computed with Algorithm 3, which computes all the \(\lambda_{i}\) with a more efficient iterative strategy. From a computational point of view, we are improving by switching from an \(O(nd)\) algorithm to an \(O(n)\) algorithm, therefore we know that for small \(d\) there is no big gain, but, as \(d\) increases, the efficient implementation can result in a lower execution time, without losing much in terms of stability. For example, Figure 1 compares the stability and the running time of the two algorithms considering \(n+1\) equidistant interpolation nodes \(x_{i}\in[0,1]\) for \(n=150\) with associated data \(y_{i}=0\), for \(i=0,\dots,n-1\) and \(y_{n}=1\). In particular, on the left we see the maximum of both relative forward errors evaluated at \(1{,}000\) random points in \([0,1]\) and on the right their execution times in seconds as a function of \(d\in\{1,2,\dots,n\}\). We observe that Algorithm 2 is slightly more precise than Algorithm 3, as expected, but the latter always wins in terms of efficiency, especially in the neighborhood of \(d=n/2\). It is worth noting that, if the stability flag is turned off, then the efficiency flag does not play any role, because Algorithm 7 takes the second barycentric form by default.

As for the weights, also in this case the user can decide to turn on the guard flag to prevent overflow and underflow errors. The basic idea is again the same as in Section 4 for the weights, that is, rescaling the numerator and the denominator properly in order to perform all floating-point operations far from the overflow and underflow regions and to avoid introducing any further rounding errors. In this case, we have to further take into account that the summations in the numerator and denominator of Equations (1) and (5) involve both subtractions and additions that may lead to cancellation errors. Moreover, from Proposition 3.1, preventing non-underflow subtractions could be tricky, especially if we are subtracting numbers of similar magnitudes. However, we recall that cancellation errors are already taken into consideration in the study of the numerical stability of the interpolant [15], so, if the method is stable, then it is proven that they cannot happen. Otherwise, one option would be to split each summation into one for the positive terms and one for the negative terms, which, in case the condition of Proposition 3.2 is satisfied, can be safely computed using Algorithm 4 and making only one final subtraction. The problem with this procedure is that, from our experiments, it seems to lose precision compared to the normal iterative algorithm when there are cancellation errors. Furthermore, since the two summations consist only of additions, if the addends have high orders of magnitude and \(n\) is large, then it is much more likely to run into overflow errors. For these reasons, we decided to safely compute only the addends of the summations in Equations (1) and (5) using Proposition 3.2 (for the special case \(M=0\)), and finally to calculate the sums with the classic iterative algorithm. Therefore, no adjustment is considered that could prevent overflow or underflow errors during the summations, but we check only afterwards if something like this has happened. In particular, considering the floating-point representation of each weight \(w_{i}=\omega_{i}\times 2^{W_{i}}\), data \(y_{i}=\upsilon_{i}\times 2^{Y_{i}}\), and difference \(x-x_{i}=\xi_{i}\times 2^{X_{i}}\), we letand we define

After that, denoting the common numerator of Equations (1) and (5), the denominator of Equation (1), and the denominator of Equation (5) by \(N\), \(D_{s}\), and \(D_{f}\), respectively, if \(U_{j}-L_{j} \lt E_{\max}-E_{\smash{\min}}\) for all \(j\in\{1,2,3\}\), then we can computewith Algorithm 4, but, since the exponents in Equation (30) do not consider additions and subtractions, at the end it is necessary to check whether overflow or underflow errors have happened. Finally, if every operation was computed safely, then we get the result with the second barycentric form asor with the first barycentric form as

Otherwise, the code stops and gives an error message to the user. Note that the guarded result of the first form involves also the rescale factor \(C_{w}\) in Equation (29), because its denominator does not depend on the weights. If \(U_{j}-L_{j}\geq E_{\max}-E_{\smash{\min}}\) for some \(j\in\{1,2,3\}\), then the algorithm tries to compute the corresponding summation as in Algorithm 1, 2, or 3. Again, if some overflow or underflow error happens, then the code stops, otherwise the corresponding rescale factor \(C_{j}\) is set to \(0\) for \(j\in\{1,2,3\}\) and the final result is computed as before.

Finally, we analysed how much the guard and the stability flags can affect the computational cost of the method. In particular, we conducted several tests to compare the running times with these flags both turned on and off and also considering variations in the input, such as \(d\), the number of nodes and data \(n\), and the number of evaluation points. Notably, if we activate only the guard flag, then we observed a maximum increase in the execution time of one order of magnitude. Instead, in the worst-case scenario where we enable both the guard and the stability flag, the increase can reach at most two orders of magnitude.

As mentioned in Section 5, the BRI class uses the functions \(\Lambda_{n}\), \(\Gamma_{d}\), and \(\operatorname{\kappa}\) internally to decide, if needed, which method should be used for the computation of \(r(x)\). Moreover, it also allows the user to evaluate them explicitly by using the functions leb, gamma, and cond with the evaluation point or vector \(x\) as input. However, for purposes of numerical stability, what really matters is the maximum of these functions in the range \([x_{0},x_{n}]\) and this is what these functions return when they have no argument as input. Let us see how they compute it.

First, we recall that \(\Lambda_{n}\), \(\Gamma_{d}\), and \(\operatorname{\kappa}\) can all be written in the common formwhere \(b_{i}(x)=w_{i}/(x-x_{i})\) for \(\Lambda_{n}\), \(b_{i}(x)=(-1)^{i}/\prod_{j=i}^{i+d}(x-x_{j})\) for \(\Gamma_{d}\), and \(b_{i}(x)=w_{i}y_{i}/(x-x_{i})\) for \(\operatorname{\kappa}\). One possible way to compute their maxima is to sample these functions very densely on the domain and search for the largest values, but, to be really accurate, we have to consider a big number of samples, thus losing efficiency. However, it is clear from the triangle inequality that a general function of the type (32) is greater than or equal to \(1\). Furthermore, in the specific case of \(\Lambda_{n}\), \(\Gamma_{d}\), and \(\operatorname{\kappa}\), the minimum with value 1 is assumed exactly at the nodes \(x_{i}\) and, from numerical experiments (see Figure 2), it appears that they are all concave in every sub-interval \([x_{i},x_{i+1}]\). This means that, considering these functions locally in every open sub-interval \((x_{i},x_{i+1})\), the point where the first derivative vanishes is a local maximum. Therefore, we apply Newton's method to find the root of \(\Lambda^{\prime}_{n}\), \(\Gamma^{\prime}_{d}\), and \(\operatorname{\kappa}^{\prime}\) in each sub-interval \((x_{i},x_{i+1})\) and we finally search for the maximum point among them. In general, in order to find a root \(t\in\mathbb{R}\) of a function \(g\colon\mathbb{R}\to\mathbb{R}\), Newton's method starts with some \(t_{0}\in\mathbb{R}\) and then generates a sequence of \(t_{1},t_{2},\dots\in\mathbb{R}\) by applying the iterative formula

Although this method usually converges quadratically, it may also fail, for example, if the initial guess \(t_{0}\) is too far from the correct solution. Consequently, the choice of the initial value \(t_{0}\) is an important step of Newton's method. Figure 2 shows that \(\Lambda_{n}\), \(\Gamma_{d}\), and \(\operatorname{\kappa}\) take their local maxima approximately in the midpoint of every sub-interval \([x_{i},x_{i+1}]\), so that we consider \(t_{0}=(x_{i}+x_{i+1})/2\) in each sub-interval \((x_{i},x_{i+1})\). After that, by Equation (32), it is easy to see that the method can be implemented straightforwardly by computing at each iteration \(j\) firstandand then using Equation (33) to set \(t_{j+1}=t_{j}-g^{\prime}(t_{j})/g^{\prime\prime}(t_{j})\), until the error \({\lvert t_{j+1}-t_{j}\rvert} \lt 10\epsilon\), where \(\epsilon\) is the machine epsilon.