In this paper,the linear space $\mathcal F$ of a special type of fractal interpolation functions (FIFs) on an interval I is considered. Each FIF in $\mathcal F$ is established from a continuous function on I. We show that, for a finite set of linearly independent continuous functions on I, we get linearly independent FIFs. Then we study a finite-dimensional reproducing kernel Hilbert space (RKHS) $\mathcal F_{\mathcal B}\subset\mathcal F$, and the reproducing kernel $\mathbf k$ for $\mathcal F_{\mathcal B}$ is defined by a basis of $\mathcal F_{\mathcal B}$. For a given data set $\mathcal D=\{(t_k, y_k) : k=0,1,\ldots,N\}$, we apply our results to curve fitting problems of minimizing the regularized empirical error based on functions of the form $f_{\mathcal V}+f_{\mathcal B}$, where $f_{\mathcal V}\in C_{\mathcal V}$ and $f_{\mathcal B}\in \mathcal F_{\mathcal B}$. Here $C_{\mathcal V}$ is another finite-dimensional RKHS of some classes of regular continuous functions with the reproducing kernel $\mathbf k^*$. We show that the solution function can be written in the form $f_{\mathcal V}+f_{\mathcal B}=\sum_{m=0}^N\gamma_m\mathbf k^*_{t_m} +\sum_{j=0}^N \alpha_j\mathbf k_{t_j}$, where ${\mathbf k}_{t_m}^\ast(\cdot)={\mathbf k}^\ast(\cdot,t_m)$ and $\mathbf k_{t_j}(\cdot)=\mathbf k(\cdot,t_j)$, and the coefficients γm and αj can be solved by a system of linear equations.

We propose a novel time-asymptotically stable, implicit–explicit, adaptive, time integration method (denoted by the $\theta $-method) for the solution of the fractional advection–diffusion-reaction (FADR) equations. The spectral analysis of the method (involving the group velocity and the phase speed) indicates a region of favourable dispersion for a limited range of Péclet number. The numerical inversion of the coefficient matrix is avoided by exploiting the sparse structure of the matrix in the iterative solver for the Poisson equation. The accuracy and the efficacy of the method is benchmarked using (a) the two-dimensional fractional diffusion equation, originally proposed by researchers earlier, and (b) the incompressible, subdiffusive dynamics of a planar viscoelastic channel flow of the Rouse chain melts (FADR equation with fractional time-derivative of order ) and the Zimm chain solution (). Numerical simulations of the viscoelastic channel flow effectively capture the nonhomogeneous regions of high viscosity at low fluid inertia (or the so-called “spatiotemporal macrostructures”), experimentally observed in the flow-instability transition of subdiffusive flows.

First, we build a computational procedure to reconstruct the convex body from a pre-given surface area measure. Nontrivially, we prove the convergence of this procedure. Then, the sufficient and necessary conditions of a Sobolev binary function to be a lightness function of a convex body are investigated. Finally, we present a computational procedure to compute the curvature function from a pre-given lightness function, and then we reconstruct the convex body from this curvature function by using the first procedure. The convergence is also proved. The main computations in both procedures are simply about the spherical harmonics.

In our previous paper (Gaster et al., 2018, arXiv:1810.11932), we showed that the theory of harmonic maps between Riemannian manifolds, especially hyperbolic surfaces, may be discretized by introducing a triangulation of the domain manifold with independent vertex and edge weights. In the present paper, we study convergence of the discrete theory back to the smooth theory when taking finer and finer triangulations, in the general Riemannian setting. We present suitable conditions on the weighted triangulations that ensure convergence of discrete harmonic maps to smooth harmonic maps, introducing the notion of (almost) asymptotically Laplacian weights, and we offer a systematic method to construct such weighted triangulations in the two-dimensional case. Our computer software Harmony successfully implements these methods to compute equivariant harmonic maps in the hyperbolic plane.

We study some properties of integro splines. Using these properties, we design an algorithm to construct splines $S_{m+1}(x)$ of neighbouring degrees to the given spline $S_m(x)$ with degree m. A local integro-sextic spline is constructed with the proposed algorithm. The local integro splines work efficiently, that is, they have low computational complexity, and they are effective for use in real time. The construction of nonlocal integro splines usually leads to solving a system of linear equations with band matrices, which yields high computational costs.

We consider the problem of numerical integration when the sampling nodes form a stationary point process on the real line. In previous papers it was argued that a naïve Riemann sum approach can cause a severe variance inflation when the sampling points are not equidistant. We show that this inflation can be avoided using a higher-order Newton–Cotes quadrature rule which exploits smoothness properties of the integrand. Under mild assumptions, the resulting estimator is unbiased and its variance asymptotically obeys a power law as a function of the mean point distance. If the Newton–Cotes rule is of sufficiently high order, the exponent of this law turns out to only depend on the point process through its mean point distance. We illustrate our findings with the stereological estimation of the volume of a compact object, suggesting alternatives to the well-established Cavalieri estimator.

The discontinuous Galerkin (DG) method provides a robust and flexible technique for the time integration of fractional diffusion problems. However, a practical implementation uses coefficients defined by integrals that are not easily evaluated. We describe specialized quadrature techniques that efficiently maintain the overall accuracy of the DG method. In addition, we observe in numerical experiments that known superconvergence properties of DG time stepping for classical diffusion problems carry over in a modified form to the fractional-order setting.

Effective and accurate high-degree spline interpolation is still a challenging task in today’s applications. Higher degree spline interpolation is not so commonly used, because it requires the knowledge of higher order derivatives at the nodes of a function on a given mesh.

In this article, our goal is to demonstrate the continuity of the piecewise polynomials and their derivatives at the connecting points, obtained with a method initially developed by Beaudoin (1998, 2003) and Beauchemin (2003). This new method, involving the discrete Fourier transform (DFT/FFT), leads to higher degree spline interpolation for equally spaced data on an interval $[0,T]$. To do this, we analyze the singularities that may occur when solving the system of equations that enables the construction of splines of any degree. We also note an important difference between the odd-degree splines and even-degree splines. These results prove that Beaudoin and Beauchemin’s method leads to spline interpolation of any degree and that this new method could eventually be used to improve the accuracy of spline interpolation in traditional problems.

The Orthogonal Least Squares (OLS) algorithm is an efficient sparse recovery algorithm that has received much attention in recent years. On one hand, this paper considers that the OLS algorithm recovers the supports of sparse signals in the noisy case. We show that the OLS algorithm exactly recovers the support of $K$-sparse signal $\boldsymbol{x}$ from $\boldsymbol{y}=\boldsymbol{\unicode[STIX]{x1D6F7}}\boldsymbol{x}+\boldsymbol{e}$ in $K$ iterations, provided that the sensing matrix $\boldsymbol{\unicode[STIX]{x1D6F7}}$ satisfies the restricted isometry property (RIP) with restricted isometry constant (RIC) $\unicode[STIX]{x1D6FF}_{K+1}<1/\sqrt{K+1}$, and the minimum magnitude of the nonzero elements of $\boldsymbol{x}$ satisfies some constraint. On the other hand, this paper demonstrates that the OLS algorithm exactly recovers the support of the best $K$-term approximation of an almost sparse signal $\boldsymbol{x}$ in the general perturbations case, which means both $\boldsymbol{y}$ and $\boldsymbol{\unicode[STIX]{x1D6F7}}$ are perturbed. We show that the support of the best $K$-term approximation of $\boldsymbol{x}$ can be recovered under reasonable conditions based on the restricted isometry property (RIP).

In the paper the correspondence between a formal multiple power series and a special type of branched continued fractions, the so-called ‘multidimensional regular C-fractions with independent variables’ is analysed providing with an algorithm based upon the classical algorithm and that enables us to compute from the coefficients of the given formal multiple power series, the coefficients of the corresponding multidimensional regular C-fraction with independent variables. A few numerical experiments show, on the one hand, the efficiency of the proposed algorithm and, on the other, the power and feasibility of the method in order to numerically approximate certain multivariable functions from their formal multiple power series.

We give bounds on the error in the asymptotic approximation of the log-Gamma function $\ln \unicode[STIX]{x1D6E4}(z)$ for complex $z$ in the right half-plane. These improve on earlier bounds by Behnke and Sommer [Theorie der analytischen Funktionen einer komplexen Veränderlichen, 2nd edn (Springer, Berlin, 1962)], Spira [‘Calculation of the Gamma function by Stirling’s formula’, Math. Comp.25 (1971), 317–322], and Hare [‘Computing the principal branch of log-Gamma’, J. Algorithms25 (1997), 221–236]. We show that $|R_{k+1}(z)/T_{k}(z)|<\sqrt{\unicode[STIX]{x1D70B}k}$ for nonzero $z$ in the right half-plane, where $T_{k}(z)$ is the $k$th term in the asymptotic series, and $R_{k+1}(z)$ is the error incurred in truncating the series after $k$ terms. We deduce similar bounds for asymptotic approximation of the Riemann–Siegel theta function $\unicode[STIX]{x1D717}(t)$. We show that the accuracy of a well-known approximation to $\unicode[STIX]{x1D717}(t)$ can be improved by including an exponentially small term in the approximation. This improves the attainable accuracy for real $t>0$ from $O(\exp (-\unicode[STIX]{x1D70B}t))$ to $O(\exp (-2\unicode[STIX]{x1D70B}t))$. We discuss a similar example due to Olver [‘Error bounds for asymptotic expansions, with an application to cylinder functions of large argument’, in: Asymptotic Solutions of Differential Equations and Their Applications (ed. C. H. Wilcox) (Wiley, New York, 1964), 16–18], and a connection with the Stokes phenomenon.

A new minimization principle for the Poisson equation using two variables – the solution and the gradient of the solution – is introduced. This principle allows us to use any conforming finite element spaces for both variables, where the finite element spaces do not need to satisfy the so-called inf–sup condition. A numerical example demonstrates the superiority of this approach.

In this paper we consider the algorithm for recovering sparse orthogonal polynomials using stochastic collocation via ℓq minimization. The main results include: 1) By using the norm inequality between ℓq and ℓ2 and the square root lifting inequality, we present several theoretical estimates regarding the recoverability for both sparse and non-sparse signals via ℓq minimization; 2) We then combine this method with the stochastic collocation to identify the coefficients of sparse orthogonal polynomial expansions, stemming from the field of uncertainty quantification. We obtain recoverability results for both sparse polynomial functions and general non-sparse functions. We also present various numerical experiments to show the performance of the ℓq algorithm. We first present some benchmark tests to demonstrate the ability of ℓq minimization to recover exactly sparse signals, and then consider three classical analytical functions to show the advantage of this method over the standard ℓ1 and reweighted ℓ1 minimization. All the numerical results indicate that the ℓq method performs better than standard ℓ1 and reweighted ℓ1 minimization.

The computation efficiency of high dimensional (3D and 4D) B-spline interpolation, constructed by classical tensor product method, is improved greatly by precomputing the B-spline function. This is due to the character of NLT code, i.e. only the linearised characteristics are needed so that the unperturbed orbit as well as values of the B-spline function at interpolation points can be precomputed at the beginning of the simulation. By integrating this fixed point interpolation algorithm into NLT code, the high dimensional gyro-kinetic Vlasov equation can be solved directly without operator splitting method which is applied in conventional semi-Lagrangian codes. In the Rosenbluth-Hinton test, NLT runs a few times faster for Vlasov solver part and converges at about one order larger time step than conventional splitting code.

In this paper we introduce high dimensional tensor product interpolation for the combination technique. In order to compute the sparse grid solution, the discrete numerical subsolutions have to be extended by interpolation. If unsuitable interpolation techniques are used, the rate of convergence is deteriorated. We derive the necessary framework to preserve the error structure of high order finite difference solutions of elliptic partial differential equations within the combination technique framework. This strategy enables us to obtain high order sparse grid solutions on the full grid. As exemplifications for the case of order four we illustrate our theoretical results by two test examples with up to four dimensions.

Iterative Filtering (IF) is an alternative technique to the Empirical Mode Decomposition (EMD) algorithm for the decomposition of non–stationary and non–linear signals. Recently in [3] IF has been proved to be convergent for any L2 signal and its stability has been also demonstrated through examples. Furthermore in [3] the so called Fokker–Planck (FP) filters have been introduced. They are smooth at every point and have compact supports. Based on those results, in this paper we introduce the Multidimensional Iterative Filtering (MIF) technique for the decomposition and time–frequency analysis of non–stationary high–dimensional signals. We present the extension of FP filters to higher dimensions. We prove convergence results under general sufficient conditions on the filter shape. Finally we illustrate the promising performance of MIF algorithm, equipped with high–dimensional FP filters, when applied to the decomposition of two dimensional signals.