Arbitrary Time Steps Research Articles

A 4-Stated DICE: Quantitatively Addressing Uncertainty Effects in Climate Change

We introduce a version of the DICE-2007 model designed for uncertainty analysis. DICE is a wide-spread deterministic integrated assessment model of climate change. Climate change, long-term economic development, and their interactions are highly uncertain. The quantitative analysis of optimal mitigation policy under uncertainty requires a recursive dynamic programming implementation of integrated assessment models. Such implementations are subject to the curse of dimensionality. Every increase in the dimension of the state space is paid for by a combination of (exponentially) increasing processor time, lower quality of the value or policy function approximations, and reductions of the uncertainty domain. The paper promotes a state-reduced, recursive dynamic programming implementation of the DICE-2007 model. We achieve the reduction by simplifying the carbon cycle and the temperature delay equations. We compare our model’s performance and that of the DICE model to the scientific AOGCM models emulated by MAGICC 6.0 and find that our simplified model performs equally well as the original DICE model. Our implementation solves the infinite planning horizon problem in an arbitrary time step. The paper is the first to carefully analyze the quality of the value function approximation using two different types of basis functions and systematically varying the dimension of the basis. We present the closed form, continuous time approximation to the exogenous (discretely and inductively defined) processes in DICE, and we present a numerically more efficient re-normalized Bellman equation that, in addition, can disentangle risk attitude from the propensity to smooth consumption over time.

A 4-Stated Dice: Quantitatively Addressing Uncertainty Effects in Climate Change

We introduce a version of the DICE-2007 model designed for uncertainty analysis. DICE is a wide-spread deterministic integrated assessment model of climate change. Climate change, long-term economic development, and their interactions are highly uncertain. The quantitative analysis of optimal mitigation policy under uncertainty requires a recursive dynamic programming implementation of integrated assessment models. Such implementations are subject to the curse of dimensionality. Every increase in the dimension of the state space is paid for by a combination of (exponentially) increasing processor time, lower quality of the value or policy function approximations, and reductions of the uncertainty domain. The paper promotes a state reduced, recursive dynamic programming implementation of the DICE-2007 model. We achieve the reduction by simplifying the carbon cycle and the temperature delay equations. We compare our model's performance and that of the DICE model to the scientific AOGCM models emulated by MAGICC 6.0 and find that our simplified model performs equally well as the original DICE model. Our implementation solves the infinite planning horizon problem in an arbitrary time step. The paper is the first to carefully analyze the quality of the value function approximation using two different types of basis functions and systematically varying the dimension of the basis. We present the closed form, continuous time approximation to the exogenous (discretely and inductively defined) processes in DICE, and we present a numerically more efficient re-normalized Bellman equation that, in addition, can disentangle risk attitude from the propensity to smooth consumption over time.

An inexact interior point method for the large-scale simulation of granular material

Non-smooth contact dynamics provides an increasingly popular simulation framework for granular material. In contrast to classical discrete element methods, this approach is stable for arbitrary time steps and produces visually acceptable results in very short computing time. Yet when it comes to the prediction of draft forces, non-smooth contact dynamics is typically not accurate enough. We therefore propose to combine the method class with an interior point algorithm for higher accuracy. Our specific algorithm is based on the so-called Jordan algebras and exploits the relation to symmetric cones in order to tackle the conical constraints that are intrinsic to frictional contact problems. In every interior point iteration a linear system has to be solved. We analyze how the interior point method behaves when it is combined with Krylov subspace solvers and incomplete factorizations. We show that efficient preconditioners and efficient linear solvers are essential for the method to be applicable to large-scale problems. Using CG as a linear solver and incomplete LU factorizations, we substantially improve the accuracy in comparison to the projected Gauß–Jacobi solver.

Unconditionally Convergent Nonlinear Solver for Multiphase Flow in Porous Media under Viscous Force, Buoyancy, and Capillarity

Numerical simulations of multiphase flow in porous media often face convergence difficulties in the nonlinear Newton solver, including erratic time stepping, large number of (Newton) iterations, and timestep cuts. Such convergence problems can lead to unacceptably large computational time and are often the main impediment to performing simulation studies of large scale problems, such as oil/gas recovery, groundwater remediation, and CO2 geological sequestration. We analyze the nonlinearity of the discrete transport (mass conservation) equation for immiscible, two-phase flow in porous media in the presence of viscous, buoyancy, and capillary forces. The critical features that cause oscillations and divergence of the Newton iterations are identified and located. Based on the analysis, we develop a nonlinear solver that guides Newton iterations safely and efficiently, such that convergence is achieved for arbitrary timestep sizes.

Generalized Picard iterations: A class of iterated Runge–Kutta methods for stiff problems

A family of iterated Runge–Kutta (RK) methods capable of handling both non-stiff and stiff problems with complex spectrum is developed. The methods are based on classical implicit RK methods, are matrix-free and have natural parallelization capability. They are equivalent to first-order explicit RK schemes with unlimited number of stages. It is shown that they converge to the solution which corresponds to the underlying implicit RK method for all asymptotically stable linear problems y′=Jy with arbitrary big time steps τ. The convergence rate is limited by the ‘stiffness ratio’ ρ(J)ρ(J−1), rather than τρ(J). The experiments show that the proposed methods outperform explicit Dormand–Prince DOP853 code in stiff case and can compete with it in non-stiff case.

Assessing the combined benefits of water recycling technologies by modelling the total urban water cycle

This study investigates the potential benefits of new technologies, modern appliances, and innovative techniques that help to improve the performance of the urban water cycle. Urbanisation is a major source of additional pressures (both qualitative and quantitative) on the environment. For example abstractions to cover the increased demands for water supply or alterations of the topographic and geomorphologic properties of the land cover result in considerable changes to the dynamics of the hydrosystem (change of average and maximum values of flows). Sustainable, water-aware technologies, like SUstainable Drainage Systems (SUDS) and rainwater harvesting schemes, can be implemented to reduce these adverse effects. These technologies introduce interactions between the components of the urban water cycle. Rainwater harvesting for example, apart from the potable water demand reduction, may have a significant influence on the generated runoff. Consequently, an integrated modelling of the urban water cycle is necessary for the simulation of the water-aware technologies and the identification of their combined benefits. In this study, two hypothetical developments implement rainwater harvesting schemes and SUDS and are simulated using the Urban Water Optioneering Tool (UWOT), which is capable of using rainfall time series of arbitrary time steps. The two hypothetical developments were studied to investigate the contribution of the water-aware technologies to the minimisation of the environmental pressures. Significantly different urban density was assigned to these developments to highlight the influence of urban density on the efficiency and reliability of the water-aware technologies. The results indicate that: (a) water-saving schemes like rainwater harvesting and greywater treatment can reduce significantly the pressures of new developments (e.g., reduction of potable water demand by 27%); (b) the reliability of the water-aware technologies decreases with urban density; and (c) if localised rainwater harvesting is implemented then the efficiency of the water appliances influences considerably the generated runoff.

Application of the constrained interpolation profile method to room acoustic problems: Examination of boundary modeling and spatial/time discretization

The constrained interpolation profile (CIP) method is proposed in the field of computational fluid dynamics, because this method is less influenced by numerical dispersions and is stable with an arbitrary time step. In this paper, the CIP method is applied to room acoustic problems, in which boundary modeling and spatial/time discretization are specifically examined. First, frequency-independent and frequency-dependent boundary conditions are introduced and validated with an acoustic tube problem. Second, in one-dimensional and three-dimensional closed sound fields, this paper shows that errors in waveforms and reverberation times using the CIP method are smaller than those using the finite difference time domain method even if the steps for time and space are larger. Third, to deal with complicated shapes, the cubic interpolation with volume/area coordinates (CIVA) method using arbitrary triangular elements is introduced; in two-dimensional sound fields, calculation of the CIVA method is stable even if the quality of triangular elements is low.

Disordered quantum walks in one lattice dimension

We study a spin-$\frac{1}{2}$12-particle moving on a one-dimensional lattice subject to disorder induced by a random, space-dependent quantum coin. The discrete time evolution is given by a family of random unitary quantum walk operators, where the shift operation is assumed to be deterministic. Each coin is an independent identically distributed random variable with values in the group of two-dimensional unitary matrices. We derive sufficient conditions on the probability distribution of the coins such that the system exhibits dynamical localization. Put differently, the tunneling probability between two lattice sites decays rapidly for almost all choices of random coins and after arbitrary many time steps with increasing distance. Our findings imply that this effect takes place if the coin is chosen at random from the Haar measure, or some measure continuous with respect to it, but also for a class of discrete probability measures which support consists of two coins, one of them being the Hadamard coin.

Selection of time step for pseudodynamic testing

Although the step degree of nonlinearity has been introduced to conduct basic analysis and error propagation analysis for the pseudodynamic testing of nonlinear systems, it cannot be reliably used to select an appropriate time step before performing a pseudodynamic test. Therefore, a novel parameter of instantaneous degree of nonlinearity is introduced to monitor the stiffness change at the end of a time step, and can thus be used to evaluate numerical and error propagation properties for nonlinear systems. Based on these properties, it is possible to select an appropriate time step to conduct a pseudodynamic test in advance. This possibility is very important in pseudodynamic testing, since the use of an arbitrary time step might lead to unreliable results or even destroy the test specimen. In this paper, guidelines are proposed for choosing an appropriate time step for accurate integration of nonlinear systems. These guidelines require estimation of the maximum instantaneous degree of nonlinearity and the solution of the initial natural frequency. The Newmark explicit method is chosen for this study. All the analytical results and the guidelines proposed are thoroughly confirmed with numerical examples.

Gradient stability of numerical algorithms in local nonequilibrium problems of critical dynamics

The critical dynamics of a spatially inhomogeneous system are analyzed with allowance for local nonequilibrium, which leads to a singular perturbation in the equations due to the appearance of a second time derivative. An extension is derived for the Eyre theorem, which holds for classical critical dynamics described by first-order equations in time and based on the local equilibrium hypothesis. It is shown that gradient-stable numerical algorithms can also be constructed for second-order equations in time by applying the decomposition of the free energy into expansive and contractive parts, which was suggested by Eyre for classical equations. These gradient-stable algorithms yield a monotonically nondecreasing free energy in simulations with an arbitrary time step. It is shown that the gradient stability conditions for the modified and classical equations of critical dynamics coincide in the case of a certain time approximation of the inertial dynamics relations introduced for describing local nonequilibrium. Model problems illustrating the extended Eyre theorem for critical dynamics problems are considered.

An energy- and charge-conserving, implicit, electrostatic particle-in-cell algorithm

This paper discusses a novel fully implicit formulation for a one-dimensional electrostatic particle-in-cell (PIC) plasma simulation approach. Unlike earlier implicit electrostatic PIC approaches (which are based on a linearized Vlasov–Poisson formulation), ours is based on a nonlinearly converged Vlasov–Ampére (VA) model. By iterating particles and fields to a tight nonlinear convergence tolerance, the approach features superior stability and accuracy properties, avoiding most of the accuracy pitfalls in earlier implicit PIC implementations. In particular, the formulation is stable against temporal (Courant–Friedrichs–Lewy) and spatial (aliasing) instabilities. It is charge- and energy-conserving to numerical round-off for arbitrary implicit time steps (unlike the earlier “energy-conserving” explicit PIC formulation, which only conserves energy in the limit of arbitrarily small time steps). While momentum is not exactly conserved, errors are kept small by an adaptive particle sub-stepping orbit integrator, which is instrumental to prevent particle tunneling (a deleterious effect for long-term accuracy). The VA model is orbit-averaged along particle orbits to enforce an energy conservation theorem with particle sub-stepping. As a result, very large time steps, constrained only by the dynamical time scale of interest, are possible without accuracy loss. Algorithmically, the approach features a Jacobian-free Newton–Krylov solver. A main development in this study is the nonlinear elimination of the new-time particle variables (positions and velocities). Such nonlinear elimination, which we term particle enslavement, results in a nonlinear formulation with memory requirements comparable to those of a fluid computation, and affords us substantial freedom in regards to the particle orbit integrator. Numerical examples are presented that demonstrate the advertised properties of the scheme. In particular, long-time ion acoustic wave simulations show that numerical accuracy does not degrade even with very large implicit time steps, and that significant CPU gains are possible.

Unconditionally gradient-stable computational schemes in problems of fast phase transitions

Equations of fast phase transitions, in which the phase boundaries move with velocities comparable with the atomic diffusion speed or with the speed of local structural relaxation, are analyzed. These equations have a singular perturbation due to the second derivative of the order parameter with respect to time, which appears due to phenomenologically introduced local nonequilibrium. To develop unconditionally stable computational schemes, the Eyre theorem [D. J. Eyre, unpublished] proved for the classical equations, based on hypotheses of local equilibrium, is used. An extension of the Eyre theorem for the case of equations for fast phase transitions is given. It is shown that the expansion of the free energy on contractive and expansive parts, suggested by Eyre for the classical equations of Cahn-Hilliard and Allen-Cahn, is also true for the equations of fast phase transitions. Grid approximations of these equations lead to gradient-stable algorithms with an arbitrary time step for numerical modeling, ensuring monotonic nonincrease of the free energy. Special examples demonstrating the extended Eyre theorem for fast phase transitions are considered.

Probabilistic initial value problem for cellular automaton rule 172

We present a method of solving of the probabilistic initial value problem for cellular automata (CA) using CA rule 172 as an example. For a disordered initial condition on an infinite lattice, we derive exact expressions for the density of ones at arbitrary time step. In order to do this, we analyze topological structure of preimage trees of finite strings of length 3. Level sets of these trees can be enumerated directly using classical combinatorial methods, yielding expressions for the number of $n$-step preimages of all strings of length 3, and, subsequently, probabilities of occurrence of these strings in a configuration obtained from the initial one after $n$ iterations of rule 172. The density of ones can be expressed in terms of Fibonacci numbers, while expressions for probabilities of other strings involve Lucas numbers. Applicability of this method to other CA rules is briefly discussed.

A multigrid method for the Cahn–Hilliard equation with obstacle potential

We present a multigrid finite element method for the deep quench obstacle Cahn–Hilliard equation. The non-smooth nature of this highly nonlinear fourth order partial differential equation make this problem particularly challenging. The method exhibits mesh-independent convergence properties in practice for arbitrary time step sizes. In addition, numerical evidence shows that this behaviour extends to small values of the interfacial parameter γ . Several numerical examples are given, including comparisons with existing alternative solution methods for the Cahn–Hilliard equation.

An advanced modelling tool for simulating complex river systems

The present paper describes MOHID River Network (MRN), a 1D hydrodynamic model for river networks as part of MOHID Water Modelling System, which is a modular system for the simulation of water bodies (hydrodynamics and water constituents). MRN is capable of simulating water quality in the aquatic and benthic phase and its development was especially focused on the reproduction of processes occurring in temporary river networks (flush events, pools formation, and transmission losses). Further, unlike many other models, it allows the quantification of settled materials at the channel bed also over periods when the river falls dry. These features are very important to secure mass conservation in highly varying flows of temporary rivers. The water quality models existing in MOHID are base on well-known ecological models, such as WASP and ERSEM, the latter allowing explicit parameterization of C, N, P, Si, and O cycles. MRN can be coupled to the basin model, MOHID Land, with computes runoff and porous media transport, allowing for the dynamic exchange of water and materials between the river and surroundings, or it can be used as a standalone model, receiving discharges at any specified nodes (ASCII files of time series with arbitrary time step). These features account for spatial gradients in precipitation which can be significant in Mediterranean-like basins. An interface has been already developed for SWAT basin model.

Two-way coupling of fluids to rigid and deformable solids and shells

We propose a novel solid/fluid coupling method that treats the coupled system in a fully implicit manner making it stable for arbitrary time steps, large density ratios, etc. In contrast to previous work in computer graphics, we derive our method using a simple back-of-the-envelope approach which lumps the solid and fluid momenta together, and which we show exactly conserves the momentum of the coupled system. Notably, our method uses the standard Cartesian fluid discretization and does not require (moving) conforming tetrahedral meshes or ALE frameworks. Furthermore, we use a standard Lagrangian framework for the solid, thus supporting arbitrary solid constitutive models, both implicit and explicit time integration, etc. The method is quite general, working for smoke, water, and multiphase fluids as well as both rigid and deformable solids, and both volumes and thin shells. Rigid shells and cloth are handled automatically without special treatment, and we support fully one-sided discretizations without leaking. Our equations are fully symmetric, allowing for the use of fast solvers, which is a natural result of properly conserving momentum. Finally, for simple explicit time integration of rigid bodies, we show that our equations reduce to form similar to previous work via a single block Gaussian elimination operation, but that this approach scales poorly, i.e. as though four spatial dimensions rather than three.

Simulating system dynamics with arbitrary time step

We suggest a dynamic simulation method that allows efficient and realistic modeling of kinetic processes, such as atomic diffusion, in which time has its actual meaning. Our method is similar in spirit to widely used kinetic Monte Carlo (KMC) techniques; however, in our approach, the time step can be chosen arbitrarily. This has an advantage in some cases, e.g., when the transition rates change with time sufficiently fast over the period of the KMC time step (e.g., due to time dependence of some external factors influencing kinetics such as moving scanning probe microscopy tip or external time-dependent field) or when the clock time is set by some external conditions, and it is convenient to use equal time steps instead of the random choice of the KMC algorithm in order to build up probability distribution functions. We show that an arbitrary choice of the time step can be afforded by building up the complete list of events including the ``residence site'' and multihop transitions. The idea of the method is illustrated in a simple ``toy'' model of a finite one-dimensional lattice of potential wells with unequal jump rates to either side, which can be studied analytically. We show that for any choice of the time step, our general kinetics method reproduces exactly the solution of the corresponding master equations for any choice of the time steps. The final kinetics also matches the standard KMC, and this allows better understanding of this algorithm, in which the time step is chosen in a certain way and the system always advances by a single hop.

Neimark–Sacker bifurcations in a non-standard numerical scheme for a class of positivity-preserving ODEs

We discuss the nature of Neimark–Sacker bifurcations occurring in a non-standard numerical scheme, for a class of positivity-preserving systems of ordinary differential equations (ODEs) which undergoes a corresponding Hopf bifurcation. Extending previous work (Alexander & Moghadas 2005 a Electron. J. Diff. Eqn. Conf . 12 , 9–19), it is shown that the type of Neimark–Sacker bifurcation (supercritical or subcritical) may be affected by the integration time-step . The general form of the scheme in the vicinity of a fixed point is given, from which the expression for the first Lyapunov coefficient for two-dimensional systems, valid for arbitrary time-step, is explicitly derived. The analysis shows that this coefficient undergoes an shift with respect to the corresponding coefficient of the original ODE. This could lead to a type of bifurcation which differs from the corresponding Hopf bifurcation in the ODE, due to changes in the sign of the first Lyapunov coefficient as varies. This is especially problematic in the vicinity of certain types of degenerate Hopf bifurcation, at which this coefficient may vanish. We also present a general method to eliminate the possible shift in the bifurcation parameter of the scheme; however, the first Lyapunov coefficient may still be subjected to an shift, leading to a possibly erroneous type of bifurcation. Examples are given to illustrate the theoretical results of the paper with applications to mathematical biology.

Implicit solution of the four-field extended-magnetohydrodynamic equations using high-order high-continuity finite elements

Here we describe a technique for solving the four-field extended-magnetohydrodynamic (MHD) equations in two dimensions. The introduction of triangular high-order finite elements with continuous first derivatives (C1 continuity) leads to a compact representation compatible with direct inversion of the associated sparse matrices. The split semi-implicit method is introduced and used to integrate the equations in time, yielding unconditional stability for arbitrary time step. The method is applied to the cylindrical tilt mode problem with the result that a nonzero value of the collisionless ion skin depth will increase the growth rate of that mode. The effect of this parameter on the reconnection rate and geometry of a Harris equilibrium and on the Taylor reconnection problem is also demonstrated. This method forms the basis for a generalization to a full extended-MHD description of the plasma with six, eight, or more scalar fields.

Ionospheric index of solar activity based on the data of measurements of the spacecraft signals characteristics

The method for estimating an ionospheric index of solar-activity (IISA) based on the processing of spacecraft radio signals is suggested. The IISA values have been obtained by comparison between the measured and calculated variations of radio-signal characteristics. To calculate the variations of radio-signal characteristics, the straight rays approximation and the solar-activity index (Wolf numbers W and/or values of F 10.7 solar flux) as a control parameter of the ionospheric model have been used. The suggested method was tested using spacecraft radio signals from the radio-navigation system “CIKADA”. The reduced differences of phases (Δ Φ) for frequencies 150 and 400 MHz were measured and the same characteristics were calculated by integration along the ray of radio-wave propagation between the receiver and the satellite. The IRI-95 has been used as a background ionospheric model. The satellite co-ordinates were determined using the orbit parameters recorded in the navigation messages. Minimization of the difference measured and calculated Δ Φ using arbitrary time steps, or during whole time intervals of observation, gives the IISA corresponding the satellite pass. Daily IISA values were obtained by averaging over all communication contacts during a day (20–30 passes). Testing this approach based on the measurements during March/April 1997, 1998, shows that on magnetically quiet days differences between IISA and the primary solar activity indices are about 5%.