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.

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.

The nonlinear Dirac equation is an important model in quantum physics with a set of conservation laws and a multi-symplectic formulation. In this paper, we propose energy-preserving and multi-symplectic wavelet algorithms for this model. Meanwhile, we evidently improve the efficiency of these algorithms in computations via splitting technique and explicit strategy. Numerical experiments are conducted during long-term simulations to show the excellent performances of the proposed algorithms and verify our theoretical analysis.

The computation of integrals in higher dimensions and on general domains, when no explicit cubature rules are known, can be ”easily” addressed by means of the quasi-Monte Carlo method. The method, simple in its formulation, becomes computationally inefficient when the space dimension is growing and the integration domain is particularly complex. In this paper we present two new approaches to the quasi-Monte Carlo method for cubature based on nonnegative least squares and approximate Fekete points. The main idea is to use less points and especially good points for solving the system of the moments. Good points are here intended as points with good interpolation properties, due to the strict connection between interpolation and cubature. Numerical experiments show that, in average, just a tenth of the points should be used mantaining the same approximation order of the quasi-Monte Carlo method. The method has been satisfactory applied to 2 and 3-dimensional problems on quite complex domains.

We present a simple, accurate method for computing singular or nearly singular integrals on a smooth, closed surface, such as layer potentials for harmonic functions evaluated at points on or near the surface. The integral is computed with a regularized kernel and corrections are added for regularization and discretization, which are found from analysis near the singular point. The surface integrals are computed from a new quadrature rule using surface points which project onto grid points in coordinate planes. The method does not require coordinate charts on the surface or special treatment of the singularity other than the corrections. The accuracy is about O(h3), where h is the spacing in the background grid, uniformly with respect to the point of evaluation, on or near the surface. Improved accuracy is obtained for points on the surface. The treecode of Duan and Krasny for Ewald summation is used to perform sums. Numerical examples are presented with a variety of surfaces.

A simple method is proposed for constructing fourth-degree cubature formulae over general product regions with no symmetric assumptions. The cubature formulae that are constructed contain at most n2 + 7n + 3 nodes and they are likely the first kind of fourth-degree cubature formulae with roughly n2 nodes for non-symmetric integrations. Moreover, two special cases are given to reduce the number of nodes further. A theoretical upper bound for minimal number of cubature nodes is also obtained.

Computing the value of a high-dimensional integral can often be reduced to the problem of finding the ratio between the measures of two sets. Monte Carlo methods are often used to approximate this ratio, but often one set will be exponentially larger than the other, which leads to an exponentially large variance. A standard method of dealing with this problem is to interpolate between the sets with a sequence of nested sets where neighboring sets have relative measures bounded above by a constant. Choosing such a well-balanced sequence can rarely be done without extensive study of a problem. Here a new approach that automatically obtains such sets is presented. These well-balanced sets allow for faster approximation algorithms for integrals and sums using fewer samples, and better tempering and annealing Markov chains for generating random samples. Applications, such as finding the partition function of the Ising model and normalizing constants for posterior distributions in Bayesian methods, are discussed.

The evolution of precipitates in stressed solids is modeled by coupling a quasi-steady diffusion equation and a linear elasticity equation with dynamic boundary conditions. The governing equations are solved numerically using a boundary integral method (BIM). A critical step in applying BIM is to develop fast algorithms to reduce the arithmetic operation count of matrix-vector multiplications. In this paper, we develop a fast adaptive treecode algorithm for the diffusion and elasticity problems in two dimensions (2D). We present a novel source dividing strategy to parallelize the treecode. Numerical results show that the speedup factor is nearly perfect up to a moderate number of processors. This approach of parallelization can be readily implemented in other treecodes using either uniform or non-uniform point distribution. We demonstrate the effectiveness of the treecode by computing the long-time evolution of a complicated microstructure in elastic media, which would be extremely difficult with a direct summation method due to CPU time constraint. The treecode speeds up computations dramatically while fulfilling the stringent precision requirement dictated by the spectrally accurate BIM.

Pyramidal elements are often used to connect tetrahedral and hexahedral elements in the finite element method. In this paper we derive three new higher order numerical cubature formulae for pyramidal elements.

This paper deals with the solvability and the convergence of a class of unsymmetric Meshless Local Petrov-Galerkin (MLPG) method with radial basis function (RBF) kernels generated trial spaces. Local weak-form testings are done with step-functions. It is proved that subject to sufficiently many appropriate testings, solvability of the unsymmetric RBF-MLPG resultant systems can be guaranteed. Moreover, an error analysis shows that this numerical approximation converges at the same rate as found in RBF interpolation. Numerical results (in double precision) give good agreement with the provided theory.

Nonlinear eigenvalue problems for fluxons in long Josephson junctions with exponentially varying width are treated. Appropriate algorithms are created and realized numerically. The results obtained concern the stability of the fluxons, the centering both magnetic field and current for the magnetic flux quanta in the Josephson junction as well as the ascertaining of the impact of the geometric and physical parameters on these quantities. Each static solution of the nonlinear boundary-value problem is identified as stable or unstable in dependence on the eigenvalues of associated Sturm-Liouville problem. The above compound problem is linearized and solved by using of the reliable Continuous analogue of Newton method.

This paper is devoted to the study of the Eulerian-Lagrangian method (ELM) for convection-diffusion equations on unstructured grids with or without accurate numerical integration. We first propose an efficient and accurate algorithm to calculate the integrals in the Eulerian-Lagrangian method. Our approach is based on an algorithm for finding the intersection of two non-matching grids. It has optimal algorithmic complexity and runs fast enough to make time-dependent velocity fields feasible. The evaluation of the integrals leads to increased precision and the unconditional stability. We demonstrate by numerical examples that the ELM with our proposed algorithm for accurate numerical integration has the following two features: first it is much more accurate and more stable than the ones with traditional numerical integration techniques and secondly the overall cost of the proposed method is comparable with the traditional ones.

In this paper we consider polynomials orthogonal with respect to the linear functional defined on the space of all algebraic polynomials by

where α,β> -1/2 are real numbers such that ℓ = |β - α| is a positive integer, and ζ∈ R{0}. We prove the existence of such orthogonal polynomials for some pairs of α and ζ and for all nonnegative integers ℓ. For such orthogonal polynomials we derive three-term recurrence relations and also some differential-difference relations. For such orthogonal polynomials the corresponding quadrature rules of Gaussian type are considered. Also, some numerical examples are included.

Using the Girsanov transformation we derive estimates for the accuracy of piecewise approximations for one-sided and two-sided boundary crossing probabilities. We demonstrate that piecewise linear approximations can be calculated using repeated numerical integration. As an illustrative example we consider the case of one-sided and two-sided square-root boundaries for which we also present analytical representations in a form of infinite power series.

The corpuscle problem of Wicksell is discussed. We give a numerical quadrature of Gauss–Chebyshev type for Wicksell's integral equation which combines a size distribution of discs on a sectional plane with that of spheres. We also give an estimation procedure of three-dimensional size distributions based on this quadrature and examine its theoretical properties. In practice, we need a smoothing technique for empirical distribution functions before applying this estimator. Simulation results are given. Our idea also is applied to the thick section case and an analysis of microscopic data is given.

We investigate convergence in a weighted L1 -norm of Hermite-Fejér and Hermite interpolation at the zeros of orthogonal polynomials associated with weights on the real line. The results are then applied to convergences of product integration rules. From the point of view of orthogonal polynomials, the new feature is that Freud and Erdös weights are treated simultaneously and that relatively few assumptions are placed on the weight. From the point of view of product integration, the rules exhibit convergence for highly oscillatory kernels (for example) and for functions of rapid growth at infinity.