Stable Numerical Methods Research Articles

Stable numerical methods for a stochastic nonlinear Schrödinger equation with linear multiplicative noise

<p style='text-indent:20px;'>This paper is concerned with fully discrete finite element approximations of a stochastic nonlinear Schrödinger (sNLS) equation with linear multiplicative noise of the Stratonovich type. The goal of studying the sNLS equation is to understand the role played by the noises for a possible delay or prevention of the collapsing and/or blow-up of the solution to the sNLS equation. In the paper we first carry out a detailed analysis of the properties of the solution which lays down a theoretical foundation and guidance for numerical analysis, we then present a family of three-parameters fully discrete finite element methods which differ mainly in their time discretizations and contains many well-known schemes (such as the explicit and implicit Euler schemes and the Crank-Nicolson scheme) with different combinations of time discetization strategies. The prototypical <inline-formula><tex-math id="M1">\begin{document}$ \theta $\end{document}</tex-math></inline-formula>-schemes are analyzed in detail and various stability properties are established for its numerical solution. An extensive numerical study and performance comparison are also presented for the proposed fully discrete finite element schemes.</p>

Second order approximation for a quasi-incompressible Navier-Stokes Cahn-Hilliard system of two-phase flows with variable density

Phase-field model has been applied extensively and successfully for simulating two-phase flows with variable density. Several stable numerical methods have been proposed for solving such a highly nonlinear and coupled system. The key issue is to design a method such that it can preserve the (exact or slightly modified) conservative/dissipative law of energy of the fluid system at the discrete level. However, most of the existing energy stable numerical methods are restricted to only the first order accuracy in time for two-phase flow with different density. The design of a temporally second order accurate numerical method for the two-phase flows with variable density still remains challenging. In this paper, we develop several second order, robust and accurate numerical schemes for a quasi-incompressible Navier-Stokes Cahn-Hilliard model of two-phase flows with variable density which is thermodynamically consistent and was originally developed in [1] for simulating two-phase flows in complex geometries. Especially, numerical schemes proposed in this paper can preserve the mass conservation or energy dissipative law. Several numerical examples are presented to validate the robustness and accuracy of our numerical schemes.

Stable and accurate numerical methods for generalized Kirchhoff–Love plates

Kirchhoff–Love plate theory is widely used in structural engineering. In this paper, efficient and accurate numerical algorithms are developed to solve a generalized Kirchhoff–Love plate model subject to three common physical boundary conditions: (i) clamped; (ii) simply supported; and (iii) free. The generalization stems from the inclusion of additional physics to the classical Kirchhoff–Love model that accounts for bending only. We solve the model equation by discretizing the spatial derivatives using second-order finite-difference schemes, and then advancing the semi-discrete problem in time with either an explicit predictor–corrector or an implicit Newmark-Beta time-stepping algorithm. Stability analysis is conducted for the schemes, and the results are used to determine stable time steps in practice. A series of carefully chosen test problems are solved to demonstrate the properties and applications of our numerical approaches. The numerical results confirm the stability and 2nd-order accuracy of the algorithms and are also comparable with experiments for similar thin plates. As an application, we illustrate a strategy to identify the natural frequencies of a plate using our numerical methods in conjunction with a fast Fourier transformation power spectrum analysis of the computed data. Then we take advantage of one of the computed natural frequencies to simulate the interesting physical phenomena known as resonance and beat for a generalized Kirchhoff–Love plate.

Dichotomy property of dispersion equation of guided waves propagating in anisotropic composite plates

Accurate prediction of dispersion curves for ultrasonic guided waves propagation in multi-layered composite laminates is crucial for the deployment of nondestructive testing procedures and structural health monitoring algorithms dedicated to aerospace composite materials. However, existing efforts mainly focused on finding ways to build complex-valued dispersion equations for guided waves (see transfer matrix method (TMM), global matrix method (GMM) and stiffness matrix method (SMM) for example) with little focus on developing efficient and stable numerical solving methods associated with the derived complex-valued equations. In this paper, the conditions under which complex-valued dispersion equations are either real- or purely imaginary-valued equations (termed as dichotomy property) are derived for both single- and multi-layered composite plates. With such a property, the complex-valued dispersion equations can be efficiently numerically solved within the real number field via the standard bisection method or the corrected phase change method. It is thus now possible to overcome numerical issues frequently reported in literature. Besides, a parallel computing technique is proposed in this paper to improve the computational efficiency of the traditional GMM. The proposed methodology provides a new standard framework to solve the dispersion equations which is stable, multipurpose, and numerically efficient.

A New Class of A Stable Summation by Parts Time Integration Schemes with Strong Initial Conditions

Since integration by parts is an important tool when deriving energy or entropy estimates for differential equations, one may conjecture that some form of summation by parts (SBP) property is involved in provably stable numerical methods. This article contributes to this topic by proposing a novel class of A stable SBP time integration methods which can also be reformulated as implicit Runge-Kutta methods. In contrast to existing SBP time integration methods using simultaneous approximation terms to impose the initial condition weakly, the new schemes use a projection method to impose the initial condition strongly without destroying the SBP property. The new class of methods includes the classical Lobatto IIIA collocation method, not previously formulated as an SBP scheme. Additionally, a related SBP scheme including the classical Lobatto IIIB collocation method is developed.

Truly multi-dimensional all-speed schemes for the Euler equations on Cartesian grids

Finite volume schemes often have difficulties to resolve the low Mach number (incompressible) limit of the Euler equations. Incompressibility is only non-trivial in multiple spatial dimensions. Low Mach fixes, however generally are applied to the one-dimensional method and the method is then used in a dimensionally split way. This often reduces its stability. Here, it is suggested to keep the one-dimensional method as it is, and only to extend it to multiple dimensions in a particular, all-speed way. This strategy is found to lead to much more stable numerical methods. Apart from the conceptually pleasing property of modifying the scheme only when it becomes necessary, the multi-dimensional all-speed extension also does not include any free parameters or arbitrary functions, which generally are difficult to choose, or might be problem dependent. The strategy is exemplified on a Lagrange Projection method and on a relaxation solver.

Vibrational Analysis of Cables Using the Non Standard Finite Difference Method

In continuous systems such as rods, cables, and shafts, the differential equation governing the behaviour of the structure can be obtained by applying Newton's second law. Only a limited group of these equations have an analytic solution; hence the provision of efficient and stable numerical methods is of particular importance. The method used in this study is a nonstandard finite difference (NSFD) scheme, which is an effective numerical method for solving the above-mentioned equations. One of the important features of this method is its relatively good high accuracy. Results obtained by nonstandard finite difference methods are more compatible with the exact behaviour of the problem than that of the standard finite difference method. In this research, the equations governing the vibrations of cables with three different boundary conditions are solved. Three solutions including the analytical solution, the standard finite difference method (SFD) and the NSFD method are compared. Results show that the error of the NSFD is significantly less than of the SFD.

Novel linear decoupled and unconditionally energy stable numerical methods for the modified phase field crystal model

In this paper, we propose a novel numerical approach to construct unconditionally energy stable schemes for the modified phase field crystal (MPFC) model. The new technique is based on the invariant energy quadratization (IEQ) method. The numerical schemes based on IEQ approach lead to time-dependent dense matrices, thus the fast Fourier transform (FFT) is difficult to be applied to solve the systems directly. By introducing a new auxiliary variable to replace the original one, we derive a novel equivalent MPFC system. A step-by-step solving approach, termed 3S-IEQ method is considered to solve this equivalent system. Some totally decoupled, linear and unconditional energy stable semi-implicit schemes based on 3S-IEQ method are very easy to construct. More importantly, the phase function ϕ and auxiliary variable η can be calculated step-by-step. Meanwhile, the proposed approach only needs to solve linear equation with constant coefficients which is easy to use FFT directly in calculation. Some numerical simulations are demonstrated to verify the accuracy and efficiency of our proposed schemes.

Stable numerical methods for determination of the molecular clusters force fields

Abstract The inverse problem of molecular force fields calculation is considered within the theory of regularization. In our strategy, we choose the stabilizing matrix F 0 F^{0} as a result of quantum mechanical calculations. The solution of the inverse problem is finding a matrix 𝐹 which is the nearest by the chosen Euclidean norm to the given ab initio F 0 F^{0} . The optimized solution is referred to as regularized quantum mechanical force field (RQMFF). Regularizing algorithms of molecular force fields calculation based on the joint treatment of experimental and ab initio quantum mechanical data have been applied to the calculations of molecular force fields (matrices of force constants) for small water clusters (H2O)𝑛 ( n = 2 , 3 n=2,3 ).

Convergence of second-order, entropy stable methods for multi-dimensional conservation laws

High-order accurate, entropy stable numerical methods for hyperbolic conservation laws have attracted much interest over the last decade, but only a few rigorous convergence results are available, particularly in multiple space dimensions. In this paper we show how the entropy stability of one such method, which is semi-discrete in time, yields a (weak) bound on oscillations. Under the assumption of L∞-boundedness of the approximations we use compensated compactness to prove convergence to a weak solution satisfying at least one entropy condition.

Chemo-mechanical coupling and material evolution in finitely deforming solids with advancing fronts of reactive fluids

A new stabilized method is presented for coupled chemo-mechanical problems involving chemically reacting fluids flowing through deformable elastic solids. A mixture theory model is employed wherein kinematics is represented via an independent set of balance laws for each of the interacting constituents. A significant feature of the mixture model is the interactive force field in the momentum balance equations that couples the constituents implicitly at the level of the governing system of equations. The constitutive relations for the constituents in the mixture model are based on maximization of the rate of entropy production. Since each constituent is not discretely modeled and the interactive effects are mathematically coupled at the local continuum level, the resulting system serves as a physics-based reduced-order model for the complex microstructure of the material system. When constitutive equations are substituted into the balance laws, they give rise to a system of coupled nonlinear PDEs. Evolving nonlinearity and coupled chemo-mechanical effects give rise to spatially localized phenomena, namely boundary layers, shear bands, and steep gradients that appear at the reaction fronts. For large reaction rates, the balance of mass of the fluid becomes a singularly perturbed equation (reaction-dominated), which may exhibit boundary and/or internal layers. Likewise, for large reaction rates and/or low diffusivity, the balance of linear momentum for the fluid constituent also becomes a singularly perturbed PDE. Presence of these features in the solution requires stable numerical methods, and we present a variational multiscale (VMS)-based stabilized finite element method for the initial-boundary value problem. Mathematical attributes of the method are investigated via a range of numerical test cases that involve diffusion of chemically reacting fluids through nonlinear elastic solids. Enhanced stabilization features and higher spatial accuracy of the models and the methods are highlighted.

Nonequilibrium Kinetics of the Initial Stage of a Phase Transition

The Kolmogorov–Feller and Einstein–Smoluchowski kinetic partial derivative equations with nonlinear coefficients are solved by new stable numerical methods. The theory of stochastic dynamic variables relates the solution of the Ito equation in the Stratanovich sense for trajectories of Wiener random processes to the transition probability density of these processes or the distribution functions of kinetic equations. The classical nucleation (formation of first-order phase transition nuclei) theory describes the nonequilibrium stage of the condensation process by a diffusion random process the phase space nuclei size, when fluctuations influence the nucleus clusterization. The model of formation of vacancy-gas defects (VGD) (pores and blisters) in a crystal lattice, which appear as a result of its irradiation by inert Xe++ gas is supplemented by the consideration of the Brownian motion of nonpoint lattice defects, which occurs under action of the superposition of pair long-range model potentials of an indirect elastic interaction of VGDs to each other and to the layer boundaries. The spacial-temporal porosity structures in a sample form for times of 10‒100 μs as a result of the Brownian motion of vacancy-gas defects, which is simulated using stable algorithms. The nonequilibrium kinetic distribution functions of VGDs on sizes and coordinates in layers of irradiated materials have been found as a result of calculations of 106 trajectories; these functions characterize the fluctuation instability of the initial stage of the phase transition and they are used to estimate local stresses and the porosity in a model volume.

Stochastic Model of Pore Nucleation when the Sample is Irradiates with Inert Gas Ions

Abstract—Non-point vacancy-gas defects (VGD) (pores, blisters) in the crystal lattice arise as a result of its irradiation with an inert gas Xe++, the process is considered as the formation of first-order phase transition nuclei by means of computational mathematics, kinetic theory, and the theory of stochastic dynamic variables and Wiener random processes. The characteristics of disordered porosity in samples consisting of layers (“dielectric/metal”) at various layer thicknesses are calculated. The formation of spherical shape defects in materials at times of the order of 10–4 s, when the porosity of the material nucleates, is preceded or accompanied by the appearance of microcracks. Irradiation of surfaces with an inert gas with ion energies of 5–10 keV refers to low-temperature blistering; in the model, the constancy of the influence of radiation fluxes stimulating a phase transition creates the condition of an “open” physical system in which it is possible to combine lattice defects into structures that cause local stresses, and the constancy of the temperature of the sample, the pressure of the implantable “monomers” of the gas, and supersaturation (similar to the parameters of vapor condensation) are characteristic of the fluctuation unstable stage of the phase transition move. Partial kinetic equations with nonlinear coefficients (Kolmogorov–Feller and Einstein–Smoluchowski) are solved by the method of stochastic molecular dynamics, which establishes a connection between the solution of Ito stochastic equations in the sense of Stratonovich and kinetic equations; stable numerical methods are used. The conditions for the formation of VGD are determined taking into account the influence of the elastic properties of the lattice, and the porosity and local stresses in thin layers of irradiated materials are calculated.

Неравновесная кинетика начальной стадии фазового перехода

Kinetic partial differential equations of Kolmogorov-Feller and Einstein-Smoluchowski equation with nonlinear coefficients are solved by a new, stable numerical methods. The theory of stochastic dynamic variables establishes the connection of the solution of Ito stochastic equations in the sense of Stratonovich for the trajectories of Wiener random processes with the transition probability density of these processes, or distribution functions of kinetic equations. The classical theory of nucleation (formation of nuclei of the first order phase transition) describes a non-equilibrium stage of the condensation process by a diffusion random process in the space of the size of the nuclei of the phase transition, when fluctuations affect the clustering of the nuclei. The model of formation of vacancy-gas defects (pores, blisters) in the crystal lattice, arising as a result of its irradiation by inert gas ions xenon, is supplemented by the consideration of Brownian motion of non-point lattice defects, occurring under the action of superposition of paired long-range potentials of indirect elastic interaction of pores between themselves and with the boundaries of the layers. Pores coordinates are changing at times of the order of 10 − 100 ms, sustainable algorithms for calculating which provide a self-consistent defnition spatial - temporal structures of porosity in the sample. According to calculations of 106 trajectories, non-equilibrium kinetic functions were found. Pores distribution in size and coordinates in the layers of the irradiated materials, they characterize the fluctuation instability the initial stage of the phase transition, they are estimated local stresses and porosity in the model volume.

Stable and efficient numerical schemes for two-dimensional Maxwell equations in lossy medium

In this paper, some stable and efficient numerical methods are developed for two-dimensional Maxwell equations. To avoid solving large scale algebraic equations, we use local one-dimensional (LOD) technique and split the original equations into several LOD equations. Then, the resulted LOD equations are discretized by Wendroff scheme which is usually applied to simulate hyperbolic-type equations. The Lie-Trotter composition and Strang composition are employed. To improve the computational efficiency and the convergent accurate, the space direction is approximated by high order compact method. Some benchmark symbols, such as stability, conservation laws of the schemes are analyzed. Some numerical examples confirm the theoretical analysis. Numerical results show that the new schemes can not only accurately simulate the electromagnetic waves, including profiles, amplitude, but also capture the energy structures exactly.

Numerical solution of the conditionally averaged equations for representing net mass flux due to convection

The representation of subgrid‐scale convection is a weak aspect of weather and climate prediction models and the assumption that no net mass is transported by convection in parametrizations is increasingly unrealistic as models enter the grey zone, partially resolving convection. The solution of conditionally averaged equations of motion (multifluid equations) is proposed in order to avoid this assumption. Separate continuity, temperature, and momentum equations are solved for inside and outside convective plumes, which interact via mass‐transfer terms, drag, and by a common pressure. This is not a convection scheme that can be used with an existing dynamical core—this requires a whole new model.This article presents stable numerical methods for solving the multifluid equations, including large transfer terms between the environment and plume fluids. Without transfer terms, the two fluids are not sufficiently coupled and solutions diverge. Two transfer terms are presented, which couple the fluids together in order to stabilize the model: diffusion of mass between the fluids (similar to turbulent entrainment) and drag between the fluids. Transfer terms are also proposed to move buoyant air into the plume fluid and vice versa as would be needed to represent initialization and termination of subgrid‐scale convection. The transfer terms are limited (clipped in size) and solved implicitly in order to achieve bounded, stable solutions.Results are presented for a well‐resolved warm bubble with rising air being transferred to the plume fluid. For stability, equations are formulated in advective rather than flux form and solved using bounded finite‐volume methods. Discretization choices are made to preserve boundedness and conservation of momentum and energy when mass is transferred between fluids.The formulation of transfer terms in order to represent subgrid convection is the subject of future work.

Light transport in homogeneous tissue with m-dependent anisotropic scattering I: Case's singular eigenfunctions solution and Chandrasekhar polynomials.

This paper is the first of two deriving the analytical solutions for light transport in infinite homogeneous tissue with an azimuth-dependent (m-dependent) anisotropic scattering kernel by two approaches, Case’s singular eigenfuncions expansion and Fourier transform, as well as proving the consistence of the two solutions. In this paper, Case’s method was applied and extended to the general m-dependent anisotropic scattering case. The explicit Green’s function of radiance distributions, which was regarded as the comparative standard for the equivalent solution via Fourier transform and inversion in our second accompanying paper, was expanded into a complete set of the discrete and continuous eigenfunctions. Considering that the two kinds of m-dependent Chandrasekhar orthogonal polynomials that play vital roles in these analytical solutions are very sensitive to the typical optical parameters of biological tissue as well as the degrees or orders, four numerical evaluation methods were benchmarked to find the stable, reliable and feasible numerical evaluation methods in high degrees and high orders.

On segmentation model for vector valued images and fast iterative solvers

In this paper, we propose a new convex variational model for segmentation of vector valued images. The data term of the proposed model is based on the coefficient of variation, which works well in vector valued images having intensity inhomogeneity. Due to convexity of the model, it is independent of the placement of initial contour. Better performance of the proposed model can be seen from experimental results qualitatively and quantitatively. Images in practice are of large sizes, which makes numerical methods more important. In this paper, we also develop fast and stable numerical methods for solution of partial differential equation arisen from the minimization of the proposed model. We have developed a novel multigrid method based on a locally supported smoother. The proposed method is compared with the existing methods in terms of iterations and CPU time for vector valued images having large sizes.

New developments of the Extended Quadrature Method of Moments to solve Population Balance Equations

Population Balance Models have a wide range of applications in many industrial fields as they allow accounting for heterogeneity among properties which are crucial for some system modelling. They actually describe the evolution of a Number Density Function (NDF) using a Population Balance Equation (PBE). For instance, they are applied to gas–liquid columns or stirred reactors, aerosol technology, crystallisation processes, fine particles or biological systems. There is a significant interest for fast, stable and accurate numerical methods in order to solve for PBEs, a class of such methods actually does not solve directly the NDF but resolves their moments. These methods of moments, and in particular quadrature-based methods of moments, have been successfully applied to a variety of systems. Point-wise values of the NDF are sometimes required but are not directly accessible from the moments. To address these issues, the Extended Quadrature Method of Moments (EQMOM) has been developed in the past few years and approximates the NDF, from its moments, as a convex mixture of Kernel Density Functions (KDFs) of the same parametric family. In the present work EQMOM is further developed on two aspects. The main one is a significant improvement of the core iterative procedure of that method, the corresponding reduction of its computational cost is estimated to range from 60% up to 95%. The second aspect is an extension of EQMOM to two new KDFs used for the approximation, the Weibull and the Laplace kernels. All MATLAB source codes used for this article are provided with this article.

Efficient energy stable schemes for isotropic and strongly anisotropic Cahn–Hilliard systems with the Willmore regularization

We develop efficient energy stable numerical methods for solving isotropic and strongly anisotropic Cahn–Hilliard systems with the Willmore regularization. The scheme, which involves adaptive mesh refinement and a nonlinear multigrid finite difference method, is constructed based on a convex splitting approach. We prove that, for the isotropic Cahn–Hilliard system with the Willmore regularization, the total free energy of the system is non-increasing for any time step and mesh sizes. A straightforward modification of the scheme is then used to solve the regularized strongly anisotropic Cahn–Hilliard system, and it is numerically verified that the discrete energy of the anisotropic system is also non-increasing, and can be efficiently solved by using the modified stable method. We present numerical results in both two and three dimensions that are in good agreement with those in earlier work on the topics. Numerical simulations are presented to demonstrate the accuracy and efficiency of the proposed methods.