Search

In this work we introduce a novel semi-implicit structure-preserving finite-volume/finite-difference scheme for the viscous and resistive equations of magnetohydrodynamics (MHD) based on an appropriate 3-split of the governing PDE system, which is decomposed into a first convective subsystem, a second subsystem involving the coupling of the velocity field with the magnetic field and a third subsystem involving the pressure-velocity coupling. The nonlinear convective terms are discretized explicitly, while the remaining two subsystems accounting for the Alfven waves and the magneto-acoustic waves are treated implicitly. The final algorithm is at least formally constrained only by a mild CFL stability condition depending on the velocity field of the pure hydrodynamic convection. To preserve the divergence-free constraint of the magnetic field exactly at the discrete level, a proper set of overlapping dual meshes is employed. The resulting linear algebraic systems are shown to be symmetric and therefore can be solved by means of an efficient standard matrix-free conjugate gradient algorithm. One of the peculiarities of the presented algorithm is that the magnetic field is defined on the edges of the main grid, while the electric field is on the faces. The final scheme can be regarded as a novel shock-capturing, conservative and structure preserving semi-implicit scheme for the nonlinear viscous and resistive MHD equations. Several numerical tests are presented to show the main features of our novel solver: linear-stability in the sense of Lyapunov is verified at a prescribed constant equilibrium solution; a 2nd-order of convergence is numerically estimated; shock-capturing capabilities are proven against a standard set of stringent MHD shock-problems; accuracy and robustness are verified against a nontrivial set of 2- and 3-dimensional MHD problems., 44 pages, 22 figures, 2 tables

We present a new method for estimating the unobservable drift of the risk-neutral spot freight rate process from Forward Freight Agreements (FFA) prices in the absence of a closed-form solution and demonstrate robustness via numerical simulations. Moreover, we conduct empirical experiments involving estimation of standard parametric models and a nonparametric model using Baltic Exchange data. We find that our nonparametric approach yields the lowest FFA pricing errors across maturities. Finally, we estimate the market price of risk, analyze its behavior in-sample and out-of-sample and observe that, when estimated using our nonparametric approach, it evolves consistently with the indices under study.

PURPOSE To demonstrate and validate electric field (E-field) calculation and peripheral nerve stimulation (PNS) prediction methods that are accurate, computationally efficient, and that could be used to inform regulatory standards. METHODS We describe a simplified method for calculating the spatial distribution of induced E-field over the volume of a body model given a gradient coil vector potential field. The method is easily programmed without finite element or finite difference software, allowing for straightforward and computationally efficient E-field evaluation. Using these E-field calculations and a range of body models, population-weighted PNS thresholds are determined using established methods and compared against published experimental PNS data for two head gradient coils and one body gradient coil. RESULTS A head-gradient-appropriate chronaxie value of 669 µs was determined by meta-analysis. Prediction errors between our calculated PNS parameters and the corresponding experimentally measured values were ~5% for the body gradient and ~20% for the symmetric head gradient. Our calculated PNS parameters matched experimental measurements to within experimental uncertainty for 73% of ∆Gmin estimates and 80% of SRmin estimates. Computation time is seconds for initial E-field maps and milliseconds for E-field updates for different gradient designs, allowing for highly efficient iterative optimization of gradient designs and enabling new dimensions in PNS-optimal gradient design. CONCLUSIONS We have developed accurate and computationally efficient methods for prospectively determining PNS limits, with specific application to head gradient coils.

Compartmental model is the classic and widely used approach in pharmacokinetics. Recently, fractional calculus is introduced to describe the time course of drugs in human body which follow the anomalous diffusion mechanism. To consider the different fractional order transmission process, the two compartmental fractional model has been proposed and studied. And, it will be of great significance to find out a simple and efficient numerical method to solve the model and estimate model parameters. This work investigates two numerical methods of the two compartmental fractional model using finite difference schemes, which are based on the shifted Grunwald-Letnikov approximate formula and L1 formula, respectively. The hybrid Nelder-Mead simplex search and particle swarm optimization, denoted as NMSS-PSO, is provided to estimate the fractional model parameters. Comparison between the numerical solution and the solution by the numerical inverse Laplace transform method establishes the validity of the developed numerical algorithms. Then, the influence of the order of fractional derivative on the drug amount in human body is also discussed. Finally, the two compartmental fractional model is applied to fit the amiodarone pharmacokinetics data based on the NMSS-PSO algorithms. This work will be of importance for the development of fractional pharmacokinetics.

The energy-conserving sampling and weighting (ECSW) method is a hyperreduction method originally developed for accelerating the performance of Galerkin projection-based reduced-order models (PROMs) associated with large-scale finite element models, when the underlying projected operators need to be frequently recomputed as in parametric and/or nonlinear problems. In this paper, this hyperreduction method is extended to Petrov-Galerkin PROMs where the underlying high-dimensional models can be associated with arbitrary finite element, finite volume, and finite difference semi-discretization methods. Its scope is also extended to cover local PROMs based on piecewise-affine approximation subspaces, such as those designed for mitigating the Kolmogorov $n$-width barrier issue associated with convection-dominated flow problems. The resulting ECSW method is shown in this paper to be robust and accurate. In particular, its offline phase is shown to be fast and parallelizable, and the potential of its online phase for large-scale applications of industrial relevance is demonstrated for turbulent flow problems with $O(10^7)$ and $O(10^8)$ degrees of freedom. For such problems, the online part of the ECSW method proposed in this paper for Petrov-Galerkin PROMs is shown to enable wall-clock time and CPU time speedup factors of several orders of magnitude while delivering exceptional accuracy.

We present a class of new explicit and stable numerical algorithms to solve the spatially discretized linear heat or diffusion equation. After discretizing the space and the time variables like conventional finite difference methods, we do not approximate the time derivatives by finite differences, but use constant neighbor and linear neighbour approximations to decouple the ordinary differential equations and solve them analytically. During this process, the timestep-size appears not in polynomial, but in exponential form with negative exponents, which guarantees stability. We compare the performance of the new methods with analytical and numerical solutions. According to our results, the methods are first and second order in time and can be much faster than the commonly used explicit or implicit methods, especially in the case of extremely large stiff systems., 21 pages

We consider the test suite for the Shallow Water (SW) equations on the sphere suggested in [27, 28]. This series of tests consists of zonally propagating wave solutions of the linearized Shallow Water (LSW) equations on the full sphere. Two series of solutions are considered. The first series [27] is referred to as "barotropic". It consists of an extension of the Rossby-Haurwitz test case in [33]. The second series [28] referred to as (Matsuno) "baroclinic", consists of a generalisation of the solution to LSW in an equatorial chanel introduced by Matsuno [17]. The Hermitian Compact Cubed Sphere (HCCS) model which is used in this paper is a Shallow Water solver on the sphere that was introduced in [4]. The spatial approximation is a center finite difference scheme based on high order differencing along great circles. The time stepping is performed by the explicit RK4 scheme or by an exponential scheme. For both test cases, barotropic and baroclinic, the results show a very good agreement of the numerical solution with the analytic one, even for long time simulations.

Reverse time migration can accurately image underground earth structures. However, for shallow seismic exploration, the seismic wave velocity is often lower than the velocity in the middle‐deep layers, which causes numerical dispersion for finite‐difference schemes and leads to poor seismic imaging quality. Suppressing numerical dispersion by grid encryption or increasing the finite‐difference order seriously reduces computational efficiency, which is not the optimal solution. To improve the imaging quality without sacrificing computational efficiency, a regularization factor is added to the acoustic wave equation to correct the phase velocity of high wave‐number components. An appropriate regularization factor can eliminate numerical dispersion that results from large grid interval schemes and can reduce the size of computational grids and improve computational efficiency. For simulations in the frequency domain, reducing the grid size also means reducing computer memory requirements. Numerical experiments indicate that the regularization factor should match the degree of numerical dispersion. Larger regularization factors can suppress serious numerical dispersion. However, excessively large regularization factors may destroy the effective wave field. Different numerical models verify the effectiveness of the improved acoustic equation for suppressing numerical dispersion and maintaining amplitudes and provide a novel means to improve the shallow seismic image quality.

The non-linear, non-isothermal, magnetohydrodynamic (MHD) laminar convection flows of\ud Buongiorno’s nanofluid past a vertical surface with Darcy-Forchheimer model is mathematically\ud investigated in the present article. Keller’s Box implicit finite difference technique is utilized to\ud solve the dimensionless conservation equations. Graphical and tabulated results are analyzed to\ud study the behavior of primary and secondary velocity, temperature, nanoparticle volume fraction,\ud shear stress rate, heat transfer rate and mass transfer rate for various emerging thermos-physical\ud parameters. The Hall current and ion slip current effects are also considered. Validations of earlier\ud solutions in the literature is also included. The study finds applications in nanomaterial fabrication\ud processing, biomedical, polymer processing, chemical engineering, crude oil purifying, etc.

SIR model is widely used for modeling the infectious diseases. This is a system of ordinary differential equations (ODEs). The numbers of susceptible, infectious, or immunized individuals are the compartments in these equations and change in time. Two parameters are the factor of differentiating these models. Here, we are not interested in solving the ODEs describing a certain SIR model. Given the observed data, we try to estimate the parameters that determine the model. For this, we propose a least squares approach using second-order centered differences for replacing the derivatives appeared in the ODEs. Then we arrive at a simple linear system that can be solved explicitly and furnish the approximations of the parameters. Numerical results over various artificial data verify the simplicity and accuracy of the new method. © 2020 John Wiley & Sons, Ltd.

This paper is a rebuttal to the claim found in the literature that the MUSCL scheme cannot be third-order accurate for nonlinear conservation laws. We provide a rigorous proof for third-order accuracy of the MUSCL scheme based on a careful and detailed truncation error analysis. Throughout the analysis, the distinction between the cell average and the point value will be strictly made for the numerical solution as well as for the target operator. It is shown that the average of the solutions reconstructed at a face by Van Leer's kappa-scheme recovers a cubic solution exactly with kappa = 1/3, the same is true for the average of the nonlinear fluxes evaluated by the reconstructed solutions, and a dissipation term is already sufficiently small with a third-order truncation error. Finally, noting that the target spatial operator is a cell-averaged flux derivative, we prove that the leading truncation error of the MUSCL finite-volume scheme is third-order with kappa = 1/3. The importance of the diffusion scheme is also discussed: third-order accuracy will be lost when the third-order MUSLC scheme is used with a wrong fourth-order diffusion scheme for convection-diffusion problems. Third-order accuracy is verified by thorough numerical experiments for both steady and unsteady problems. This paper is intended to serve as a reference to clarify confusions about third-order accuracy of the MUSCL scheme, as a guide to correctly analyze and verify the MUSCL scheme for nonlinear equations, and eventually as the basis for clarifying third-order unstructured-grid schemes in a subsequent paper.

We present the fourth-order compact finite difference (4cFD) discretizations for the long time dynamics of the nonlinear Klein-Gordon equation (NKGE), while the nonlinearity strength is characterized by $\varepsilon^p$ with a constant $p \in \mathbb{N}^+$ and a dimensionless parameter $\varepsilon \in (0, 1]$. Based on analytical results of the life-span of the solution, rigorous error bounds of the 4cFD methods are carried out up to the time at $O(\varepsilon^{-p})$. We pay particular attention to how error bounds depend explicitly on the mesh size $h$ and time step $\tau$ as well as the small parameter $\varepsilon \in (0, 1]$, which indicate that, in order to obtain `correct' numerical solutions up to the time at $O(\varepsilon^{-p})$, the $\varepsilon$-scalability (or meshing strategy requirement) of the 4cFD methods should be taken as: $h = O(\varepsilon^{p/4})$ and $\tau = O(\varepsilon^{p/2})$. It has better spatial resolution capacity than the classical second order central difference methods. By a rescaling in time, it is equivalent to an oscillatory NKGE whose solution propagates waves with wavelength at $O(1)$ in space and $O(\varepsilon^p)$ in time. It is straightforward to get the error bounds of the oscillatory NKGE in the fixed time. Finally, numerical results are provided to confirm our theoretical analysis.

Solute transport is usually modeled by the advection-dispersion-reaction equation. In the standard approach, mechanical dispersion is a tensor with principal directions parallel and perpendicular to the flow vector. Since realistic scenarios include nonuniform and unsteady flow fields, the governing equation has full tensor mechanical dispersion. When conventional grid-based numerical methods are used, approximation of the cross terms arising from the off-diagonal terms cause nonphysical solution with oscillations. As an example, for the common scenario of contaminant input into a domain with zero initial concentration, the cross-dispersion terms can result in negative concentrations that can wreak havoc in reactive transport applications. To address this issue, we use the well-known flux-corrected-transport (FCT) technique for a standard finite volume method. Although FCT has most often been used to eliminate oscillations resulting from discretization of the advection term for explicit time stepping, we show that it can be adapted for full-tensor dispersion and implicit time stepping. Unlike other approaches based on new discretization techniques (e.g., mimetic finite difference, nonlinear finite volume), FCT has the advantage of being flexible and widely applicable. Implementation of FCT requires solving an additional system of equations at each time step, using a modified "low order" matrix and a modified right-hand-side vector. To demonstrate the flexibility of FCT, we have modified the well-known and widely used groundwater solute transport simulator, MT3DMS. We apply the new simulator, MT3DMS-FCT, to several benchmark problems that suffer from negative concentrations when using MT3DMS. The new results are mass conservative and strictly nonnegative.

In this paper, a complete simulation of a trombone using finite difference time-domain (FDTD) methods is proposed. In particular, we propose the use of a novel method to dynamically vary the number of grid points associated to the FDTD method, to simulate the fact that the physical dimension of the trombone’s resonator dynamically varies over time. We describe the different elements of the model and present the results of a real-time simulation.

Explicit numerical finite difference schemes for partial differential equations are well known to be easy to implement but they are particularly problematic for solving equations whose solutions admit shocks, blowups and discontinuities. Here we present an explicit numerical scheme for solving non-linear advection-diffusion equations admitting shock solutions that is both easy to implement and stable. The numerical scheme is obtained by considering the continuum limit of a discrete time and space stochastic process for non-linear advection-diffusion. The stochastic process is well posed and this guarantees the stability of the scheme. Several examples are provided to highlight the importance of the formulation of the stochastic process in obtaining a stable and accurate numerical scheme.

Two real‐domain numerical approximation methods for accurate computation of finite strain elasticity moduli are developed and their accuracy and computational efficiency investigated, with reference to hyperelastic constitutive models with known analytical solutions. The methods are higher‐order and higher floating‐point precision numerical approximation, the latter being novel in this context. A general formula for higher‐order approximation finite difference schemes is derived and a new procedure is proposed to implement increased floating‐point precision. The accuracy of the approximated elasticity moduli is investigated numerically using higher‐order approximations in standard double precision and increased quadruple precision. It is found that as the order of the approximation increases, the elasticity moduli tend towards the analytical solution. Using higher floating‐point precision, the approximated elasticity moduli for all orders of approximation are found to be more accurate than the standard double precision evaluation of the analytical moduli. Application of the techniques to a finite element problem shows that the numerically approximated methods obtain convergence equivalent to the analytical method but require greater computational effort. It is concluded that numerical approximation of elasticity moduli is a powerful and effective means of implementing advanced constitutive models in the finite element method without prior derivation of difficult analytical solutions.

Mesh-free discretization, flexibly distributing nodes without computationally expen-sive meshing process, is able to deal with staircase problem, oversampling and under-sampling problems and saves plenty of nodes through distributing nodes suitably with respect to irregular boundaries and model parameters. However, the time-domain mesh-free discretization usually exhibits poorer stability than that in regular grid discretization. In order to reach unconditional stability and easy implementation in parallel computing, we develop the frequency-domain finite-difference method in a mesh-free discretization, incorporated with two perfectly matched layer boundary conditions. Furthermore, to maintain the flexibility of mesh-free discretization, the nodes are still irregularly distributed in the absorbing zone, which complicates the situation of artificial boundary reflections. In this paper, we implement frequency-domain acoustic wave modelling in a mesh-free system. First, we present the perfectly matched layer boundary condition to suppress spurious reflections. Moreover, we develop the complex frequency shifted–perfectly matched layer boundary condition to improve the attenuation of grazing waves. In addition, we employ the radial-basis-function-generated finite difference method in the mesh-free discretization to calculate spatial derivatives. The numerical experiment on a rectangle homogeneous model shows the effectiveness of the perfectly matched layer boundary condition and the complex frequency shifted–perfectly matched layer boundary condition, and the latter one is better than the former one when absorbing large angle incident waves. The experiment on the Marmousi model suggests that the complex frequency shifted–perfectly matched layer boundary condition works well for complicated models.

Implicit-solvent models are widely used to study the electrostatics in dissolved biomolecules, which are parameterized using force fields. Standard force fields treat the charge distribution with point charges; however, other force fields have emerged which offer a more realistic description by considering polarizability. In this work, we present the implementation of the polarizable and multipolar force field atomic multipole optimized energetics for biomolecular applications (AMOEBA), in the boundary integral Poisson-Boltzmann solver PyGBe. Previous work from other researchers coupled AMOEBA with the finite-difference solver APBS, and found difficulties to effectively transfer the multipolar charge description to the mesh. A boundary integral formulation treats the charge distribution analytically, overlooking such limitations. This becomes particularly important in simulations that need high accuracy, for example, when the quantity of interest is the difference between solvation energies obtained from separate calculations, like happens for binding energy. We present verification and validation results of our software, compare it with the implementation on APBS, and assess the efficiency of AMOEBA and classical point-charge force fields in a Poisson-Boltzmann solver. We found that a boundary integral approach performs similarly to a volumetric method on CPU. Also, we present a GPU implementation of our solver. Moreover, with a boundary element method, the mesh density to correctly resolve the electrostatic potential is the same for standard point-charge and multipolar force fields. Finally, we saw that for binding energy calculations, a boundary integral approach presents more consistent results than a finite difference approximation for multipolar force fields. © 2019 Wiley Periodicals, Inc.

Generalized solitary waves with exponentially small non-decaying far field oscillations have been studied in a range of singularly-perturbed differential equations, including higher-order Korteweg-de Vries (KdV) equations. Many of these studies used exponential asymptotics to compute the behaviour of the oscillations, revealing that they appear in the solution as special curves known as Stokes lines are crossed. Recent studies have identified similar behaviour in solutions to difference equations. Motivated by these studies, the seventh-order KdV and a hierarchy of higher-order KdV equations are investigated, identifying conditions which produce generalized solitary wave solutions. These results form a foundation for the study of infinite-order differential equations, which are used as a model for studying lattice equations. Finally, a lattice KdV equation is generated using finite-difference discretization, in which a lattice generalized solitary wave solution is found.

We derived, implemented, and thoroughly tested the complete analytic expression for atomic forces, consisting of the Hellmann-Feynman term and the Pulay correction, for the Fermi-Lowdin orbital self-interaction correction (FLO-SIC) method. Analytic forces are shown to be numerically accurate through an extensive comparison to forces obtained from finite differences. Using the analytic forces, equilibrium structures for a small set of molecules were obtained. This work opens the possibility of routine self-interaction free geometrical relaxations of molecules using the FLO-SIC method. © 2018 Wiley Periodicals, Inc.

We present a numerical method for simulating 2D flow through a channel with deformable walls. The fluid is assumed to be incompressible and viscous. We consider the highly viscous regime, where fluid dynamics are described by the Stokes equations, and the less viscous regime described by the Navier–Stokes equations. The model is formulated as an immersed boundary problem, with the channel defined by compliant walls that are immersed in a larger computational fluid domain. The channel traverses through the computational domain, and the walls do not form a closed region. When the walls deviate from their equilibrium position, they exert singular forces on the underlying fluid. We compute the numerical solution to the model equations using the immersed interface method, which preserves sharp jumps in the solution and its derivatives. The immersed interface method typically requires a closed immersed interface, a condition that is not met by the present configuration. Thus, a contribution of the present work is the extension of the immersed interface method to immersed boundary problems with open interfaces. Numerical results indicate that this new method converges with second‐order accuracy in both space and time, and can sharply capture discontinuities in the fluid solution.

Some mathematical models of applied problems lead to the need of solving boundary value problems with a fractional power of an elliptic operator. In a number of works, approximations of such a nonlocal operator are constructed on the basis of an integral representation with a singular integrand. In the present paper, new and more convenient integral representations are proposed for operators with fractional powers. Approximations are based on the classical quadrature formulas. The results of numerical experiments on the accuracy of quadrature formulas are presented. The proposed approximations are used for numerical solving a model two-dimensional boundary value problem for fractional diffusion.