Finite Time Step Research Articles

Efficient event-driven simulations shed new light on microtubule organization in the plant cortical array

The dynamics of the plant microtubule cytoskeleton is a paradigmatic example of the complex spatiotemporal processes characterising life at the cellular scale. This system is composed of large numbers of spatially extended particles, each endowed with its own intrinsic stochastic dynamics, and is capable of non-equilibrium self-organisation through collisional interactions of these particles. To elucidate the behaviour of such a complex system requires not only conceptual advances, but also the development of appropriate computational tools to simulate it. As the number of parameters involved is large and the behaviour is stochastic, it is essential that these simulations be fast enough to allow for an exploration of the phase space and the gathering of sufficient statistics to accurately pin down the average behaviour as well as the magnitude of fluctuations around it. Here we describe a simulation approach that meets this requirement by adopting an event-driven methodology that encompasses both the spontaneous stochastic changes in microtubule state as well as the deterministic collisions. In contrast with finite time step simulations this technique is intrinsically exact, as well as several orders of magnitude faster, which enables ordinary PC hardware to simulate systems of $\sim 10^3$ microtubules on a time scale $\sim 10^{3}$ faster than real time. In addition we present new tools for the analysis of microtubule trajectories on curved surfaces. We illustrate the use of these methods by addressing a number of outstanding issues regarding the importance of various parameters on the transition from an isotropic to an aligned and oriented state.

Numerical Langevin Simulations: Equilibrium Dynamics and Nonequilibrium Thermodynamics

Stochastic differential equations of motion (e.g., Langevin dynamics) provide a popular framework for simulating molecular systems. Any computational algorithm must discretize these equations, yet the resulting finite time step integration schemes suffer from several practical shortcomings: for example, equilibrium dynamical properties diverge from those of the continuous equations of motion; path-sampling applications are impeded by a path action that is either unknown or cumbersome; and thermodynamic statistics of driven systems are skewed, leading to biased inference using recently-developed nonequilibrium fluctuation theorems. We show how any finite time step Langevin integrator, even with a time-independent Hamiltonian, can be thought of as a driven, nonequilibrium physical process. Once an appropriate work-like quantity (the shadow work) is defined, nonequilibrium fluctuation theorems can characterize or correct for the errors introduced by the use of finite time steps. In particular, we demonstrate that amending estimators based on nonequilibrium work theorems to include this shadow work removes the time step dependent error from estimates of free energies. We also quantify, for the first time, the magnitude of deviations between the sampled stationary distribution and the desired equilibrium distribution for equilibrium Langevin simulations of solvated systems of varying size. What is more, we show that the incorporation of a novel time step rescaling in the deterministic updates of position and velocity can correct a number of dynamical defects in these integrators. Finally, we identify a particular splitting that has essentially universally appropriate properties for the simulation of Langevin dynamics for molecular systems in equilibrium, nonequilibrium, and path sampling contexts.

Dynamic impact of a particle

In this work, we consider a moving particle which drops down onto a stationary rigid foundation and bounces off after its contact. The equation of its motion is formulated by a second-order ordinary differential equation. The particle satisfies the Signorini contact conditions which can be interpreted in terms of complementarity conditions. The existence of weak solutions is shown by using a finite time step and the necessary a priori estimates which allow us to pass to the limit. The uniqueness of the solutions can be proved under some additional assumptions. Conservation of energy is also investigated theoretically and numerically. Numerical solutions are computed via both finite- and infinite-dimensional approaches.

Dynamic Stability of Plane Free Surface of Liquid in Axisymmetric Tanks

When liquid filled containers are excited vertically, it is known that, for some combinations of frequency and amplitude, the free surface undergoes unbounded motion leading to instability, called parametric instability or parametric resonance, while for other combinations the free surface remains plane. In this paper, the stability of the plane free surface is investigated theoretically when the vessel is a vertical axisymmetric container. The effect of coupled horizontal excitation on the stability is examined. The dynamics of sloshing flows under specified excitations are simulated numerically using fully nonlinear finite element method based on non-linear potential flow theory. A mixed Eulerian-Lagrangian technique combined with 4th-order Runge-Kutta method is employed to advance the solution in time. A regridding technique based on cubic spline is applied to the free surface for every finite time step to avoid possible numerical instabilities.

The intrinsic dependence structure of peak, volume, duration, and average intensity of hyetographs and hydrographs.

[1] The information contained in hyetographs and hydrographs is often synthesized by using key properties such as the peak or maximum value Xp, volume V, duration D, and average intensity I. These variables play a fundamental role in hydrologic engineering as they are used, for instance, to define design hyetographs and hydrographs as well as to model and simulate the rainfall and streamflow processes. Given their inherent variability and the empirical evidence of the presence of a significant degree of association, such quantities have been studied as correlated random variables suitable to be modeled by multivariate joint distribution functions. The advent of copulas in geosciences simplified the inference procedures allowing for splitting the analysis of the marginal distributions and the study of the so-called dependence structure or copula. However, the attention paid to the modeling task has overlooked a more thorough study of the true nature and origin of the relationships that link , and I. In this study, we apply a set of ad hoc bootstrap algorithms to investigate these aspects by analyzing the hyetographs and hydrographs extracted from 282 daily rainfall series from central eastern Europe, three 5 min rainfall series from central Italy, 80 daily streamflow series from the continental United States, and two sets of 200 simulated universal multifractal time series. Our results show that all the pairwise dependence structures between , and I exhibit some key properties that can be reproduced by simple bootstrap algorithms that rely on a standard univariate resampling without resort to multivariate techniques. Therefore, the strong similarities between the observed dependence structures and the agreement between the observed and bootstrap samples suggest the existence of a numerical generating mechanism based on the superposition of the effects of sampling data at finite time steps and the process of summing realizations of independent random variables over random durations. We also show that the pairwise dependence structures are weakly dependent on the internal patterns of the hyetographs and hydrographs, meaning that the temporal evolution of the rainfall and runoff events marginally influences the mutual relationships of , and I. Finally, our findings point out that subtle and often overlooked deterministic relationships between the properties of the event hyetographs and hydrographs exist. Confusing these relationships with genuine stochastic relationships can lead to an incorrect application of multivariate distributions and copulas and to misleading results.

A Theory in Finite Differences for Computation of Interchanging Principal Stresses and Principal Axes Rotation

Suppose consecutive ordinates of three arbitrary and nonproportional stress-time functions of plane state of stress are entered into a computer by a little finite time step. The theory proposed solves the following problem: correct computation of the ordinates of the principal stress-time functions and the angle of principal axes rotation. This problem is not as simple as researchers approached it prior to the computer era. First of all, the correct solution for the principal stresses and the principal axes rotation require correct interchange of the principal stresses while computing them, i.e., correct interchange of the plus/minus signs in the well-known equations for them. For the interchange analysis, an ellipse of stress transformation in the three-dimensional stress-coordinate space is revealed. By changing a coordinate scale, the ellipse turns into a circumference that is an analog to, but different from, a Mohr circle. The correct solution also requires treating the principal axes rotation in little finite differences per little time differences during which little finite elements appear as building the stressing path in the three-dimensional stress-coordinate space. Based on the ellipse/circumference mentioned, three interchange conditions are revealed. The third one is the most important. And, a necessity is also revealed for dividing some stressing path's elements into two subelements. Based on all the findings, the main commands of an algorithm for computing the ordinates of the principal stress-time functions and the angle of principal axes rotation are presented. The correct solution of the problem has been achieved thanks to new notions taken from the so-called integration of damage differentials (IDD) theory. In fact, the paper presents a new contribution to the variable plane stress state analysis.

A domain decomposition method for pseudo-spectral electromagnetic simulations of plasmas

Pseudo-spectral electromagnetic solvers (i.e. representing the fields in Fourier space) have extraordinary precision. In particular, Haber et al. presented in 1973 a pseudo-spectral solver that integrates analytically the solution over a finite time step, under the usual assumption that the source is constant over that time step. Yet, pseudo-spectral solvers have not been widely used, due in part to the difficulty for efficient parallelization owing to global communications associated with global FFTs on the entire computational domains.A method for the parallelization of electromagnetic pseudo-spectral solvers is proposed and tested on single electromagnetic pulses, and on Particle-In-Cell simulations of the wakefield formation in a laser plasma accelerator.The method takes advantage of the properties of the Discrete Fourier Transform, the linearity of Maxwell’s equations and the finite speed of light for limiting the communications of data within guard regions between neighboring computational domains.Although this requires a small approximation, test results show that no significant error is made on the test cases that have been presented.The proposed method opens the way to solvers combining the favorable parallel scaling of standard finite-difference methods with the accuracy advantages of pseudo-spectral methods.

DESIGN OF A LOAD FOLLOWING CONTROLLER FOR APR+ NUCLEAR PLANTS

A load-following operation in APR+ nuclear plants is necessary to reduce the need to adjust the boric acid concentration and to efficiently control the control rods for flexible operation. In particular, a disproportion in the axial flux distribution, which is normally caused by a load-following operation in a reactor core, causes xenon oscillation because the absorption cross-section of xenon is extremely large and its effects in a reactor are delayed by the iodine precursor. A model predictive control (MPC) method was used to design an automatic load-following controller for the integrated thermal power level and axial shape index (ASI) control for APR+ nuclear plants. Some tracking controllers employ the current tracking command only. On the other hand, the MPC can achieve better tracking performance because it considers future commands in addition to the current tracking command. The basic concept of the MPC is to solve an optimization problem for generating finite future control inputs at the current time and to implement as the current control input only the first control input among the solutions of the finite time steps. At the next time step, the procedure to solve the optimization problem is then repeated. The support vector regression (SVR) model that is used widely for function approximation problems is used to predict the future outputs based on previous inputs and outputs. In addition, a genetic algorithm is employed to minimize the objective function of a MPC control algorithm with multiple constraints. The power level and ASI are controlled by regulating the control banks and part-strength control banks together with an automatic adjustment of the boric acid concentration. The 3-dimensional MASTER code, which models APR+ nuclear plants, is interfaced to the proposed controller to confirm the performance of the controlling reactor power level and ASI. Numerical simulations showed that the proposed controller exhibits very fast tracking responses.

A viscoelastic Timoshenko beam with Coulomb law of friction

This work focuses on the analysis and numerical simulations of dynamic frictional contact between a Timoshenko beam and a stationary rigid obstacle. The beam is assumed to be viscoelastic, and is clamped at its left end while it is free at its right end. When the beam moves horizontally and its right end contacts the rigid obstacle, contact forces arise and then the vertical motion of its right end is considered with friction. Thus the right end of the beam satisfies two contact conditions: the Signorini unilateral contact condition and Coulomb law of dry friction. Moreover, the slip rate dependence of the coefficient of friction is taken into account. The existence of weak solutions to the problem is shown by using finite time stepping and the necessary a priori estimates that allow us to vanish the time step in the limiting process. The energy balance in the system is investigated both theoretically and numerically. A fully discrete numerical scheme for the variational formulation of the problem is implemented and numerical results are presented.

An incremental variational formulation of dissipative magnetostriction at the macroscopic continuum level

This paper outlines a new variational-based modeling and computational implementation of macroscopic continuum magneto-mechanics involving non-linear, inelastic material behavior, with a special focus on dissipative magnetostriction. It is based on a constitutive variational principle that optimizes a generalized incremental work function with respect to the internal state variables. In an incremental setting at finite time steps, this variational problem defines a quasi-hyper-magnetoelastic potential for the stresses and the magnetic induction, and incorporates energy storage as well as dissipative mechanisms. The existence of this potential further allows the incremental boundary-value problem of quasi-static inelastic magneto-mechanics to be recast into a principle of stationary incremental energy. The second focus of this paper is on the careful construction of the energy storage and dissipation functions for the model problem of hysteretic magnetostriction at the macroscopic level. It is then demonstrated that the proposed model is capable of predicting the ferromagnetic and field-induced strain hysteresis curves characteristic of magnetostrictive material response in good agreement with experiments. The numerical solution of the coupled non-linear boundary-value problem is based on a monolithic multi-field finite element implementation. As a consequence of the proposed incremental variational principle, the discretization of the multi-field problem appears in a compact symmetric format. In this sense, the proposed formulation provides a canonical framework for the simulation of boundary-value-problems in dissipative magnetostriction at the macro-level. The performance of the proposed algorithm is tested by application to relevant numerical examples.

On evolving deformation microstructures in non-convex partially damaged solids

The paper outlines a relaxation method based on a particular isotropic microstructure evolution and applies it to the model problem of rate independent, partially damaged solids. The method uses an incremental variational formulation for standard dissipative materials. In an incremental setting at finite time steps, the formulation defines a quasi-hyperelastic stress potential. The existence of this potential allows a typical incremental boundary value problem of damage mechanics to be expressed in terms of a principle of minimum incremental work. Mathematical existence theorems of minimizers then induce a definition of the material stability in terms of the sequential weak lower semicontinuity of the incremental functional. As a consequence, the incremental material stability of standard dissipative solids may be defined in terms of weak convexity notions of the stress potential. Furthermore, the variational setting opens up the possibility to analyze the development of deformation microstructures in the post-critical range of unstable inelastic materials based on energy relaxation methods. In partially damaged solids, accumulated damage may yield non-convex stress potentials which indicate instability and formation of fine-scale microstructures. These microstructures can be resolved by use of relaxation techniques associated with the construction of convex hulls. We propose a particular relaxation method for partially damaged solids and investigate it in one- and multi-dimensional settings. To this end, we introduce a new isotropic microstructure which provides a simple approximation of the multi-dimensional rank-one convex hull. The development of those isotropic microstructures is investigated for homogeneous and inhomogeneous numerical simulations.

A Model to Simulate Pumping Tests in Wells with High Loss and Backflow

Traditional methods of analyzing pumping tests in single wells fail when the well loss is very high due to a low transmissivity skin. Because of the restricted rate at which water can enter a high loss well from the aquifer, well casing storage becomes a significant factor. Additionally, if a slug of water enters the well from the pump column immediately after the pump is switched off, it has a long-lasting significant effect on the recovering water level in the well because it cannot be absorbed rapidly by the aquifer. A theoretical model is derived here that simulates the water level in a well in these circumstances. In the model, the continuously changing rate of water inflow from the aquifer to the well is approximated by a step function with a finite difference time step. It is demonstrated by a real example that the model can be applied easily to analyze pumping tests, including tests with a varying pumping rate. The analysis confirms suspected high well loss, calculates the unknown rate of backflow, and determines the aquifer's transmissivity.

Finite difference based vibration simulation analysis of a segmented distributedpiezoelectric structronic plate system

Electrical modeling of piezoelectric structronic systems by analog circuits has thedisadvantages of huge circuit structure and low precision. However, studies of electricalsimulation of segmented distributed piezoelectric structronic plate systems (PSPSs) byusing output voltage signals of high-speed digital circuits to evaluate the real-time dynamicdisplacements are scarce in the literature. Therefore, an equivalent dynamic modelbased on the finite difference method (FDM) is presented to simulate the actualphysical model of the segmented distributed PSPS with simply supported boundaryconditions. By means of the FDM, the four-ordered dynamic partial differentialequations (PDEs) of the main structure/segmented distributed sensor signals/controlmoments of the segmented distributed actuator of the PSPS are transformed tofinite difference equations. A dynamics matrix model based on the Newmark-β integration method is established. The output voltage signal characteristics of the lower modes (m ≤ 3,n ≤ 3) with different finite difference mesh dimensions and different integration time steps areanalyzed by digital signal processing (DSP) circuit simulation software. The control effectsof segmented distributed actuators with different effective areas are consistent with theresults of the analytical model in relevant references. Therefore, the method of digitalsimulation for vibration analysis of segmented distributed PSPSs presented in this papercan provide a reference for further research into the electrical simulation of PSPSs.

Particle-in-cell modeling of relativistic laser–plasma interaction with the adjustable-damping, direct implicit method

Implicit particle-in-cell codes offer advantages over their explicit counterparts in that they suffer weaker stability constraints on the need to resolve the higher frequency modes of the system. This feature may prove particularly valuable for modeling the interaction of high-intensity laser pulses with overcritical plasmas, in the case where the electrostatic modes in the denser regions are of negligible influence on the physical processes under study. To this goal, we have developed the new two-dimensional electromagnetic code ELIXIRS (standing for ELectromagnetic Implicit X-dimensional Iterative Relativistic Solver) based on the relativistic extension of the so-called Direct Implicit Method [D. Hewett, A.B. Langdon, Electromagnetic direct implicit plasma simulation, J. Comput. Phys. 72 (1987) 121–155]. Dissipation-free propagation of light waves into vacuum is achieved by an adjustable-damping electromagnetic solver. In the high-density case where the Debye length is not resolved, satisfactory energy conservation is ensured by the use of high-order weight factors. In this paper, we first derive the electromagnetic direct implicit method as a simplified Newton scheme. Its linear properties are then investigated through numerically solving the relation dispersions obtained for both light and plasma waves, accounting for finite space and time steps. Finally, our code is successfully benchmarked against explicit particle-in-cell simulations for two kinds of physical problems: plasma expansion into vacuum and relativistic laser–plasma interaction. In both cases, we will demonstrate the robustness of the implicit solver for crude discretizations, as well as the gains in efficiency which can be realized over standard explicit simulations.

Computed chaos or numerical errors

Discrete numerical methods with finite time-steps represent a practical technique to solve initial-value problems involving nonlinear differential equations. These methods seem particularly useful to the study of chaos since no analytical chaotic solution is currently available. Using the well-known Lorenz equations as an example, it is demonstrated that numerically computed results and their associated statistical properties are time-step dependent. There are two reasons for this behavior. First, chaotic differential equations are unstable so that any small error is amplified exponentially near an unstable manifold. The more serious and lesser-known reason is that stable and unstable manifolds of singular points associated with differential equations can form virtual separatrices. The existence of a virtual separatrix presents the possibility of a computed trajectory actually “jumping” through it due to the finite time-steps of discrete numerical methods. Such behavior violates the uniqueness theory of differential equations and amplifies the numerical errors explosively. These reasons imply that, even if computed results are bounded, their independence on time-step should be established before accepting them as useful numerical approximations to the true solution of the differential equations. However, due to these exponential and explosive amplifications of numerical errors, no computed chaotic solutions of differential equations independent of integration-time step have been found. Thus, reports of computed non-periodic solutions of chaotic differential equations are simply consequences of unstably amplified truncation errors, and are not approximate solutions of the associated differential equations.

A framework for multiple robust explicit MPC controllers for linear systems

This article presents a multi-mode explicit controller for constrained linear systems with bounded disturbances using a switching strategy based on Model Predictive Control (MPC). In the proposed approach, the system switches among several MPC controllers having different performance levels. The switching is done so as to achieve increasing levels of performance as time evolves, reaching the desired controller in finite time steps. The conditions needed for switching and robust convergence of the multi-mode MPC controllers are provided. Compared with standard robust explicit MPC implementations, the proposed approach has the flexibility of having a large domain of attraction, a good asymptotic behaviour and a small number of partitions.

An explicit dynamics extended finite element method. Part 2: Element-by-element stable-explicit/explicit dynamic scheme

This paper presents a general explicit time integration scheme for dynamics simulations using the eXtended Finite Element Method with standard critical time step. We use the generalized mass lumping technique proposed in Part I of this paper. This technique allows us to consider arbitrary enrichment functions in the X-FEM for explicit dynamics simulations. In this second part, the proposed approach allows the use of standard finite elements critical time step estimates. For that purpose, we develop a classical element-by-element strategy that couples the standard central difference scheme with the unconditionally stable-explicit scheme proposed by Chang [S.Y. Chang, An explicit method with improved stability property, Int. J. Numer. Methods Engrg. 77 (8) (2008) 1100–1120]. This scheme coupled with X-FEM allows us to recover FE critical time steps independently of the enrichment functions considered. Furthermore, a study of the stability property of this new explicit scheme is proposed in the X-FEM framework, both with fixed enrichments and evolving enrichments with time. Some examples illustrate the good properties of the new explicit numerical time scheme and some applications to dynamics crack growth are given.

Using genetic algorithms to characterize ferrofluid topographies in externally applied magnetic fields

Both ferrofluidics and genetic algorithms are relatively new fields. Due to complex physical interactions, ferrofluidic topographies and assemblies have only been solved using finite time step, Lattice Boltzmann, and finite-element methods in very simple magnetic field configurations. In this paper, we show that it is possible (and highly advantageous) to employ genetic algorithms to solve for the fluid topographies, which can be extended to include more complex magnetic fields.

A symplectic integrator with arbitrary vector and scalar potentials

We study a new class of symplectic integrators for particles in arbitrary, time-dependent vector and scalar potentials. The methods were introduced in [Y.K. Wu, E. Forest, D.S. Robin, Phys. Rev. E 68 (2003) 046502] and are based on the ability to integrate Hamiltonians of the form ( p i − a i ( q ) ) 2 exactly for a finite time-step. We show that the integrators are symplectic in the non-relativistic case but not symplectic in the full six-dimensional phase space for relativistic Hamiltonians.

A Closed-Form Approximation of Likelihood Functions for Discretely Sampled Diffusions: The Exponent Expansion

In this paper we discuss a closed-form approximation of the likelihood functions of an arbitrary diffusion process. The approximation is based on an exponential ansatz of the transition probability for a finite time step $\Delta t$, and a series expansion of the deviation of its logarithm from that of a Gaussian distribution. Through this procedure, dubbed {\em exponent expansion}, the transition probability is obtained as a power series in $\Delta t$. This becomes asymptotically exact if an increasing number of terms is included, and provides remarkably accurate results even when truncated to the first few (say 3) terms. The coefficients of such expansion can be determined straightforwardly through a recursion, and involve simple one-dimensional integrals. We present several examples of financial interest, and we compare our results with the state-of-the-art approximation of discretely sampled diffusions [Sahalia, Journal of Finance, Vol. 54, p. 1361 (1999)]. We find that the exponent expansion provides a similar accuracy in most of the cases, but a better behavior in the low-volatility regime. Furthermore the implementation of the present approach turns out to be simpler. Within the functional integration framework the exponent expansion allows one to obtain remarkably good approximations of the pricing kernels of financial derivatives. This is illustrated with the application to simple path-dependent interest rate derivatives. Finally we discuss how these results can also be used to increase the efficiency of numerical (both deterministic and stochastic) approaches to derivative pricing.