Exponential B-Spline Approximation

• Exponential B-spline approximation is a technique that extends classical B-spline methods by incorporating exponential terms and a tunable tension parameter for flexible, localized basis functions.
• It generates C2 continuous, piecewise-exponential basis functions that interpolate between cubic B-splines and piecewise-linear hat functions, enhancing numerical accuracy.
• Its efficient sparse-matrix implementations in collocation, Galerkin, and finite element frameworks provide robust solutions for PDEs with sharp gradients and boundary layers.

Exponential B-spline approximation generalizes polynomial B-spline methods by incorporating exponential terms with a tunable tension parameter, yielding localized, $C^2$ continuous, piecewise-exponential basis functions that encompass both the classical cubic B-spline and piecewise-linear “hat” functions as limiting cases. This framework offers enhanced flexibility for the high-order numerical approximation of functions and partial differential equations, especially those featuring boundary layers or sharp gradients. The methods admit efficient sparse-matrix representations and are widely used in collocation, Galerkin, and finite element solvers for time-dependent and steady PDEs, including Burgers’, Gardner, fractional sub-diffusion, Klein-Gordon, and advection-diffusion equations (Zhu et al., 2016, Görgülü et al., 2015, Hepsona et al., 2017, Hepson et al., 2016, Gorgulu et al., 2016).

1. Basis Function Construction and Properties

Exponential B-splines on a uniform mesh $a = x_0 < x_1 < \cdots < x_M = b$, $h = (b-a)/M$, are defined with a real tension parameter $p > 0$. Auxiliary quantities $s = \sinh(ph)$ and $c = \cosh(ph)$ enter closed-form expressions for the basis. For each index $j$, the basis function $B_j(x)$ is supported on $[x_{j-2}, x_{j+2}]$ and defined piecewise as

$$B_j(x) = \begin{cases} e\,(x_{j-2} - x) - \frac{e}{p} \sinh(p(x_{j-2} - x)), & x \in [x_{j-2}, x_{j-1}] \ a_0 + b_0 (x_j - x) + c_0 e^{p(x_j - x)} + d_0 e^{-p(x_j - x)}, & x \in [x_{j-1}, x_j] \ a_0 + b_0 (x - x_j) + c_0 e^{p(x - x_j)} + d_0 e^{-p(x - x_j)}, & x \in [x_j, x_{j+1}] \ e\,(x - x_{j+2}) - \frac{e}{p} \sinh(p(x - x_{j+2})), & x \in [x_{j+1}, x_{j+2}] \ 0 & \text{otherwise} \end{cases}$$

with explicit formulas for $e, a_0, b_0, c_0, d_0$ in terms of $p, h, s, c$ (Zhu et al., 2016).

At mesh nodes,

$$B_j(x_k) = \begin{cases} 1 & k = j \ \frac{s - ph}{2(phc - s)} & k = j \pm 1 \ 0 & |k-j| \geq 2 \end{cases}$$

and the partition of unity property holds: $\sum_j B_j(x) = 1$.

The key parameter $p$ controls basis behavior:

• $p \to 0$: $B_j$ reduces to the classical cubic B-spline.
• $p \to \infty$: $B_j$ becomes a piecewise-linear "hat" function.
• Intermediate $p$: $B_j$ interpolates between cubic and linear, reducing spurious oscillations near steep gradients (Zhu et al., 2016, Pool et al., 2022).

2. Derivatives and Evaluations

On each subinterval, $B_j$ and its derivatives admit closed forms. For example, on $[x_{j-2}, x_{j-1}]$, setting $\xi = x_{j-2} - x$,

$$B_j(x) = e \xi - \frac{e}{p} \sinh(p \xi), \quad B_j'(x) = -e + e \cosh(p \xi), \quad B_j''(x) = -p e \sinh(p \xi)$$

On central intervals, expressions involve combinations of exponentials, polynomials, and coefficients ($a_0$, $b_0$, etc.) (Zhu et al., 2016).

At the knots, explicit derivative values are: $$B_j'(x_{j\pm1}) = \frac{\mp p(1-c)}{2(phc-s)}, \quad B_j''(x_j) = -\frac{p^2 s}{phc-s}, \quad B_j''(x_{j\pm1}) = \frac{p^2 s}{2(phc-s)}$$ These support efficient computation of spatial derivatives in discretized PDEs.

3. Collocation and Galerkin Schemes

Approximation of a function $u(x,t)$ proceeds via expansion: $u_N(x,t) = \sum_j \alpha_j(t) B_j(x)$ Collocation methods enforce the PDE at node points using the above basis and its derivatives. The resulting linear or nonlinear systems inherit banded (often tridiagonal) structure due to the compact support of $B_j$ (Zhu et al., 2016, Hepsona et al., 2017, Pool et al., 2022).

A typical time-stepping scheme combines spatial exponential B-spline collocation (or Galerkin projection) with implicit time discretization, e.g., Crank–Nicolson, generating a matrix equation: $A \boldsymbol{\alpha}^n = \text{RHS}$ where $A$ is tridiagonal, enabling $O(M)$ solvers such as Thomas algorithm.

In the Galerkin framework, the mass, stiffness, and convection matrices are built from integrals over basis products and their derivatives (locally analytic or quadrature-computable), and the solution advances via time-stepping schemes. For nonlinear terms, inner iterations or linearization schemes (e.g., Rubin–Graves, predictor–corrector) are employed (Görgülü et al., 2015, Gorgulu et al., 2016).

4. Error Analysis and Approximation Order

Exponential B-splines of order $N$ reproduce exponential polynomials up to degree $N-1$, offering exponential polynomial reproduction beyond classical B-splines. On uniform grids, the interpolation error for $C^N$ functions is $\mathcal{O}(h^N)$ in the $L^\infty$ norm (Massopust, 2019). For cubic order ($N=4$), observed convergence is typically $O(h^4)$ for smooth problems (Zhu et al., 2016, Ersoy et al., 2016).

The tension parameter $p$ mildly influences the error constants; $p \approx 1$ is often optimal for sharp layers or shocks, while smaller $p$ is appropriate for globally smooth profiles (Zhu et al., 2016, Gorgulu et al., 2016, Pool et al., 2022).

5. Stability, Conservation, and Implementation

Schemes using exponential B-splines inherit stability properties from the underlying time-discretization (e.g., unconditional stability of Crank–Nicolson for linear problems). Von Neumann-type analysis shows that for the considered PDEs, the amplification factor satisfies $|\zeta^n| \leq |\zeta^0|$ regardless of $h$ or $\tau$ (Zhu et al., 2016, Hepson et al., 2016).

Conservation properties, such as invariants for the Klein–Gordon or Gardner equations, are preserved up to $10^{-6}$ accuracy in long simulations, attributed to the $C^2$ structure and symmetric formulation (Hepson et al., 2016, Hepsona et al., 2017).

Practical implementation comprises precomputation of basis constants, assembly of sparse system matrices (mass, stiffness), elimination or modification for boundary conditions, and time-stepping. Tuning of $p$ for performance can be automated by residual minimization.

Example pseudocode skeleton for 1D fractional subdiffusion (paraphrased from (Zhu et al., 2016)):

1. Setup: mesh, p, basis constants. 2. Assemble tridiagonal matrices for mass and stiffness using basis values and derivatives at nodes. 3. For each timestep: a. Update right-hand side (incorporating previous solutions and forcing). b. Solve tridiagonal system for new coefficients. 4. Recover approximate solution as linear combination of basis functions.

6. Generalizations: Complex Order and Fractal B-splines

Integral-order exponential B-splines generalize to complex orders $z$ with $E_{z,a}(x)$ defined via Fourier transform: $$\widehat{E}_{z,a}(\omega) = \left(\frac{1 - e^{-(a + i\omega)}}{a + i\omega}\right)^z$$ In the time-domain these can be written as infinite sums involving truncated powers and exponentials (Massopust, 2013, Massopust, 2019).

These complex-order splines admit fractional differential characterizations, generate multiresolution analyses, and support wavelet constructions. They interpolate smoothly between polynomial and exponential B-splines and allow tunable smoothness and decay via $\Re(z)$ and $a$ respectively (Massopust, 2013).

Fractal (self-referential) exponential B-splines are produced by contractive Read–Bajraktarević operators, generating uncountable families that retain global approximation or interpolation properties but introduce local self-similarity. Geometric convergence and error estimates are governed by the contraction parameters (Massopust, 2019).

7. Applications and Comparative Performance

Exponential B-spline methods have been extensively applied to:

• Fractional sub-diffusion equations via collocation and GMMP time discretization: yielding unconditional stability, rapid tridiagonal solvers, and observed $O(h^4),\ O(\tau)$ accuracy (Zhu et al., 2016).
• Nonlinear PDEs: Burgers’, RLW, Gardner, advection-diffusion, nonlinear Schrödinger, and Klein-Gordon equations, achieving high accuracy, robustness for traveling fronts and shocks, and improved resolution over polynomial B-splines in the presence of sharp features (Görgülü et al., 2015, Ersoy et al., 2016, Pool et al., 2022, Hepsona et al., 2017, Hepson et al., 2016, Ersoy et al., 2016).
• Conservation laws: Monitoring of invariants demonstrates numerical conservation to near-machine precision for integrable systems (Hepson et al., 2016, Hepsona et al., 2017).

Comparison studies indicate that exponential B-spline methods with optimized tension parameter generally outperform polynomial B-spline and standard Galerkin approaches in maximum error, particularly in resolving boundary layers and sharp gradients without introducing spurious oscillations (Gorgulu et al., 2016, Pool et al., 2022, Ersoy et al., 2016).

---

References:

• Exponential B-spline collocation for fractional sub-diffusion (Zhu et al., 2016)
• Exponential B-spline Galerkin for RLW (Görgülü et al., 2015); advection-diffusion (Gorgulu et al., 2016); Burgers' (Gorgulu et al., 2016)
• Gardner (Hepsona et al., 2017), Klein-Gordon (Hepson et al., 2016), Schrödinger (Ersoy et al., 2016)
• Generalizations and self-referential forms (Massopust, 2019, Massopust, 2013)
• Comparative review of spline methods (Pool et al., 2022)