Midpoint Scheme Research Articles

On the properties of numerical solutions of dynamical systems obtained using the midpoint method

The article considers the midpoint scheme as a finite-difference scheme for a dynamical system of the form ̇ = (). This scheme is remarkable because according to Cooper’s theorem, it preserves all quadratic integrals of motion, moreover, it is the simplest scheme among symplectic Runge-Kutta schemes possessing this property. The properties of approximate solutions were studied in the framework of numerical experiments with linear and nonlinear oscillators, as well as with a system of several coupled oscillators. It is shown that in addition to the conservation of all integrals of motion, approximate solutions inherit the periodicity of motion. At the same time, attention is paid to the discussion of introducing the concept of periodicity of an approximate solution found by the difference scheme. In the case of a nonlinear oscillator, each step requires solving a system of nonlinear algebraic equations. The issues of organizing computations using such schemes are discussed. Comparison with other schemes, including those symmetric with respect to permutation of and .̂

On the properties of numerical solutions of dynamical systems obtained using the midpoint method

The article considers the midpoint scheme as a finite-difference scheme for a dynamical system of the form

Linear second-order IMEX-type integrator for the (eddy current) Landau–Lifshitz–Gilbert equation

Abstract Combining ideas from Alouges et al. (2014, A convergent and precise finite element scheme for Landau–Lifschitz–Gilbert equation. Numer. Math., 128, 407–430) and Praetorius et al. (2018, Convergence of an implicit-explicit midpoint scheme for computational micromagnetics. Comput. Math. Appl., 75, 1719–1738) we propose a numerical algorithm for the integration of the nonlinear and time-dependent Landau–Lifshitz–Gilbert (LLG) equation, which is unconditionally convergent, formally (almost) second-order in time, and requires the solution of only one linear system per time step. Only the exchange contribution is integrated implicitly in time, while the lower-order contributions like the computationally expensive stray field are treated explicitly in time. Then we extend the scheme to the coupled system of the LLG equation with the eddy current approximation of Maxwell equations. Unlike existing schemes for this system, the new integrator is unconditionally convergent, (almost) second-order in time, and requires the solution of only two linear systems per time step.

Strange Dynamic Behavior Induced by Time Delay for a Rijke Tube System with Periodic Excitation

Thermoacoustic instability problems widely exist in many real-world applications such as gas turbines and rocket motors. A Rijke tube is a typical thermoacoustic system, and it is difficult to analyze such a system due to the nonlinearity and time delay. In this paper, a set of nonlinear ordinary differential equations with time delay which represent a Rijke tube system will be studied. The state space of such a tube system consists of velocity and pressure, and the periodic motion can be discretized based on an implicit midpoint scheme. Through Newton–Raphson method, the node points on the periodic motion will be solved, and the analytical solution of such a periodic motion for Rijke tube system can be recovered using a set of Fourier representations. According to the theory of discrete maps, the stability of the periodic motion will be obtained. With such a proposed technique, the stable and unstable branches of period-1 motions with one mode and multiple modes will be presented. The Neimark–Sacker bifurcation will be observed, and the periodic and quasi-periodic motions related to Neimark–Sacker bifurcation will be discussed. The strange dynamic behavior for stable period-1 motion at the neighborhood of Neimark–Sacker bifurcation will be observed, and the motion can only be attracted to the stable period-1 limit cycle when the initial conditions are given exactly from the analytic solution of such a stable period-1 motion. Otherwise, it will be attracted to another coexisting quasi-periodic motion.

Computational micromagnetics with Commics

We present our open-source Python module Commics for the study of the magnetization dynamics in ferromagnetic materials via micromagnetic simulations. It implements state-of-the-art unconditionally convergent finite element methods for the numerical integration of the Landau–Lifshitz–Gilbert equation. The implementation is based on the multiphysics finite element software Netgen/NGSolve. The simulation scripts are written in Python, which leads to very readable code and direct access to extensive post-processing. Together with documentation and example scripts, the code is freely available on GitLab. Program summaryProgram title: CommicsProgram Files doi:http://dx.doi.org/10.17632/29wv9h78h7.1Licensing provisions: GPLv3Programming language: Python3Nature of problem: Numerical integration of the Landau–Lifshitz–Gilbert equation in three space dimensionsSolution method: Tangent plane scheme [1]: original first-order version, projection-free version, second-order version, efficient second-order IMEX version; Midpoint scheme [2]: original version, IMEX version;Magnetostatic Maxwell equations are treated by the hybrid FEM–BEM method [3]Additional comments including restrictions and unusual features: An installation of the finite element software Netgen/NGSolve and an installation of the boundary element library BEM++ are required.

The numerical influence of additional parameters of inertia representations for quaternion-based rigid body dynamics

Different inertia representations can lead to different formulations of the differential-algebraic equations for the quaternion-based rigid body dynamics. In this paper, the inertia representations are classified into $\alpha $-type and $\gamma $-type, according to the additional parameters in the kinetic energy. These two types of representations and the corresponding parameters $\alpha $ and $\gamma $ are theoretically equivalent if the constraint $\boldsymbol{q}^{T}\boldsymbol{q} = 1$ is satisfied exactly. Nevertheless, the error estimation demonstrates that they can present entirely different numerical features in simulation and suggests that the parameter $\gamma $ can be used to optimize the numerical performance of the integrations in simulation. To further verify the numerical difference between the inertia representations of $\alpha $-type and $\gamma $-type, the corresponding modified Hamilton’s equations are discretized by the IMS (implicit midpoint scheme), EMS (energy–momentum preserving scheme) and Gauss–Lobatto SPARK methods. Numerical performance for the examples of the spinning symmetrical top is shown to result from the comprehensive effect of the discretization schemes including the distribution of discretized points and the convergence order, the inertia representations and their combinations. Numerical results further suggest that the integrations of $\gamma $-type are superior to those of $\alpha $-type and the optimized values of $\gamma $ can be used to achieve better numerical accuracy, convergence speed and stability.

Symplectic scheme for the Schrödinger equation with fractional Laplacian

In this paper, the symplectic scheme is presented for solving the space fractional Schrödinger equation (SFSE) with one dimension. First, the symplectic conservation laws are investigated for space semi-discretization systems of the SFSE based on the existing second-order central difference scheme and the existing fourth-order compact scheme. Then, a fourth-order central difference scheme is developed in space discretization, and the resulting semi-discretization system is shown to be a finite dimension Hamiltonian system of ordinary differential equations. Moreover, we get the full discretization scheme for the Hamiltonian system by symplectic midpoint scheme in time direction. In particular, the space semi-discretization and the full discretization are shown to preserve some properties of the SFSE. At last, numerical experiments are given to verify the efficiency of the scheme, and show that these symplectic difference schemes can be applied to long time simulation for one-dimensional SFSEs.

Improvement of numerical simulation for GMAW based on a new model with thermocapillary effect

In this paper, we present a numerical investigation about the metal transfer of GMAW with a new model based on the phase field model. Unlike most published work, the thermocapillary effect and mixture energy resulting from the new research of multiphase fluids are taken into model of phase transfer and interface change which is different from volume of fluid method (VOF) that has been widely used in welding problem. We discretize the whole model including continuity condition, Navier–Stokes equation, phase field model, energy equation and Maxwell’s equations with a continuous finite element method in space and a midpoint scheme in time and a penalty formulation is applied to the continuity condition due to the stability in the pressure. Metal transfer of GMAW with pulse current is computed as a numerical example where the numerical result have been verified by the experimental data obtained by using a high-speed camera. Moreover, we also compare the numerical results with and without considering the thermocapillary effects to show that the new model has a higher precision in predicting the droplet in welding. • We present a new model with thermocapillary, mix-energy, phase field model for GMAW. • The accuracy of the new model is examined with the data of high-speed photography. • The new model has a higher precision in predicting the metal transfer than before.

Numerical Analysis of the Midpoint Scheme for the Generalized Benjamin-Bona-Mahony Equation with White Noise Dispersion

Numerical Analysis of the Midpoint Scheme for the Generalized Benjamin-Bona-Mahony Equation with White Noise Dispersion

2D continuum viscodamage-embedded discontinuity model with second order mid-point scheme

This paper deals with numerical modeling of dynamic failure phenomena in rate-sensitive brittle and/or ductile materials. To this end, a two-dimensional continuum viscodamage-embedded discontinuity model, which is based on our previous work (see Do et al. 2017), is developed. More specifically, the pre-peak nonlinear and rate-sensitive hardening response of the material behavior, representing the fracture-process zone creation, is described by a rate-dependent continuum damage model. Meanwhile, an embedded displacement discontinuity model is used to formulate the post-peak response, involving the macro-crack creation accompanied by exponential softening. The numerical implementation in the context of the finite element method exploiting the second-order mid-point scheme is discussed in detail. In order to show the performance of the model several numerical examples are included.

A novel kind of efficient symplectic scheme for Klein–Gordon–Schrödinger equation

In this paper, we construct a family of high order compact symplectic (S-HOC) schemes for the Klein–Gordon–Schrödinger (KGS) equation. The KGS can be cast into a Hamiltonian form. At first, we discretize the Hamiltonian system in space by a high order compact method which has higher convergent rate than general finite difference methods. Then the semi-discretized system is approximated in time by the Euler midpoint scheme which preserves the symplectic structure of the original system. The conserved quantities of the scheme, including symplectic structure conservation law, charge conservation law and energy conservation law, are discussed. The local truncation error and global error of the numerical solvers are investigated. Finally, some numerical verifications are presented to numerically validate the theoretical analysis. The numerical results are persuasive and illustrate the theoretical analysis.

Generalized, Energy-conserving Numerical Simulations of Particles in General Relativity. I. Time-like and Null Geodesics

The numerical integration of particle trajectories in curved spacetimes is fundamental for obtaining realistic models of the particle dynamics around massive compact objects such as black holes and neutron stars. Generalized algorithms capable of handling generic metrics are required for studies of both standard spacetimes (Schwarzschild and Kerr metrics) and non-standard spacetimes (e.g. Schwarzschild metric plus non-classical perturbations or multiple black hole metrics). The most commonly employed explicit numerical schemes (e.g. Runge-Kutta) are incapable of producing highly accurate results at critical points, e.g. in the regions close to the event horizon where gravity causes extreme curvature of the spacetime, at an acceptable computational cost. Here, we describe a generalized algorithm for the numerical integration of time-like (massive particles) and null (photons) geodesics in any given 3+1 split spacetime. We introduce a new, exactly energy-conserving implicit integration scheme based on the preservation of the underlying Hamiltonian, and we compare its properties with a standard fourth-order Runge-Kutta explicit scheme and an implicit midpoint scheme. We test the numerical performance of the three schemes against analytic solutions of particle and photon orbits in Schwarzschild and Kerr spacetimes. We also prove the versatility of our framework in handling more exotic metrics such as Morris-Thorne wormholes and quantum-perturbed Schwarzschild black holes. The generalized approach is also discussed in the perspective of future extensions to more complex particle dynamics, e.g. the addition of the Lorentz force acting on charged particles, which allows for test particle diagnostics in GRMHD simulations.

Numerical simulation and convergence analysis for a system of nonlinear singularly perturbed differential equations arising in population dynamics

In this article, we consider a system of nonlinear singularly perturbed differential equations with two different parameters. To solve this system, we develop a weighted monotone hybrid scheme on a nonuniform mesh. The proposed scheme is a combination of the midpoint scheme and the upwind scheme involving the weight parameters. The weight parameters enable the method to switch automatically from the midpoint scheme to the upwind scheme as the nodal points start moving from the inner region to the outer region. The nonuniform mesh in particular the adaptive grid is constructed using the idea of equidistributing a positive monitor function involving the solution gradient. The method is shown to be second order convergent with respect to the small parameters. Numerical experiments are presented to show the robustness of the proposed scheme and indicate that the estimate is optimal.

Effects of Brownian motion and thermophoresis on nanofluids in a rotating circular groove: A numerical simulation

We present a numerical study of the movement of nanoparticles (Cu) and heat transfer of nanofluids with Newtonian and non-Newtonian base fluids in this paper. The computing domain is a circular groove with a fixed slow speed of revolution. Two kinds of models with different definitions of thermal conductivity and concentrations of nanoparticles are considered in the numerical calculation process, as are the influences of thermophoresis and Brownian motion. A continuous finite element scheme in space and a modified midpoint scheme in time are used, and numerical solutions are obtained by using the finite element method. Compared with base fluids only, the heat transfer enhancement of nanofluids can be found easily and is more intuitive. The correlation of the heat transfer enhancement and power law index of base fluids used to distinguish different power law fluids is investigated. This contributes to determining the reasons for heat transfer enhancement in nanofluids. It is found that the Maxwell Model has a higher enhancement of heat transfer than the Traditional Model and the thermophoresis effect breaks the dynamic equilibrium of nanoparticles in the rotating circular groove.

Symplectic Runge–Kutta methods for Hamiltonian systems driven by Gaussian rough paths

We consider Hamiltonian systems driven by multi-dimensional Gaussian processes in rough path sense, which include fractional Brownian motions with Hurst parameter H∈(1/4,1/2]. We prove that the phase flow preserves the symplectic structure almost surely and this property could be inherited by symplectic Runge–Kutta methods, which are implicit methods in general. If the vector fields satisfy some smoothness and boundedness conditions, we obtain the pathwise convergence rates of Runge–Kutta methods. When vector fields are linear, we get the solvability of the midpoint scheme for skew symmetric cases, and obtain its pathwise convergence rate. Numerical experiments verify our theoretical analysis.

Trajectory errors of different numerical integration schemes diagnosed with the MPTRAC advection module driven by ECMWF operational analyses

Abstract. The accuracy of trajectory calculations performed by Lagrangian particle dispersion models (LPDMs) depends on various factors. The optimization of numerical integration schemes used to solve the trajectory equation helps to maximize the computational efficiency of large-scale LPDM simulations. We analyzed global truncation errors of six explicit integration schemes of the Runge–Kutta family, which we implemented in the Massive-Parallel Trajectory Calculations (MPTRAC) advection module. The simulations were driven by wind fields from operational analysis and forecasts of the European Centre for Medium-Range Weather Forecasts (ECMWF) at T1279L137 spatial resolution and 3 h temporal sampling. We defined separate test cases for 15 distinct regions of the atmosphere, covering the polar regions, the midlatitudes, and the tropics in the free troposphere, in the upper troposphere and lower stratosphere (UT/LS) region, and in the middle stratosphere. In total, more than 5000 different transport simulations were performed, covering the months of January, April, July, and October for the years 2014 and 2015. We quantified the accuracy of the trajectories by calculating transport deviations with respect to reference simulations using a fourth-order Runge–Kutta integration scheme with a sufficiently fine time step. Transport deviations were assessed with respect to error limits based on turbulent diffusion. Independent of the numerical scheme, the global truncation errors vary significantly between the different regions. Horizontal transport deviations in the stratosphere are typically an order of magnitude smaller compared with the free troposphere. We found that the truncation errors of the six numerical schemes fall into three distinct groups, which mostly depend on the numerical order of the scheme. Schemes of the same order differ little in accuracy, but some methods need less computational time, which gives them an advantage in efficiency. The selection of the integration scheme and the appropriate time step should possibly take into account the typical altitude ranges as well as the total length of the simulations to achieve the most efficient simulations. However, trying to summarize, we recommend the third-order Runge–Kutta method with a time step of 170 s or the midpoint scheme with a time step of 100 s for efficient simulations of up to 10 days of simulation time for the specific ECMWF high-resolution data set considered in this study. Purely stratospheric simulations can use significantly larger time steps of 800 and 1100 s for the midpoint scheme and the third-order Runge–Kutta method, respectively.

Rapid sampling of stochastic displacements in Brownian dynamics simulations with stresslet constraints.

Brownian Dynamics simulations are an important tool for modeling the dynamics of soft matter. However, accurate and rapid computations of the hydrodynamic interactions between suspended, microscopic components in a soft material are a significant computational challenge. Here, we present a new method for Brownian dynamics simulations of suspended colloidal scale particles such as colloids, polymers, surfactants, and proteins subject to a particular and important class of hydrodynamic constraints. The total computational cost of the algorithm is practically linear with the number of particles modeled and can be further optimized when the characteristic mass fractal dimension of the suspended particles is known. Specifically, we consider the so-called "stresslet" constraint for which suspended particles resist local deformation. This acts to produce a symmetric force dipole in the fluid and imparts rigidity to the particles. The presented method is an extension of the recently reported positively split formulation for Ewald summation of the Rotne-Prager-Yamakawa mobility tensor to higher order terms in the hydrodynamic scattering series accounting for force dipoles [A. M. Fiore et al., J. Chem. Phys. 146(12), 124116 (2017)]. The hydrodynamic mobility tensor, which is proportional to the covariance of particle Brownian displacements, is constructed as an Ewald sum in a novel way which guarantees that the real-space and wave-space contributions to the sum are independently symmetric and positive-definite for all possible particle configurations. This property of the Ewald sum is leveraged to rapidly sample the Brownian displacements from a superposition of statistically independent processes with the wave-space and real-space contributions as respective covariances. The cost of computing the Brownian displacements in this way is comparable to the cost of computing the deterministic displacements. The addition of a stresslet constraint to the over-damped particle equations of motion leads to a stochastic differential algebraic equation (SDAE) of index 1, which is integrated forward in time using a mid-point integration scheme that implicitly produces stochastic displacements consistent with the fluctuation-dissipation theorem for the constrained system. Calculations for hard sphere dispersions are illustrated and used to explore the performance of the algorithm. An open source, high-performance implementation on graphics processing units capable of dynamic simulations of millions of particles and integrated with the software package HOOMD-blue is used for benchmarking and made freely available in the supplementary material.

Large scale Brownian dynamics of confined suspensions of rigid particles.

We introduce methods for large-scale Brownian Dynamics (BD) simulation of many rigid particles of arbitrary shape suspended in a fluctuating fluid. Our method adds Brownian motion to the rigid multiblob method [F. Balboa Usabiaga et al., Commun. Appl. Math. Comput. Sci. 11(2), 217-296 (2016)] at a cost comparable to the cost of deterministic simulations. We demonstrate that we can efficiently generate deterministic and random displacements for many particles using preconditioned Krylov iterative methods, if kernel methods to efficiently compute the action of the Rotne-Prager-Yamakawa (RPY) mobility matrix and its "square" root are available for the given boundary conditions. These kernel operations can be computed with near linear scaling for periodic domains using the positively split Ewald method. Here we study particles partially confined by gravity above a no-slip bottom wall using a graphical processing unit implementation of the mobility matrix-vector product, combined with a preconditioned Lanczos iteration for generating Brownian displacements. We address a major challenge in large-scale BD simulations, capturing the stochastic drift term that arises because of the configuration-dependent mobility. Unlike the widely used Fixman midpoint scheme, our methods utilize random finite differences and do not require the solution of resistance problems or the computation of the action of the inverse square root of the RPY mobility matrix. We construct two temporal schemes which are viable for large-scale simulations, an Euler-Maruyama traction scheme and a trapezoidal slip scheme, which minimize the number of mobility problems to be solved per time step while capturing the required stochastic drift terms. We validate and compare these schemes numerically by modeling suspensions of boomerang-shaped particles sedimented near a bottom wall. Using the trapezoidal scheme, we investigate the steady-state active motion in dense suspensions of confined microrollers, whose height above the wall is set by a combination of thermal noise and active flows. We find the existence of two populations of active particles, slower ones closer to the bottom and faster ones above them, and demonstrate that our method provides quantitative accuracy even with relatively coarse resolutions of the particle geometry.

A New Model for Simulating Gas Metal Arc Welding based on Phase Field Model

Lots of physical process, such as metal melting, multiphase fluids flow, heat and mass transfer and thermocapillary effect (Marangoni) and so on, will occur in gas metal arc welding (GMAW) which should be considered as a mixture system. In this paper, based on the previous work, we propose a new model to simulate GMAW including Navier-Stokes equation, the phase field model and energy equation. Unlike most previous work, we take the thermocapillary effect into the phase field model considering mixture energy which is different of volume of fluid method (VOF) widely used in GMAW before. We also consider gravity, electromagnetic force, surface tension, buoyancy effect and arc pressure in momentum equation. The spray transfer especially the projected transfer in GMAW is computed as numerical examples with a continuous finite element method and a modified midpoint scheme. Pulse current is set as welding current as the numerical example to show the numerical simulation of metal transfer which fits the theory of GMAW well. From the result compared with the data of high-speed photography and VOF model, the accuracy and stability of the model and scheme are easily validated and also the new model has the higher precieion.

An enhanced SDOF model to predict the behaviour of a steel column impacted by a rigid body

The transient dynamic response of a steel beam-column subjected to impact loading is a complex phenomenon involving large localized plastic deformations and non-smooth contact interactions. Exposed to high intensity of the contact force generated from impact, the beam-column may undergo large displacement and inelastic deformation. Previous research has shown that a calibrated elasto-plastic single degree of freedom system is able to reproduce both the displacement and the force time-history of a steel beam subjected to non-impulsive loading or low-velocity impact. In these models, the static force-displacement curve is derived from either experiments or detailed 3D nonlinear analysis. Tri-linear resistance function has been extensively used to reproduce the different stages of the response including catenary effects. A rigorous treatment of such a complex problem calls for the use of non-smooth analysis tools to handle the impulsive nature of the impact force, the unilateral constraint, the impenetrability condition and the discontinuity of the velocity in a rigorous manner. In this paper, we present a non-smooth elasto-plastic single degree of freedom model under impact loading that permits the use of arbitrary resistance function. Adopting the non-smooth framework offers tools such as differential measures and convex analysis concepts to deal with unilateral contact incorporating Newton’s impact law. The mid-point scheme is adopted to avoid numerical unrealistic energy decay or blowup. Furthermore, the non-penetration condition is numerically satisfied by imposing the constraint at only the velocity level to guarantee energy-momentum conservation [1]. The explicit expression of resistance functions of the beam that are used in the SDOF model are obtained from a simplified nonlinear static analysis of two beam-column models. In the analysis, a linear relation between normal force and bending moment is assumed for the plastification of the hinges. Two proposals to simplify the explicit expressions of the model’s response behavior are given. Performing an energy-based analysis, we predict maximum displacement that is needed to absorb the kinetic energy arising from the impact for different coefficient of restitution. The numerical examples underline the validity of the model by showing good agreement with the predictions of reference models.