Coarse-Grained Descriptions of Dynamics for Networks with both Intrinsic and Structural Heterogeneities

TL;DR

This paper introduces a systematic approach to identify coarse-grained variables that effectively summarize the dynamics of large, heterogeneous networks, enabling reduced models for analysis and simulation.

Contribution

It proposes a novel method linking node descriptors to collective network behavior, incorporating heterogeneity quantification for complex coupled systems.

Findings

Successfully identified collective variables for heterogeneous oscillator networks

Demonstrated the approach's applicability to networks with intrinsic and structural heterogeneities

Discussed computational implementation and potential extensions

Abstract

Finding accurate reduced descriptions for large, complex, dynamically evolving networks is a crucial enabler to their simulation, analysis, and, ultimately, design. Here we propose and illustrate a systematic and powerful approach to obtaining good collective coarse-grained observables-- variables successfully summarizing the detailed state of such networks. Finding such variables can naturally lead to successful reduced dynamic models for the networks. The main premise enabling our approach is the assumption that the behavior of a node in the network depends (after a short initial transient) on the node identity: a set of descriptors that quantify the node properties, whether intrinsic (e.g. parameters in the node evolution equations) or structural (imparted to the node by its connectivity in the particular network structure). The approach creates a natural link with modeling and "computational enabling technology" developed in the context of Uncertainty Quantification. In our case, however, we will not focus on ensembles of different realizations of a problem, each with parameters randomly selected from a distribution. We will instead study many coupled heterogeneous units, each characterized by randomly assigned (heterogeneous) parameter value(s). One could then coin the term Heterogeneity Quantification for this approach, which we illustrate through a model dynamic network consisting of coupled oscillators with one intrinsic heterogeneity (oscillator individual frequency) and one structural heterogeneity (oscillator degree in the undirected network). The computational implementation of the approach, its shortcomings and possible extensions are also discussed.

Figures

Tables

Table 1. Table 1 : Frequently encountered probability distributions and the corresponding weighted orthogonal polynomial families.

Discrete		Continuous
Distribution	Polynomials	Distribution	Polynomials
Binomial	Kravchuk	Gaussian	Hermite
Poisson	Charlier	Gamma	Laguerre
Negative binomial	Meixner	Beta	Jacobi
Hypergeometric	Hahn	Uniform	Legendre

Equations

$$\frac{\mathrm{d}{\phi }_i(t)}{\mathrm{d}t}={\hat{\omega }}_i+\frac{K}{N}\sum _{j=1}^N{A}_{i,j}\sin ({\phi }_j(t)-{\phi }_i(t)),$$

$${\theta }_i(t)\stackrel{\mathrm{def}}{=}\left({\phi }_i(t)-\frac{1}{N}\sum _{j=1}^N{\phi }_j(t)\right),i=1,2,\dots ,N-1,$$

$$\frac{\mathrm{d}{\theta }_i}{\mathrm{d}t}={\hat{\omega }}_i-\frac{1}{N}\left[\sum _{j=1}^N{\hat{\omega }}_j\right]+\frac{K}{N}\left[\sum _{j=1}^N{A}_{i,j}\sin ({\theta }_j-{\theta }_i)\right],$$

$$\begin{array}[]{rccl}\theta_{i}(t)&=&f(t;\bm{p}_{i})\approx&\sum_{k=1\ldots{M}}\alpha_{k}(t)\psi^{(k)}(\bm{p}_{i}).\end{array}$$

$$w_i=Np(1-q(i-1)/N{)}^r,i=1,2,\dots ,N$$

$$P_{ij}=P_{ji}=\min \left(\frac{w_i{w}_j}{{\sum }_k{w}_k},1\right).$$

$${\psi }^{(\mathbf{\gamma })}(\omega ,\kappa )={\psi }_{\omega }^{({\gamma }_{k,\omega })}(\omega ){\psi }_{\kappa }^{({\gamma }_{k,\kappa })}(\kappa ).$$

$$\begin{array}[]{rcl}\theta(t)&\approx&\sum_{k=1}^{M}\alpha_{k}(t)\psi^{(k)}(\omega,\kappa)\equiv\sum_{k=1}^{M}\alpha_{k}(t)\xi^{(\gamma_{k,\omega})}(\omega)\zeta^{(\gamma_{k,\kappa})}(\kappa)\\ G&=&\{\bm{\gamma}_{k}=(\gamma_{k,\omega},\gamma_{k,\kappa})\}:0\leq\gamma_{k,\omega},\gamma_{k,\kappa}\in\mathbb{Z},\quad\gamma_{k,\omega}+\gamma_{k,\kappa}\leq p_{\mathrm{max}}\},\\ {M}&=&||G||=(1+p_{\text{max}})(2+p_{\text{max}})/2\end{array}$$

$${\alpha }_k=\frac{{\int }_{D}f(\mathbf{x}){\psi }^{(k)}(\mathbf{x})\mathrm{d}\Gamma (\mathbf{x})}{{\int }_D({\psi }^{(k)}(\mathbf{x}){)}^2\mathrm{d}\Gamma (\mathbf{x})},$$

$$\begin{array}[]{rcl}\sigma&=&||f(\bm{x})-\sum_{{k}=1}^{{M}}a_{k}\psi^{(k)}(\bm{x})||_{2}^{2}\\ &\approx&\hat{\sigma}=\sum_{i=1}^{n_{\text{samp}}}w(\bm{x}_{i})\left(f(\bm{x}_{i})-\sum_{{k}=1}^{{M}}\alpha_{k}\psi^{(k)}(\bm{x}_{i})\right)^{2},\end{array}$$

$${\alpha }_k=\frac{{\sum }_{i=1}^{n_{\text{samp}}}w({\mathbf{x}}_i)f({\mathbf{x}}_i){\psi }^{(k)}({\mathbf{x}}_i)}{{\sum }_{i=1}^{n_{\text{samp}}}w({\mathbf{x}}_i)({\psi }^{(k)}({\mathbf{x}}_i){)}^2}.$$

$$\sum _{i=1}^{n_{\text{samp}}}w({\mathbf{x}}_i)g({\mathbf{x}}_i)$$

$${\int }_{D}g(\mathbf{x})\mathrm{d}\Gamma (\mathbf{x}),$$

$${\int }_{D}g(\mathbf{x})\rho (\mathbf{x})\mathrm{d}\mathbf{x}.$$

$$\underset{n_{\text{samp}}\to \infty }{\lim }\frac{1}{n_{\text{samp}}}\sum _{i=1}^{n_{\text{samp}}}g({\mathbf{x}}_i)={\int }_{D}g(\mathbf{x})\rho (\mathbf{x})\mathrm{d}\mathbf{x}.$$

$$\begin{array}[]{rcl}\Phi_{{\tau},C}&:&\mathbb{R}^{M}\rightarrow\mathbb{R}^{M}\\ \Phi_{{\tau},F}&:&\mathbb{R}^{{N}-1}\rightarrow\mathbb{R}^{{N}-1}\\ \Phi_{{\tau},F}[\bm{\theta}(t)]&=&\bm{\theta}(t+{\tau})\\ &=&\int_{s=t}^{s=t+{\tau}}\frac{\mathrm{d}\bm{\theta}(s)}{\mathrm{d}t}\mathrm{d}s\\ \Phi_{{\tau},C}[\bm{\alpha}(t)]&=&\bm{\alpha}(t+{\tau})\\ &\equiv&\left(R\circ\Phi_{{\tau},F}\circ L\right)[\bm{\alpha}(t)]\end{array}$$

$$\mathbf{\alpha }(t+h)=\mathbf{\alpha }(t)+h\frac{{\mathbf{\Phi }}_{\tau ,C}[\mathbf{\alpha }(t)]-\mathbf{\alpha }(t)}{\tau }$$

$${\mathbf{F}}_{\tau }[\mathbf{\alpha }(t)]={\mathbf{\Phi }}_{\tau ,C}[\mathbf{\alpha }(t))]-\mathbf{\alpha }(t).$$

$${\lambda }_i\approx {\hat{\lambda }}_i=\ln ({\mu }_i+1)/\tau$$

$$⟨{\psi }_{i,j},{\psi }_{m,n}{⟩}_{\rho \sigma }={\delta }_{im}{\delta }_{jn}$$

$$⟨f,g{⟩}_{\rho \sigma }=\iint f(x,y)g(x,y)\rho (x)\sigma (y)\mathrm{d}x\mathrm{d}y$$

$$⟨{\varphi }_i,{\varphi }_j{⟩}_{\rho }={\delta }_{ij}$$

$$⟨{\theta }_i,{\theta }_j{⟩}_{\sigma }={\delta }_{ij}.$$

$${\psi }_{i,j}(x,y)=\sum _{p\le i,q\le j}{A}_{ijpq}{\varphi }_p(x){\theta }_q(y)$$

$$\begin{array}[]{rcl}\delta_{im}\delta_{jn}&=&\sum_{p\leq i,q\leq j}\sum_{r\leq m,s\leq n}A_{ijpq}A_{mnrs}\langle\phi_{p}(x),\phi_{r}(x)\rangle_{\rho}\langle\theta_{q}(y),\theta_{s}(y)\rangle_{\sigma}\\ &=&\sum_{p\leq i,q\leq j}\sum_{r\leq m,s\leq n}A_{ijpq}A_{mnrs}\delta_{pr}\delta_{qs}\\ &=&\sum_{p\leq\min(i,m),q\leq\min(j,n)}A_{ijpq}A_{mnpq}\end{array}$$

$$A_{ijmn}={\delta }_{im}{\delta }_{jn}.$$

$$\begin{array}[]{rcl}&&A\xrightleftharpoons[b]{a}B\\ \frac{\mathrm{d}}{\mathrm{d}t}\bm{x}&=&{\mathbf{J}}\cdot\bm{x}+V\left[\begin{array}[]{c}x_{A,\text{in}}\\ x_{B,\text{in}}\end{array}\right],\end{array}$$

Full text

1 Keywords:

nonlinear dynamics, heterogeneity, instability, uncertainty, bifurcation, data mining, machine learning, networks, model reduction.

2 Introduction

Model reduction for dynamical systems has been an important research direction for decades; accurate reduced models are very useful, and often indispensable for the understanding, analysis, and ultimately for the design of large/complex dynamical systems. The relevant tools and techniques range from center manifold reduction close to bifurcation points [15] to singular perturbation techniques (analytical [3] or computational [22]) and (more recently) to data-driven reduction methods (like PCA [17], or nonlinear manifold learning techniques [7, 10]). While many such tools are well established for ODEs and PDEs, the dynamics of networked dynamical systems (e.g. [38]) pose additional challenges.

We are interested here in large sets of dynamic units (agents, oscillators, cells) linked in a prescribed (and, for this paper, fixed) coupling pattern. Every unit consists here of a (relatively small) set of ordinary differential equations. The units are intrinsically heterogeneous, meaning that the parameters of this set of ODEs are sampled from a probability distribution. Once the overall system of ODEs modeling a large network is assembled, any generic dynamical system model reduction technique can be tried. For all-to-all coupled heterogeneous units, in particular, there has been extensive work taking advantage of the overall model structure, leading to the systematic reduction of such intrinsically heterogeneous assemblies [29, 40, 41].

Our illustrative example is a simulation of coupled phase oscillators whose dynamics are governed by the equations

[TABLE]

where $i\in 1,\ldots,{N}$, and $A_{i,j}\in\{0,1\}$ is the adjacency matrix for a network with identical edges. This model was originally formulated by Yoshiki Kuramoto with all-to-all coupling (i.e., $A_{i,j}=1\forall i,j$) [27, 28]. While we work with the simplified Kuramoto oscillator system, the methods presented here have been shown to work for more realistic coupled-oscillator systems as well, such as Hodgkin-Huxley-like neurons [5], metabolizing cells [4], gene-expression oscillations in circadian rhythms (ongoing work), or other candidate systems [2]. In cases such as the Hodgkin-Huxley, where each unit is described by multiple dynamic variables, (e.g, membrane potential and gating variables), the analysis used here is repeated for each per-cell variable.

In order to construct a frame that moves with the average phase angle, we use states $\bm{\theta}\in\mathbb{R}^{{N}-1}$, where

[TABLE]

and $\theta_{N}(t)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}-\sum_{j=1}^{{N}-1}\theta_{j}(t)$. (Hereafter, explicit time-dependence of $\theta(t)$ and $\varphi(t)$ is usually not indicated, but can be assumed.) Since the vector $\bm{\hat{}}\omega$ of natural frequencies is not time-dependent, the transformation is particular to this problem and not really relevant to the model reduction technique discussed in this paper, though it ensures the existence of a steady $\bm{\theta}$ state for sufficiently high values of ${K}$. This transformation is used to generate a new dynamical system

[TABLE]

$\forall i\in[1,{N}-1]$.

The main idea involves mathematical “technology transfer” from the field of Uncertainty Quantification (UQ) [14, 52]. We assume that the long-time behavior of each unit in the assembly is characterized by (is a function of) its identity–the value(s) of the heterogeneous parameter(s). The problem can then be formulated as being distributed in (heterogeneous) parameter space ($\bm{p}$-space) in a manner analogous to spatiotemporal processes “distributed” over physical space. The state $\theta_{i}(t)$ for any unit $i$ with identity $\bm{p}_{i}$ can be approximated in terms of appropriate basis functions not in physical space, but rather in “identity” space: heterogeneous parameter space.

[TABLE]

The basis functions $\psi^{(k)}(\bm{p})$ are constructed as orthogonal polynomials in §4.2. Because we are modeling cases in which the behavior of each unit is assumed to be a smooth function of identity (that is, nodes with similar identities are expected to behave similarly) a relatively short truncation of such a series may well be accurate enough if the right basis functions in $\bm{p}$ space are chosen. In such a case, the number of ODEs to be solved reduces from the number of units $\mathcal{O}({N})$, to the number of terms in the series $\mathcal{O}({M})$. This approach, and its links to UQ modeling/computational developments (like the use of Smolyak grids for nonintrusive collocation-based simulation) has been explored in [31]. Yet these developments were only applicable for all-to-all coupled, intrinsically heterogeneous assemblies of units.

The purpose of this paper is to generalize this approach by introducing a simple, yet nontrivial, extension. Namely, we consider networks with non-trivial coupling structure, i.e. not all-to-all coupled. In this work all connections have the same strength; further extension to weighted connections is nontrivial, and the subject of current research. Each unit is now also characterized, beyond its ODE parameter values, by its connectivity in the network—the nature of its coupling with other units in the network, which in turn is quantified by features such as the unit’s degree (its count of undirected connections). Different nodes have different connectivity features (imposed by the network structure)–we can therefore think of connectivity as a type of heterogeneity of our building block units: structural heterogeneity rather than intrinsic heterogeneity.

It should be noted that the steady state of similar dynamical systems have been predicted analytically [16, 46]. However, [16] considers the $N\rightarrow\infty$ limit, whereas we deal with finite $N$. Additionally, [46] requires that all the phases be known exactly to use a self-consistency argument. We are trying to obtain a reduced description, so necessarily we do not know (nor want to know) all the phases of all oscillators, but rather to approximate them.

In the system Eq. (1) to which we apply our coarse-graining strategy, we will see that the only connectivity feature that appreciably affects unit dynamics is the unit degree $\hat{\kappa}_{i}=\left(\sum_{j=1}^{N}A_{ij}\right)\in[0,N]$. That is, though we have previously shown that intrinsically similar nodes in this system [35] (and others [36]) have similar dynamics with all-to-all coupling, we show here that, with a nontrivial coupling topology network, nodes which additionally have the same (similar) degrees also have similar dynamics (possibly after a short initial transient). The degree of a node can be treated as another heterogeneous node parameter, whose probability distribution is the network degree distribution. We argue that the same approach which uses the distribution of an intrinsic heterogeneity, and led to the reduction of all-to-all unit assemblies in [35] can be naturally extended to include a distribution over structural heterogeneity that leads to reduction of unit assemblies coupled in networks.

We demonstrate this in the simplest nontrivial representative setting we can put together: a set of coupled phase oscillators, characterized by heterogeneous frequencies $\hat{\omega}_{i}$ sampled from a prescribed distribution (here a truncated Gaussian)—but now not all-to-all coupled. Instead, the coupling ${\mathbf{A}}$ is in the form of a complex network, generated by a Chung-Lu process similar to that described in [30, 6]. The ideas presented here work also for general networks, if the degrees of the nodes are large enough; the Chung-Lu network is used as a convenient example. Likewise, the particular degree distribution tested here (shown in Fig. 1) is not itself important, and networks with other degree distributions, such as power-law, or even nonmonotonic degree distributions, can also be used. Here, we first generate a weight sequence $w_{i}$, as

[TABLE]

with parameters $p=0.50$, $q=0.90$, and $r=0.50$. With $\bm{w}$, we generate the connection probabilities ${\mathbf{P}}$ via

[TABLE]

We concretely generate an adjacency matrix ${\mathbf{A}}$ from these probabilities by inverse transform sampling only above the diagonal of ${\mathbf{A}}$, and copying to the bottom triangle, to ensure that the network is undirected, with no self-loops.

Instead of following the behavior of each individual oscillator, we exploit the observation that similar oscillators have similar behavior and can be tracked together. For all to all coupling, “similar oscillators” is taken to imply similar natural frequencies, and we write the oscillator state as a function of natural frequency [35] and time. However, when a non-trivially structured coupling exists, “similar oscillators” implies not only intrinsic similarity but also structural similarity. In addition to the intrinsic explanatory parameter of the natural frequencies, the degree of each node appears to be an explanatory parameter which suffices (in our simple model) to capture the influence of the coupling structure on the behavior of each oscillator; yet other features such as in-degree or local clustering coefficient are also worth considering. If two oscillators have similar $\hat{\kappa}$ values and similar $\hat{\omega}$ values, then their time-dependent behavior is observed here to be similar, possibly after a short transient. Finding the relationship between oscillator characteristics (intrinsic and structural) and oscillator states generates a coarse-grained description, whereby the system state can be encoded in fewer independent variables.

We illustrate a number of coarse-grained modeling tasks facilitated by this reduction: accelerated simulation (via coarse projective integration), accelerated fixed point computation, continuation, and coarse-grained stability analysis (via time-stepper based coarse Newton-Krylov GMRES [19, 20] and Arnoldi algorithms [47]). This creates a natural link between the reduction approach we present here, and our so-called Equation-Free framework for complex systems modeling [25, 26].

In the end, what makes it all possible is the fundamental assumption about how heterogeneity (intrinsic as well as structural) affects the solution: “nearby” parameter values and “nearby” connectivities imply “nearby” dynamics. This is not always the case for any network, and so testing that this assumption holds must be performed on a case-by-case basis. For our networks, such a check is demonstrated in Fig. 2. However, whether such a parameterization is possible is linked to the question of whether frequency-synchronization emerges, a subject for which an extensive literature exists [9].

The remainder of the paper is organized as follows. First, we describe our illustrative example, a network of heterogeneous phase oscillators. We then give the form of our low-dimensional representation of the system state. We use our ability to transform back and forth between the two state representations to perform several computational tasks in the coarse-grained space, including the solution of initial value problems, the computation of fixed points, their stabilities and bifurcations. Appendices include an analysis of the validity of using higher-order coarse-grained integration schemes.

3 An illustrative example of heterogeneous coupled oscillator networks

Our illustrative example is a network of coupled Kuramoto oscillators with heterogeneous natural frequencies $\hat{\omega}_{i}$, coupled in a stochastically generated network (here, a Chung-Lu network [6] with parameters $p=0.50$, $q=0.90$, and $r=0.50$, an example instance of which is shown in Fig. 1). This type of model system was used in some previous reduction studies [35, 43]. The number of oscillators in the network is also a parameter we will vary; our base case is $N=196$. A basic premise, which is corroborated by Fig. 6, is that the network is large enough (the number $N$ of nodes is large enough) for the single realization to be representative of the expectation over all consistent network realizations. The fine dynamics are governed by the system of coupled ordinary differential equations (ODEs) Eq. (1), where the natural frequencies $\hat{\omega}_{i}$ and the node degrees (numbers of neighbors) $\hat{\kappa}_{i}$ are heterogeneous across the oscillators constituting the network. We remind the reader of our assumption that, of all structural node features that may affect the dynamics, it will be the node degree that matters here–so that the identity of node $i$ is sufficiently described by its intrinsic parameter $\hat{\omega}_{i}$ and its structural parameter $\hat{\kappa}_{i}$. This assumption is supported first by Fig. 2, and later, as we will see, more quantitatively by Fig. 6.

We further define rescaled versions of the two heterogeneous parameters, $x_{i}=(\hat{x}_{i}-\text{mean}(\{\hat{x}_{j}\})/\text{stddev}(\{\hat{x}_{j}\})$ for $x=\omega,\kappa$ and $i,j\in 1,\ldots,{N}-1$. These two transformations do not affect the fine dynamics of Eq. (1) or Eq. (3); only the numerics of the implementation of the restriction $R$ to a coarse-grained state representation, to be developed below. Without axis markings, Fig. 2, for instance, would look the same whether $\omega\times\kappa$ or $\hat{\omega}\times\hat{\kappa}$ were used for plotting.

The emergent functional dependence of the $\theta_{i}$ on the $\omega_{i}$ (the intrinsic heterogeneity only) was discussed in the all-to-all coupling context in [34, 35]. There, we used a one-dimensional polynomial chaos expansion (PCE) to describe the reduced problem for $A_{ij}=1\,\forall i,j\in[1,{N}],\,i\neq j$, so that all nodes have degree $\hat{\kappa}_{i}={N}-1$. In this paper, we again expect the oscillator states to quickly become smooth (and time-dependent!) functions of their identities, but node $i$’s identity now includes both $\hat{\omega}_{i}$ and $\hat{\kappa}_{i}$. For ${K}$ sufficiently large so that a steady state of Eq. (1) exists, we indeed observe that the states of randomly initialized oscillators quickly approach an apparently smooth surface in $\omega\times\kappa$ space (see Fig. 2) suggesting that a low-order series truncation of the type described in Eq. (8) may constitute a good description. This motivates the use of a functional fit of the coefficients (the few $\alpha_{k}(t)$ in Eq. (8)) to the data (the many $\theta_{i}(t)$) as a coarse representation.

In previous work [43] we have used a projection onto the eigenvectors of the discrete Laplacian on the graph to describe the dependence of oscillator state on structural heterogeneity, while using a one-term/linear fit to account for dependence on intrinsic frequency.

4 Low-dimensional representation

4.1 Polynomial chaos

Given our observation that oscillator behavior quickly becomes a function of oscillator identity, we want to describe the long-term dynamics of the oscillator phase angles as a smooth function $\theta=\theta(t;\omega,\kappa)$. The phase angle of the $i$-th oscillator is then given by $\theta_{i}(t)\equiv\theta(t;\omega_{i},\kappa_{i})$. Since our two heterogeneities (the intrinsic and the structural) are here independent, the basis functions are a tensor product of two independent polynomial bases

[TABLE]

This is a special case; the formulation will still in principle be applicable for parameters with correlated joint probability distributions if one constructs an appropriate set of basis functions [37, 8].

We now express the network dynamics in the form of a series expansion in a truncation of this tensor product basis as

[TABLE]

where the $\alpha_{k}(t)$ are time-dependent coefficients, $\xi^{(\gamma_{k,\omega})}(\omega)$ are basis functions arising from the intrinsic heterogeneity dependence and $\zeta^{(\gamma_{k,\kappa})}(\kappa)$ are basis functions arising from the structural heterogeneity dependence. Within the truncation of the set of functions $G$ included in the basis, the ordering of the basis can be chosen arbitrarily, and so we substitute the vector index $\bm{\gamma}_{k}=(\gamma_{k,\omega},\gamma_{k,\kappa})$ with a scalar index $1\leq k\leq{M}$.

The analogy with UQ now manifests itself in our choice of the two independent basis sets: each one of them is chosen to be a polynomial chaos basis in the corresponding heterogeneous (in analogy to random) parameter. Each set of polynomials is orthogonal with respect to the probability density of the corresponding heterogeneous parameter, and the joint heterogeneity probability density is just the product of the two unidimensional, independent heterogeneity probability densities, so that Eq. (7) is satisfied, as shown in Appendix A.

Note that in Eq. (8), we specify that $\gamma_{k,\omega}+\gamma_{k,\kappa}\leq p_{\mathrm{max}}$. This allows us to say that our two-dimensional polynomials are of total degree $\leq p_{\mathrm{max}}$. An alternative truncation rule would be to require that $\gamma_{k,\omega}\leq p_{\mathrm{max}}$ and $\gamma_{k,\kappa}\leq p_{\mathrm{max}}$ (or even to place separate bounds on $\gamma_{k,\omega}$ and $\gamma_{k,\kappa}$, allowing for some anisotropy in the details, a topic for separate investigation). Both approaches can be found in the literature.

This allows, per Eq. (8), an approximation of the behavior as a time-dependent two-dimensional surface in one intrinsic dimension (here, the (normalized) natural frequencies $\omega$), and one structural dimension (here, the (normalized) node degrees $\kappa$). We repeat that this tensor product basis, limited to those polynomials of total order less than some desired maximum, is a truncated orthogonal basis for the 2D space weighted by probability densities that are products of two marginal distributions. Some standard distributions and their corresponding families of orthogonal polynomials [52] are given in Table 1.

In classical Galerkin methods an inner product is taken between the governing evolution equations and each basis function, producing ODEs for the dynamics of the expansion coefficients (the $\alpha_{k}(t)$ in Eq. (8)) by exploiting orthogonality. A similar approach could be taken here (through analytical computation of the inner product integrals if possible, else through numerical quadrature). Instead, we do not directly calculate the temporal rates-of-change of the expansion coefficients, but infer features of the coefficient dynamics from brief bursts of simulation of the dynamics of the full system (see §5). This equation-free approach relies on our ability to go back and forth between fine descriptions of the system state (the $\theta_{i}$ values), and coarse ones (the $\alpha_{k}$ values). This is analogous to a non-intrusive (black-box, input-output) approach to using polynomial chaos in UQ.

4.2 Equation-free numerics

The choice of the polynomial basis sets follows the selection of the appropriate heterogeneity distribution. For several frequently encountered distributions, the bases have been tabulated (e.g. Table 1) from original generalized polynomial chaos references (see e.g. [52]). If not already available in such tables, one can construct the basis polynomials e.g. through Gram-Schmidt orthogonalization with inner products in the space weighted by the distribution function. A couple of nontrivial considerations arising in our case are that (a) our structural heterogeneity parameter (the node degree) takes integer values, and so the degree distribution is discrete; and (b) often we may encounter problems for which the heterogeneity distribution is not explicitly known, but has to be estimated from specific system realizations (so, from large enough samples). For the case of explicitly unknown but sampled distributions (whether discrete or continuous) we have used here the moments of the sampling of the heterogeneity parameters for our particular network realization to extract the corresponding polynomials (using SVD-based pseudo-inverses) [39].

Results of the moment-based polynomial generation method, which we used to generate all 1D polynomials used in this paper, are shown in Fig. 3. Here, the marginal samplings of degrees and natural frequencies were used separately to generate 1D polynomials via moments, and then a 2D basis was defined from the tensor product of these 1D bases, with the restriction $\mathcal{O}(\omega)+\mathcal{O}(\kappa)\leq p_{\mathrm{max}}$ placed on the total polynomial degree of the 2D basis functions used.

For the expansion in Eq. (8), we used a set of 2D basis functions.

Once the relevant polynomials have been constructed, the appropriately defined inner product also allows us (whether for continuous or for discrete distributions, explicitly known or not) to find the coefficients $\alpha_{k}$ in Eq. (8) for a given observation $\theta(t)$ of the system states for a particular system realization (sampling of the distribution). This can be accomplished directly (via numerical approximations of the relevant inner products), using the orthogonality of the $\psi^{(k)}$,

[TABLE]

where we use the Lebesgue integral ${\int_{D}}{}g(\bm{x}){\mathrm{d}\Gamma(\bm{x})}=\mathbb{E}[g(\bm{x})]$. For our problem, with one continuous and one discrete variable, this can be written concretely as ${\int_{D}}{}g(\bm{x}){\mathrm{d}\Gamma(\bm{x})}=\int\sum_{\kappa=0}^{{N}}g(\kappa,\omega)\rho_{\kappa}(\kappa)\rho_{\omega}(\omega)\mathrm{d}\omega$, $(\omega,\kappa)\in D$, where $\rho_{\kappa}$ is the (discrete) probability mass function for the degrees, and $\rho_{\omega}$ is the (continuous) probability distribution function for the natural frequencies. In separate work [5], we examine the computation of this integral for the case when the problem can be recast as a PDE, and so the coupling sum in Eq. (1) and Eq. (3) can also be written as an Lebesgue integral. There, we consider standard Monte Carlo integration in addition to Gaussian quadrature and a repurposing of anchored ANOVA. These latter methods have the benefit of allowing integrals to be computed using only a few key virtual oscillators, with anchored ANOVA having the additional benefit of decreased scaling sensitivity to the number of random dimensions (2 in this paper). However, these benefits require that the original model be recast as PDEs continuous in both time and the random dimensions.

Here, we take the alternate approach of finding the $\alpha_{k}$ indirectly through least squares fitting, minimizing the squared residual norm $\sigma$ with respect to the coefficients $\alpha_{k}$ (here using a QR algorithm).

[TABLE]

where $\lim_{{n_{\text{samp}}}\rightarrow\infty}\hat{\sigma}=\sigma$ according to the law of large numbers, and the weights $w(\bm{x}_{i})$ are still to be decided.

For a (large) finite sample $\bm{x}_{i}$, $i=1,\,\ldots,\,{n_{\text{samp}}}$, if we take the partial derivative $\partial\hat{\sigma}/\partial\alpha_{k}$, and use the fact that $\sum_{i=1}^{n_{\text{samp}}}w(\bm{x}_{i})\psi^{(l)}(\bm{x}_{i})\psi^{(k)}(\bm{x}_{i})=0$ if ${l}\neq{k}$ (orthogonality) to remove some terms, then we find that

[TABLE]

So, as long as we accept that

[TABLE]

is a good approximation to

[TABLE]

we obtain the same formulas for the $\alpha_{k}$.

Suppose $\bm{x}$ has a density $\rho(\bm{x}\in D)$, so Eq. (13) can be written as

[TABLE]

If the $\bm{x}_{i}$ are chosen randomly in accordance with $\rho(\bm{x})$, then Eq. (12), where $w(\bm{x}_{i})=1/{n_{\text{samp}}}$, is a good approximation to Eq. (14). This is just Monte Carlo integration, and the law of large numbers gives

[TABLE]

5 Coarse Computational Modeling Tasks

Beyond their conceptual simplification value, collective (coarse) variables can be valuable in facilitating the computer-assisted study of complex dynamical systems by accelerating tasks such as direct simulation, continuation, stability and bifurcation analysis for different types of solutions. To accomplish this acceleration, the Equation-Free approach [25, 26] is predicated on the ability to map between corresponding fine and coarse descriptions of the same system.

This is accomplished through the definition of a restriction operator $R:\mathbb{R}^{{N}-1}\rightarrow\mathbb{R}^{M}$ which maps from fine states $\bm{\theta}(t)$ to corresponding coarse states $\bm{\alpha}(t)$ by minimizing the residual $\sigma(\bm{\alpha}(t))$ from Eq. (10). We also need to define the counterpart of restriction: a lifting operator $L:\mathbb{R}^{M}\rightarrow\mathbb{R}^{{N}-1}$ which maps $\bm{\alpha}$ vectors to $\bm{\theta}$ vectors by setting the $\theta_{i}$ values equal to the right-hand-side of the approximant in Eq. (8) evaluated at the corresponding $(\omega_{i},\kappa_{i})$.

One more important thing to note before proceeding to demonstrating the approach is that the lifting operator is, in general, a one-to-many relation; there are many fine realizations of the process that are mapped to the same coarse representation—coarse-graining (e.g. averaging) loses information. If the problem can be usefully coarse-grained, any of these consistent fine realizations, or the average of several of them, can be used practically in the definition of the coarse time-stepper below; we may think of the coarse-timestepper as the expected value over all such consistent realizations. In singularly perturbed multiscale problems one can clearly see how the memory of the details of the lifting are quickly forgotten, suggesting that any consistent fine realization is “good enough” to estimate this expectation [12, 23].

The $L$ and $R$ operators combine to define a coarse timestepper $\Phi_{{\tau},C}$, in

[TABLE]

This is the timestepper for the (unavailable) coarse-grained dynamical system, approximated through observing the results of short bursts of appropriately initialized fine-grained simulations. A single evaluation of this coarse time-stepper, by itself, does not provide any computational savings; it is the way we design, and process the results of, several such coarse time-steps that leads to computational benefits. Using traditional numerical analysis codes (initial value solvers, fixed point solvers) as templates for wrapper codes around the coarse timestepper, tasks like accelerated simulation, coarse-grained stability and bifurcation analysis, optimization, and controller design, can be performed. This wrapper technology is described in detail (and fruitfully used to explore model coarse-graining across disciplines) in a series of publications [50, 24]. What is important here is not the established wrapper algorithms technology; it is the selection of coarse observables, leading to the appropriate definition of the coarse time-stepper, that makes the entire program feasible and useful.

5.1 Coarse Initial Value Problems

We can use the coarse timestepper to accelerate the computation of dynamic trajectories of the system, through Coarse Projective Integration (CPI) [11, 32]. Given a coarse initial condition $\bm{\alpha}(t=0)$ we lift to a consistent fine scale state $L[\bm{\alpha}(t=0)]$ and use it to initialize a fine scale numerical integrator. We run for a short time $\tau$ (the inner step) and we record the final coarse state by restricting the corresponding fine state, $R[\bm{\theta}(t+\tau)]$. We use these two coarse states to estimate the coarse time derivative, which we then use in the forward Euler formula to project forward in time the coarse state for a (large, coarse) time step $h$ (the outer step). This constitutes the simplest coarse projective forward Euler integration scheme:

[TABLE]

A slightly more sophisticated approach would also take two points separated by an inner step size $\tau$ to approximate the rate-of-change of $\bm{\alpha}$, but only after first performing a healing integration in the fine equations [12]. If the lifted representation of the state $\bm{\alpha}(t)$, as projected from the previous timestep $\bm{\alpha}(t-h)$, is slightly off the hypothetical slow manifold in the fine space, this short healing trajectory dampens the fast components which are not captured in the coarse representation [51, 1, 13].

It is important to note here that in all our CPI computations we used a single network (with a single $\omega$ vector) to generate the polynomial basis functions and to lift to at every projective step. This can be thought of as a “single instance” CPI; one may also consider CPI for the expected behavior over all networks that share the same degree as well as $\omega$ distributions – in which case one should lift to many consistent network realizations and average over them. This issue will be examined more closely below.

Results of applying the simpler scheme to the coupled oscillator problem with coarse variables obtained by a $2D$ PCE fit are shown in Fig. 4. In general, the coarse time-stepper results can be used to estimate the coarse right-hand-side function $\dot{\bm{\alpha}}$, which is not available in closed form. On demand numerical estimates of this right-hand-side through short bursts of appropriately initialized fine simulation allow us to use other existing integrators, such as MATLAB’s ode45, to computationally approximate coarse trajectories (also shown in Fig. 4). Appendix B contains a quick illustration of the useful properties of such projective initial value solvers. It is shown there that, under reasonable conditions, the order of a projective integrator templated on a two-step Runge-Kutta initial value solver (including the additional estimation step for the coarse time derivatives) is the same as the order of the actual Runge-Kutta initial value solver. Fig. 5) confirms this for projective integration of the fine equations for our model.

This approximation of the coarse right-hand-side function can be used for other computational tasks besides the computation of dynamic trajectories.

5.2 Coarse Fixed Point Computation

The coarse time-stepper can be used to define a coarse difference

[TABLE]

Steady states of the fine timestepper are clearly zeroes of this difference; one expects that the zeroes of the coarse difference correspond to coarse steady states of the original problem. $F_{\tau}$ can therefore be used to find coarse steady states involving only ${M}$ variables. Iterative, matrix-free linear algebra lends itself to finding zeroes of such a problem in the absence of explicit equations for the dynamics of the coarse variables $\alpha_{j}$, We used a Krylov-type matrix free technique (Newton-Krylov GMRES) to converge to such coarse steady states.

A Newton-Krylov iteration to find such a state is depicted in Fig. 6. In Newton-Krylov GMRES (Generalized Minimal RESidual), the inner linear problem of an outer (nonlinear) Newton-type solver is solved by GMRES, in which the solution to ${\mathbf{B}}\bm{x}=\bm{b}$ is assembled in a space derived from the $n^{th}$ Krylov subspace $\{{\mathbf{B}}^{j}\bm{r}_{0}\}_{j=0,\ldots,n-1}$, where $\bm{r}_{0}$ is the residual of the initial iterate (and ${\mathbf{B}}$, which is not computed, is the Jacobian of Eq. (18)).

Here we work again with ${N}=196$ node networks. The basis polynomials are computed from a single network realization (a large, $10000$ node Chung-Lu network); because the support of the degree distributions for a ${N}=10000$ and for a ${N}=196$ network are not the same, the degrees and frequencies are normalized as described in §3. But now we construct 32 realizations of networks consistent with the chosen degree distribution, and perform our fixed point computation for each one of them. We do not regenerate polynomials for each of these realizations; we observe them on the “large sample” polynomials; this is justified in Fig. 6, where we recover the Probabilist’s Hermite polynomials as ${N}$ increases, by generating with sample moments. The effect of regenerating polynomials within a larger computation is considered further in [44, 45].

We report the sensitivity of the results to the basis size ${M}$ in Fig. 6, where the error bars are indicative of the variation across our $32$ network samples. We can also see in Fig. 6 that the dimension of the Krylov subspace (at the final Newton iteration before convergence) initially closely follows the size of the 2D basis we use (it uses ”all of ${M}$”) but later on plateaus.

5.3 Coarse stability computations: Eigenvalues and eigenvectors

(Coarse) eigenvalues of the Jacobian of the (coarse) difference map Eq. (18) upon convergence to its (coarse) fixed points can be used to establish the stability of said fixed points and help determine the nature of their potential (coarse) bifurcations. These eigenvalues ${\mu}_{i}$ are related by

[TABLE]

to the corresponding eigenvalues ${{\hat{\lambda}}}_{i}$ of the Jacobian of the (unavailable) coarse differential evolution equations; in turn, these should coincide with the leading eigenvalues $\lambda_{i}$ of the actual problem (the leading eigenvalues of the fine differential equations). For ${M}>3$ coarse eigenvalues $\mu_{i}$ were obtained through the Jacobian-free implicitly restarted Arnoldi Method (IRAM) [19, 33] applied to the coarse difference operator Eq. (18). For ${M}\leq 3$, a forward finite-difference Jacobian with a fixed step size of $0.001$ was computed and the eigenpairs calculated directly with the QZ algorithm implemented in MATLAB’s eig.

As the number of coarse variables is increased, and therefore the quality of the coarse approximation improves, one expects these coarse eigenvalue estimates to approach the leading eigenvalues of the analytical fine Jacobian of Eq. (3), located through any eigensolver. In Fig. 7, we demonstrate this convergence of the approximate eigenvalues ${\hat{\lambda}}_{i}$ to the leading fine eigenvalues $\lambda_{i}$ with increasing ${M}$. Fig. 7 shows the corresponding convergence of a coarse eigenvector (the one corresponding to the smallest absolute value of ${\hat{\lambda}}_{i}$) to the fine eigenvector corresponding to smallest absolute value of $\lambda_{i}$. Note that lifting is necessary to make a comparison between $\theta$ and $\alpha$ eigenvectors. On the other hand, the transformation Eq. (19) is necessary not because of our coarse and fine spaces, but because the ${\mu}$ eigenvalues come from the Jacobian of a finite-time flow map while the $\lambda$ eigenvalues come from the Jacobian of a vector of infinitesimal differential equations.

In performing computations involving finite differences, we used a value of $\sqrt{\text{machine precision}}$ (according to [42]), which is approximately $10^{-7}$ for IEEE standard double-precision floating point variables in MATLAB.

5.4 Coarse Continuation/Bifurcation Diagrams

To build a coarse bifurcation diagram (See Fig. 8), we performed pseudo-arclength continuation [18, 21] for the coarse fixed points. We computed branches of coarse solutions to $\bm{0}=F_{\tau}(\bm{\alpha};{K})$ as the global parameter ${K}$ is varied. To trace out these solution branches, steps were taken in (pseudo-)arclength along the branch rather than in ${K}$. This allows the continuation to extend naturally beyond turning points. At some point along this continued branch of solutions, one of the computed eigenvalues becomes positive. At this point, the line color is changed to indicate that the new branch comprises unstable solutions.

Beyond these branches (to the left) we know that a limit cycle solution arises: a periodic orbit characterized by one free rogue oscillator, which performs full rotations and only momentarily slows down as it passes through the remaining pack of clustered oscillators [34]. In dynamical systems terminology this is a “SNIPER” (saddle-node infinite period) bifurcation [49]. The insets in Fig. 8 show transient dynamics in terms of the synchronization index $r$ in $re^{i{\bar{\theta}}}=\frac{1}{{N}}\sum_{j=1}^{N}e^{i\theta_{j}}$ (the real magnitude of the complex Kuramoto order parameter; see e.g. [48]). The presence of this rogue oscillator means that the coarse representation of Eq. (8) is not particularly accurate/informative to the left of ${K}_{c}$, without explicitly including the rogue’s value of $\theta$ in the set of coarse variables, as was done, for example, in [35].

6 Discussion

In this paper we have demonstrated that a general network of coupled, intrinsically heterogeneous, oscillators can be usefully described using a small number of collective dynamic variables. These variables are the time-dependent coefficients of an expansion of the complete state of the network in terms of a set of orthogonal polynomials. The polynomials are products of univariate polynomials in the parameters describing the intrinsic heterogeneity of a given oscillator, and a structural heterogeneous property (here, the degree) indicative of the connectivity of the oscillators in the network. Our results extend previous work which only considered all-to-all coupled networks, in which the state of an oscillator was a function of only its intrinsic heterogeneity [35, 29]. Our expansion (and subsequent truncation of the expansion) in this form is motivated by the large body of work in the field of uncertainty quantification; the difference being that here we have heterogeneous parameters characterizing a single network, rather than many realizations of a dynamical system, each with different (uncertain) parameters. We anticipate that this new link between the two fields (network dynamics and UQ) may provide many more fruitful opportunities for mathematical/computational technology transfer that can enhance our understanding and ability to usefully describe and analyze dynamics on complex networks.

Although we have only considered Kuramoto-type oscillators in a specific Chung-Lu network, our methods do not rely on either the type of oscillator used or the specific network (as long as the mean degree is not small). Thus they should be widely applicable to many non-trivial networks of neurons which exhibit synchrony for some range of parameters

Using this reduced description of a network, we demonstrated a number of standard computational tasks using the equation-free framework, in which differential equations describing the evolution of the expansion coefficients are not explicitly derived, but rather estimated on-the-fly. Specifically, we demonstrated coarse projective integration, the computation of coarse fixed points and their stability, as well as parametric analysis through continuation.

The success of our method relied on the rapid development of correlations between the state of an oscillator and its heterogeneous identifying parameters, in this case, its intrinsic frequency and its degree. A potential shortcoming of the method would arise when such a strong dependence does not develop–that is, when “similar” oscillators do not behave “similarly” (e.g. when the initial conditions, or something more than just the degree, like the clustering coefficient of every node, matters). This implies that additional “heterogeneity dimensions” must be introduced, in analogy to when, say, a two-dimensional flow loses stability and becomes three-dimensional. One such case we have encountered [36] is when the oscillators in a network (an all-to-all network of Hodgkin-Huxley neurons) split in two subsets and in each subset a distinct relation of state to identity was established. Knowing the identity of the oscillator was not enough, in that case, to characterize dynamics; one also needed to know in which cluster the oscillator belonged. We could regard both our case as well as this other one as special cases where, for every oscillator identity, there is a distribution of behaviors–a strongly peaked unimodal distribution in this paper, and a strongly peaked bimodal distribution in [36]; for that matter, in our study of breakup in multiple communities/clusters, one obtained a multimodal distribution. We anticipate that, in the spirit of stochastic PDEs in physical space, our approach might be extended to evolve state distributions in heterogeneity space (as opposed to state functions in heterogeneity space).

In other networks it may be that the state of a node depends on more than just these two properties: networks with weighted edges provide an obvious context in which this may occur. We believe (and are actively pursuing this research direction) that the approach introduced here can also be usefully extended to help in determining reduced descriptions for such networks.

The tensor product basis used here relied on a lack of correlation across the heterogeneities. As mentioned in §4.1, this reliance can be overcome by generating the full multidimensional basis all at once, and our current work addresses this possibility. This is likely to be the case in dynamical systems in which it is is useful to retain multiple structural heterogeneities. Degree is one of several structural parameters–others include a node’s participation in motifs like triangles (complete graphs on three nodes), cherries (triangles with one edge removed) or its local clustering coefficient. This progression can be continued to higher-order statistics of the node connectivity by noting that using the degree of each node as the representative structural heterogeneity is equivalent to considering the per node counts of the two-node one-edge motif. As more structural heterogeneities are considered, it is reasonable to expect that these heterogeneities will not be statistically independent.

Author Contributions

TB, IGK, CL and CWG planned the work, and TB and YW performed the computations. All the authors contributed to the writing of the manuscript.

Funding

This work was partially supported by the US National Science Foundation and by NIH grant 1U01EB021956-01.

Acknowledgements

IGK, TB and CL gratefully acknowledge the hospitality of the Institute of Advanced Studies of the Technical University of Münich at Garching where IGK is a Hans Fischer Senior Fellow. TB and IGK thank Thomas Thiem for useful comments/discussions.

Appendix A Orthogonality of the tensor product basis for independent weightings

Here, we demonstrate that if we have polynomials in two variables that are orthogonal with respect to a two-dimensional weight function that is the product of two one-dimensional weight functions, then the orthogonal polynomials are the product of the respective one-dimensional orthogonal sets.

A.1 Problem setting

We suppose that we are given two one-dimensional orthogonal function sets defined on (possibly infinite) intervals $I_{x}$ and $I_{y}$ based on the weight functions $\rho(x)$, $x\in I_{x}$, and $\sigma(y)$, $y\in I_{y}$ and we want to construct a set of two-variable orthonormal polynomials, $\psi_{i,j}(x,y)$ of degrees $i$ in $x$ and $j$ in $y$ such that

[TABLE]

where

[TABLE]

and $\{\phi_{i}(x)\}$ and $\{\theta_{j}(y)\}$ are the orthonormal sets of polynomials corresponding to the one-dimensional systems based on the weights $\rho(x)$ and $\sigma(y)$ respectively, that is,

[TABLE]

and

[TABLE]

We will see that these restrictions lead to the product polynomials $\psi_{i,j}(x,y)=\phi_{i}(x)\theta_{j}(y)$.

A.2 Proof

Since $\psi_{i,j}(x,y)$ has maximum degree of $i$ in $x$ and $j$ in $y$, it can be written as

[TABLE]

Substituting Eq. (22) in Eq. (20) we get

[TABLE]

We prove the result by induction on $k=m+n$ showing that

[TABLE]

For $k=0$ from Eq. (23) we have immediately that $A^{2}_{0000}=1$ and $A_{ij00}A_{0000}=0$ for $i+j>0$. Since the signs of the one-dimensional orthogonal polynomials are arbitrary, we can take $A_{0000}=1$. The second relations implies that $A_{ij00}=0$ for $i+j>0$. This establishes the result for $k=0$.

Assuming that it is true for $k-1$, we examine Eq. (23) with $m+n=k$. Setting $i=m$ and $j=n$ in Eq. (23) we have $A^{2}_{mnmn}=1$ allowing us to choose $A_{mnmn}=1$. Then, for any $(i,j)\neq(m,n)$ we have $A_{ijmn}A_{mnmn}=0$ implying that $A_{ijmn}=0$, thus establishing the result.

Appendix B Accuracy of Integrators using a Chord Slope rather than the True Derivative.

In this paper we used fine-scale simulation to estimate the derivatives of the coarse variables from their chord slopes. We then used the chord slope in various numerical tasks such as finding steady states and integration. In this appendix we address the impact of the errors due to using chord slopes on integrators. Note that this is unrelated to any errors due to switching back and forth between fine and coarse variables which must be analyzed separately.

We have used both fixed step-size methods and off-the-shelf codes that automatically pick the step size to control an error estimate. We show below the following two lemmas that are valid for Runge-Kutta (RK) and multi-step methods (and probably valid for other classes of methods, but each needs to be examined individually): (1) When a fixed step-size method of order $p$ is used, if the chord is of length ${\tau}=\mathcal{O}({h}^{p})$ then the resulting method continues to be order $p$; (2) When an automatic method is used, if the error in the chord estimate of the derivative is bounded by the error control amount, then the one-step error of the resulting method is generally increased by no more than a constant factor. We discuss these results for the equation $\dot{y}=f(y)$.

Fixed step-size methods. Any method—whether RK or multi-step—advances from one step to the next by computing an approximation to $y(t+{h})-y(t)$ from a linear combination of approximations to the derivatives at various points—all multiplied by ${h}$ (and also a linear combination of past values in the case of multi-step methods). If the method has order $p$ it means that the approximation to the change from $t$ to $t+{h}$ has an error of ${\rm O}({h}^{p+1})$. Using a chord slope approximation to a derivative over distance ${\tau}$ gives an error in the derivative estimate of ${\rm O}({\tau})$, so the error in the integration formula is ${\rm O}({\tau}{h})$. Hence, if ${\tau}={\rm O}({h}^{p})$ then the additional error is ${\rm O}({h}^{p+1})$, so the formula remains of order $p$.

If the method handles stiff equations and hence needs a Jacobian to solve the nonlinear implicit equation, as long as the Jacobian is calculated numerically (which is the case in most automatic codes), this result does not change, as now the non-linear equation to be solved involves the chord slope rather than the derivative.

Automatic methods. Automatic method adjust the step size (and possibly also the order) from step to step to control an error estimate. In the early days of automatic integrators, there was a lot of discussion about whether to control the error per step, or the error per unit step. In error-per-unit-step, the error estimate for a step is controlled to be less than ${h}\epsilon$ where ${h}$ was the step size used and $\epsilon$ was the desired error. The reasoning behind this approach was that, if the errors accumulated more or less linearly from step to step (that depends very much on the stability or otherwise of the differential equation) then the error at the integration end point should be roughly proportional to the integration interval times $\epsilon$. However, if the system is stiff, errors are strongly damped from step to step, so the dominant error is the error in the last step, around $\epsilon$. Many automatic codes, having gone to the effort to estimate the error in a step, correct the solution by the error estimate (although there is now no direct estimate of the higher-order error, an assumption that higher derivatives behave similarly to lower derivatives can provide some comfort). An interesting aspect of this process is that if the original control was on error per step, the error control now becomes per-unit-step because of the increased order.

Regardless of the actual control mechanism used, all methods add a linear combination of derivative estimates multiplied by ${h}$ (and possibly a combination of past solution values) so that the influence of the slope estimates is of the form ${h}\sum_{i=0}^{k}\beta_{i}s_{i}$. If each $s_{i}$ is controlled to have an error of no greater than $\epsilon$, the contribution to the error is no worse than $A{h}\epsilon$, where $A=\sum_{i=0}^{k}|\beta_{i}|$. Thus the one step error is no worse than $A+1$ larger than before.

Most automatic codes provide both a relative error tolerance and an absolute error tolerance. Generally we are interested in errors relative to the size of the solution. However, when the solution is almost zero it may be impossible to get an error which is small relative to the solution, so the reason for an absolute error control is to defuse that problem. If either relative or absolute error is bounded, the step proceeds.

We demonstrate the convergence of error as ${h}$ is decreased for various ${\tau}$ on a two-ODE continuously-stirred tank reactor (CSTR) problem with in- and out-flow in Fig. 9. “Full” integration is by ode45 with the default tolerances, and the projective integration is a second-order Runge-Kutta scheme. For this example, we model the reaction

[TABLE]

where \bm{x}=\left[\begin{array}[]{c}x_{A}\\ x_{B}\end{array}\right] is the system state, and {\mathbf{J}}=\left[\begin{array}[]{cc}-a-V&b\\ a&-b-V\end{array}\right] is the Jacobian. The system parameters are the reaction rates $a=2$ and $b=1$, the input flow rate $V=0.4$, and the input concentrations $x_{A,\text{in}}=1$ and $x_{B,\text{in}}=0$. The eigenvalues of ${\mathbf{J}}$ are $\lambda_{1}=-(a+b+V)$ and $\lambda_{2}=-V$. We first take the initial condition to be \bm{x}=\left[\begin{array}[]{c}1\\ 0\end{array}\right], integrate until $t=-1/\lambda_{1}/2$, half of the first timescale. We use this state as the common initial condition for the integrators in Fig. 9. As can be seen, the error of the method decreases with the inner step size $\tau$ until a bound governed by the outer step size $h$ is reached. As predicted by the analysis above, this governing function is close to the curve $\rm O(h^{p})$.