Search

For a large class of orthogonal basis functions, there has been a recent identification of expansion methods for computing accurate, stable approximations of a quantity of interest. This paper presents, within the context of uncertainty quantification, a practical implementation using basis adaptation, and coherence motivated sampling, which under assumptions has satisfying guarantees. This implementation is referred to as Basis Adaptive Sample Efficient Polynomial Chaos (BASE-PC). A key component of this is the use of anisotropic polynomial order which admits evolving global bases for approximation in an efficient manner, leading to consistently stable approximation for a practical class of smooth functionals. This fully adaptive, non-intrusive method, requires no a priori information of the solution, and has satisfying theoretical guarantees of recovery. A key contribution to stability is the use of a presented correction sampling for coherence-optimal sampling in order to improve stability and accuracy within the adaptive basis scheme. Theoretically, the method may dramatically reduce the impact of dimensionality in function approximation, and numerically the method is demonstrated to perform well on problems with dimension up to 1000. [ABSTRACT FROM AUTHOR]

Abstract: In this paper, we approximate the semi-classical Schrödinger equation in the presence of electromagnetic field by the Hagedorn wave packets approach. By operator splitting, the Hamiltonian is divided into the modified part and the residual part. The modified Hamiltonian, which is the main new idea of this paper, is chosen by the fact that Hagedorn wave packets are localized both in space and momentum so that a crucial correction term is added to the truncated Hamiltonian, and is treated by evolving the parameters associated with the Hagedorn wave packets. The residual part is treated by a Galerkin approximation. We prove that, with the modified Hamiltonian only, the Hagedorn wave packets dynamics give the asymptotic solution with error , where ε is the scaled Planck constant. We also prove that, the Galerkin approximation for the residual Hamiltonian can reduce the approximation error to , where k depends on the number of Hagedorn wave packets added to the dynamics. This approach is easy to implement, and can be naturally extended to the multidimensional cases. Unlike the high order Gaussian beam method, in which the non-constant cut-off function is necessary and some extra error is introduced, the Hagedorn wave packets approach gives a practical way to improve accuracy even when ε is not very small. [Copyright &y& Elsevier]

This paper presents two discrete maximum principles (DMP) for the numerical solution of second order elliptic equations arising from the weak Galerkin finite element method. The results are established by assuming an h -acute angle condition for the underlying finite element triangulations. The mathematical theory is based on the well-known De Giorgi technique adapted in the finite element context. Some numerical results are reported to validate the theory of DMP. [ABSTRACT FROM AUTHOR]

In this paper, we propose a Newton–Krylov solver and a Picard–Krylov solver for finite element discrete problem of stationary incompressible magnetohydrodynamic equations in three dimensions. Using a mixed finite element method, we discretize the velocity and the pressure by H 1 ( Ω ) -conforming finite elements and discretize the magnetic field by H ( curl , Ω ) -conforming edge elements. An efficient preconditioner is proposed to accelerate the convergence of GMRES method for solving linearized discrete problems. By extensive numerical experiments, we demonstrate the robustness of the Newton–Krylov solver for relatively large physical parameters and the optimality with respect to the number of degrees of freedom. Moreover, the numerical experiments show that the Newton–Krylov solver is more robust than the Picard–Krylov solver for large Reynolds number. [ABSTRACT FROM AUTHOR]

Multilevel Monte Carlo (MLMC) is a recently proposed variation of Monte Carlo (MC) simulation that achieves variance reduction by simulating the governing equations on a series of spatial (or temporal) grids with increasing resolution. Instead of directly employing the fine grid solutions, MLMC estimates the expectation of the quantity of interest from the coarsest grid solutions as well as differences between each two consecutive grid solutions. When the differences corresponding to finer grids become smaller, hence less variable, fewer MC realizations of finer grid solutions are needed to compute the difference expectations, thus leading to a reduction in the overall work. This paper presents an extension of MLMC, referred to as multilevel control variates (MLCV) , where a low-rank approximation to the solution on each grid, obtained primarily based on coarser grid solutions, is used as a control variate for estimating the expectations involved in MLMC. Cost estimates as well as numerical examples are presented to demonstrate the advantage of this new MLCV approach over the standard MLMC when the solution of interest admits a low-rank approximation and the cost of simulating finer grids grows fast. [ABSTRACT FROM AUTHOR]

Abstract: This paper presents a new computational framework for the simulation of solid mechanics on general overlapping grids with adaptive mesh refinement (AMR). The approach, described here for time-dependent linear elasticity in two and three space dimensions, is motivated by considerations of accuracy, efficiency and flexibility. We consider two approaches for the numerical solution of the equations of linear elasticity on overlapping grids. In the first approach we solve the governing equations numerically as a second-order system (SOS) using a conservative finite-difference approximation. The second approach considers the equations written as a first-order system (FOS) and approximates them using a second-order characteristic-based (Godunov) finite-volume method. A principal aim of the paper is to present the first careful assessment of the accuracy and stability of these two representative schemes for the equations of linear elasticity on overlapping grids. This is done by first performing a stability analysis of analogous schemes for the first-order and second-order scalar wave equations on an overlapping grid. The analysis shows that non-dissipative approximations can have unstable modes with growth rates proportional to the inverse of the mesh spacing. This new result, which is relevant for the numerical solution of any type of wave propagation problem on overlapping grids, dictates the form of dissipation that is needed to stabilize the scheme. Numerical experiments show that the addition of the indicated form of dissipation and/or a separate filter step can be used to stabilize the SOS scheme. They also demonstrate that the upwinding inherent in the Godunov scheme, which provides dissipation of the appropriate form, stabilizes the FOS scheme. We then verify and compare the accuracy of the two schemes using the method of analytic solutions and using problems with known solutions. These latter problems provide useful benchmark solutions for time dependent elasticity. We also consider two problems in which exact solutions are not available, and use a posterior error estimates to assess the accuracy of the schemes. One of these two problems is additionally employed to demonstrate the use of dynamic AMR and its effectiveness for resolving elastic “shock” waves. Finally, results are presented that compare the computational performance of the two schemes. These demonstrate the speed and memory efficiency achieved by the use of structured overlapping grids and optimizations for Cartesian grids. [Copyright &y& Elsevier]

Abstract: In this follow up paper to our previous study in Bayona et al. (2011) , we present a new technique to compute the solution of PDEs with the multiquadric based RBF finite difference method (RBF-FD) using an optimal node dependent variable value of the shape parameter. This optimal value is chosen so that, to leading order, the local approximation error of the RBF-FD formulas is zero. In our previous paper (Bayona et al., 2011) we considered the case of an optimal (constant) value of the shape parameter for all the nodes. Our new results show that, if one allows the shape parameter to be different at each grid point of the domain, one may obtain very significant accuracy improvements with a simple and inexpensive numerical technique. We analyze the same examples studied in Bayona et al. (2011) , both with structured and unstructured grids, and compare our new results with those obtained previously. We also find that, if there are a significant number of nodes for which no optimal value of the shape parameter exists, then the improvement in accuracy deteriorates significantly. In those cases, we use generalized multiquadrics as RBFs and choose the exponent of the multiquadric at each node to assure the existence of an optimal variable shape parameter. [Copyright &y& Elsevier]

Abstract: A key idea in finite difference weighted essentially non-oscillatory (WENO) schemes is a combination of lower order fluxes to obtain a higher order approximation. The choice of the weight to each candidate stencil, which is a nonlinear function of the grid values, is crucial to the success of WENO schemes. For the system case, WENO schemes are based on local characteristic decompositions and flux splitting to avoid spurious oscillation. But the cost of computation of nonlinear weights and local characteristic decompositions is very high. In this paper, we investigate hybrid schemes of WENO schemes with high order up-wind linear schemes using different discontinuity indicators and explore the possibility in avoiding the local characteristic decompositions and the nonlinear weights for part of the procedure, hence reducing the cost but still maintaining non-oscillatory properties for problems with strong shocks. The idea is to identify discontinuity by an discontinuity indicator, then reconstruct numerical flux by WENO approximation in discontinuous regions and up-wind linear approximation in smooth regions. These indicators are mainly based on the troubled-cell indicators for discontinuous Galerkin (DG) method which are listed in the paper by Qiu and Shu (J. Qiu, C.-W. Shu, A comparison of troubled-cell indicators for Runge–Kutta discontinuous Galerkin methods using weighted essentially non-oscillatory limiters, SIAM Journal of Scientific Computing 27 (2005) 995–1013). The emphasis of the paper is on comparison of the performance of hybrid scheme using different indicators, with an objective of obtaining efficient and reliable indicators to obtain better performance of hybrid scheme to save computational cost. Detail numerical studies in one- and two-dimensional cases are performed, addressing the issues of efficiency (less CPU time and more accurate numerical solution), non-oscillatory property. [ABSTRACT FROM AUTHOR]

Abstract: The recently developed Flexible Local Approximation MEthod (FLAME) produces accurate difference schemes by replacing the usual Taylor expansion with Trefftz functions – local solutions of the underlying differential equation. This paper advances and casts in a general form a significant modification of FLAME proposed recently by Pinheiro and Webb: a least-squares fit instead of the exact match of the approximate solution at the stencil nodes. As a consequence of that, FLAME schemes can now be generated on irregular stencils with the number of nodes substantially greater than the number of approximating functions. The accuracy of the method is preserved but its robustness is improved. For demonstration, the paper presents a number of numerical examples in 2D and 3D: electrostatic (magnetostatic) particle interactions, scattering of electromagnetic (acoustic) waves, and wave propagation in a photonic crystal. The examples explore the role of the grid and stencil size, of the number of approximating functions, and of the irregularity of the stencils. [Copyright &y& Elsevier]

Abstract: A set of optimised boundary closure schemes is presented for use with compact central finite difference schemes in computational aeroacoustics (CAA) involving non-trivial boundaries. The boundary schemes are given in a form of non-central compact finite differences. They maintain fourth-order accuracy, a pentadiagonal matrix system and seven-point stencil which the main interior scheme employs. This paper introduces a new strategy to optimise the boundary schemes in the spectral domain and achieve the best resolution characteristics given a strict tolerance for the dispersion and dissipation errors. The boundary schemes are derived from sophisticated extrapolation of solutions outside the domain. The extrapolation functions are devised by combining polynomials and trigonometric series which contain extra control variables used to optimise the resolution characteristics. The differencing coefficients of the boundary schemes are determined in association with the existing coefficients of the interior scheme which is also optimised through an improved procedure in this paper. The accuracy of the proposed schemes is demonstrated by their application to CAA benchmark problems. [Copyright &y& Elsevier]

In this work, we apply a Matrix version of the so-called Discrete Empirical Interpolation (MDEIM) for the efficient reduction of nonaffine parametrized systems arising from the discretization of linear partial differential equations. Dealing with affinely parametrized operators is crucial in order to enhance the online solution of reduced-order models (ROMs). However, in many cases such an affine decomposition is not readily available, and must be recovered through (often) intrusive procedures, such as the empirical interpolation method (EIM) and its discrete variant DEIM. In this paper we show that MDEIM represents a very efficient approach to deal with complex physical and geometrical parametrizations in a non-intrusive, efficient and purely algebraic way. We propose different strategies to combine MDEIM with a state approximation resulting either from a reduced basis greedy approach or Proper Orthogonal Decomposition. A posteriori error estimates accounting for the MDEIM error are also developed in the case of parametrized elliptic and parabolic equations. Finally, the capability of MDEIM to generate accurate and efficient ROMs is demonstrated on the solution of two computationally-intensive classes of problems occurring in engineering contexts, namely PDE-constrained shape optimization and parametrized coupled problems. [ABSTRACT FROM AUTHOR]

This paper proposes an approach for high-order time integration within a multi-domain setting for time-fractional differential equations. Since the kernel is singular or nearly singular, two main difficulties arise after the domain decomposition: how to properly account for the history/memory part and how to perform the integration accurately. To address these issues, we propose a novel hybrid approach for the numerical integration based on the combination of three-term-recurrence relations of Jacobi polynomials and high-order Gauss quadrature. The different approximations used in the hybrid approach are justified theoretically and through numerical examples. Based on this, we propose a new multi-domain spectral method for high-order accurate time integrations and study its stability properties by identifying the method as a generalized linear method. Numerical experiments confirm hp -convergence for both time-fractional differential equations and time-fractional partial differential equations. [ABSTRACT FROM AUTHOR]

In this paper, a tensor decomposition approach combined with Chebyshev spectral differentiation is presented to solve the high dimensional transient Fokker–Planck equations (FPE) arising in the simulation of polymeric fluids via multi-bead-spring (MBS) model. Generalizing the authors' previous work on the stationary FPE, the transient solution is obtained in a single CANDECOMP/PARAFAC decomposition (CPD) form for all times via the alternating least squares algorithm. This is accomplished by treating the temporal dimension in the same manner as all other spatial dimensions, thereby decoupling it from them. As a result, the transient solution is obtained without resorting to expensive time stepping schemes. A new, relaxed approach for imposing the vanishing boundary conditions is proposed, improving the quality of the approximation. The asymptotic behavior of the temporal basis functions is studied. The proposed solver scales very well with the dimensionality of the MBS model. Numerical results for systems up to 14 dimensional state space are successfully obtained on a regular personal computer and compared with the corresponding matrix Riccati differential equation (for linear models) or Monte Carlo simulations (for nonlinear models). [ABSTRACT FROM AUTHOR]

In this paper, we apply positivity-preserving local discontinuous Galerkin (LDG) methods to solve parabolic equations with blow-up solutions. This model is commonly used in combustion problems. However, previous numerical methods are mainly based on a second order finite difference method. This is because the positivity-preserving property can hardly be satisfied for high-order ones, leading to incorrect blow-up time and blow-up sets. Recently, we have applied discontinuous Galerkin methods to linear hyperbolic equations involving δ -singularities and obtained good approximations. For nonlinear problems, some special limiters are constructed to capture the singularities precisely. We will continue this approach and study parabolic equations with blow-up solutions. We will construct special limiters to keep the positivity of the numerical approximations. Due to the Dirichlet boundary conditions, we have to modify the numerical fluxes and the limiters used in the schemes. Numerical experiments demonstrate that our schemes can capture the blow-up sets, and high-order approximations yield better numerical blow-up time. [ABSTRACT FROM AUTHOR]

In this paper we construct a new difference analog of the Caputo fractional derivative (called the L 2 - 1 σ formula). The basic properties of this difference operator are investigated and on its basis some difference schemes generating approximations of the second and fourth order in space and the second order in time for the time fractional diffusion equation with variable coefficients are considered. Stability of the suggested schemes and also their convergence in the grid L 2 -norm with the rate equal to the order of the approximation error are proved. The obtained results are supported by the numerical calculations carried out for some test problems. [ABSTRACT FROM AUTHOR]

Fourier extensions have been shown to be an effective means for the approximation of smooth, nonperiodic functions on bounded intervals given their values on an equispaced, or in general, scattered grid. Related to this method are two parameters. These are the extension parameter T (the ratio of the size of the extended domain to the physical domain) and the oversampling ratio η (the number of sampling nodes per Fourier mode). The purpose of this paper is to investigate how the choice of these parameters affects the accuracy and stability of the approximation. Our main contribution is to document the following interesting phenomenon: namely, if the desired condition number of the algorithm is fixed in advance, then the particular choice of such parameters makes little difference to the algorithm's accuracy. As a result, one is free to choose T without concern that it is suboptimal. In particular, one may use the value T=2 - which corresponds to the case where the extended domain is precisely twice the size of the physical domain - for which there is known to be a fast algorithm for computing the approximation. In addition, we determine the resolution power (points-per-wavelength) of the approximation to be equal to Tη, and address the trade-off between resolution power and stability. [ABSTRACT FROM AUTHOR]

Abstract: In this paper we consider the numerical solution of the Stokes and Darcy coupled equations, which frequently appears in porous media modeling. The main contribution of this work is as follows: First, we introduce a new formulation for the Stokes/Darcy coupled equations, subject respectively to the Beavers–Joseph–Saffman interface condition and an alternative matching interface condition. Secondly, we prove the well-posedness of these weak problems by using the classical saddle point theory. Thirdly, some spectral approximations to the weak problems are proposed and analyzed, and some error estimates are provided. It is found that the new formulations significantly simplify the error analysis and numerical implementation. Finally, some two-dimensional spectral and spectral element numerical examples are provided to demonstrate the efficiency of our methods. [Copyright &y& Elsevier]

Abstract: In this paper, we extend a class of the Hermite weighted essentially non-oscillatory (HWENO) schemes for solving the Hamilton–Jacobi equations by Qiu and Shu (2005) [24] to two dimensional unstructured meshes. The idea of the reconstruction in the HWENO schemes comes from the original WENO schemes, however both the function and its first two derivative values are evolved via time advancing and used in the reconstructions, while only the function values are evolved and used in the original WENO schemes which are nodal based approximations. The third and fourth order HWENO schemes using the combinations of second order approximations with nonlinear weights and TVD Runge–Kutta time discretization method are used here. Comparing with the original WENO schemes for Hamilton–Jacobi equations, one major advantage of HWENO schemes presented here is its compactness in the reconstructions. Extensive numerical tests are performed to illustrate the capability and high order accuracy of the methodologies. [Copyright &y& Elsevier]

Abstract: In this paper, we propose efficient numerical methods for computing ground states of spin-1 Bose–Einstein condensates (BECs) with/without the Ioffe–Pritchard magnetic field . When , a numerical method is introduced to compute the ground states and it is also applied to study properties of ground states. Numerical results suggest that the densities of components in ground states are identical for any nonzero . In particular, if is a constant, the ground states satisfy the single-mode approximation. When , efficient and simpler numerical methods are presented to solve the ground states of spin-1 BECs based on their ferromagnetic/antiferromagnetic characterizations. Numerical simulations show that our methods are more efficient than those in the literature. In addition, some conjectures are made from our numerical observations. [Copyright &y& Elsevier]

Abstract: In the field of uncertainty quantification (UQ), epistemic uncertainty often refers to the kind of uncertainty whose complete probabilistic description is not available, largely due to our lack of knowledge about the uncertainty. Quantification of the impacts of epistemic uncertainty is naturally difficult, because most of the existing stochastic tools rely on the specification of the probability distributions and thus do not readily apply to epistemic uncertainty. And there have been few studies and methods to deal with epistemic uncertainty. A recent work can be found in [J. Jakeman, M. Eldred, D. Xiu, Numerical approach for quantification of epistemic uncertainty, J. Comput. Phys. 229 (2010) 4648–4663], where a framework for numerical treatment of epistemic uncertainty was proposed. The method is based on solving an encapsulation problem, without using any probability information, in a hypercube that encapsulates the unknown epistemic probability space. If more probabilistic information about the epistemic variables is known a posteriori, the solution statistics can then be evaluated at post-process steps. In this paper, we present a new method, similar to that of Jakeman et al. but significantly extending its capabilities. Most notably, the new method (1) does not require the encapsulation problem to be in a bounded domain such as a hypercube; (2) does not require the solution of the encapsulation problem to converge point-wise. In the current formulation, the encapsulation problem could reside in an unbounded domain, and more importantly, its numerical approximation could be sought in norm. These features thus make the new approach more flexible and amicable to practical implementation. Both the mathematical framework and numerical analysis are presented to demonstrate the effectiveness of the new approach. [Copyright &y& Elsevier]

Abstract: The steady-state simplified () approximations to the linear Boltzmann equation have been proven to be asymptotically higher-order corrections to the diffusion equation in certain physical systems. In this paper, we present an asymptotic analysis for the time-dependent simplified equations up to . Additionally, equations of arbitrary order are derived in an ad hoc way. The resulting equations are hyperbolic and differ from those investigated in a previous work by some of the authors. In two space dimensions, numerical calculations for the and equations are performed. We simulate neutron distributions of a moving rod and present results for a benchmark problem, known as the checkerboard problem. The equations are demonstrated to yield significantly more accurate results than diffusion approximations. In addition, for sufficiently low values of N, they are shown to be more efficient than models of comparable cost. [Copyright &y& Elsevier]

Abstract: A continuum model for self-organized dynamics is numerically investigated. The model describes systems of particles subject to alignment interaction and short-range repulsion. It consists of a non-conservative hyperbolic system for the density and velocity orientation. Short-range repulsion is included through a singular pressure which becomes infinite at the jamming density. The singular limit of infinite pressure stiffness leads to phase transitions from compressible to incompressible dynamics. The paper proposes an Asymptotic-Preserving scheme which takes care of the singular pressure while preventing the breakdown of the Courant–Friedrichs–Lewy (CFL) stability condition near congestion. It relies on a relaxation approximation of the system and an elliptic formulation of the pressure equation. Numerical simulations of impinging clusters show the efficiency of the scheme to treat congestions. A two-fluid variant of the model provides a model of path formation in crowds. [Copyright &y& Elsevier]

Abstract: In this paper we are aimed to solve the non-dissipative Degasperis–Procesi equation based on the formulation. To resolve the computational difficulty at the wave crest where the first-order derivative may diverge and the shockpeakon solution may form, the first-order spatial derivative term in the two-step equations will be approximated in a conservative form. The resulting equations will be approximated by the symplecticity-preserving time-stepping scheme and the spatial discretization scheme that can optimize the numerical wavenumber for the first-order spatial derivative term. This scheme will be developed in a three-point grid stencil with the accuracy order of seventh within the combined compact finite difference framework. Besides the validation of numerical accuracy, we will in particular address the discrete conservation of Hamiltonians even when peakon collides with antipeakon and generates, as a result, a shockpeakon. We will also demonstrate the capability of applying the proposed numerical method to sharply resolve some important features of the third-order dispersive DP equation. [Copyright &y& Elsevier]

Abstract: In this paper, we present a fast and efficient numerical method to solve a two-dimensional fractional evolution equation on a finite domain. This numerical method combines the alternating direction implicit (ADI) approach with the second-order difference quotient in space, the backward Euler in time and order one convolution quadrature approximating the integral term. By using the discrete energy method, we prove that the ADI scheme is unconditionally stable and the numerical solution converges to the exact one with order , where k is the temporal grid size and are spatial grid sizes in the x and y directions, respectively. Two numerical examples with known exact solution are also presented, and the behavior of the error is analyzed to verify the order of convergence of the ADI-Euler method. [Copyright &y& Elsevier]

Abstract: Numerical schemes to approximate the Cahn–Hilliard equation have been widely studied in recent times due to its connection with many physically motivated problems. In this work we propose two type of linear schemes based on different ways to approximate the double-well potential term. The first idea developed in the paper allows us to design a linear numerical scheme which is optimal from the numerical dissipation point of view meanwhile the second one allows us to design unconditionally energy-stable linear schemes (for a modified energy). We present first and second order in time linear schemes to approximate the CH problem, detailing their advantages over other linear schemes that have been previously introduced in the literature. Furthermore, we compare all the schemes through several computational experiments. [Copyright &y& Elsevier]

Abstract: With the ultimate goal of devising effective absorbing boundary conditions (ABCs) for general elastic media, we investigate the accuracy aspects of local ABCs designed for untilted non-elliptic anisotropy in the frequency domain (time-harmonic analysis). While simple space–time transformations are available to treat the wavemodes with opposing phase and group velocities present in elliptic anisotropic media, no such transformations are known to exist for the case of non-elliptic anisotropy. In this paper, we use the concept of layer groupings along with an unconventional stretching of the finite element mesh to guarantee the accuracy of local ABCs designed to treat all propagating wavemodes, even those with opposing phase and group velocities. The local ABC used here is the perfectly matched discrete layer (PMDL) which is a simple variant of perfectly matched layers (PMLs) that is also equivalent to rational approximation-based local ABCs (rational ABCs); it inherits the straightforward approximation properties of rational ABCs along with the versatility of PML. The approximation properties of PMDL quantified through its reflection matrix allow us to (a) show that it is impossible to design an accurate PMDL with wavenumber-independent parameters, (b) theoretically demonstrate the ability of wavenumber-dependent parameters to ensure accuracy, and finally (c) design a practical though unconventional stretching of the finite element PMDL mesh that facilitates the implementation of wavenumber-dependent parameters. The validity of this work is demonstrated through a series of numerical experiments. [Copyright &y& Elsevier]

Abstract: The mixed boundary value problem of the Laplace equation is considered. The method of fundamental solutions (MFS) approximates the exact solution to the Laplace equation by a linear combination of independent fundamental solutions with different source points. The accuracy of the numerical solution depends on the distribution of source points. In this paper, a weighted greedy QR decomposition (GQRD) is proposed to choose significant source points by introducing a weighting parameter. An index called an average degree of approximation is defined to show the efficiency of the proposed method. From numerical experiments, it is concluded that the numerical solution tends to be more accurate when the average degree of approximation is larger, and that the proposed method can yield more accurate solutions with a less number of source points than the conventional GQRD. [Copyright &y& Elsevier]

Abstract: In this paper, we present the development of new explicit group relaxation methods which solve the two dimensional second order hyperbolic telegraph equation subject to specific initial and Dirichlet boundary conditions. The explicit group methods use small fixed group formulations derived from a combination of the rotated five-point finite difference approximation together with the centered five-point centered difference approximation on different grid spacings. The resulting schemes involve three levels finite difference approximations with second order accuracies. Analyses are presented to confirm the unconditional stability of the difference schemes. Numerical experimentations are also conducted to compare the new methods with some existing schemes. [Copyright &y& Elsevier]

Abstract: Data assimilation is an essential tool for predicting the behavior of real physical systems given approximate simulation models and limited observations. For many complex systems, there may exist several models, each with different properties and predictive capabilities. It is desirable to incorporate multiple models into the assimilation procedure in order to obtain a more accurate prediction of the physics than any model alone can provide. In this paper, we propose a framework for conducting sequential data assimilation with multiple models and sources of data. The assimilated solution is a linear combination of all model predictions and data. One notable feature is that the combination takes the most general form with matrix weights. By doing so the method can readily utilize different weights in different sections of the solution state vectors, allow the models and data to have different dimensions, and deal with the case of a singular state covariance. We prove that the proposed assimilation method, termed direct assimilation, minimizes a variational functional, a generalized version of the one used in the classical Kalman filter. We also propose an efficient iterative assimilation method that assimilates two models at a time until all models and data are assimilated. The mathematical equivalence of the iterative method and the direct method is established. Numerical examples are presented to demonstrate the effectiveness of the new method. [Copyright &y& Elsevier]

Abstract: Given local frequency domain wave data, the Numerical Micro-Local Analysis (NMLA) method (Benamou et al., 2004) and its recent improved version (Benamou et al., 2011) gives a pointwise numerical approximation of the number of rays, their slowness vectors and corresponding wavefront curvatures. With time domain wave data and assuming the source wavelet is given, the method also estimates the travel-time. The paper provides a non technical presentation of the improved NMLA algorithm and presents a numerical application which can be interpreted as a high frequency asymptotic version of the classical time reversal method (Borcea et al., 2003) . A detailed technical presentation of the algorithm is available in Benamou et al. (2011) and more numerical experiments can be found in Collino and Marmorat (2011) . [Copyright &y& Elsevier]

Abstract: The aim of this paper is the numerical simulation of compressible hydrodynamic strong implosions, which take place for instance in Inertial Confinement Fusion. It focuses in particular on two crucial issues, for such calculations: very low CFL number at the implosion center and approximation error on the initial geometry. Thus, we propose an exceptional points method, which is translation invariant and suitable to curved meshes. This method is designed for cell-centered Lagrangian schemes (GLACE , EUCCLHYD ). Several numerical examples on significant test cases are presented to show the relevance of our approach. [Copyright &y& Elsevier]

Abstract: In this paper, we suggest a very simple and effective post-processing procedure to increase the order of accuracy in time for numerical results obtained by the trapezoidal rule. We first derive a new exact, closed-form, a-priori error estimator for time integration of linear elastodynamics equations by the trapezoidal rule with non-uniform time increments. Based on this error estimator, we suggest a new post-processing procedure (containing additional time integration of elastodynamics equations by the trapezoidal rule with few time increments) that systematically improves the order of accuracy of numerical results, with the increase in the number of additional time increments used for post-processing. For example, the use of just one additional time increment for post-processing after time integration with any number of uniform time increments, renders the order of accuracy of numerical results equal to 10/3. Numerical examples of the application of the new techniques to a system with a single degree of freedom and to a multi-degree system confirm the corresponding increase in the order of convergence of numerical results after post-processing. Because the same trapezoidal rule is used for basic computations and post-processing, the new technique retains all of the properties of the trapezoidal rule, requires no writing of a new computer program for its implementation, and can be easily used with any existing commercial and research codes for elastodynamics. [Copyright &y& Elsevier]

Abstract: High-order compact finite difference scheme with operator splitting technique for solving two-dimensional time fractional diffusion equation is considered in this paper. A Grünwald–Letnikov approximation is used for the Riemann–Liouville time derivative, and the second order spatial derivatives are approximated by the compact finite differences to obtain a fully discrete implicit scheme. Alternating direction implicit (ADI) method is used to split the original problem into two separate one-dimensional problems. The local truncation error is analyzed and the stability is discussed by the Fourier method. The proposed scheme is suitable when the order of the time fractional derivative γ lies in the interval . A correction term is added to maintain high accuracy when . Numerical results are provided to verify the accuracy and efficiency of the proposed algorithm. [Copyright &y& Elsevier]

Abstract: In this paper, we introduce a fictitious dynamics for describing the only fast relaxation of a stiff ordinary differential equation (ODE) system towards a stable low-dimensional invariant manifold in the phase-space (slow invariant manifold – SIM). As a result, the demanding problem of constructing SIM of any dimensions is recast into the remarkably simpler task of solving a properly devised ODE system by stiff numerical schemes available in the literature. In the same spirit, a set of equations is elaborated for local construction of the fast subspace, and possible initialization procedures for the above equations are discussed. The implementation to a detailed mechanism for combustion of hydrogen and air has been carried out, while a model with the exact Chapman–Enskog solution of the invariance equation is utilized as a benchmark. [Copyright &y& Elsevier]

Abstract: Accurate computation of time-dependent well bore pressure is important in well test analysis – a branch of petroleum engineering where reservoir properties are estimated by comparing measured pressure responses at an oil well to results from a mathematical model. Similar methods are also used in groundwater engineering. In this paper we present the new approach of decoupled overlapping grids for accurately computing time-dependent pressure at the oil well. Our method is implemented in two stages: a global stage with a simple point or line source well approximation, and a local post-process stage with the well modeled correctly as an internal boundary. We investigate the accuracy of our method for a representative 2D problem in both homogeneous and heterogeneous isotropic domains, and compare our results with the widely used Peaceman well index solution (in the homogeneous case), and the approximate solution on locally refined grids. We also present a theoretical analysis that explains the observed behavior of the error in our method for the homogeneous case. [Copyright &y& Elsevier]

Abstract: The bidomain model is widely used in electro-cardiology to simulate spreading of excitation in the myocardium and electrocardiograms. It consists of a system of two parabolic reaction diffusion equations coupled with an ODE system. Its discretisation displays an ill-conditioned system matrix to be inverted at each time step: simulations based on the bidomain model therefore are associated with high computational costs. In this paper we propose a preconditioning for the bidomain model either for an isolated heart or in an extended framework including a coupling with the surrounding tissues (the torso). The preconditioning is based on a formulation of the discrete problem that is shown to be symmetric positive semi-definite. A block LU decomposition of the system together with a heuristic approximation (referred to as the monodomain approximation) are the key ingredients for the preconditioning definition. Numerical results are provided for two test cases: a 2D test case on a realistic slice of the thorax based on a segmented heart medical image geometry, a 3D test case involving a small cubic slab of tissue with orthotropic anisotropy. The analysis of the resulting computational cost (both in terms of CPU time and of iteration number) shows an almost linear complexity with the problem size, i.e. of type nlog α (n) (for some constant α) which is optimal complexity for such problems. [Copyright &y& Elsevier]

Abstract: Numerical simulation of metamaterials has attracted more and more attention since 2000, after the first metamaterial with negative refraction index was successfully constructed. In this paper we construct a fully-discrete leap-frog type finite element scheme to solve the three-dimensional time-dependent Maxwell’s equations when metamaterials are involved. First, we obtain some superclose results between the interpolations of the analytical solutions and finite element solutions obtained using arbitrary orders of Raviart–Thomas–Nédélec mixed spaces on regular cubic meshes. Then we prove the superconvergence result in the discrete l 2 norm achieved for the lowest-order Raviart–Thomas–Nédélec space. To our best knowledge, such superconvergence results have never been obtained elsewhere. Finally, we implement the leap-frog scheme and present numerical results justifying our theoretical analysis. [Copyright &y& Elsevier]

Abstract: In this paper, we propose a composite preconditioner for the electromagnetic scattering from a large cavity. The electromagnetic cavity problem is described by the Helmholtz equation with a nonlocal boundary condition on the aperture of the cavity and Dirichlet (or Neumann) boundary conditions on the walls of the cavity. The preconditioner proposed here combines the optimal sine transform based approximation with a layered medium model. Using fast Fourier transforms, the computational cost of every iteration is on an N × N uniform partition of the unit square. Numerical results for a model problem show that the new preconditioner is more efficient than those recently considered in the literature. For the cavity with a small portion of non-layered media, we propose a sparse preconditioned conjugate orthogonal conjugate gradient solver combined with the new preconditioner. Numerical results for a model problem are reported to demonstrate the efficiency of the sparse solver. [Copyright &y& Elsevier]

Abstract: We present in this paper a new kernel-independent fast multipole method (FMM), named as FKI-FMM, for pairwise particle interactions with translation-invariant kernel functions. FKI-FMM creates, using numerical techniques, sufficiently accurate and compressive representations of a given kernel function over multi-scale interaction regions in the form of a truncated Fourier series. It provides also economic operators for the multipole-to-multipole, multipole-to-local, and local-to-local translations that are typical and essential in the FMM algorithms. The multipole-to-local translation operator, in particular, is readily diagonal and does not dominate in arithmetic operations. FKI-FMM provides an alternative and competitive option, among other kernel-independent FMM algorithms, for an efficient application of the FMM, especially for applications where the kernel function consists of multi-physics and multi-scale components as those arising in recent studies of biological systems. We present the complexity analysis and demonstrate with experimental results the FKI-FMM performance in accuracy and efficiency. [Copyright &y& Elsevier]

Abstract: This paper presents a new compact approximation method for the discretisation of second-order elliptic equations in one and two dimensions. The problem domain, which can be rectangular or non-rectangular, is represented by a Cartesian grid. On stencils, which are three nodal points for one-dimensional problems and nine nodal points for two-dimensional problems, the approximations for the field variable and its derivatives are constructed using integrated radial basis functions (IRBFs). Several pieces of information about the governing differential equation on the stencil are incorporated into the IRBF approximations by means of the constants of integration. Numerical examples indicate that the proposed technique yields a very high rate of convergence with grid refinement. [Copyright &y& Elsevier]

Abstract: In this paper, a compact finite difference scheme for the fractional sub-diffusion equations is derived. After a transformation of the original problem, the L1 discretization is applied for the time-fractional part and fourth-order accuracy compact approximation for the second-order space derivative. The unique solvability of the difference solution is discussed. The stability and convergence of the finite difference scheme in maximum norm are proved using the energy method, where a new inner product is introduced for the theoretical analysis. The technique is quite novel and different from previous analytical methods. Finally, a numerical example is provided to show the effectiveness and accuracy of the method. [ABSTRACT FROM AUTHOR]

Abstract: Techniques for rapidly computing approximate solutions to elliptic PDEs such as Laplace’s equation are well established. For problems involving general domains, and operators with constant coefficients, a highly efficient approach is to rewrite the boundary value problem as a Boundary Integral Equation (BIE), and then solve the BIE using fast methods such as, e.g., the Fast Multipole Method (FMM). The current paper demonstrates that this procedure can be extended to elliptic difference equations defined on infinite lattices, or on finite lattice with boundary conditions of either Dirichlet or Neumann type. As a representative model problem, a lattice equivalent of Laplace’s equation on a square lattice in two dimensions is considered: discrete analogs of BIEs are derived and fast solvers analogous to the FMM are constructed. Fast techniques are also constructed for problems involving lattices with inclusions and local deviations from perfect periodicity. The complexity of the methods described is O(N boundary + N source + N inc) where N boundary is the number of nodes on the boundary of the domain, N source is the number of nodes subjected to body loads, and N inc is the number of nodes that deviate from perfect periodicity. This estimate should be compared to the O(N domainlogN domain) estimate for FFT based methods, where N domain is the total number of nodes in the lattice (so that in two dimensions, ). Several numerical examples are presented. [ABSTRACT FROM AUTHOR]

Abstract: In the following paper, we discuss the exhaustive use and implementation of stabilization finite element methods for the resolution of the 3D time-dependent incompressible Navier–Stokes equations. The proposed method starts by the use of a finite element variational multiscale (VMS) method, which consists in here of a decomposition for both the velocity and the pressure fields into coarse/resolved scales and fine/unresolved scales. This choice of decomposition is shown to be favorable for simulating flows at high Reynolds number. We explore the behaviour and accuracy of the proposed approximation on three test cases. First, the lid-driven square cavity at Reynolds number up to 50,000 is compared with the highly resolved numerical simulations and second, the lid-driven cubic cavity up to Re=12,000 is compared with the experimental data. Finally, we study the flow over a 2D backward-facing step at Re=42,000. Results show that the present implementation is able to exhibit good stability and accuracy properties for high Reynolds number flows with unstructured meshes. [ABSTRACT FROM AUTHOR]

Abstract: Fractional diffusion equations model phenomena exhibiting anomalous diffusion that can not be modeled accurately by the second-order diffusion equations. Because of the nonlocal property of fractional differential operators, the numerical methods have full coefficient matrices which require storage of O(N 2) and computational cost of O(N 3) where N is the number of grid points. In this paper we develop a fast finite difference method for fractional diffusion equations, which only requires storage of O(N) and computational cost of O(N log2 N) while retaining the same accuracy and approximation property as the regular finite difference method. Numerical experiments are presented to show the utility of the method. [ABSTRACT FROM AUTHOR]

Abstract: This paper investigates a novel approximate Bayesian inference procedure for numerically solving inverse problems. A hierarchical formulation which determines automatically the regularization parameter and the noise level together with the inverse solution is adopted. The framework is of variational type, and it can deliver the inverse solution and regularization parameter together with their uncertainties calibrated. It approximates the posteriori probability distribution by separable distributions based on Kullback–Leibler divergence. Two approximations are derived within the framework, and some theoretical properties, e.g. variance estimate and consistency, are also provided. Algorithms for their efficient numerical realization are described, and their convergence properties are also discussed. Extensions to nonquadratic regularization/nonlinear forward models are also briefly studied. Numerical results for linear and nonlinear Cauchy-type problems arising in heat conduction with both smooth and nonsmooth solutions are presented for the proposed method, and compared with that by Markov chain Monte Carlo. The results illustrate that the variational method can faithfully capture the posteriori distribution in a computationally efficient way. [ABSTRACT FROM AUTHOR]

Abstract: We propose a new interpolation technique for the CIP method applied to curvilinear coordinates. The CIP method can hardly maintain third-order accuracy on curvilinear coordinates. The reason for the degeneracy in accuracies has not been discussed in detail. This paper reveals the problems of the CIP method on curvilinear coordinates and presents an improved CIP method to solve the advection equation accurately. The features of the presented method are: (1) the metric computation on the upwind stencil is defined in the same manner as in the advection phase of the CIP method; and (2) gradient values in the physical domain in the computation on the curvilinear coordinates are used. Various test problems show that the improved CIP method has approximate third-order accuracy. [Copyright &y& Elsevier]

Abstract: In this paper, we construct spatially consistent explicit second order discretizations for time dependent hyperbolic problems, starting from a given residual distribution (RD) discrete approximation of the steady operator. We review the existing knowledge on consistent RD mass matrices and highlight the relations between different definitions. We then introduce our explicit approach which is based on three main ingredients: first recast the RD discretization as a stabilized Galerkin scheme, then use a shifted time discretization in the stabilization operator, and lastly apply high order mass lumping on the Galerkin component of the discretization. The discussion is particularly relevant for schemes of the residual distribution type which we will use for all our numerical experiments. However, similar ideas can be used in the context of residual-based finite volume discretizations such as the ones proposed in . The schemes are tested on a wide variety of classical problems confirming the theoretical expectations. [Copyright &y& Elsevier]

Abstract: A new PDE solver was introduced recently, in Part I of this two-paper sequence, on the basis of two main concepts: the well-known Alternating Direction Implicit (ADI) approach, on one hand, and a certain “Fourier Continuation” (FC) method for the resolution of the Gibbs phenomenon, on the other. Unlike previous alternating direction methods of order higher than one, which only deliver unconditional stability for rectangular domains, the new high-order FC-AD (Fourier-Continuation Alternating-Direction) algorithm yields unconditional stability for general domains—at an cost per time-step for an N point spatial discretization grid. In the present contribution we provide an overall theoretical discussion of the FC-AD approach and we extend the FC-AD methodology to linear hyperbolic PDEs. In particular, we study the convergence properties of the newly introduced FC(Gram) Fourier Continuation method for both approximation of general functions and solution of the alternating-direction ODEs. We also present (for parabolic PDEs on general domains, and, thus, for our associated elliptic solvers) a stability criterion which, when satisfied, ensures unconditional stability of the FC-AD algorithm. Use of this criterion in conjunction with numerical evaluation of a series of singular values (of the alternating-direction discrete one-dimensional operators) suggests clearly that the fifth-order accurate class of parabolic and elliptic FC-AD solvers we propose is indeed unconditionally stable for all smooth spatial domains and for arbitrarily fine discretizations. To illustrate the FC-AD methodology in the hyperbolic PDE context, finally, we present an example concerning the Wave Equation—demonstrating sixth-order spatial and fourth-order temporal accuracy, as well as a complete absence of the debilitating “dispersion error”, also known as “pollution error”, that arises as finite-difference and finite-element solvers are applied to solution of wave propagation problems. [Copyright &y& Elsevier]

Abstract: In this paper, we present a method of decomposing a highly oscillatory wave field into a sparse superposition of Gaussian beams. The goal is to extract the necessary parameters for a Gaussian beam superposition from this wave field, so that further evolution of the high frequency waves can be computed by the method of Gaussian beams. The methodology is described for with numerical examples for . In the first example, a field generated by an interface reflection of Gaussian beams is decomposed into a superposition of Gaussian beams. The beam parameters are reconstructed to a very high accuracy. The data in the second example is not a superposition of a finite number of Gaussian beams. The wave field to be approximated is generated by a finite difference method for a geometry with two slits. The accuracy in the decomposition increases monotonically with the number of beams. [Copyright &y& Elsevier]

Abstract: High-order compact finite difference scheme for solving one-dimensional fractional diffusion equation is considered in this paper. After approximating the second-order derivative with respect to space by the compact finite difference, we use the Grünwald–Letnikov discretization of the Riemann–Liouville derivative to obtain a fully discrete implicit scheme. We analyze the local truncation error and discuss the stability using the Fourier method, then we prove that the compact finite difference scheme converges with the spatial accuracy of fourth order using matrix analysis. Numerical results are provided to verify the accuracy and efficiency of the proposed algorithm. [Copyright &y& Elsevier]