High-accuracy approximation of piecewise smooth functions using the Truncation and Encode approach

Uncertainty Quantification (UQ) is the science that studies the quantification and the reduction of uncertainties in real applications with an intensive computational component. An example would be the computation of the fuel consumption of a car. Suppose we know how to compute the consumption as a function of some parameters (car and wind velocity, wheel condition, car weight, etc.), but we do not know the exact value of these parameters because they are random variables. Hence, the fuel consumption is also a random variable. UQ studies how to infer the model of the fuel consumption random variable by assuming models for the parameters. To carry out such study, many simulations of the consumption, changing the parameters among simulations, maybe required. These computations tend to be very expensive, because they may involve numerically solving systems of partial differential equations.

One of the aims of UQ is to find the statistical moments of the random variable that provides the solution to a given problem. In [1. 11], the authors propose a new multi-scale technique, the Truncation and Encode approach (TE henceforth) to be applied to lower the cost of the computation (by quadrature rules) of the moments of the solution of the following stochastic PDE

∂tu(x,t,ξ)+∂xf(x,t,ξ,u(x,t,ξ))=0,$$ \begin{equation*}% \label{eq:PDE} \partial_t u(x,t,\xi) + \partial_x f(x,t,\xi,u(x,t,\xi)) = 0, \end{equation*} $$

where ξ ∈ ℝ^N is a random vector, with known probabilistic distribution. Notice that for each ξ, a PDE must be solved to compute u_ξ(x, t) = u(x, t, ξ), namely

∂tuξ(x,t)+∂xfξ(x,t,uξ(x,t))=0.$$ \begin{equation*}%{eq:PDE_xi} \partial_t u_\xi(x,t) + \partial_x f_\xi(x,t,u_\xi(x,t)) = 0. \end{equation*} $$

Consequently, the calculation of these moments is a computationally intensive task, since function evaluations for given parameter values usually require a complex numerical simulation. Many numerical methods have been specially developed for this kind of UQ problems (for instance, methods generalizing or improving the classical Monte Carlo method, or those based on the Polynomial Chaos representation [14]) and all of them try to use as few function evaluations as possible.

The TE technique is specifically designed to meet a target accuracy in the computation of the desired integral, by ensuring a predetermined accuracy in the required function evaluations. It is specially suited for the integration of functions which are only piecewise smooth, since in this case other alternatives such as the Polynomial Chaos method may be highly inefficient.

In this paper we analyze a basic ingredient of the TE technique designed in [1, 11]: the influence of the approximation method for piecewise smooth functions in the global efficiency of the TE technique. For this, we consider a simplified framework in which we seek to compute integrals of the type

E=∫01f(ξ)dξ$$ \begin{equation*}\nonumber E=\int_0^1 f(\xi) d\xi \end{equation*} $$

by means of a quadrature rule using point values on an equally spaced grid. For the sake of simplicity, let us assume that we use the trapezoidal rule and a uniform mesh of grid-size h_K = 2^−K, i.e.

EK=2−K(12v0K+∑i=02KviK+12v2KK),viK=f(i2−K).$$ \begin{equation*}%{eq:integration_rule1} E_K = 2^{-K}(\frac12 v_0^K + \sum_{i=0}^{2^K} v_i^K + \frac12 v_{2^K}^K), \qquad v_i^K=f(i2^{-K}). \end{equation*} $$

By standard results, we know that when f is sufficiently smooth

|E−EK|=O(2−2K)$$ |E-E_K|= O(2^{-2K}) $$

while

|E−EK|=O(2−K)$$ |E- E_K| = O(2^{-K}) $$

when f is only piecewise smooth.

The TE algorithm seeks to replace as many as possible viK $v_i^K$ values, which require function evaluations of f, with modifiedv^iK $\hat v_i^K$ values, so that

E^K=2−K(12v^0K+∑i=02Kv^iK+12v^2KK)$$ \begin{equation*}%{eq:integration_rule2} \hat E_K = 2^{-K}(\frac12 \hat v_0^K + \sum_{i=0}^{2^K} \hat v_i^K + \frac12 \hat v_{2^K}^K) \end{equation*} $$

satisfies

|E−E^K|≤ε$$ |E-\hat E_K| \leq \epsilon $$

where ε is a user-dependent predetermined accuracy, even when f is only piecewise smooth.

Notice that

|EK−E^K|=2−K|(12v0K+∑i=12K−1v^iK+12v2KK)−(12v^0K+∑i=12K−1v^iK+12v^2KK)|≤‖vK−v^K‖1≤‖vK−v^K‖∞,$$ \begin{align*} |E_K - \hat E_K| &= 2^{-K}\left |(\frac12 v_0^K + \sum_{i=1}^{2^K-1} \hat v_i^K + \frac12 v_{2^K}^K) - (\frac12 \hat v_0^K + \sum_{i=1}^{2^K-1} \hat v_i^K + \frac12 \hat v_{2^K}^K)\right | \leq \|v^K - \hat v^K\|_1 \leq \|v^K - \hat v^K\|_\infty, \end{align*} $$

where

‖vK−v^K‖1=2−K∑i=02K|viK−v^iK|,‖vK−v^K‖∞=sup0≤i≤2K|viK−v^iK|.$$ \begin{equation*}\nonumber \|v^K - \hat v^K\|_1 = 2^{-K}\sum_{i=0}^{2^K} |v_i^K - \hat v_i^K|, \qquad \|v^K-\hat v^K\|_\infty = \sup_{0\leq i \leq 2^K} |v_i^K - \hat v_i^K|. \end{equation*} $$

Hence,

|E−E^K|≤|E−EK|+||v^K−vK||1.$$ \begin{equation} |E - \hat E_K| \leq |E - E_K| + ||\hat v^K - v^K ||_1. \end{equation} $$(1)

Thus, for a given (piecewise smooth) f, it is possible to ensure a target accuracy by choosing appropriately the mesh (i.e. K) and the target accuracy ε so that ||v^K−vK||∞≤ε $||\hat v^K - v^K ||_\infty \leq \epsilon$ .

The TE algorithm in [1, 11] is based on Harten’s Multiresolution Framework (MRF), which provides a set of tools to manage data in a multi-scale setting. In the last two decades, the MRF has been successfully applied in various contexts, and in particular in the design of adaptive schemes for the numerical solution of conservation laws and systems [6[7]–8, 13]. The TE algorithm follows a strategy similar to that used in the cost-effective schemes described in [7] and used in [6, 8, 13]. Both algorithms use the ideas in Harten’s MRF to compute a sequence v^K $\hat v^K$ such that ||v^K−vK||≤ε $||\hat v^K - v^K|| \leq \epsilon$ for a specified norm, with vK=(f(ξiK))i=02K $v^K=\left ( f(\xi_i^K)\right )_{i=0}^{2^K}$ .

These strategies proceed from coarse to fine resolution levels as follows: At each resolution level k ≤ K, (associated to a uniform grid with spacing 2^−k), a vector v^k $\hat v^k$ is obtained, whose components either coincide with or approximate the values of f on the corresponding grid, i.e. v^ik≈f(i2−k) $\hat v^k_i\approx f(i2^{-k})$ . Approximations are computed in [1, 7, 11] via an interpolatory reconstruction. The TE authors also consider a multilevel thresholding strategy (εk)k=0K $(\epsilon^k)_{k=0}^K$ , which should be properly chosen to guarantee certain precision in the output v^K $\hat v^K$ :

‖vK−v^K‖∞<≈ε.$$ \begin{equation*}%{eq:error_control0} \| v^K - \hat v^K\|_\infty \lessapprox \epsilon. \end{equation*} $$

In this work we show that the simple strategy ε^k = ε, also considered in [7], is more effective than the one proposed in [1, 11] for the current purpose, ε^k = 2^K−kε. Another advantage of taking ε^k = ε is that K can be taken as large as needed without increasing the number of function evaluations. Thus, |E − E_K| can be made arbitrary small in (1), which implies |E−E^K|≲ε $|E - \hat E_K|\lessapprox \epsilon$ .

In this work we also analyze the use of high accuracy interpolation techniques, and we establish the growth rate of the number of evaluations (n_eval) respect to the precision (‖vK−v^K‖∞) $(\|v^K - \hat v^K\|_\infty)$ :

neval=O(‖vK−v^K‖∞−1/s)≡‖vK−v^K‖∞=O(neval−s),$$ \begin{equation*}\nonumber n_{eval} = O(\|v^K - \hat v^K\|_\infty^{-1/s}) \quad \equiv \quad \|v^K - \hat v^K\|_\infty = O(n_{eval}^{-s}), \end{equation*} $$

where s is the approximation order of the interpolation technique. This result is independent of the smoothness of f. That is, the number of function evaluations increases slowly even if f has some discontinuities. Hence the TE algorithm is very convenient in such situations.

A comparison between linear and non-linear interpolation techniques is also carried out. Our numerical results will show that non-linear techniques are dispensable in the TE strategy. Since linear ones are faster, computationally speaking, they are preferable. In particular, we make the comparison between piecewise polynomial Lagrange interpolation of degrees 1 and 3, and the nonlinear PCHIP interpolation developed in [5]. In [1, 11], the ENO interpolation was chosen instead of PCHIP, but PCHIP interpolation is monotone, non-oscillatory, stable as a subdivision process and faster than the ENO technique.

The paper is organized as follows: In section 2 we recall Harten’s MRF and define the specific prediction schemes to be considered for our study. In section 3 we recall the Truncation and Encode strategy described in [1, 11] and in section 4 we perform a series of numerical experiments that confirm our observations. We close the paper with some conclusions and perspectives for future work.

In this section, we describe the TE method restricted to the approximation of functions.

Given a target error ε > 0, a piecewise smooth function f: [0, 1]→ ℝ and a suitable finest resolution level K, the goal is to obtain some v^K∈ℝnK+1 $\hat v^K\in\mathbb{R}^{n_K+1}$ such that

‖vK−v^K‖∞≲ε,vK=DKf,$$ \begin{equation} \|v^K-\hat v^K\|_\infty\lessapprox \epsilon, \qquad v^K = \mathscr{D}_K f, \end{equation} $$(8)

using as few evaluations of f as we can.

The algorithm computes recursively v^k $\hat v^k$ , from the coarsest level k = 0 to the finest one k = K. The construction of (v^k)k=0K $(\hat v^k)_{k=0}^K$ begins by evaluating f at ξ⁰ and ξ¹: v^0=D0f $\hat v^0 = \mathscr{D}_0 f$ and v^1=D1f $\hat v^1 = \mathscr{D}_1 f$ . The details are computed using (5):

d^j0=v^2j+11−(P01v^0)2j+1=dj0=v2j+11−(P01v0)2j+1.$$ \begin{equation*} \hat d^0_j = \hat v_{2j+1}^{1} - (P_{0}^{1}\hat v^{0})_{2j+1} = d^0_j = v_{2j+1}^{1} - (P_{0}^{1} v^{0})_{2j+1}. \end{equation*} $$

We compute v^k $\hat v^k$ , k > 1, iteratively: For each 0 ≤ i ≤ n_k, and for each 1 ≤ k ≤ K(K pre-fixed), we set v^2ik+1=v^ik $\hat v^{k+1}_{2i}=\hat v^k_i$ and

v^<!-- ^ -->2i+1k+1=(Pkk+1v^<!-- ^ -->k)2i+1,if |d^<!-- ^ -->jk−<!-- - -->1|<<!-- < -->ε<!-- ? -->kf(ξ<!-- ? -->2i+1k+1),if |d^<!-- ^ -->jk−<!-- - -->1|≥<!-- = -->ε<!-- ? -->k,∀<!-- ? -->i∈<!-- ? -->{2j,2j+1},$$ \begin{equation} \label{eq:recurrence} \hat v_{2i+1}^{k+1}= \begin{cases} (P_k^{k+1}\hat v^k)_{2i+1}, &\text{if }|\hat d^{k-1}_j|\lt\epsilon^k \\ f(\xi_{2i+1}^{k+1}), &\text{if }|\hat d^{k-1}_j|\geq\epsilon^k \end{cases}, \qquad \forall i\in\{2j,2j+1\}, \end{equation} $$(9)

where

d^jk−1=v^2j+1k−(Pk−1kv^k−1)2j+1.$$ \begin{equation} \hat d_j^{k-1}=\hat v_{2j+1}^{k} - (P_{k-1}^{k}\hat v^{k-1})_{2j+1}. \end{equation} $$(10)

Figure 1 shows a sketch of of recurrence (9). It has a simple interpretation: On one hand, if |d^jk−1|<εk $|\hat d^{k-1}_j|\lt\epsilon^k$ , then a 'small' prediction error indicates that Pk−1k $P_{k-1}^k$ approximates well the function at ξ2j+1k $\xi^k_{2j+1}$ ; thus Pkk+1 $P_k^{k+1}$ may also provide a good approximation at ξ2i+1k+1 $\xi^{k+1}_{2i+1}$ , i ∈ {2j, 2j + 1}, and we do not need to evaluate f at these points.

Fig. 1

On the other hand, |d^jk−1|≥εk $|\hat d^{k-1}_j|\geq\epsilon^k$ means that Pk−1k $P_{k-1}^k$ did not provide a sufficiently 'good' approximation to the function values, hence we prefer to evaluate f to obtain v^2i+1k+1 $\hat v^{k+1}_{2i+1}$ as the exact value of v2i+1k+1 $v^{k+1}_{2i+1}$ . It should be noted that if |d^jk−1|<εk $|\hat d^{k-1}_j|\lt\epsilon^k$ then d^ik=0 $\hat d^{k}_i=0$ , and if |d^jk−1|≥εk $|\hat d^{k-1}_j|\geq \epsilon^k$ then d^ik=f(ξ2i+1k+1)−(Pkk+1v^k)2i+1 $\hat d^k_i = f(\xi^{k+1}_{2i+1}) - (P_k^{k+1}\hat v^k)_{2i+1}$ , i ∈ {2j, 2j + 1}. Also, |d^jk−1|≥εk $|\hat d^{k-1}_j|\geq \epsilon^k$ implies that v4j+1k+1 $v^{k+1}_{4j+1}$ (i = 2j) and v4j+3k+1 $v^{k+1}_{4j+3}$ (i = 2j + 1) must be computed, so f is evaluated twice.

In the end of the recurrence, we obtain the vectors v^K $\hat v^K$ and (d^k)k=0K−1 $(\hat d^k)_{k=0}^{K-1}$ . Moreover, we know the MR decomposition of v^k $\hat v^k$ because of (10):

Mkv^k=(v^0;d^0;…;d^k−1)=(v0;d0;d^1;…;d^k−1).$$ \begin{equation*}\nonumber M^k \hat v^k=(\hat v^0;\hat d^0;\ldots;\hat d^{k-1})=(v^0;d^0;\hat d^1;\ldots;\hat d^{k-1}). \end{equation*} $$

In [1] the authors arrive to the conclusion that (8) if fulfilled if ε^k = 2^K−kε. However, we notice in the numerical experiments of [1, 11] that this thresholding strategy leads to

‖vK−v^K‖∞→K→∞0.$$ \begin{equation*} \| v^K - \hat v^K \|_\infty \overset{K\rightarrow\infty}{\longrightarrow}0. \end{equation*} $$

Thus, the real precision of v^K $\hat v^K$ depends on ε and K. Moreover, if a large K is needed, because a very fine mesh is demanded, then it turns out that ‖vK−v^K‖∞ $\| v^K - \hat v^K \|_\infty$ is unnecessarily small and the technique looses efficiency. We will see that the thresholding criterion ε^k = 2^K−kε proposed in [1, 11] can be relaxed to ε^k = ε, to ensure (8) with a precision that stays close to ε. Our strategy reduces the number of evaluations. In addition, the non-dependency of K allows us to take this value as large as necessary to ensure the desired accuracy in the numerical integration of f.

Another advantage of using ε^k = ε is that it is not mandatory to prefix the finest level K. So we can select another criterion to stop the recurrence, for example, that certain amount of evaluations has been reached. On the opposite, taking ε = 2^K−kε demands knowing K before the execution starts.

The purpose of the present section is to examine the efficiency of the TE algorithm under the various strategies considered in the previous sections. In particular, we show that the threshold strategy ε^k = ε, k = 0, leads to the desired target accuracy in a much more efficient manner than ε^k = 2^K−kε, as proposed in [1, 11].

In this paper, the efficiency is a measure of the number of evaluations required to reach a target error in the computation of v^K $\hat v^K$ . We show in this section that the efficiency of the TE strategy relies on the interpolation technique used in the prediction operator. In addition, we show numerically that the strategy ε^k = ε ensures

2−sε≲‖vk−v^k‖∞≲ε,$$ \begin{equation*} 2^{-s}\epsilon \lessapprox\|v^k - \hat v^k\|_\infty\lessapprox \epsilon, \end{equation*} $$

where s is the order of approximation of the local rule, but also we show that the number of evaluations increases slowly respect to the precision:

neval=O(‖vk−v^k‖∞−1/s)=O(ε−1/s).$$ \begin{equation*} n_{eval} = O(\|v^k - \hat v^k\|_\infty^{-1/s}) = O(\epsilon^{-1/s}). \end{equation*} $$

Hence, highly accurate interpolation techniques improves the efficiency.

On the other hand, we will see that non-linear interpolation techniques does not provide, in general, any advantage respect to the linear ones. Thus, linear prediction operators are preferable, because they are faster to compute.

The TE authors consider the following efficiency measures. The first one is

μcr=nK+1nw+n0+1,nw=#{|d^ik|≥ε},$$ \begin{equation*} \mu_{cr} = \frac{n_K +1}{n_w + n_0 + 1}, \qquad n_w=\#\{|\hat d^k_i| \geq \epsilon\}, \end{equation*} $$

where n_w is the number of detail coefficients greater than ε. The second one is

τ=nK+1neval,$$ \begin{equation*} \tau = \frac{n_K+1}{n_{eval}}, \end{equation*} $$

where n_eval is the total amount of function evaluations used to obtain v^K $\hat v^K$ . τ = 1 if the function was evaluated everywhere in ξ^K, and the larger τ is, the more vik $v_i^k$ were interpolated instead of evaluated.

In the TE method, each |d^ik−1|≥ε $|\hat d^{k-1}_i| \geq \epsilon$ implies two new function evaluations at the resolution level k + 1. So n_eval can be written as

neval = 2nw+n0.$$ \begin{equation} n_{eval} = 2 n_w + n_0. \end{equation} $$(11)

Using (11), τ can be expressed as

τ=nK+12nw+n0⇒μcr≈2τ,$$ \begin{equation*}\nonumber \tau = \frac{n_K+1}{2 n_w + n_0} \quad \Longrightarrow \quad \mu_{cr}\approx 2 \tau, \end{equation*} $$

as shown in Tables 1 and 2 of [1, 11]. In this paper, we will only consider τ.

Table 2

Parameters as specified in Table 1 for f₁, ε = ε^k = 10^-1 and the polygonal rule.

In Tables 1 to 2, we carry out the same experiment as in section 5.1 of [11]. Let us consider the parameters of Table 1 and the following functions, which are shown in Figure 2:

Fig. 2

f1(x)=sin⁡<!-- ? -->(2π<!-- p -->x2),f2(x)=sin⁡<!-- ? -->(2π<!-- p -->x2),x≤<!-- = -->11/20sin⁡<!-- ? -->(2π<!-- p -->x2)+1,x><!-- > -->11/20,x∈<!-- ? -->[0,1].$$ \begin{equation} \label{eq:f1_f2} f_1(x)=\sin(2\pi x^2), \qquad f_2(x)=\begin{cases} \sin(2\pi x^2), & x\leq 11/20 \\ \sin(2\pi x^2) + 1, & x\gt 11/20 \end{cases}, \qquad x\in[0,1]. \end{equation} $$(12)

In Tables 2 to 4 we apply the TE strategy with ε^k = 10⁻¹, k = 0, to the smooth f₁ function in Figure 2 and display the parameters of Table 1. We use the prediction operators, Pkk+1 $P_k^{k+1}$ described in Remark 2, i.e. the polygonal, cubic and PCHIP rules, while in [11] it was used the polygonal, cubic and ENO rules. As mentioned previously, we consider that the PCHIP prediction is more convenient because of its stability. As in the original papers, we set n₀ = 1, so ξ⁰ = {0, 1}. From these experiments we extract the following conclusions.

Table 3

Parameters as specified in Table 1 for f₁, ε = ε^k = 10^-1 and the cubic rule.

Table 4

Parameters as specified in Table 1 for f₁, ε = ε^k = 10^-1 and the PCHIP rule.

As observed in Table 2, the polygonal rule uses 13 evaluations for K = 3, but 17 for K = 4. The cubic and PCHIP rules only uses 13 evaluation for any K = 3, as shown in Tables 3 and 4. From this fact, we deduce that there is a level k₀ such that for all k = k₀ the TE algorithm does not evaluate f₁ anymore. Actually, this property may hold for a general smooth function and a given threshold ε > 0.

Notice also that eK=‖vK−v^K‖∞ $e_K = \|v^K-\hat v^K\|_\infty$ is always smaller than ε, but it stays close to ε, a clear advantage of using ε^k = ε. Indeed, it can be observed in Tables 1 and 2 of [1, 11] that e_K tends to zero when K increases, when the strategy ε^k = 2^K−kε was selected. Since our target accuracy is e_K⪅ε with as few evaluations as possible, if the TE strategy uses just the necessary evaluations, then e_K should be close to ε. Thus, the fact that e_K tends to zero with the criterion ε^k = 2^K−kε means that the function is being evaluated more than it is actually needed.

As exposed in Section 1, the next condition is fulfilled:

|E−E^K|≤|E−EK|+‖vK−v^K‖1≤|E−EK|+‖vK−v^K‖∞.$$ \begin{equation} |E-\hat E_K| \leq |E - E_K| + \|v^K - \hat v^K \|_1 \leq |E - E_K| + \|v^K - \hat v^K \|_\infty. \end{equation} $$(13)

Since f₁ is smooth, then |E − E_K| = O (2^−2K). Thus, for large values of K, (13) becomes

|E−E^K|≲‖vK−v^K‖1≤‖vK−v^K‖∞≲ε,$$ \begin{equation} |E-\hat E_K| \lessapprox \|v^K - \hat v^K \|_1 \leq \|v^K - \hat v^K \|_\infty \lessapprox \epsilon, \end{equation} $$(14)

which can also be checked in Tables 2 to 4 for any K > 3.

The smallest e_K is achieved with the cubic rule, and the largest with the polygonal rule. Between them is the PCHIP rule, with less e_K and less evaluations than the polygonal rule, but larger e_K with the same evaluations than the cubic rule. Hence the cubic rule here seems to be better than PCHIP rule for smooth functions, which is probably a consequence of the better approximation properties of the cubic rule at all smooth regions, while the PCHIP formula its as accurate as the cubic rule only at monotone regions.

In Tables 5 to 7 we repeat the same experiment for the function f₂, which has a discontinuity jump in ξ = 11/20. Once again, we use ε^k = ε = 10⁻¹, and we explore the effect of increasing K.

Table 5

Parameters as specified in Table 1 for f₂, ε = ε^k = 10^-1 and the polygonal rule.

Table 6

Parameters as specified in Table 1 for f₂, ε = ε^k = 10^-1 and the cubic rule.

Table 7

Parameters as specified in Table 1 for f₂, ε = ε^k = 10^-1and the PCHIP rule

For K > 4 in Table 5 and K > 3 in Tables 6 and 7, the TE algorithm evaluates f₂ only two more times in each new level. This happens because there is only one detail coefficient satisfying |d^jk−1|≥ε $|\hat d_j^{k-1}|\geq \epsilon$ , which implies two new evaluations at level k + 1. We refer to Figure 1 for the sake of clarity. This djk−1 $d_j^{k-1}$ is located around the discontinuity, which is easily seen in Fig. 5 of [11].

In spite of the discontinuity of f₂, ‖vK−v^K‖∞ $\|v^K-\hat v^K\|_\infty$ is smaller than ε, and also is near ε, which supports the criterion ε^k = ε. Now |E − E_K| = O (2^−K), thus (14) is satisfied for K > 5.

Note that the polygonal rule uses more f₂ evaluations, but also its e_K are smaller than the other rules. e_K has similar magnitudes in the PCHIP and cubic rules and both use the same number of evaluations.

Once again, e_K converges to some value close to ε for K → ∞, in contrast to observed in Tables 1 and 2 of [1, 11], where e_K → 0, which confirms that the thresholding strategy ε^k = ε, instead of ε^k = 2^K−kε, is the most appropriate to attain a predetermined target accuracy in the computation.

In order to validate the error control property

2−sε≲‖vK−v^K‖∞≲ε,$$ \begin{equation} 2^{-s}\epsilon \lessapprox\|v^K - \hat v^K\|_\infty\lessapprox \epsilon, \end{equation} $$(15)

we repeat the same experiment for a fixed value of K, K = 15, and several ε values. We do not consider τ in Tables 8 to 13, because if K is fixed, then τ=O(neval−1) $\tau = O(n_{eval}^{-1})$ . In the following, we study the approximation of the smooth function f₁, which is shown in Tables 8 to 10.

Table 8

Parameters as specified in Table 1 for f₁, K = 15 and the polygonal rule. ε^k = ε varies between 10^-1 and 10^-10.

Table 9

Parameters as specified in Table 1 for f₁, K = 15 and the polygonal rule. ε^k = ε varies between 10^-1 and 10^-10.

Table 10

Parameters as specified in Table 1 for f₁, K = 15 and the PCHIP rule. ε^k = ε varies between 10^-1 and 10^-10.

Table 11

Parameters as specified in Table 1 for f₂, K = 15 and the polygonal rule. ε^k = ε varies between 10^-1 and 10^-10.

Table 12

Parameters as specified in Table 1 for f₂, K = 15 and the polygonal rule. ε^k = ε varies between 10^-1 and 10^-10.

Table 13

Parameters as specified in Table 1 for f₂, K = 15 and the polygonal rule. ε^k = ε varies between 10^-1 and 10^-10.

The first feature we observe is that (15) is satisfied, which provides a strong control on the error eK=‖vK−v^K‖∞ $e_K = \|v^K - \hat v^K\|_\infty$ . We also notice that the growth rate of the number of evaluations, n_eval, is larger in the polygonal rule. We can easily see it plotting the results of the tables, shown in Figure 3 and 4. From such figures is deduced that the growth rate is neval=O(eK−1/2) $n_{eval}=O(e_K^{-1/2})$ for the polygonal rule and neval=O(eK−1/4) $n_{eval}=O(e_K^{-1/4})$ for the cubic and PCHIP rules. It can be equivalently formulated as e_K = O ((n_eval)⁻²) and e_K = O ((n_eval)⁻⁴). Remember that 2 and 4 are the highest orders of approximation of the considered local rules: PCHIP has only second order of approximation, in general, but forth order in monotone regions. Despite this, PCHIP verifies e_K = O ((n_eval)⁻⁴) because, in the TE algorithm, the rule is applied in small monotone regions along the whole interval [0,1], where PCHIP has order 4. In addition, we clearly observe that cubic is a bit better than PCHIP, but much better than polygonal.

Fig. 3

Fig. 4

Another issue we can observe in Tables 8 to 10 is that |E − Ê_K| decrees together with ε, if ε≥10⁻⁸, but for smaller ε, the precision of the integral gets stuck around 10⁻¹¹. It can be explained using the triangle inequality:

|E−E^K|≥|E−EK|−|EK−E^K|≥|E−EK|−‖vK−vK‖1≥|E−EK|−ε.$$ \begin{equation*} |E - \hat E_K| \geq |E- E_K| - |E_K - \hat E_K| \geq |E - E_K| - \|v^K - v^K\|_1 \geq |E - E_K| - \epsilon. \end{equation*} $$

Since f₁ is an smooth function, then

|E−E15|=2−3012|f″(η)|≅7.8e-11|f″(η)|,η∈[0,1].$$ \begin{equation} |E-E_{15}| = \frac{2^{-30}}{12} |f''(\eta)| \cong \text{7.8e-11} |f''(\eta)|, \qquad \eta\in[0,1]. \end{equation} $$(16)

Thus,

|E−E^K|≥7.8e-11|f″(η)|−ε.$$ \begin{equation*} |E - \hat E_K| \geq \text{7.8e-11} |f''(\eta)|- \epsilon. \end{equation*} $$

This says that, if ε → 0, then |E − Ê_K| is lower bounded by |f″(η)|, which should be close to the number 1.5e-11, shown in Tables 8 to 11.

Comparing Figure 4 in this paper with Fig. 4 of [11], we note that our plot displays less oscillations. This may indicate that the dependency between e_K and the number of evaluations is stronger and smoother when the criterion ε^k = ε is selected.

We repeat the experiment for the function f₂, which has a discontinuity jump in ξ = 11/20, in Tables 5 to 7.

Notice that (15) is satisfied and the growth of n_eval is faster in the polygonal rule. Observing the Figures 5 and 6, corresponding to these tables, we observe that the growth rates are neval=O(eK−1/2) $n_{eval}=O(e_K^{-1/2})$ for polygonal rule and neval=O(eK−1/4) $n_{eval}=O(e_K^{-1/4})$ for cubic and PCHIP rules. As shown in Figure 5, cubic rule is better than PCHIP rule to approximate v^K. However, PCHIP rule is now better than cubic rule to compute the integral. At least, it is true for ε > 10⁻⁴, before |E − Ê_K| gets stuck.

Fig. 5

Fig. 6

A similar argument to (16) can be used to explain that |E − Ê_K| → 1.1e-06 in Tables 5 to 7. Since f₂ is not continuous, the error of the trapezoidal rule decreases as O (2^−K). Thus

|E−E^K|≥O(2−K)−ε.$$ \begin{equation*}\nonumber |E - \hat E_K| \geq O(2^{-K}) - \epsilon. \end{equation*} $$

Again, a comparison between our Figure 4 and Fig. 6 of [11] seems to indicate that the strategy ε^k = ε makes the TE approach more robust than ε^k = 2^K−kε.

We conclude that PCHIP only provides a better performance respect to the usual cubic interpolation when the integral of a discontinuous function is computed. However, the cubic rule was better in all the other cases: integration of a continuous functions and approximation of continuous or even discontinuous functions. The TE authors also arrived to the conclusion that no advantages were seen using the ENO interpolation for f₁, but ENO was better for f₂ in terms of accuracy and amount of evaluations.

The results included in the present paper indicate that non-linear interpolatory reconstructions may not lead to substantial advantages over linear ones within the TE framework. We have observed that, in general, the cubic interpolation lead to better results than the PCHIP interpolation, while both always fulfill neval=O(eK−1/4) $n_{eval}=O(e_K^{-1/4})$ . On the other hand, the usual reason of using ENO or PCHIP interpolatory techniques is the capability of approximating jumps in the data without spurious oscillations. But in TE strategy this is not needed, because the algorithm itself has a mechanism to detect the problematic zones of the function. In addition, linear interpolations are always faster than nonlinear ones.

We wish to remark again that the choice ε^k = ε confers a high control over the error. In contrast, the original selection ε^k = 2^K−kε done by the TE authors only guarantees that ε is an upper bound of the error. Observe that in [1, 11] all the numerical experiments were done with ε = 10⁻¹, but the error converges to zero as soon as K → ∞. In our case, the error converges to some number under ε, but greater than 2^−sε, where s is the order of the interpolation technique. This is a practical advantage, because K determines the resolution of the approximation v^K $\hat v^K$ , and we can set it as large as we need, but the number of evaluations and the error depends exclusively on ε.

In this paper we have examined some features of the Truncation and Encode strategy described in [1, 11], in particular its ability to compute integrals by quadrature rules to a given target accuracy.

In the TE strategy, the values used in the quadrature rule are computed by combining function evaluations and suitably interpolated values in a multilevel fashion, following the basic guidelines of Harten’s MRF.

Our theoretical observations and numerical experiments allow us to conclude that the TE technique becomes more efficient if the thresholding strategy proposed in [1, 11] is relaxed to ε^k = ε for all resolution levels. In addition, higher order prediction schemes lead to improved the performance, but it does not seem necessary to resort to nonlinear prediction schemes.

More research should be done to support these conclusions at a theoretical level.

The authors acknowledge support from Project MTM2014-54388 (MINECO, Spain) and the FPU14/02216 grant (MECD, Spain).