Low-dimensional Ordinary Differential Equations Research Articles

The art of solving a large number of non-stiff, low-dimensional ordinary differential equation systems on GPUs and CPUs

This paper discusses the main performance barriers for solving a large number of independent ordinary differential equation systems on processors (CPU) and graphics cards (GPU). With a naïve approach, for instance, the utilisation of a CPU can be as low as 4% of its theoretical peak processing power. The main barriers identified by the detailed analysing of the hardware architectures and profiling using hardware performance monitoring units are as follows. First, exploitation of the SIMD capabilities of the CPU via vector registers. The solution is to implement/enforce explicit vectorisation. Second, hiding instruction latencies on both CPUs and GPUs that can be achieved with increasing (instruction-level) parallelism. Third, the efficient handling of large timescale differences or event handling using the massively parallel architecture of GPUs. A viable option to overcome this difficulty is asynchronous time stepping. The above optimisation techniques and their implementation possibilities are discussed and tested on three program packages: MPGOS written in C++ and specialised only for GPUs; ODEINT implemented in C++, which supports execution on both CPUs and GPUs; finally, DifferentialEquations.jl written in Julia that also supports execution on both CPUs and GPUs. The tested systems (Lorenz equation, Keller–Miksis equation and a pressure relief valve model) are non-stiff and have low dimension. Thus, the performance of the codes are not limited by memory bandwidth, and Runge–Kutta type solvers are efficient and suitable choices. The employed hardware are an Intel Core i7-4820K CPU with 30.4 GFLOPS peak double-precision performance per cores and an Nvidia GeForce Titan Black GPU that has a total of 1707 GFLOPS peak double-precision performance.

Linear Quadratic Mean Field Games: Asymptotic Solvability and Relation to the Fixed Point Approach

Mean field game theory has been developed largely following two routes. One of them, called the direct approach, starts by solving a large-scale game and next derives a set of limiting equations as the population size tends to infinity. The second route is to apply mean field approximations and formalize a fixed point problem by analyzing the best response of a representative player. This paper addresses the connection and difference of the two approaches in a linear quadratic (LQ) setting. We first introduce an asymptotic solvability notion for the direct approach, which means for all sufficiently large population sizes, the corresponding game has a set of feedback Nash strategies in addition to a mild regularity requirement. We provide a necessary and sufficient condition for asymptotic solvability and show that in this case the solution converges to a mean field limit. This is accomplished by developing a re-scaling method to derive a low-dimensional ordinary differential equation (ODE) system, where a non-symmetric Riccati ODE has a central role. We next compare with the fixed point approach which determines a two-point boundary value (TPBV) problem, and show that asymptotic solvability implies feasibility of the fixed point approach, but the converse is not true. We further address non-uniqueness in the fixed point approach and examine the long time behavior of the non-symmetric Riccati ODE in the asymptotic solvability problem.

Independence screening for high dimensional nonlinear additive ODE models with applications to dynamic gene regulatory networks.

Mechanism-driven low-dimensional ordinary differential equation (ODE) models are often used to model viral dynamics at cellular levels and epidemics of infectious diseases. However, low-dimensional mechanism-based ODE models are limited for modeling infectious diseases at molecular levels such as transcriptomic or proteomic levels, which is critical to understand pathogenesis of diseases. Although linear ODE models have been proposed for gene regulatory networks (GRNs), nonlinear regulations are common in GRNs. The reconstruction of large-scale nonlinear networks from time-course gene expression data remains an unresolved issue. Here, we use high-dimensional nonlinear additive ODEs to model GRNs and propose a 4-step procedure to efficiently perform variable selection for nonlinear ODEs. To tackle the challenge of high dimensionality, we couple the 2-stage smoothing-based estimation method for ODEs and a nonlinear independence screening method to perform variable selection for the nonlinear ODE models. We have shown that our method possesses the sure screening property and it can handle problems with non-polynomialdimensionality. Numerical performance of the proposed method is illustrated with simulated data and a real data example for identifying the dynamic GRN of Saccharomyces cerevisiae.

Low-dimensional spike rate models derived from networks of adaptive integrate-and-fire neurons: Comparison and implementation

The spiking activity of single neurons can be well described by a nonlinear integrate-and-fire model that includes somatic adaptation. When exposed to fluctuating inputs sparsely coupled populations of these model neurons exhibit stochastic collective dynamics that can be effectively characterized using the Fokker-Planck equation. This approach, however, leads to a model with an infinite-dimensional state space and non-standard boundary conditions. Here we derive from that description four simple models for the spike rate dynamics in terms of low-dimensional ordinary differential equations using two different reduction techniques: one uses the spectral decomposition of the Fokker-Planck operator, the other is based on a cascade of two linear filters and a nonlinearity, which are determined from the Fokker-Planck equation and semi-analytically approximated. We evaluate the reduced models for a wide range of biologically plausible input statistics and find that both approximation approaches lead to spike rate models that accurately reproduce the spiking behavior of the underlying adaptive integrate-and-fire population. Particularly the cascade-based models are overall most accurate and robust, especially in the sensitive region of rapidly changing input. For the mean-driven regime, when input fluctuations are not too strong and fast, however, the best performing model is based on the spectral decomposition. The low-dimensional models also well reproduce stable oscillatory spike rate dynamics that are generated either by recurrent synaptic excitation and neuronal adaptation or through delayed inhibitory synaptic feedback. The computational demands of the reduced models are very low but the implementation complexity differs between the different model variants. Therefore we have made available implementations that allow to numerically integrate the low-dimensional spike rate models as well as the Fokker-Planck partial differential equation in efficient ways for arbitrary model parametrizations as open source software. The derived spike rate descriptions retain a direct link to the properties of single neurons, allow for convenient mathematical analyses of network states, and are well suited for application in neural mass/mean-field based brain network models.

Low-dimensional spike rate dynamics of coupled adaptive model neurons

The spiking activity of single neurons can be well described by a two-dimensional integrate-and-fire model that includes neuronal adaptation [1] caused by slowly decaying potassium currents [2]. For fluctuating inputs sparsely coupled spiking model neurons exhibit stochastic population dynamics which can be effectively characterized using the Fokker-Planck equation (see, e.g., [3-5]). This approach leads to a model with an infinite-dimensional state space and non-standard boundary conditions. However, the spike rate dynamics can be approximated by a low-dimensional ordinary differential equation in different ways [4,6,7]. Although these approximation techniques are interrelated it is not clear which reduced model best reproduces the spike rate of the original spiking network, depending on the statistics of the input. Here we first extend each of these reduction methods to account for neuronal adaptation and then evaluate the reduced models in terms of spike rate reproduction accuracy for a range of biologically plausible input statistics, computational demand and implementation complexity (see, e.g., Figure ​Figure1).1). These reduced descriptions are well suited for (i) application in neural mass/mean-field based brain network models, having a link to single neuron properties retained and being computationally efficient, and (ii) mathematical analyses of, e.g., asynchronous and rhythmic network states. Figure 1 Simulation of a large population of adaptive exponential integrate-and-fire (aEIF) neurons driven by a stochastic current with time-varying moments. Instantaneous spike rate and adaptation current averaged over 200,000 neurons are shown in gray. Overlaid ...

APOD‐based control of linear distributed parameter systems under sensor/controller communication bandwidth limitations

The synthesis of a model‐based control structure for general linear dissipative distributed parameter systems (DPSs) is explored in this manuscript. Discrete‐time distributed state measurements (called process snapshots) are used by a continuous‐time regulator to stabilize the process. The main objective of this article is to identify a criterion to minimize the communication bandwidth between sensors and controller (snapshots acquisition frequency) using linear systems analysis and still achieve closed‐loop stability. This objective is addressed by adding a modeling layer to the regulator. Theoretically, DPSs can be well described by low dimensional ordinary differential equation models when represented in functional spaces; practically, the model accuracy hinges on finding basis functions for these spaces. Adaptive proper orthogonal decomposition is used to identify statistically important basis functions and establish locally accurate reduced order models which are then used in controller design. The proposed approach is successfully applied toward thermal regulation in a tubular chemical reactor. © 2014 American Institute of Chemical Engineers AIChE J, 61: 434–447, 2015

Cocirculation of infectious diseases on networks.

We consider multiple diseases spreading in a static configuration model network. We make standard assumptions that infection transmits from neighbor to neighbor at a disease-specific rate and infected individuals recover at a disease-specific rate. Infection by one disease confers immediate and permanent immunity to infection by any disease. Under these assumptions, we find a simple, low-dimensional ordinary differential equations model which captures the global dynamics of the infection. The dynamics depend strongly on initial conditions. Although we motivate this Rapid Communication with infectious disease, the model may be adapted to the spread of other infectious agents such as competing political beliefs, or adoption of new technologies if these are influenced by contacts. As an example, we demonstrate how to model an infectious disease which can be prevented by a behavior change.

Low-Dimensional Modelling of Turbulence Using the Proper Orthogonal Decomposition: A Tutorial

The proper orthogonal decomposition identifies basis functions or modes which optimally capture the average energy content from numerical or experimental data. By projecting the Navier–Stokes equations onto these modes and truncating, one can obtain low-dimensional ordinary differential equation models for fluid flows. In this paper we present a tutorial on the construction of such models. In addition to providing a general overview of the procedure, we describe two different ways to numerically calculate the modes, show how symmetry considerations can be exploited to simplify and understand them, comment on how parameter variations are captured naturally in such models, and describe a generalization of the procedure involving projection onto uncoupled modes that allow streamwise and cross-stream components to evolve independently. We illustrate for the example of plane Couette flow in a minimal flow unit – a domain whose spanwise and streamwise extent is just sufficient to maintain turbulence.

Parametrically forced envelope Rossby solitons: Soliton-oscillation

In this paper, the non-dimensional equivalent barotropic vorticity equation with dissipation and an external vorticity source induced by diabatic heating was reduced to the generalized periodically forced nonlinear Schrödinger (NLS) equation similar to the parametrically excited NLS equation derived by Miles (1984) by means of the multiplescale method. This forced NLS equation can be reduced to a low-dimensional ordinary differential equation describing the amplitude and phase of the forced envelope Rossby soliton by the inverse-scattering perturbation method. The numerical calculations indicate that for weaker dissipation the envelope Rossby soliton, under a constant amplitude travelling diabatic heating, exhibits multi-periodic oscillations, and its phase portrait is found to be a limit cycle. However, for the forcing of one-year modulated diabatic heating, the forced envelope Rossby soliton can execute chaotic motion when the controlling parameter reaches a certain extent. In addition, the oscillatory behaviour of the forced envelope Rossby soliton is found to be able to explain the establishment and decay of vortex pair blocking usually observed in two oceans.

FEIGENBAUM’S UNIVERSALITY IN A LOW-DIMENSIONAL FLUID MODEL

We present a low-dimensional truncated model for a viscous fluid contained in a two-dimensional square box, obtained by truncating a dynamical system of amplitudes for the velocity field. This low-dimensional model exhibits a route to chaos via a period doubling cascade (Feigenbaum’s Scenario). In order to compute with high accuracy the period doubling, a numerical method based on the first order variational equations and a Poincaré map has been developed. This methodology can also be applied to the analysis of bifurcations of periodic orbits in low-dimensional ordinary differential equations. This method allows to detect not only the presence of bifurcations but also the computation of stable and unstable periodic orbits. On the other hand, the chaotic dynamics of the system is analyzed in detail by the computation of the Liapunov exponents for long-time integrations. For this purpose, a numerical scheme based on renormalization techniques has been constructed.