A spectral collocation algorithm for two-point boundary value problem in fiber Raman amplifier equations

A novel algorithm implementing Chebyshev spectral collocation (pseudospectral) method in combination with Newton’s method is proposed for the nonlinear two-point boundary value problem (BVP) arising in solving propagation equations in fiber Raman amplifier. Moreover, an algorithm to train the known linear solution for use as a starting solution for the Newton iteration is proposed and successfully implemented. The exponential accuracy obtained by the proposed Chebyshev pseudospectral method is demonstrated on a case of the Raman propagation equations with strong nonlinearities. This is in contrast to algebraic accuracy obtained by typical solvers used in the literature. The resolving power and the efficiency of the underlying Chebyshev grid are demonstrated in comparison to a known BVP solver.     

A particle–particle hybrid method for kinetic and continuum equations

We present a coupling procedure for two different types of particle methods for the Boltzmann and the Navier–Stokes equations. A variant of the DSMC method is applied to simulate the Boltzmann equation, whereas a meshfree Lagrangian particle method, similar to the SPH method, is used for simulations of the Navier–Stokes equations. An automatic domain decomposition approach is used with the help of a continuum breakdown criterion. We apply adaptive spatial and time meshes. The classical Sod’s 1D shock tube problem is solved for a large range of Knudsen numbers. Results from Boltzmann, Navier–Stokes and hybrid solvers are compared. The CPU time for the hybrid solver is 3–4 times faster than for the Boltzmann solver.     

A Direct Spectral Collocation Poisson Solver in Polar and Cylindrical Coordinates

In this paper, we present a direct spectral collocation method for the solution of the Poisson equation in polar and cylindrical coordinates. The solver is applied to the Poisson equations for several different domains including a part of a disk, an annulus, a unit disk, and a cylinder. Unlike other Poisson solvers for geometries such as unit disks and cylinders, no pole condition is involved for the present solver. The method is easy to implement, fast, and gives spectral accuracy. We also use the weighted interpolation technique and nonclassical collocation points to improve the convergence.     

ADI spectral collocation methods for parabolic problems

We discuss the Crank–Nicolson and Laplace modified alternating direction implicit Legendre and Chebyshev spectral collocation methods for a linear, variable coefficient, parabolic initial-boundary value problem on a rectangular domain with the solution subject to non-zero Dirichlet boundary conditions. The discretization of the problems by the above methods yields matrices which possess banded structures. This along with the use of fast Fourier transforms makes the cost of one step of each of the Chebyshev spectral collocation methods proportional, except for a logarithmic term, to the number of the unknowns. We present the convergence analysis for the Legendre spectral collocation methods in the special case of the heat equation. Using numerical tests, we demonstrate the second order accuracy in time of the Chebyshev spectral collocation methods for general linear variable coefficient parabolic problems.     

Numerical methods for stochastic partial differential equations with multiple scales

A new method for solving numerically stochastic partial differential equations (SPDEs) with multiple scales is presented. The method combines a spectral method with the heterogeneous multiscale method (HMM) presented in [W. E, D. Liu, E. Vanden-Eijnden, Analysis of multiscale methods for stochastic differential equations, Commun. Pure Appl. Math., 58(11) (2005) 1544–1585]. The class of problems that we consider are SPDEs with quadratic nonlinearities that were studied in [D. Blömker, M. Hairer, G.A. Pavliotis, Multiscale analysis for stochastic partial differential equations with quadratic nonlinearities, Nonlinearity, 20(7) (2007) 1721–1744]. For such SPDEs an amplitude equation which describes the effective dynamics at long time scales can be rigorously derived for both advective and diffusive time scales. Our method, based on micro and macro solvers, allows to capture numerically the amplitude equation accurately at a cost independent of the small scales in the problem. Numerical experiments illustrate the behavior of the proposed method.     

Numerical approximations for fractional diffusion equations via a Chebyshev spectral-tau method

In this paper, a class of fractional diffusion equations with variable coefficients is considered. An accurate and efficient spectral tau technique for solving the fractional diffusion equations numerically is proposed. This method is based upon Chebyshev tau approximation together with Chebyshev operational matrix of Caputo fractional differentiation. Such approach has the advantage of reducing the problem to the solution of a system of algebraic equations, which may then be solved by any standard numerical technique. We apply this general method to solve four specific examples. In each of the examples considered, the numerical results show that the proposed method is of high accuracy and is efficient for solving the time-dependent fractional diffusion equations.     

Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods

A robust methodology is presented for efficiently solving partial differential equations using Chebyshev spectral techniques. It is well known that differential equations in one dimension can be solved efficiently with Chebyshev discretizations, O(N) operations for N unknowns, however this efficiency is lost in higher dimensions due to the coupling between modes. This paper presents the “quasi-inverse“ technique (QIT), which combines optimizations of one-dimensional spectral differentiation matrices with Kronecker matrix products to build efficient multi-dimensional operators. This strategy results in O(N^2D?1) operations for N^D unknowns, independent of the form of the differential operators. QIT is compared to the matrix diagonalization technique (MDT) of Haidvogel and Zang [D.B. Haidvogel, T. Zang, The accurate solution of Poisson’s equation by expansion in Chebyshev polynomials, J. Comput. Phys. 30 (1979) 167–180] and Shen [J. Shen, Efficient spectral-Galerkin method. II. Direct solvers of second- and fourth-order equations using Chebyshev polynomials, SIAM J. Sci. Comp. 16 (1) (1995) 74–87]. While the cost for MDT and QIT are the same in two dimensions, there are significant differences. MDT utilizes an eigenvalue/eigenvector decomposition and can only be used for relatively simple differential equations. QIT is based upon intrinsic properties of the Chebyshev polynomials and is adaptable to linear PDEs with constant coefficients in simple domains. We present results for a standard suite of test problems, and discuss of the adaptability of QIT to more complicated problems.     

磁流体方腔槽道流整体线性稳定性数值方法研究

将Chebyshev谱配置法和基于非均匀网格的高阶FD-q差分格式运用于磁流体方腔槽道流整体线性稳定性研究,比较两类数值方法的优缺点.Chebyshev谱配置法收敛快且精度高,但需要构造非常庞大的满矩阵,存储量和计算开销巨大;高阶FD-q差分格式采用了基于Kosloff-Tal-Ezer变换的Chebyshev谱配置点作为离散网格,在保持较高网格收敛精度的同时,差分格式可以采用稀疏矩阵进行存储,显著降低了存储量和计算开销.相比传统的谱配置法,基于非均匀网格的高阶FD-q差分格式计算效率得到显著的提升,将高阶FD-q差分格式运用于非正则模线性最优瞬态增长的计算,计算效果良好.     

Three-dimensional elliptic solvers for interface problems and applications

Second-order accurate elliptic solvers using Cartesian grids are presented for three-dimensional interface problems in which the coefficients, the source term, the solution and its normal flux may be discontinuous across an interface. One of our methods is designed for general interface problems with variable but discontinuous coefficient. The scheme preserves the discrete maximum principle using constrained optimization techniques. An algebraic multigrid solver is applied to solve the discrete system. The second method is designed for interface problems with piecewise constant coefficient. The method is based on the fast immersed interface method and a fast 3D Poisson solver. The second method has been modified to solve Helmholtz/Poisson equations on irregular domains. An application of our method to an inverse interface problem of shape identification is also presented. In this application, the level set method is applied to find the unknown surface iteratively.     

An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations

In recent years, there has been a growing interest in analyzing and quantifying the effects of random inputs in the solution of ordinary/partial differential equations. To this end, the spectral stochastic finite element method (SSFEM) is the most popular method due to its fast convergence rate. Recently, the stochastic sparse grid collocation method has emerged as an attractive alternative to SSFEM. It approximates the solution in the stochastic space using Lagrange polynomial interpolation. The collocation method requires only repetitive calls to an existing deterministic solver, similar to the Monte Carlo method. However, both the SSFEM and current sparse grid collocation methods utilize global polynomials in the stochastic space. Thus when there are steep gradients or finite discontinuities in the stochastic space, these methods converge very slowly or even fail to converge. In this paper, we develop an adaptive sparse grid collocation strategy using piecewise multi-linear hierarchical basis functions. Hierarchical surplus is used as an error indicator to automatically detect the discontinuity region in the stochastic space and adaptively refine the collocation points in this region. Numerical examples, especially for problems related to long-term integration and stochastic discontinuity, are presented. Comparisons with Monte Carlo and multi-element based random domain decomposition methods are also given to show the efficiency and accuracy of the proposed method.     

Finite element formulation of the first- and second-order discrete ordinates equations for radiative heat transfer calculation in three-dimensional participating media

Two finite element methods (FEMs), FEDOM1 and FEDOM2 (standing for the first and the second finite element discrete ordinates methods, respectively), are formulated and numerically tested. The reference second-order discrete equation is modified in its scattering terms and is applied to the problems of absorbing/emitting and anisotropically scattering media by using the FEM. Numerical features of the developed FEMs are compared with one of the discrete ordinates interpolation method (DOIM), which uses a finite difference scheme. Prediction results of radiative heat transfer by these two FEMs are compared with reference solutions and verified in three-dimensional enclosures containing participating media. The results of FEDOM1 and FEDOM2 agree well with exact solutions for the problem of absorbing/emitting medium with various range of optical thickness. Generally, the two FEMs show more accurate results than DOIM. And FEDOM1 shows more accurate results than FEDOM2 in most of the test problems. Both of the developed FEMs show reasonable results compared with published Monte Carlo solutions for the tested absorbing/emitting and anisotropically scattering media. Although the FEDOM2 is not as accurate as the FEDOM1, it shows its own advantages that it reduces CPU time and memory space of dependent variable to half.     

Collocation Methods for Hyperbolic Partial Differential Equations with Singular Sources

A numerical study is given on the spectral methods and the high order WENO finite difference scheme for the solution of linear and nonlinear hyperbolic partial differential equations with stationary and non-stationary singular sources. The singular source term is represented by the $δ$-function. For the approximation of the $δ$-function, the direct projection method is used that was proposed in [6]. The $δ$-function is constructed in a consistent way to the derivative operator. Nonlinear sine-Gordon equation with a stationary singular source was solved with the Chebyshev collocation method. The $δ$-function with the spectral method is highly oscillatory but yields good results with small number of collocation points. The results are compared with those computed by the second order finite difference method. In modeling general hyperbolic equations with a non-stationary singular source, however, the solution of the linear scalar wave equation with the non-stationary singular source using the direct projection method yields non-physical oscillations for both the spectral method and the WENO scheme. The numerical artifacts arising when the non-stationary singular source term is considered on the discrete grids are explained.     

Spectrally-accurate algorithm for moving boundary problems for the Navier–Stokes equations

A spectral algorithm based on the immersed boundary conditions (IBC) concept is developed for simulations of viscous flows with moving boundaries. The algorithm uses a fixed computational domain with flow domain immersed inside the computational domain. Boundary conditions along the edges of the time-dependent flow domain enter the algorithm in the form of internal constraints. Spectral spatial discretization uses Fourier expansions in the stream-wise direction and Chebyshev expansions in the normal-to-the-wall direction. Up to fourth-order implicit temporal discretization methods have been implemented. It has been demonstrated that the algorithm delivers the theoretically predicted accuracy in both time and space. Performances of various linear solvers employed in the solution process have been evaluated and a new class of solver that takes advantage of the structure of the coefficient matrix has been proposed. The new solver results in a significant acceleration of computations as well as in a substantial reduction in memory requirements.     

Efficient algebraic multigrid for migration–diffusion–convection–reaction systems arising in electrochemical simulations

The article discusses components and performance of an algebraic multigrid (AMG) preconditioner for the fully coupled multi-ion transport and reaction model (MITReM) with nonlinear boundary conditions, important for electrochemical modeling. The governing partial differential equations (PDEs) are discretized in space by a combined finite element and residual distribution method. Solution of the discrete system is obtained by means of a Newton-based nonlinear solver, and an AMG-preconditioned BICGSTAB Krylov linear solver. The presented AMG preconditioner is based on so-called point-based classical AMG. The linear solver is compared to a standard direct and several one-level iterative solvers for a range of geometries and chemical systems with scientific and industrial relevance. The results indicate that point-based AMG methods, carefully designed, are an attractive alternative to more commonly employed numerical methods for the simulation of complex electrochemical processes.     

Some theoretical aspects of elastic wave modeling with a recently developed spectral element method

A spectral element method has been recently developed for solving elastodynamic problems. The numerical solutions are obtained by using the weak formulation of the elastodynamic equation for heterogeneous media, based on the Galerkin approach applied to a partition, in small subdomains, of the original physical domain. In this work, some mathematical aspects of the method and the associated algorithm implementation are systematically investigated. Two kinds of orthogonal basis functions, constructed with Legendre and Chebyshev polynomials, and their related Gauss-Lobatto collocation points are introduced. The related integration formulas are obtained. The standard error estimations and expansion convergence are discussed. An element-by-element pre-conditioned conjugate gradient linear solver in the space domain and a staggered predictor/multi-corrector algorithm in the time integration are used for strong heterogeneous elastic media. As a consequence, neither the global matrices nor the effective force vector is assembled. When analytical formulas are used for the element quadrature, there is even no need for forming element matrix in order to further save memory without losing much in computational efficiency. The element-by-element algorithm uses an optimal tensor product scheme which makes this method much more efficient than finite-element methods from the point of view of both memory storage and computational time requirements. This work is divided into two parts. The first part mainly focuses on theoretical studies with a simple numerical result for the Che-byshev spectral element, and the second part, mainly with the Legendre spectral element, will give the algorithm implementation, numerical accuracy and efficiency analyses, and then the detailed modeling example comparisons of the proposed spectral element method with a pseudo-spectral method, which will be seen in another work by Lin, Wang and Zhang.     

Time-Domain Numerical Solutions of Maxwell Interface Problems with Discontinuous Electromagnetic Waves

This paper is devoted to time domain numerical solutions of two-dimensional (2D) material interface problems governed by the transverse magnetic (TM) and transverse electric (TE) Maxwell's equations with discontinuous electromagnetic solutions. Due to the discontinuity in wave solutions across the interface, the usual numerical methods will converge slowly or even fail to converge. This calls for the development of advanced interface treatments for popular Maxwell solvers. We will investigate such interface treatments by considering two typical Maxwell solvers – one based on collocation formulation and the other based on Galerkin formulation. To restore the accuracy reduction of the collocation finite-difference time-domain (FDTD) algorithm near an interface, the physical jump conditions relating discontinuous wave solutions on both sides of the interface must be rigorously enforced. For this purpose, a novel matched interface and boundary (MIB) scheme is proposed in this work, in which new jump conditions are derived so that the discontinuous and staggered features of electric and magnetic field components can be accommodated. The resulting MIB time-domain (MIBTD) scheme satisfies the jump conditions locally and suppresses the staircase approximation errors completely over the Yee lattices. In the discontinuous Galerkin time-domain (DGTD) algorithm – a popular Galerkin Maxwell solver, a proper numerical flux can be designed to accurately capture the jumps in the electromagnetic waves across the interface and automatically preserves the discontinuity in the explicit time integration. The DGTD solution to Maxwell interface problems is explored in this work, by considering a nodal based high order discontinuous Galerkin method. In benchmark TM and TE tests with analytical solutions, both MIBTD and DGTD schemes achieve the second order of accuracy in solving circular interfaces. In comparison, the numerical convergence of the MIBTD method is slightly more uniform, while the DGTD method is more flexible and robust.     

Automatic spectral collocation for integral,integro-differential,and integrally reformulated differential equations

Automatic Chebyshev spectral collocation methods for Fredholm and Volterra integral and integro-differential equations have been implemented as part of the chebfun software system. This system enables a symbolic syntax to be applied to numerical objects in order to pose and solve problems without explicit references to discretization. The same objects can be used in matrix-free iterative methods in linear algebra, in order to avoid very large dense matrices or allow application to problems with nonsmooth coefficients. As a further application of the ability to implement operator equations, a method of Greengard [1] for the recasting of differential equations as integral equations is generalized to mth order boundary value and generalized eigenvalue problems. In the integral form, large condition numbers associated with differentiation matrices in high-order problems are avoided. The ability to implement the recasting process generally follows from implementation of the operator expressions in chebfun. The integral method also can be extended to first-order systems, although chebfun syntax does not currently allow easy implementation in this case.     

Navier–Stokes spectral solver in a sphere

The paper presents the first implementation of a primitive variable spectral method for calculating viscous flows inside a sphere. A variational formulation of the Navier–Stokes equations is adopted using a fractional-step time discretization with the classical second-order backward difference scheme combined with explicit extrapolation of the nonlinear term. The resulting scalar and vector elliptic equations are solved by means of the direct spectral solvers developed recently by the authors. The spectral matrices for radial operators are characterized by a minimal sparsity – diagonal stiffness and tridiagonal mass matrix. Closed-form expressions of their nonzero elements are provided here for the first time, showing that the condition number of the relevant matrices grows as the second power of the truncation order. A new spectral elliptic solver for the velocity unknown in spherical coordinates is also described that includes implicitly the Coriolis force in a rotating frame, but requires a minimal coupling between the modal velocity components in the Fourier space. The numerical tests confirm that the proposed method achieves spectral accuracy and ensures infinite differentiability to all orders of the numerical solution, by construction. These results indicate that the new primitive variable spectral solver is an effective alternative to the spectral method recently proposed by Kida and Nakayama, where the velocity field is represented in terms of poloidal and toroidal functions.     

风扇/压气机数值稳定性模型

本文给出了一种用于预测轴流压气机旋转失速发生的数值稳定性模型,并描述了模型的理论基础以及求解模型方程的方法。模型的建立基于线性稳定性理论的基础之上,通过特征值的虚部来判断系统的稳定性。一种全局性方法被用于求解离散系统的所有特征值。离散的方法包括二阶精度的有限差分方法以及切比雪夫谱配置的方法。计算结果表明,用这两种方法均可以得到准确的特征值,但是收敛的速度有所不同。另外,将文中的分析与实验进行了比较,结果表明实验中的失速点可以用这种模型进行合理的预测。     

Efficient stochastic Galerkin methods for random diffusion equations

We discuss in this paper efficient solvers for stochastic diffusion equations in random media. We employ generalized polynomial chaos (gPC) expansion to express the solution in a convergent series and obtain a set of deterministic equations for the expansion coefficients by Galerkin projection. Although the resulting system of diffusion equations are coupled, we show that one can construct fast numerical methods to solve them in a decoupled fashion. The methods are based on separation of the diagonal terms and off-diagonal terms in the matrix of the Galerkin system. We examine properties of this matrix and show that the proposed method is unconditionally stable for unsteady problems and convergent for steady problems with a convergent rate independent of discretization parameters. Numerical examples are provided, for both steady and unsteady random diffusions, to support the analysis.