In this paper we present a method to generate independent samples for a general random variable, either continuous or discrete. The algorithm is an extension of the Acceptance-Rejection method, and it is particularly useful for kinetic simulation in which the rates are fluctuating in time and have singular limits, as occurs for example in simulation of recombination interactions in a plasma. Although it depends on some additional requirements, the new method is easy to implement and rejects less samples than the Acceptance-Rejection method.

We present new large time step methods for the shallow water flows in the low Froude number limit. In order to take into account multiscale phenomena that typically appear in geophysical flows nonlinear fluxes are split into a linear part governing the gravitational waves and the nonlinear advection. We propose to approximate fast linear waves implicitly in time and in space by means of a genuinely multidimensional evolution operator. On the other hand, we approximate nonlinear advection part explicitly in time and in space by means of the method of characteristics or some standard numerical flux function. Time integration is realized by the implicit-explicit (IMEX) method. We apply the IMEX Euler scheme, two step Runge Kutta Cranck Nicolson scheme, as well as the semi-implicit BDF scheme and prove their asymptotic preserving property in the low Froude number limit. Numerical experiments demonstrate stability, accuracy and robustness of these new large time step finite volume schemes with respect to small Froude number.

The formation of singularities in relativistic flows is not well understood. Smooth solutions to the relativistic Euler equations are known to have a finite lifespan; the possible breakdown mechanisms are shock formation, violation of the subluminal conditions and mass concentration. We propose a new hybrid Glimm/central-upwind scheme for relativistic flows. The scheme is used to numerically investigate, for a family of problems, which of the above mechanisms is involved.

In this work, we concern with the numerical comparison between different kinds of design points in least square (LS) approach on polynomial spaces. Such a topic is motivated by uncertainty quantification (UQ). Three kinds of design points are considered, which are the Sparse Grid (SG) points, the Monte Carlo (MC) points and the Quasi Monte Carlo (QMC) points. We focus on three aspects during the comparison: (i) the convergence properties; (ii) the stability, i.e. the properties of the resulting condition number of the design matrix; (iii) the robustness when numerical noises are present in function values. Several classical high dimensional functions together with a random ODE model are tested. It is shown numerically that (i) neither the MC sampling nor the QMC sampling introduces the low convergence rate, namely, the approach achieves high order convergence rate for all cases provided that the underlying functions admit certain regularity and enough design points are used; (ii)The use of SG points admits better convergence properties only for very low dimensional problems (say d ≤ 2); (iii)The QMC points, being deterministic, seem to be a good choice for higher dimensional problems not only for better convergence properties but also in the stability point of view.

Nonlinear filter problems arise in many applications such as communications and signal processing. Commonly used numerical simulation methods include Kalman filter method, particle filter method, etc. In this paper a novel numerical algorithm is constructed based on samples of the current state obtained by solving the state equation implicitly. Numerical experiments demonstrate that our algorithm is more accurate than the Kalman filter and more stable than the particle filter.

In this article we detail the methodology developed to construct an efficient interface description technique — the robust conservative level set (RCLS) — to simulate multiphase flows on mixed-element unstructured meshes while conserving mass to machine accuracy. The approach is tailored specifically for industry as the three-dimensional unstructured approach allows for the treatment of very complex geometries. In addition, special care has been taken to optimise the trade-off between accuracy and computational cost while maintaining the robustness of the numerical method. This was achieved by solving the transport equations for the liquid volume fraction using a WENO scheme for polyhedral meshes and by adding a flux-limiter algorithm. The performance of the resulting method has been compared against established multiphase numerical methods and its ability to capture the physics of multiphase flows is demonstrated on a range of relevant test cases. Finally, the RCLS method has been applied to the simulation of the primary breakup of a flat liquid sheet of kerosene in co-flowing high-pressure gas. This quasi-DNS/LES computation was performed at relevant aero-engine conditions on a three-dimensional mixed-element unstructured mesh. The numerical results have been validated qualitatively against theoretical predictions and experimental data. In particular, the expected breakup regime was observed in the simulation results. Finally, the computation reproduced faithfully the breakup length predicted by a correlation based on experimental data. This constitutes a first step towards a quantitative validation.

In this work, we are concerned with a time-splitting Fourier pseudospectral (TSFP) discretization for the Klein-Gordon (KG) equation, involving a dimensionless parameter ε∈(0,1]. In the nonrelativistic limit regime, the small ε produces high oscillations in exact solutions with wavelength of O(ε²) in time. The key idea behind the TSFP is to apply a time-splitting integrator to an equivalent first-order system in time, with both the nonlinear and linear subproblems exactly integrable in time and, respectively, Fourier frequency spaces. The method is fully explicit and time reversible. Moreover, we establish rigorously the optimal error bounds of a second-order TSFP for fixed ε = O(1), thanks to an observation that the scheme coincides with a type of trigonometric integrator. As the second task, numerical studies are carried out, with special efforts made to applying the TSFP in the nonrelativistic limit regime, which are geared towards understanding its temporal resolution capacity and meshing strategy for O(ε²)-oscillatory solutions when 0 < ε ≪ 1. It suggests that the method has uniform spectral accuracy in space, and an asymptotic O(ε⁻²∆t²) temporal discretization error bound (∆t refers to time step). On the other hand, the temporal error bounds for most trigonometric integrators, such as the well-established Gautschi-type integrator in [6], are O(ε⁻⁴∆t²). Thus, our method offers much better approximations than the Gautschi-type integrator in the highly oscillatory regime. These results, either rigorous or numerical, are valid for a splitting scheme applied to the classical relativistic NLS reformulation as well.

In Markov Chain Monte Carlo (MCMC) simulations, thermal equilibria quantities are estimated by ensemble average over a sample set containing a large number of correlated samples. These samples are selected in accordance with the probability distribution function, known from the partition function of equilibrium state. As the stochastic error of the simulation results is significant, it is desirable to understand the variance of the estimation by ensemble average, which depends on the sample size (i.e., the total number of samples in the set) and the sampling interval (i.e., cycle number between two consecutive samples). Although large sample sizes reduce the variance, they increase the computational cost of the simulation. For a given CPU time, the sample size can be reduced greatly by increasing the sampling interval, while having the corresponding increase in variance be negligible if the original sampling interval is very small. In this work, we report a few general rules that relate the variance with the sample size and the sampling interval. These results are observed and confirmed numerically. These variance rules are derived for the MCMC method but are also valid for the correlated samples obtained using other Monte Carlo methods. The main contribution of this work includes the theoretical proof of these numerical observations and the set of assumptions that lead to them.

An iterative discontinuous Galerkin (DG) method is proposed to solve the nonlinear Poisson Boltzmann (PB) equation. We first identify a function space in which the solution of the nonlinear PB equation is iteratively approximated through a series of linear PB equations, while an appropriate initial guess and a suitable iterative parameter are selected so that the solutions of linear PB equations are monotone within the identified solution space. For the spatial discretization we apply the direct discontinuous Galerkin method to those linear PB equations. More precisely, we use one initial guess when the Debye parameter λ=O(1), and a special initial guess for λ≪1 to ensure convergence. The iterative parameter is carefully chosen to guarantee the existence, uniqueness, and convergence of the iteration. In particular, iteration steps can be reduced for a variable iterative parameter. Both one and two-dimensional numerical results are carried out to demonstrate both accuracy and capacity of the iterative DG method for both cases of λ=O(1) and λ≪1. The (m+1)th order of accuracy for L² and mth order of accuracy for H¹ for P^m elements are numerically obtained.

High order discretization schemes play more important role in fractional operators than classical ones. This is because usually for classical derivatives the stencil for high order discretization schemes is wider than low order ones; but for fractional operators the stencils for high order schemes and low order ones are the same. Then using high order schemes to solve fractional equations leads to almost the same computational cost as first order schemes but the accuracy is greatly improved. Using the fractional linear multistep methods, Lubich obtains the ν-th order (ν≤6) approximations of the α-th derivative (α>0) or integral (α<0) [Lubich, SIAM J. Math. Anal., 17, 704-719, 1986], because of the stability issue the obtained scheme can not be directly applied to the space fractional operator with α∈(1,2) for time dependent problem. By weighting and shifting Lubich's 2nd order discretization scheme, in [Chen & Deng, SINUM, arXiv:1304.7425] we derive a series of effective high order discretizations for space fractional derivative, called WSLD operators there. As the sequel of the previous work, we further provide new high order schemes for space fractional derivatives by weighting and shifting Lubich's 3rd and 4th order discretizations. In particular, we prove that the obtained 4th order approximations are effective for space fractional derivatives. And the corresponding schemes are used to solve the space fractional diffusion equation with variable coefficients.

A new weak boundary procedure for hyperbolic problems is presented. We consider high order finite difference operators of summation-by-parts form with weak boundary conditions and generalize that technique. The new boundary procedure is applied near boundaries in an extended domain where data is known. We show how to raise the order of accuracy of the scheme, how to modify the spectrum of the resulting operator and how to construct non-reflecting properties at the boundaries. The new boundary procedure is cheap, easy to implement and suitable for all numerical methods, not only finite difference methods, that employ weak boundary conditions. Numerical results that corroborate the analysis are presented.