Fast numerical derivatives based on multi-interval Fourier extension

TL;DR

This paper introduces a fast, stable numerical differentiation algorithm that combines multi-interval Fourier extension with adaptive domain partitioning, effectively handling noisy data and complex functions.

Contribution

The paper proposes a novel method that integrates multi-interval Fourier extension with adaptive partitioning for efficient and stable numerical differentiation from noisy data.

Findings

Significant accuracy improvements over existing methods.

Robustness to noise and function oscillations.

Efficient computation through precomputed matrices and SVDs.

Abstract

We present a computationally efficient algorithm for stable numerical differentiation from noisy, uniformly-sampled data on a bounded interval. The method combines multi-interval Fourier extension approximations with an adaptive domain partitioning strategy: a global precomputation of local Fourier sampling matrices and their thin SVDs is reused throughout a recursive bisection procedure that selects locally-resolved Fourier fits. Each accepted subinterval stores a compact set of Fourier coefficients that are subsequently used to reconstruct the derivative via a precomputed differentiation operator. The stopping criterion balances fitting error and an explicit noise-level bound, and the algorithm automatically refines the partition where the function exhibits rapid oscillations or boundary activity. Numerical experiments demonstrate significant improvements over existing methods, achieving accurate derivative reconstruction for challenging functions. The approach provides a robust framework for ill-posed differentiation problems while maintaining computational efficiency.

Equations

$$\Phi (\mathbf{v})=\Vert \mathcal{F}\mathbf{v}-f^{\delta }{\Vert }_{L^2(-1,1)}^2+\alpha \Vert \mathcal{R}\mathbf{v}{\Vert }_{l^2}^2.$$

$$\mathcal{R}\mathbf{v}=(c_0,e^{\frac{\pi }{2}}{c}_1,e^{\frac{\pi }{2}}{s}_1,\dots ,e^{\frac{n\pi }{2}}{c}_n,e^{\frac{n\pi }{2}}{s}_n,\dots {)}^T,$$

$$\Vert f-f^{\delta }{\Vert }_{L^2(-1,1)}\le \delta ,$$

$$\Vert D^k(f^{\delta }-f){\Vert }_{L^2(-1,1)}=\mathcal{O}({\delta }^{\frac{p-k}{p}})$$

$$\Phi =\{{\varphi }_{\ell }{\}}_{\ell \in \mathbb{Z}},{\varphi }_0(t)=\frac{1}{\sqrt{2}},{\varphi }_{\ell }(t)=e^{\mathrm{i}\ell t},\ell =\pm 1,\pm 2,\dots .$$

$$\mathcal{F}:{\ell }^2\to L^2(\Lambda ),c=\{{\hat{c}}_{\ell }{\}}_{\ell =-\infty }^{\infty }\mapsto \sum _{\ell =-\infty }^{\infty }{\hat{c}}_{\ell }{\varphi }_{\ell },$$

$${\mathcal{P}}_N:{\ell }^2\to {\ell }^2,c\mapsto (\dots ,0,{\hat{c}}_{-N},\dots ,{\hat{c}}_N,0,\dots {)}^{\top },$$

$$B:\mathcal{D}(B)\subset {\ell }^2\to {\ell }^2,c\mapsto {\left\{\sqrt{1+|\ell {|}^{2p}}{\hat{c}}_{\ell }\right\}}_{\ell =-\infty }^{\infty },$$

$$W:\mathcal{D}(W)\subset {\ell }^2\to {\ell }^2,c\mapsto {\left\{e^{|\ell |}{\hat{c}}_{\ell }\right\}}_{\ell =-\infty }^{\infty }.$$

$$\mathcal{F}c=g,$$

$$\mathcal{F}c=g^{\delta },$$

$$\Vert g-g^{\delta }{\Vert }_{L^2(\Lambda )}\le \delta .$$

$$\Vert \mathcal{F}{c}^{\star }-g{\Vert }_{L^2(\Lambda )}=\underset{c\in {\ell }^2}{\inf }\Vert \mathcal{F}c-g{\Vert }_{L^2(\Lambda )}.$$

$$\Vert c_B^{†}{\Vert }_{{\ell }^2}=\inf \left\{\Vert B{c}^{\star }{\Vert }_{{\ell }^2}|c^{\star }\text{ is a least-squares solution of eqrefcompactoperatoreq}\right\}.$$

$${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_k\ge \cdots >0.$$

$${\mathcal{T}}_{k}g=\sum _{i=1}^k\frac{1}{{\sigma }_i}⟨g,u_i⟩v_i.$$

$$\Vert {\mathcal{T}}_k\Vert \le \frac{1}{{\sigma }_k}.$$

$$\Vert I-G{\mathcal{T}}_k\Vert \le 1,\Vert (I-G{\mathcal{T}}_{k-1})G\Vert \le {\sigma }_k.$$

$$g=\sum _{i=1}^{\infty }⟨g,u_i⟩u_i.$$

$$(I-G{\mathcal{T}}_{k-1})g=\sum _{i=k}^{\infty }⟨g,u_i⟩u_i,$$

$$\Vert (I-G{\mathcal{T}}_{k-1})g{\Vert }^2=\sum _{i=k}^{\infty }|⟨g,u_i⟩{|}^2\le \Vert g{\Vert }^2.$$

$$\Vert I-G{\mathcal{T}}_{k-1}\Vert \le 1.$$

$$c=\sum _{i=1}^{\infty }⟨c,v_i⟩v_i,\text{so that}Gc=\sum _{i=1}^{\infty }{\sigma }_i⟨c,v_i⟩u_i.$$

$$(I-G{\mathcal{T}}_{k-1})Gc=\sum _{i=k}^{\infty }{\sigma }_i⟨c,v_i⟩u_i,$$

$$\Vert (I-G{\mathcal{T}}_{k-1})Gc{\Vert }^2=\sum _{i=k}^{\infty }{\sigma }_i^2|⟨c,v_i⟩{|}^2\le {\sigma }_k^2\sum _{i=k}^{\infty }|⟨c,v_i⟩{|}^2\le {\sigma }_k^2\Vert c{\Vert }^2.$$

$$\Vert (I-G{\mathcal{T}}_{k-1})G\Vert \le {\sigma }_k.$$

$$c^{\delta ,\eta }=W^{-1}{\mathcal{T}}_{\eta }{g}^{\delta },$$

$$\Vert \mathcal{F}{c}^{\delta ,\eta }-g^{\delta }{\Vert }_{L^2(\Lambda )}\le \tau \delta <\Vert \mathcal{F}{c}^{\delta ,\eta -1}-g^{\delta }{\Vert }_{L^2(\Lambda )},\tau >1.$$

$${\mathbf{g}}_{\eta }^{\delta }={\mathcal{Q}}_{\eta }^{\delta }{g}^{\delta }:=\mathcal{F}{c}^{\delta ,\eta }$$

$$c_{B,N}={\mathcal{P}}_N{c}_B^{†},c^{\eta }=W^{-1}{\mathcal{T}}_{\eta }g,c_N^{\eta }=W^{-1}{\mathcal{T}}_{\eta }\mathcal{F}{c}_{B,N}.$$

$$\Vert B{c}_B^{†}{\Vert }_{{\ell }^2}\le E.$$

$$\Vert c_B^{†}-c_{B,N}{\Vert }_{{\ell }^2}\le N^{-p}E,\Vert W{c}_{B,N}{\Vert }_{{\ell }^2}\le C_{N}E,$$

Full text

1 Introduction

Numerical differentiation, the process of approximating derivatives from discrete and often noisy function values, constitutes a fundamental yet inherently ill-posed problem in scientific computing. Its significance stems from pervasive applications across diverse fields, such as heat flux estimation[3], volatility modeling in finance [9, 32], solving some differential and integral equations [12, 17, 10], and discontinuity reconstruction and edge detection[26, 25]. The core challenge lies in the intrinsic instability: minute perturbations or noise in the input function data inevitably amplify, leading to large, unbounded errors in the computed derivatives. This ill-posedness necessitates the application of specialized regularization techniques to recover stable and meaningful approximations.

A wide spectrum of computational strategies has been developed to address this challenge. Traditional approaches include finite difference schemes and polynomial interpolation [22, 16, 24], which, while simple, exhibit pronounced sensitivity to noise and grid non-uniformity. Spline-based methods [31, 13, 8, 18] offer greater robustness by leveraging piecewise polynomial smoothing, though controlling the inherent smoothness regularization remains critical. Mollification techniques [19, 11] stabilize differentiation through convolution with smooth kernels, effectively damping high-frequency noise. Methods utilizing reproducing kernel Hilbert spaces (RKHS) or radial basis functions (RBF) [6, 30, 7, 31] provide flexibility for scattered data approximation but depend heavily on kernel selection and regularization parameter tuning. Tikhonov regularization, in its standard and modified forms [34, 5, 20, 10, 4], remains a dominant paradigm, penalizing undesirable solution norms within a least-squares minimization of the data misfit. Spectral methods, particularly those employing truncated Fourier or Legendre series [3, 35, 36, 23, 33, 34], achieve high accuracy for smooth periodic functions, benefiting from convergence rates that adapt to the underlying smoothness. Integration-based approaches [14, 29, 28] exploit the inherent stability of the inverse operation (integration), while partial differential equation (PDE)-based methods [21, 27] reformulate differentiation as an inverse source problem, leveraging PDE theory for regularization, albeit often at the cost of increased complexity.

Among these, Fourier-based spectral methods are particularly attractive for their potential high-order convergence on periodic smooth data. To overcome their limitation to periodic functions, Fourier extension methods have emerged as a powerful technique [15, 2, 38, 39, 37]. The core idea involves extending the target function $f$, defined on a limited interval $[-1,1]$ and belonging to a Sobolev space $H^{p}(-1,1)$, to a function $f_{e}$ that is periodic on a larger interval $[-T,T](T>1)$ , such that $f_{e}(x)=f(x)$ on $[-1,1]$. This extension mitigates the Gibbs phenomenon by effectively transferring the non-periodicity boundary issues to the larger domain. Building upon this foundation, one of the authors of this paper and his collaborators introduced a significant advancement by incorporating a novel super-order regularization strategy within the Tikhonov framework [5]. Their method seeks an approximation $g=\mathscr{F}\mathbf{v}$ (where $\mathscr{F}$ maps the coefficient vector $\mathbf{v}$ to a trigonometric series) to the noisy data $f^{\delta}$ on $[-1,1]$ by minimizing the functional:

[TABLE]

The key innovation lies in the regularization operator $\mathscr{R}$, defined as

[TABLE]

which corresponds to penalizing the $L^{2}(-2,2)$ norm of $\sum_{l=0}^{\infty}D^{l}g/l!$. This super-order penalty imposes an exponentially weighted constraint on all derivatives simultaneously. Suppose that the noisy data satisfy the condition

[TABLE]

then the regularization parameter $\alpha$ is selected adaptively using the discrepancy principle ($\|\mathscr{F}\mathbf{v}_{\alpha}^{\delta}-f^{\delta}\|_{L^{2}(-1,1)}=C\delta,C>1$). Theoretically, this approach achieves order-optimal convergence rates

[TABLE]

for $k\leq p$ without requiring explicit knowledge of the smoothness index $p$.

While the Fourier extension method with super-order regularization demonstrates strong potential–particularly in effectively addressing non-periodic functions and achieving adaptive convergence—several limitations remain. Both numerical experiments and theoretical analysis indicate the need for further improvement in the following aspects:

• •

Difficulty with high-frequency functions: The regularization penalty suppresses high-frequency components, and floating-point precision imposes additional constraints. As a result, the method fails to accurately approximate derivatives of high-frequency functions such as $\cos(100x^{2})$, even when data noise is negligible.

• •

Inhibited accuracy for low-frequency functions: Solving the minimization problem requires forming the corresponding normal equation. According to [2], the best attainable approximation accuracy is $\mathcal{O}(\sqrt{\epsilon})$, where $\epsilon$ is machine precision. This error is further magnified during differentiation, making it difficult to achieve higher accuracy even for smooth, low-frequency functions.

• •

Lack of fast algorithms: To date, this class of methods lacks efficient implementations, which limits its practical applicability due to high computational costs.

To overcome these challenges, the authors recently proposed a local Fourier extension algorithm [37]. This approach avoids solving normal equations by truncating the singular values of the operator $\mathscr{F}$ directly, thereby improving the approximation accuracy to $\mathcal{O}(\epsilon)$. Furthermore, it leverages the rapid convergence of Fourier approximations to achieve high accuracy with few nodes for low-frequency functions. For high-frequency functions, interval partitioning effectively reduces local frequency, and the reuse of a unified discrete matrix across subdomains enhances computational efficiency. In this paper, we apply this method to numerical differentiation, addressing the three key limitations outlined above.

The remainder of this paper is organized as follows. In Section 2 we introduce the multi-interval Fourier extension formulation of the numerical differentiation problem and present the main regularization ideas that underpin our approach. Section 3 is devoted to the numerical implementation: we describe the uniform sampling discretization, the precomputation of the local Fourier bases and their singular value decompositions, and the adaptive domain partitioning strategy used to produce stable local fits. In Section 4 we summarize a comprehensive set of numerical experiments that illustrate the method’s accuracy and robustness to noise for a variety of smooth and highly oscillatory test functions. Finally, Section 5 collects concluding remarks and possible directions for future work.

2 Numerical Differentiation Method Based on Multi-Interval Fourier Extension

2.1 Single-Interval Generalized Truncated Singular Value Method

Let $\Omega=[0,2\pi]$, $\Lambda=\left[0,\frac{2\pi}{T}\right]$ with $T>1$. Define the system

[TABLE]

As shown in [15], this system forms an orthonormal basis for $L^{2}(\Omega)$ and a frame for $L^{2}(\Lambda)$.

Define the following operators:

[TABLE]

Consider the operator equation

[TABLE]

and its perturbed form

[TABLE]

with the noise assumption

[TABLE]

By the Sobolev extension theorem, for $g\in H^{p}(\Lambda)$, equation (4) admits multiple solutions. Among them, we define:

• The least-squares solution ${c}^{\star}\in\ell^{2}$ satisfies

[TABLE]

• The $B$-weighted best approximation ${c}_{B}^{\dagger}\in\ell^{2}$ satisfies

[TABLE]

Let $G=\mathcal{F}W^{-1}$, which is clearly a compact operator. Let $(\sigma_{i},v_{i},u_{i})$ be a singular system of $G$ with

[TABLE]

For any $g\in L^{2}(\Lambda)$, define the truncated inverse operator:

[TABLE]

It holds that

[TABLE]

We further have the following estimate:

Lemma 1**.**

The following bounds hold:

[TABLE]

Proof.

For any $g\in L^{2}(\Lambda)$, write

[TABLE]

Then we can get

[TABLE]

and

[TABLE]

Therefore,

[TABLE]

Moreover, for any unite vector ${c}\in\mathcal{D}(G)$, write

[TABLE]

It can be deduced that

[TABLE]

and hence,

[TABLE]

We conclude that

[TABLE]

∎

We define the generalized truncated SVD (GTSVD) solution to (5) as

[TABLE]

where the truncation parameter $\eta$ is chosen according to the discrepancy principle:

[TABLE]

Then we set

[TABLE]

as the approximation to $g$.

We also define

[TABLE]

Since $g\in H^{p}(\Lambda)$, there exists a constant $E>0$ such that

[TABLE]

And we have the following result:

Lemma 2**.**

The approximants ${c}_{B}^{\dagger}$ and ${c}_{B,N}$ satisfy

[TABLE]

where

[TABLE]

Proof.

First,

[TABLE]

Moreover,

[TABLE]

∎

Lemma 3** ([5]).**

If the vector sequence $\{{c}^{\delta}\}$ satisfies

[TABLE]

for some nonnegative constants $k_{1},k_{2},k_{3}$, then there exists a constant $E_{0}$ such that

[TABLE]

Lemma 4** ([1]).**

Let $\Theta=(a,b)\subset\mathbb{R}$ and suppose $f\in H^{s}(\Theta)$. Then for any $\epsilon\in(0,\epsilon_{0}]$ and $0\leq j\leq p$, there exists a constant $C$ (depending on $\epsilon_{0}$, $j$, and $p$) such that

[TABLE]

Theorem 5**.**

Suppose $g\in H^{p}(\Lambda)$ and let $\mathbf{g}_{\eta}^{\delta}$ be defined as in (12), with $\delta$ satisfying (6). Then for any $0\leq j\leq p$, the following convergence rate holds:

[TABLE]

Proof.

Using (6), (11), (14), and the triangle inequality, we obtain

[TABLE]

Using (8), (9), and (14), we estimate

[TABLE]

Let $\mathcal{U}=I-G\mathcal{T}_{\eta-1}$. Note that

[TABLE]

so

[TABLE]

Choose $N$ such that

[TABLE]

then

[TABLE]

Combining (20)–(23) and using the bound (15), we conclude that there exist constants $k_{1},k_{2},k_{3}$ such that

[TABLE]

Then, by Lemma 3, there exists a constant $E_{0}$ such that

[TABLE]

Therefore,

[TABLE]

Moreover,

[TABLE]

Finally, applying the interpolation inequality (Lemma 4) completes the proof of the theorem. ∎

Remark 2.1**.**

According to the regularization theory, the result in (19) is of optimal order, and the design of the method is independent of the smoothness $p$ and is adaptive. However, we need to note that the above method has serious numerical defects. Limited by the machine precision (set as $\epsilon$), we cannot obtain the exact value of the singular value of $\sigma_{i}<\epsilon$. According to the results in [39], even if $g$ is accurate, we can only obtain the approximate accuracy of

[TABLE]

For high-frequency functions, the norm $\|W{c}_{B}^{{\dagger}}\|_{\ell^{2}}$ is very large. We know that $e^{40}\approx 2.35e+17$, while $\epsilon=1e-16$ in general double-precision calculations. That is to say, when the function g contains high-frequency components that cannot be ignored, the above method cannot obtain satisfactory results even from the perspective of approximation, not to mention unstable calculations such as derivatives.

2.2 Multi-interval processing

Let $f\in H^{p}(a,b)$ be a non-periodic function defined on the interval $I=[a,b]$ and $f^{\delta}$ be its perturbation measurement data. Points $a=a_{0}<a_{1}<\cdots<a_{J}=b$ divides $I$ into $J$ subintervals $I_{1},I_{2},\ldots,I_{J}$. Let $f_{j}(x),f_{j}^{\delta}(x)$ be functions defined on interval $I_{j}$, satisfying:

[TABLE]

and

[TABLE]

Our goal is to obtain an effective approximation of $f$ through the perturbed data $f^{\delta}$.

Let $\{g_{k}(t)\}_{j=1}^{J}$ and $\{g^{\delta}_{k}(t)\}_{j=1}^{J}$ are two families of functions defined on $\Lambda$, such that

[TABLE]

where

[TABLE]

Then we have

[TABLE]

Now substitute $g^{\delta}_{j}$, $j=1,2,...J$ into the right side of equation (5), we can get a set of solution vectors from (10) and (11):

[TABLE]

and the corresponding approximations of the function $g_{j}$:

[TABLE]

Now we use

[TABLE]

as the approximation of $f$ and we have the following estimation.

Theorem 6**.**

Suppose that $f\in H^{p}(a,b)$ and the condition (29) holds, ${\bf f}^{\bm{\delta}}$ is defined by (35), and let

[TABLE]

then we have

[TABLE]

Proof.

[TABLE]

∎

Remark 2.2**.**

• •

Comparing the theoretical results of (19) and (26), it seems that there is no substantial advantage of the multi-interval method, but we need to note that the function $f(x)=e(\text{i}\pi\omega x)$ corresponds to

[TABLE]

That is to say, for high-frequency functions, we can make each $g_{j}$ a low-frequency function by dividing the region, and its corresponding Fourier coefficients converge quickly, and the corresponding $\|W{c}_{B}^{{\dagger}}\|_{\ell^{2}}$ can be maintained within the range allowed by the machine precision, so that the function $g_{j}$ can be effectively approximated at a smaller discrete scale.

• •

Later, we will show in the numerical implementation part that we can use the same discrete matrix to solve in each subinterval $I_{j}$, thereby greatly reducing the computational complexity of the problem.

• •

The results of (26) also show that over-division of the interval $[a,b]$ will reduce the accuracy. The scale of division should be related to the frequency of the function f. How to determine the scale of division, we will give a posterior method combining errors in the numerical implementation part.

3 Numerical implement based on uniform sampling

3.1 Discretization of equation (5)

From the previous analysis, we can see that we need to repeatedly solve equation (5) with different right-hand terms. For fixed $n$, let

[TABLE]

Then take the sampling ratio $\gamma\geq 1$ and let

[TABLE]

where $\lceil\cdot\rceil$ is the rounding symbol. When we select parameters $n,\gamma,T$ and take

[TABLE]

Then we can obtain the discrete form of equation (5):

[TABLE]

where

[TABLE]

Remark 3.1**.**

The parameters $n,\gamma$ and $T$ have a substantial impact on the properties of the discrete matrix, which in turn affects the final calculation results. Reference [37] conducted a large number of empirical tests on the relevant parameters, which we will not repeat here. In the following, we use the following parameters for calculation:

[TABLE]

Take

[TABLE]

where ${\mathbf{W}}_{n}$ is a diagonal matrix with

[TABLE]

as the diagonal vector. Now we can obtain the SVD of matrix $\mathbf{G}_{\gamma,n}^{T}$ as

[TABLE]

where ${\bm{\sigma}}_{j}$ are the singular values of $\mathbf{G}_{\gamma,n}^{T}$ such that

[TABLE]

while ${\bf u}_{i}$ and ${\bf v}_{i}$ are the left and right singular vectors of $\mathbf{G}_{\gamma,n}^{T}$, respectively.

Now if we assume that the exact data ${\bf g}$ and the perturbation data ${\bf g}^{\delta}$ satisfy:

[TABLE]

Then the generalized truncated singular value decomposition solution of (40) can be given as:

[TABLE]

where ${\bm{\eta}}$ is determined by the discrepancy principle

[TABLE]

with ${\bm{\tau}}>1$.

3.2 Adaptive domain partitioning strategy

The adaptive domain partitioning strategy forms the computational backbone of our multi-interval approach. The algorithm performs local Fourier-extension fits on adaptively refined subintervals and accepts a fit when its reconstruction error, appropriately measured and weighted by the assumed noise level, falls below a relaxed tolerance. To allow unlimited bisection, the number of global sampling nodes is chosen as

[TABLE]

where $m$ is the minimal number of nodes required in a single interval and $r$ is a nonnegative integer chosen so that the global grid aligns with dyadic subdivisions; after downsampling, local fits always use $m$ reference samples so that linear algebra quantities are consistent with the same $m$-point local fitting rule.

Next, we present a detailed description of the multi-interval (recursive bisection) implementation and give pseudocode that corresponds precisely to the MATLAB implementation. The overall procedure consists of the following four main stages:

Precomputation: Given $n$, $m$, $L$, $T$, compute:

• •

Sub-Fourier matrix $\mathbf{G}\in\mathbb{C}^{m\times(2n+1)}$

• •

Extended matrix $\mathbf{TG}\in\mathbb{C}^{L\times(2n+1)}$ for high-resolution reconstruction

• •

Derivative matrices $\mathbf{dG}=\mathbf{G}\cdot\operatorname{diag}(ik\cdot 2\pi/T)$, $\mathbf{dTG}=\mathbf{TG}\cdot\operatorname{diag}(ik\cdot 2\pi/T)$

• •

SVD of $\mathbf{G}=\mathbf{U}\mathbf{S}\mathbf{V}^{T}$

These are computed once and reused throughout. 2. 2.

Local Fourier Extension Fitting (Algorithm 1): For subinterval $[a_{\text{sub}},b_{\text{sub}}]$ with noisy data $\mathbf{y}_{d}\in\mathbb{R}^{nn}$:

• •

Compute downsampling ratio $r_{\text{ds}}=(nn-1)/(m-1)$

• •

Extract reduced vector $\mathbf{y}_{d,\text{red}}=\mathbf{y}_{d}[1:r_{\text{ds}}:nn]$

• •

Solve via truncated SVD: $\mathbf{Coef}=\text{LocalFourierFit}(\mathbf{y}_{d,\text{red}},\delta_{\text{loc}},m,\mathbf{U},\mathbf{S},\mathbf{V},\mathbf{G})$

• •

Reconstruct: $\mathbf{FE}=\Re(\mathbf{TG}\cdot\mathbf{Coef})$

• •

Interpolate to original resolution: $\mathbf{refinedFE}=\operatorname{interpft}(\mathbf{FE},r_{\text{ds}}\cdot L)$ 3. 3.

A Posteriori Error Criterion: Compute residual

[TABLE]

Given noise bound $\delta_{1}$, define $\epsilon_{\text{target}}=\sqrt{nn/3}\delta_{1}$. If

[TABLE]

accept the subinterval; otherwise, bisect and recurse. 4. 4.

Derivative Reconstruction (Algorithm 3): For each accepted subinterval with coefficients $\mathbf{Coef}$:

• •

Compute derivative approximation: $\mathbf{dFE}=\mathbf{dTG}\cdot\mathbf{Coef}\cdot Scale$ where $Scale=1/(b_{\text{sub}}-a_{\text{sub}})$

• •

Interpolate: $\mathbf{dFE}_{\text{interp}}=\operatorname{interpft}(\mathbf{dFE},r_{\text{ds}}\cdot L)$

• •

Concatenate results across subintervals

Remark 3.2**.**

Implementation Notes

• •

Global SVD reuse*: The SVD of $\mathbf{G}$ is computed once, enabling fast projections without refactorization.*

• •

Consistent downsampling*: $M=2^{r}(m-1)+1$ ensures $r_{\text{ds}}$ is integer-valued at all recursion levels.*

• •

Termination condition*: Recursion stops when $nn\leq m$ to prevent undersampling.*

• •

Noise adaptation*: The local tolerance $\delta_{\text{loc}}=\delta_{1}\sqrt{m/3}$ assumes uniform noise variance $\delta_{1}^{2}/3$.*

4 Numerical tests

In this section, we assess the performance of the proposed algorithm through a series of numerical experiments. All computations were performed on a Windows 10 system with 16 GB of memory and an Intel(R) Core(TM) i7-8500U CPU @ 1.80 GHz, using MATLAB 2016b. We compare the algorithm (M1) introduced herein with the full-data algorithm (M2) of [5]. For M2, we fix the parameters $T=2$ and $\gamma=2$ to facilitate implementation. Unless otherwise specified, the maximum recursion depth is set to $r=6$, corresponding to

[TABLE]

sampling nodes for the reconstruction of derivatives from perturbed data. For lower-frequency functions, smaller values of $r$ may be used in practice to reduce the computational cost.

The relative derivative error is defined by

[TABLE]

We consider the following three classes of test functions:

• •

Low-frequency functions:

[TABLE]

• •

High-frequency functions within the interval:

[TABLE]

• •

High-frequency functions near the interval boundaries:

[TABLE]

Figures 1–6 present the numerical results. For low-frequency functions (Figs. 1(b) and 2(b)), the two methods exhibit similar performance when the error level is relatively high. However, once the error level drops below a certain threshold, the accuracy of M2 saturates and no longer improves with decreasing data error, whereas the approximation error of M1 continues to decrease. This behavior arises because M2 solves the normal equations; see the analysis in [39] for further details.

For high-frequency functions, whether the oscillations occur within the interval or are concentrated near the boundaries (Figs. 3(b)–6(b)), M2 fails to produce accurate derivative estimates, while M1 consistently delivers high-accuracy approximations in both cases. Figures 1(d)–6(d) depict the distribution of the points $\{a_{j}\}$ for different test functions. The distribution is largely consistent with the frequency characteristics of the underlying function, with slight deviations attributable to data noise.

Overall, these experiments demonstrate that the proposed method M1 effectively overcomes the limitations of M2, providing a stable and efficient approach for numerical differentiation.

5 Conclusions and remarks

We conclude with the following remarks:

• •

We have proposed a numerical differentiation method based on multi-interval Fourier extension. By partitioning the domain and reusing precomputed matrices, the method substantially reduces computational complexity. Moreover, the frequency-reduction effect of interval partitioning mitigates the machine-precision limitations inherent in high-order regularization methods, thereby enabling accurate reconstruction of high-frequency functions. This offers valuable insights for the treatment of ill-posed problems.

• •

Although numerical experiments for higher-order derivatives are not presented, Theorem 2 ensures that the proposed approach applies to derivatives of arbitrary order.

• •

The reconstructed function produced by our method is, in general, discontinuous. This discontinuity can be removed, if desired, by applying a piecewise spline post-processing step to enforce smoothness. As this is a straightforward procedure, it is omitted here.