1-d and n-d Functions
Mathwrist White Paper Series

Mathematically modelling often uses piecewise functions for simplicity. FunctPiecewise is an intermediate class derived from Function1D base class to handle piecewise lookup. Concrete piecewise functions are further derived from this intermediate class with both eval() and integral() functions being implemented.

Given $n$ number of knot points $(x_0,\mathrm{\cdots },x_{n-1})$ in ascending order, these concrete piecewise functions are strictly defined within the range $[x_0,x_{n-1}]$. When numerically equal comparison is involved, we use to test if $a$ and $b$ are equal with numerical tolerance $ϵ$. When $a\le b$ comparison is involved, it should be interpretted as $a$ is numerically not greater than $b$.

Chebyshev curve is a 1-d function strictly defined within a fixed domain $x\in [a,b]$. It is a direct subtype of the Function1D class.

A Chebyshev curve of order $k$ is defined as $f(x)={\sum }_{i=0}^k{\theta }_i{T}_i(t(x))$, where $T_i()$ is the $i$-th order Chebyshev basis polynomial of the first kind. ${\theta }_i$ is the coefficient of $T_i()$. Chebyshev basis polynomial $T_i(t)$ is defined by the following recurrence relation over the canonical domain $t\in [-1,1]$,

Mathwrist NPL mapps the physical variable $x\in [a,b]$ to canonical variable $t\in [-1,1]$ of the Chebyshev basis functions by $t(x)=2\frac{x-a}{b-a}-1$. Once the range $[a,b]$ and the order $k$ of Chebyshev curve is determined, the function is fully parametrized by the coefficients ${\mathrm{\Theta }}^T=({\theta }_0,\mathrm{\cdots },{\theta }_k)$. Again, we can write $f(x)={\mathrm{\Theta }}^T𝐓(t(x))$, where $𝐓(t(x))$ is the vector of Chebyshev basis functionals. Our main references include [1], [6] and [7].

Class FunctPiecewise2D is an intermediate class, derived from the general FunctionND base class. The 2-d piecewise functions strictly work within a $m\times n$ rectangular grid, defined by $m$ knot points $(x_0,\mathrm{\cdots },x_{m-1})$ along $x$ axis and $n$ knot points $(y_0,\mathrm{\cdots },y_{n-1})$ along $y$ axis. In other words, its working domain is $[x_0,x_{m-1}]\times [y_0,y_{n-1}]$.

Chebyshev surface is a 2-d function defined strictly within a fixed domain $(x,y)\in [a,b]\times [c,d]$. It is a direct subtype of FunctionND class.

Similar to B-spline surface, Chebyshev surface is defined as the tensor product of two vectors of Chebyshev basis functionals (first kind).

, where $𝚯$ is the basis coefficient matrix. Function $s(x)$ mapps $x\in [a,b]$ to $s\in [-1,1]$, $t(y)$ mapps $y\in [c,d]$ to $t\in [-1,1]$.

For a unknown function $f(x)$ that is continuous over domain $[a,b]$, we observe a list of discrete data points $\{(x_k,f(x_k))\}$, $k=0,\mathrm{\cdots },n,x_k\in [a,b]$. One way to approximate this unknown function $f(x)$ is to construct an interpolation function that “exactly” matches all those data points $\{(x_k,f(x_k))\}$.

By Weierstrass Approximation theorem, we know that for any $ϵ>0$, there exists a polynomial $p(x)$ such that . It is then a natural choice to use polynomials as the approximation functions due to algebraic simplicity. It can be shown that there exists a unique polynomial $p(x)$ of degree at most $n$ such that $p(x_k)=f(x_k)$, $k=0,\mathrm{\cdots },n$. Further, this polynomial $p(x)$ can be directly constructed by

This representation form is called the $n$-th Lagrange interpolation polynomial. The interpolation polynomial $p(x)$ has an error term

, where $f^{(n+1)}(\xi (x))$ is the $(n+1)$-th derivative for some $\xi (x)\in (x_0,x_n)$.

It is possible to construct an interpolation polynomial to match both function values and derivatives. In particular, a Hermite polynomial $H(x)$ of at most $2n+1$ degree agrees with $f(x)$ and $f^{\prime }(x)$ at $(x_0,\mathrm{\cdots },x_n)$. $H(x)$ can be constructed as

The Hermite polynomial has an error term

The formulation of Lagrange polynomials, equations (7) and (8), is very convenient and useful in theoretical analysis. However, matching an additional data point by Lagrange or Hermite interpolation both involves increasing the degree of the interpolation polynomial. High degree polynomials with roughly equally spaced data points tend to have very large oscillating error when $x$ moves away from an interpolation data point. The error can go extremely wild near the end points, known as “Runge phenomenon”.

Instead of using a single high degree polynomial, cubic spline interpolation uses piecewise cubic polynomials that jointly form a smooth curve to approximate the unknown function $f(x)$. Let $S_i(x)$ be the cubic polynomial over the interval $(x_{i-1},x_i)$ for data point $(x_0,\mathrm{\cdots },x_n)$. The set of piecewise polynomials satisfy the following conditions,

There are a few ways to parameterize a cubic spline interpolation. Suppose we write $S_i(x)$ as

We have $n$ number of piecewise polynomials $S_i(x)$ hence $4n$ number of coefficients $(a_i,b_i,c_i,d_i)$ to solve. Conditions (14) contribute $4n-2$ equations. Two additional conditions are usually applied to the end points $x_0$ and $x_n$. Natural boundary conditions require $S_1^{\prime \prime }(x_0)=0$ and $S_n^{\prime \prime }(x_n)=0$. A clamped cubic spline has $S_1^{\prime }(x_0)=f^{\prime }(x_0)$ and $S_n^{\prime }(x_n)=f^{\prime }(x_n)$. All together, we have $4n$ equations to solve the $4n$ number of coefficients.

In Mathwrist NPL, we generalize the idea of cubic spline interpolation to be able to match one, two or all three quantities of $(f(x_k),f^{\prime }(x_k),f^{\prime \prime }(x_k))$ at a data point $x_k$. Effectively, users can control the level, slope and curvature of the spline curve at any data point, not just the boundary points in natural spline and clamped spline.

Let $p(x)$ be the cubic spline interpolation function. We choose a formulation of using B-spline basis of order 4, $p(x)={\sum }_{i=0}^{M-1}{\theta }_i{\varphi }_{i,4}(x)$, where $M$ is the total number of B-spline basis functions. ${\varphi }_{i,4}(x)$ is the $i$-th B-spline basis of order 4. ${\theta }_i$ is the coefficient of the $i$-th basis. Here $M$ is determined from the input number of unique data points $x_k$ and the values we want to match at $x_k$. As a minimum, we need at least $M=4$ input data point and value combinations. We then place $M-2$ spline knot points internally using an automatic knot-point placement scheme. Lastly, ${\theta }_i$ is computed by solving a $M\times M$ linear system.

In order to use the generalized cubic spline interpolation, users need create objects from the InterpCubicSplineCurve class, which is a subtype derived from base class Function1D. Next, we need prepare a buffer of interpolation Point objects, where the Point data type is defined as a nested struct in class InterpCubicSplineCurve,

A data point includes coordinate $x$ and the value from one of $f(x)$, $f^{\prime }(x)$ or $f^{\prime \prime }(x)$ to match. If an application needs to match multiple values at the same $x$, users will need populate multiple points with different value types. It is not necessary to sort data points by $x$. But the min and max of $x$ in supplied data points are interpreted as the domain range $[a,b]$ that $p(x)$ is defined. Once the cubic spline $p(x)$ is fit to interpolation data points, it can be used as a regular 1-d function to compute $p(x)$, $p^{\prime }(x)$ and $p^{\prime \prime }(x)$ for any $x\in [a,b]$. A function evaluation at $x\notin [a,b]$ will trigger an exception.

The code below illustrates an example of interpolating $f(x)=\sin (x)$ over $[0,\frac{\pi }{2}]$ with natural boundary condition $f^{\prime \prime }(0)=0$, slope control $f^{\prime }(\frac{\pi }{2})=0$ and three values $f(0)$, $f(\frac{\pi }{4})$ and $f(\frac{\pi }{2})$ to match.

The idea of numerically approximating a finite integral ${\int }_a^{b}f(x)𝑑x$ starts from approximating the integrand function $f(x)$ by an interpolation polynomial passing through a list of data points $\{x_0,\mathrm{\cdots },x_n\}$. This interpolation polynomial can be constructed from (7). Consequently,

This approximation is called a quadrature rule with quadrature weights

The degree of accuracy of a quadrature rule is measured by the highest integer $n$ such that the rule produces exact result for all $f(x)=x^k,k=0,\mathrm{\cdots },n$. By construction, the interpolation polynomial $p(x)$ has degree $n$. The error term (9) then is 0 for $f(x)=x^k,k\le n$. Therefore regardless how the data points are placed, the qudrature rule constructed from $n+1$ points $\{x_0.\mathrm{\cdots },x_n\}$ has at least $n$ degree of accuracy.

The family of Newton-Cotes quadrature formulas, i.e. the classic trapezoidal, Simpson’s rules, use equally spaced data points. Large oscillating error occurs when we increase the number of data points with a single interpolation polynomial. To avoid this, low degree polynomials are used to approximate $f(x)$ in small subintervals of $[a,b]$. Then we can sum up the quadrature results of subintervals and derive composite quadrature rules.

Romberg integration is a further improvement of composite quadrature rules. When equally spaced data points are used, it is possible to cancel out the lower order error term by Richardson’s extrapolation. Practically, Romberg integration can take accuracy order or error tolerance as input. The method starts iterations from low order approximation and terminates until the desired order of accuracy is reached, or the approximation error becomes smaller than the tolerance.

Unlike Newton-Cotes quadrature formulas, Gaussian quadrature chooses data points $\{x_0,\mathrm{\cdots },x_n\}$ in an optimal way such that the degree of accuracy becomes to $2n+1$. It can be shown that $2n+1$ is the highest accuracy that one can possibly reach with $n+1$ data points.

It turns out that the optimal data points $\{x_0,\mathrm{\cdots },x_n\}$ over interval $[-1,1]$ are the roots of the $(n+1)$-th Legendre polynomial. Gaussian quadrature weights computed from (15) are strictly positive. This is a very nice property for stable numerical integration. One way to prove the $2n+1$ degree of accuracy can be found in [1], section 4.7.

In Mathwrist NPL, we provide the Integrate1D class that implements an adaptive 3-point Gaussian quadrature method. Generally speaking, the idea of adaptive methods applies to equally spaced composite rules as well. When the integrand function $f(x)$ has dramatic curvature changes, the adaptive methods tend to use less number of data points over smooth area and increase the number of data points if curvature changes a lot. The details of adaptive quadrature can be found in [1], section 4.6.

A client application first needs define the integrand function $f(x)$ by inheriting the Function1D interface class. The quadrature rule only requires to evaluate function values $f(x)$, not derivatives. Next we need pass a constant reference of the integrand function to the quadrature constructor. Meanwhile, the constructor also accepts an integration precision parameter $ϵ$ that controls if the adaptive method needs refine integral intervals. Let $I(a,b)$ denote the quadrature result over an interval $[a,b]$.

At the $i$-th iteration of the adaptive method, assume we are working on sub-interval $[a_i,b_i]$, which is further splitted to two smaller subintervals $[a_l,b_l]$ and $[a_r,b_r]$. If

, then we consider the quadrature result over $[a_i,b_i]$ is accurate enough and doesn’t need to be refined with more data points. The code block below illustrates an example of numerically integrating ${\int }_0^1{x}^{-x}𝑑x$.

The root finding problem of a continuous univariate funciton $f(x)$ is to find the root $x^{*}$ such that $f(x^{*})=0$. Root finding methods can be classified into two categories. The first class methods only use the function values $f(x)$. Some well-known first class methods include bisection, brent, secant, etc. The second class methods use both $f(x)$ and its derivatives, typically only $f^{\prime }(x)$ in practice. The most powerful method in the second class probably is Newton-Raphson method.

Let $\{x_n,n=0,\mathrm{\cdots },\mathrm{\infty }\}$ be a sequence converging to root $x^{*}$. If there exist $\lambda >0$ and $\alpha >0$ such that

, then we say the sequence converges to $x^{*}$ at the order of $\alpha$ with asymptotic error constant $\lambda$. If $\alpha =1$ with , the sequence has linear convergence rate. If $\alpha =2$, the sequence has second order convergence rate.

Bisection method works within a bracketing interval, hence guarantees convergence. It has linear convergence rate. Brent method combines bisection with quadratic interpolation and has the advantage of both guaranteed convergence and superlinear convergence rate, . Secant method also has superlinear convergence rate but does not always converge. A secant method with false position correction (so-called “Regula Falsi”) can on-the-fly establishes a bracket interval hence ensures convergence. However, it involves more function evaluation and only has linear convergence rate.

For a smooth function $f(x)$, $f^{\prime }(x)\ne 0$, given an initial guess $x_0$, the Newton-Raphson method generates a sequnce of updates by

It can be shown that the Newton-Raphson method is equivalent to a fixed point iteration problem, i.e. [1] section 2.2, 2.3. By fixed point iteration convergence theorems, when $x_0$ is within a nearby area of $x^{*}$, the method converges at quadratic rate. On the other hand, when $x_0$ is far away from $x^{*}$, the method could diverge for some functions $f(x)$.

In Mathwrist NPL, we provide a root-finding solver class RootSolver1D that can either uses $f(x)$ only or both $f(x)$ and $f^{\prime }(x)$ depending on user control. This solver uses a safeguarded search strategy that combines bisection, rational interpolation and Newton method together to benefit from the advantage of each individual method.

The solver takes 3 parameters as input. Parameter $a$ and $b$ suggests a bracketing interval $[a,b]$. If this interval actually does not bracket, the solver may extend this interval until it brackets or terminates when the max number of iterations is reached. For the best performance of this algorithm, it is recommended to provide a decent guess of $[a,b]$ based on user’s domain knowledge.

A client application will need define the objective function $f(x)$ by inheritance from the base class Function1D. When $f^{\prime }(x)$ is not available or very expensive to compute, the client application may disable using $f^{\prime }(x)$ on the solver and then ignore the calculaiton of $f^{\prime }(x)$ and $f^{\prime \prime }(x)$ inside the Function1D::eval() implementation. By default, RootSolver1D enables the usage of $f^{\prime }(x)$ and expects the objective function to provide the calculation of $f^{\prime }(x)$.

The following example demonstrates a use case of solving $f(x)=x-2\sin (x)=0$.

The RootSolver1D class is derived from the IterativeMethod base class hence shares all the termination control implemented in IterativeMethod class. In particular, the IterativeMethod::converge_tolerance() method defines a tolerance quantity $ϵ$ for $f(x)$. When , the root finding solver terminates and returns $x_n$ as the solution.

For a univariate function $f(x)$, we want to solve the following minimization problem

, where $x_l$ and $x_u$ are the lower and upper bound of the domain that $f(x)$ is defined. Generally speaking, there could be multiple local optimals in $[x_l,x_u]$. Practically, we are only interested at finding a local optimal $x^{*}$ within a resonable range.

Given a bracketing interval $[a,b]$, bisection, Golden section and Fibonacci search use interval reduction strategy to locate a local minimum $x^{*}\in [a,b]$ until $a$ and $b$ converge. These methods guarantee convergence but have only linear converge rate.

Superlinear or higher order convergence methods can be constructed by using a quadratic model function $\widehat{f}(x)$ or cubic polynomials to approximate the objective function $f(x)$. For example, if we choose $\widehat{f}(x+p)$ to be the second order Taylor expansion at point $x$,

, when $f^{\prime \prime }(x)>0$, it suggests a step move $p=-\frac{f^{\prime }(x)}{f^{\prime \prime }(x)}$ to reach the minimum of $\widehat{f}(x)$. This is just the Newton step in the equivalent root finding problem $f^{\prime }(x)=0$. The same limitation and advantage of root-finding Newton method applies in the 1-d minimization context as well. Hence, it suggests a safeguarded approach similar to what we discussed in root finding problems.

The classic Brent method is one of such safeguarded algorithms that uses a quadratic form $\widehat{f}(x)$ to predict the minimum point. This minimum point is safeguarded by the bracketing interval and interval reduction strategy. It does careful bookkeeping on recent data points and use only function values of these points to fit the quadratic form. This method has superlinear convergence rate.

The MiniSolver1D class in NPL combines quadratic polynomial approximation $\widehat{f}(x)$,

, together with a Golden section method in a safeguarded strategy to solve the 1-d minimization problem (17). The minimization solver can be customized to use $f(x)$ only or both $f(x)$ and $f^{\prime }(x)$. Users can also set a control of max step length to avoid the step move being too large at each iteration.

The solver has an initial bracketing stage and a subsequent safeguarded searching stage. Three parameters are given to the solver as input. $x_0$ is the initial guess of $x^{*}$. $[a,b]\in (x_l,x_u)$ is a good guess of the bracketing interval. Regardless whether $f^{\prime }(x)$ is used, it is possible that we are able to move in a monotonically decreasing sequence all the way down to the boundary point $x_l$ or $x_u$. In this situation, the boundary point is returned as the solution by the end of the bracketing stage.

Users need supply an objective function $f(x)$ exactly the same way as the root-finding case in section 5.3. By default, the MiniSolver1D solver uses derivatives $f^{\prime }(x)$. This behavior can be explicitly enabled or disabled by calling the solver class member function void enable_derivative(bool use_df). Users can also put a control on the max step length that $x$ is allowed to move between iterations by calling void set_max_step(double step) function.

The following example demonstrates a use case of solving the minimum of function $f(x)=-(16x^2-24x+5)e^{-x}$.