On the estimation of the Mori-Zwanzig memory integral

Zhu, Yuanran; Dominy, Jason M.; Venturi, Daniele

doi:10.1063/1.5003467

We develop a thorough mathematical analysis to deduce conditions for the accuracy and convergence of different approximations of the memory integral in the Mori-Zwanzig (MZ) equation. In particular, we derive error bounds and sufficient convergence conditions for short-memory approximations, the t-model, and hierarchical (finite-memory) approximations. In addition, we derive useful upper bounds for the MZ memory integral, which allow us to estimate a priori the contribution of the MZ memory to the dynamics. Such upper bounds are easily computable for systems with finite-rank projections. Numerical examples are presented and discussed for linear and nonlinear dynamical systems evolving from random initial states.

I. INTRODUCTION

The Mori-Zwanzig (MZ) formulation is a technique from irreversible statistical mechanics that allows the development of formally exact evolution equations for the quantities of interest such as macroscopic observables in high-dimensional dynamical systems.^45,10,41,42 One of the main advantages of developing such exact equations is that they provide a theoretical starting point to avoid integrating the full (possibly high-dimensional) dynamical system and instead solve directly for the quantities of interest, thus reducing the computational cost significantly. Computing the solution to the Mori-Zwanzig equation, however, is a challenging task that relies on approximations and appropriate numerical schemes. One of the main difficulties lies in the approximation of the MZ memory integral (convolution term), which encodes the effects of the so-called orthogonal dynamics in the time evolution of the quantity of interest. The orthogonal dynamics are essentially a high-dimensional flow satisfying a complex integro-differential equation. Over the years, many techniques have been proposed for approximating the MZ memory integral, the most efficient ones being problem-dependent.^31,41 For example, in applications to statistical mechanics, Mori’s continued fraction method^27,18 has been quite successful in determining exact solutions to several prototype problems, such as the dynamics of the auto-correlation function of a tagged oscillator in a harmonic chain.^17,23 Other effective approaches to approximate the MZ memory integral rely on perturbation methods,^6,42,8 mode coupling theories,^2,31 functional approximation techniques,^{40,19,43,21,29} and data-driven methods.^5,25 In a parallel effort, the applied mathematics community has, in recent years, attempted to derive general easy-to-compute representations of the memory integral. In particular, various approximations such as the t-model,^10,12,34 the modified t-model,⁷ and, more recently, renormalized perturbation methods³⁶ were proposed to address the approximation of the memory integral in situations where there is no clear separation of scales between resolved and unresolved dynamics.

The main objective of this paper is to develop rigorous estimates of the memory integral and provide convergence analysis of different approximation models of the Mori-Zwanzig equation, such as the short-memory approximation,³³ the t-model,¹² and hierarchical methods.³⁵ In particular, we study the MZ equation corresponding to two broad classes of projection operators: (i) infinite-rank projections (e.g., Chorin’s projection¹⁰) and (ii) finite-rank projections (e.g., Mori’s projection²⁸). We develop our analysis for both state-space and probability density function (PDF) space formulations of the MZ equation. These two descriptions are connected by the same duality principle that pairs the Koopman and Frobenious-Perron operators.¹⁵

This paper is organized as follows. In Sec. II, we outline the general procedure to derive the MZ equation in the phase space and discuss the common choices of projection operators. In Sec. III, we derive error bounds for the MZ memory integral and provide the convergence analysis of different memory approximation methods, including the t-model,^34,10,12 the short-memory approximation,³³ and the hierarchical memory approximation technique.³⁵ Such estimates are built upon semigroup estimation methods we present in Appendix A. In Sec. IV, we present numerical examples demonstrating the accuracy of the memory approximation/estimation methods we develop throughout the paper. The main findings are summarized in Sec. V. Convergence analysis of the MZ memory term in the probability density function space formulation is presented in Appendix B.

II. THE MORI-ZWANZIG FORMULATION

Consider the nonlinear dynamical system

\frac{d x}{d t} = F (x), x (0) = x_{0}

(1)

evolving on a smooth manifold $S$ ⁠. For simplicity, let us assume that $S = R^{n}$ ⁠. We will consider the dynamics of scalar-valued observables $g : S \to C$ ⁠, and for concreteness, it will be desirable to identify structured spaces of such observable functions. In Ref. 15, it was argued that C^*-algebras of observables such as $L^{\infty} (S, C)$ (the space of all measureable, essentially bounded functions on $S$ ⁠) and $C_{0} (S, C)$ (the space of all continuous functions on $S$ ⁠, vanishing at infinity) make natural choices. In what follows, we do not require the observables to comprise a C^*-algebra, but we will want them to comprise a Banach space as the estimation theorems of Sec. III make an extensive use of the norm of this space. Having the structure of a Banach space of observables also gives greater context to the meaning of the linear operators $L$ ⁠, $K$ ⁠, $P$ ⁠, and $Q$ to be defined hereafter.

The dynamics of any scalar-valued observable g(x) (quantity of interest) can be expressed in terms of a semi-group $K (t, s)$ of operators acting on the Banach space of observables. This is the Koopman operator²⁴ which acts on the function g as

g (x (t)) = [K (t, s) g] (x (s)),

(2)

where

K (t, s) = e^{(t - s) L}, L g (x) = F (x) \cdot \nabla g (x) .

(3)

Rather than computing the Koopman operator applicable to all observables, it is often more tractable to compute the evolution only of a (closed) subspace of the quantities of interest. This subspace can be described conveniently by means of a projection operator $P$ with the subspace as its image. Both $P$ and the complementary projection $Q = I - P$ act on the space of observables. The nature, mathematical properties, and connections between $P$ and the observable g are discussed in detail in Ref. 15 and summarized in Sec. II A. For now, it suffices to assume that $P$ is a bounded linear operator and that $P^{2} = P$ ⁠. The MZ formalism describes the evolution of observables initially in the image of $P$ ⁠. Because the evolution of observables is governed by the semi-group $K (t, s)$ ⁠, we seek an evolution equation for $K (t, s) P$ ⁠. By using the definition of the Koopman operator (3) and the well-known Dyson identity

e^{t L} = e^{t Q L} + \int_{0}^{t} e^{s L} P L e^{(t - s) Q L} d s,

we obtain the operator equation

\frac{d}{d t} e^{t L} = e^{t L} P L + e^{t Q L} Q L + \int_{0}^{t} e^{s L} P L e^{(t - s) Q L} Q L d s .

(4)

By applying this equation to an observable function u₀, we obtain the well-known MZ equation in phase space

\frac{\partial}{\partial t} e^{t L} u_{0} = e^{t L} P L u_{0} + e^{t Q L} Q L u_{0} + \int_{0}^{t} e^{s L} P L e^{(t - s) Q L} Q L u_{0} d s .

(5)

Acting on the left with $P$ ⁠, we obtain the evolution equation for projected dynamics

\frac{\partial}{\partial t} P e^{t L} u_{0} = P e^{t L} P L u_{0} + \int_{0}^{t} P e^{s L} P L e^{(t - s) Q L} Q L u_{0} d s .

(6)

Note, in fact, that the projection of the second term in (5), i.e., $P e^{t Q L} Q L x_{0} = 0$ ⁠, vanishes since $P Q = 0$ ⁠.

A. Projection operators

In this section, we make a summary on the commonly used projection operators $P$ in the Mori-Zwanzig framework. To make our definition mathematically sound, we begin by assuming that the Liouville operator (3) acts on observable functions in a C^*-algebra $A$ ⁠, for instance, $L^{\infty} (M, Σ, μ)$ ⁠, where $M$ is a space such as $R^{N}$ ⁠, Σ is a σ-algebra on $M$ ⁠, and μ is a measure on Σ. Let $σ \in A_{*}$ be a positive linear functional on $A$ ⁠. We define the weighted pre-inner product

{⟨ f, g ⟩}_{σ} = σ (f^{*} g) .

This can be used to define a Hilbert space $H = L^{2} (M, σ)$ ⁠, which is the completion of the quotient space

H^{'} = {f \in A : σ (f^{*} f) < \infty} / {f \in A : σ (f^{*} f) = 0}

endowed with the inner product ⟨·,·⟩_σ. The L² norm induced by the inner product is denoted as ∥·∥_σ. In the rest of the paper, the positive linear functional σ is always induced by a probability distribution $\tilde{σ}$ through

σ (u) = \int_{M} \tilde{σ} (ω) u (ω) d ω,

where $\tilde{σ}$ is typically chosen to be the probability density of the initial condition ρ₀ or the equilibrium distribution ρ_eq in statistical physics. To conform to the literature, we also use notation ${⟨ \cdot, \cdot ⟩}_{ρ_{0}}$ ⁠, ${⟨ \cdot, \cdot ⟩}_{ρ_{e q}}, {⟨ \cdot, \cdot ⟩}_{e q}$ to represent the weighted inner product corresponding to different probability measures $\tilde{σ} (ω) d ω$ ⁠. With the Hilbert space determined, we now focus on the following two broad classes of orthogonal projections on $H$ ⁠.

1. Infinite-rank projections

The first class of projection operators to consider in this setting are the conditional expectations $P$ such that $P_{*} σ = σ$ ⁠. In this case, the properties of conditional expectations (in particular, that $P [P (f) g P (h)] = P (f) P (g) P (h)$ ³⁹) and the fact that $P_{*} σ = σ$ imply that

\begin{matrix} {⟨ P f, g ⟩}_{σ} & = & σ [{(P f)}^{*} g] = P_{*} (σ) [{(P f)}^{*} g] = σ [P ({(P f)}^{*} g)] = σ [{(P f)}^{*} (P g)] \\ {⟨ f, P g ⟩}_{σ} & = & σ [f^{*} P g] = P_{*} (σ) [f^{*} P g] = σ [P (f^{*} P g)] = σ [(P f^{*}) (P g)] = σ [{(P f)}^{*} (P g)] \end{matrix}

so that

{⟨ P f, g ⟩}_{σ} = {⟨ f, P g ⟩}_{σ}

for all $f, g \in H$ ⁠. It follows that

{⟨ Q f, g ⟩}_{σ} = {⟨ f, g ⟩}_{σ} - {⟨ P f, g ⟩}_{σ} = {⟨ f, g ⟩}_{σ} - {⟨ f, P g ⟩}_{σ} = {⟨ f, Q g ⟩}_{σ} .

Therefore both $P$ and $Q$ are self-adjoint (i.e., orthogonal) projections onto closed subspaces of $H$ and hence contractions $‖ P ‖_{σ} \leq 1$ ⁠, $‖ Q ‖_{σ} \leq 1$ ⁠. Chorin’s projection^12,10 is one of this class and is defined as

(P g) ({\hat{x}}_{0}) = \frac{\int_{- \infty}^{+ \infty} g (\hat{x} (t; {\hat{x}}_{0}, {\tilde{x}}_{0}), \tilde{x} (t; {\hat{x}}_{0}, {\tilde{x}}_{0})) ρ_{0} ({\hat{x}}_{0}, {\tilde{x}}_{0}) d {\tilde{x}}_{0}}{\int_{- \infty}^{+ \infty} ρ_{0} ({\hat{x}}_{0}, {\tilde{x}}_{0}) d {\tilde{x}}_{0}} = E_{ρ_{0}} [g | {\hat{x}}_{0}] .

(7)

Here $x (t; x_{0}) = (\hat{x} (t; {\hat{x}}_{0}, {\tilde{x}}_{0}), \tilde{x} (t; {\hat{x}}_{0}, {\tilde{x}}_{0}))$ is the flow map generated by (1) split into resolved (⁠ $\hat{x}$ ⁠) and unresoved (⁠ $\tilde{x}$ ⁠) variables, and $g (x) = g (\hat{x}, \tilde{x})$ is the quantity of interest. For Chorin’s projection, the positive functional σ defining the Hilbert space $H$ may be taken to be integration with respect to the probability measure $ρ_{0} ({\hat{x}}_{0}, {\tilde{x}}_{0})$ ⁠. Clearly, if x₀ is deterministic, then $ρ_{0} ({\hat{x}}_{0}, {\tilde{x}}_{0})$ is a product of Dirac delta functions. On the other hand, if $\hat{x} (0)$ and $\tilde{x} (0)$ are statistically independent, i.e., $ρ_{0} ({\hat{x}}_{0}, {\tilde{x}}_{0}) = {\hat{ρ}}_{0} ({\hat{x}}_{0}) {\tilde{ρ}}_{0} ({\tilde{x}}_{0})$ ⁠, then the conditional expectation $P$ simplifies to

(P u) ({\hat{x}}_{0}) = \int_{- \infty}^{+ \infty} u (\hat{x} (t; {\hat{x}}_{0}, {\tilde{x}}_{0}), \tilde{x} (t; {\hat{x}}_{0}, {\tilde{x}}_{0})) {\tilde{ρ}}_{0} ({\tilde{x}}_{0}) d {\tilde{x}}_{0} .

(8)

In the special case where $u (\hat{x}, \tilde{x}) = \hat{x} (t; {\hat{x}}_{0}, {\tilde{x}}_{0})$ ⁠, we have

(P \hat{x}) ({\hat{x}}_{0}) = \int_{- \infty}^{+ \infty} \hat{x} (t; {\hat{x}}_{0}, {\tilde{x}}_{0}) {\tilde{ρ}}_{0} ({\tilde{x}}_{0}) d {\tilde{x}}_{0},

(9)

i.e., the conditional expectation of the resolved variables $\hat{x} (t)$ given the initial condition ${\hat{x}}_{0}$ ⁠. This means that an integration of (9) with respect to ${\hat{ρ}}_{0} ({\hat{x}}_{0})$ yields the mean of the resolved variables, i.e.,

E_{ρ_{0}} [\hat{x} (t)] = \int_{- \infty}^{\infty} (P \hat{x}) ({\hat{x}}_{0}) {\hat{ρ}}_{0} ({\hat{x}}_{0}) d {\hat{x}}_{0} = \int_{- \infty}^{\infty} \hat{x} (t, x_{0}) ρ_{0} (x_{0}) d x_{0} .

(10)

Obviously, if the resolved variables $\hat{x} (t)$ evolve from a deterministic initial state ${\hat{x}}_{0}$ ⁠, then the conditional expectation (9) represents the the average of the reduced-order flow map $\hat{x} (t; {\hat{x}}_{0}, {\tilde{x}}_{0})$ with respect to the PDF of ${\tilde{x}}_{0}$ ⁠, i.e., the flow map

P e^{t L} \hat{x} (0) = X_{0} (t; {\hat{x}}_{0}) = \int_{- \infty}^{+ \infty} \hat{x} (t; {\hat{x}}_{0}, {\tilde{x}}_{0}) {\tilde{ρ}}_{0} ({\tilde{x}}_{0}) d {\tilde{x}}_{0} .

(11)

In this case, the MZ equation (6) is an exact (unclosed) partial differential equation (PDE) for the multivariate field $X_{0} (t, {\hat{x}}_{0})$ ⁠. In order to close such an equation, a mean field approximation of the type $P f (\hat{x}) = f (P \hat{x})$ was introduced by Chorin et al. in Refs. 10, 12, and 11, together with the assumption that the probability distribution of x₀ is invariant under the flow generated by (1).

2. Finite-rank projections

Another class of projections is defined by choosing a closed (typically finite-dimensional) linear subspace $V \subset H = L^{2} (M, σ)$ and letting $P$ to be the orthogonal projection onto V in the σ inner product. An example of such projection is Mori’s projection,⁴⁶ widely used in statistical physics. For finite-dimensional V, given a linearly independent set {u₁, …, u_M} ⊂ V that spans V, $P$ can be defined by first constructing the positive definite Gram matrix $G_{i j} = {⟨ u_{i}, u_{j} ⟩}_{σ}$ ⁠. Then

P f = \sum_{i, j = 1}^{M} {(G^{- 1})}_{i j} {⟨ u_{i}, f ⟩}_{σ} u_{j} .

(12)

This projection is orthogonal with respect to the $L_{σ}^{2}$ inner product. In statistical physics, a common choice for the positive functional σ that generates $H$ is integration with respect to the Gibbs canonical distribution $ρ_{e q} = e^{- β H} / Z$ ⁠, for the Hamiltonian $H = H (p, q)$ and the associated partition function Z. Here q are generalized coordinates while p are kinetic momenta.

III. ANALYSIS OF THE MEMORY INTEGRAL

In this section, we develop a thorough mathematical analysis of the MZ memory integral

\int_{0}^{t} P e^{s L} P L e^{(t - s) Q L} Q L u_{0} d s = \int_{0}^{t} P e^{s L} P e^{(t - s) L Q} L Q L u_{0} d s

(13)

and its approximation. We begin by describing the behavior of the semigroup norms $‖ e^{t L} ‖$ ⁠, $‖ e^{t Q L Q} ‖$ ⁠, and $‖ e^{t L Q} ‖$ as functions of time, for different choices of projection $P$ and different norms. As we will see, the analysis will give clear computable bounds only in some circumstances, illustrating the difficulty of this problem and the need for further development and insight.

A. Semigroup estimates

For any Liouville operator $L$ of the form (3) acting on $A = L^{\infty} (M, Σ, μ)$ and for any σ identified with an element of $L^{1} (M, Σ, μ)$ ⁠, the functional $L_{*} σ$ (assuming that σ lies in the domain of $L_{*}$ ⁠) is absolutely continuous with respect to σ (essentially because $L$ acts locally) and the Radon-Nikodym derivative²⁰ may be identified with the negative of the divergence of the vector field F with respect to the measure induced by σ, i.e.,

\frac{d L_{*} σ}{d σ} = - {div}_{σ} F, i.e., (L_{*} σ) (u) = - σ (u {div}_{σ} F),

(14)

where

{div}_{σ} F = \frac{\nabla \cdot (\tilde{σ} (x) F (x))}{\tilde{σ} (x)} .

(15)

When $M = R^{N}$ and σ has the form $σ (u) = \int_{R^{N}} \tilde{σ} (x) u (x) d x$ ⁠, this can be shown more directly using integration by parts. By assuming that $\tilde{σ} (x)$ or F_i decays to 0 at ∞, we have

\begin{matrix} L_{*} σ (u) = σ (L u) & = & \int_{R^{N}} \tilde{σ} \sum_{i = 1}^{N} F_{i} \frac{\partial u}{\partial x_{i}} d x = - \int_{R^{N}} u \sum_{i = 1}^{N} \frac{\partial}{\partial x_{i}} (\tilde{σ} F_{i}) d x = - \int_{R^{N}} u [\frac{1}{\tilde{σ}} \sum_{i = 1}^{N} \frac{\partial}{\partial x_{i}} (\tilde{σ} F_{i})] \tilde{σ} d x \end{matrix}

(16a)

= σ (u [- \frac{1}{\tilde{σ}} \sum_{i = 1}^{N} \frac{\partial}{\partial x_{i}} (\tilde{σ} F_{i})]) = σ (u [- {div}_{σ} F])

(16b)

from which we see

{div}_{σ} F = \frac{\nabla \cdot (\tilde{σ} F)}{\tilde{σ}} = \frac{1}{\tilde{σ}} \sum_{i = 1}^{N} \frac{\partial}{\partial x_{i}} (\tilde{σ} F_{i}) = \nabla \cdot F + F \cdot \nabla (\ln \tilde{σ}) .

Therefore,

{⟨ v, - u {div}_{σ} F ⟩}_{σ} = σ (- v^{*} u {div}_{σ} F) = (L_{*} σ) (v^{*} u) = σ (L (v^{*} u)) = σ (L {(v)}^{*} u + v^{*} L (u)) = {⟨ L v, u ⟩}_{σ} + {⟨ v, L u ⟩}_{σ}

and the following are equivalent: (i) $L_{*} σ = 0$ (i.e., σ is invariant); (ii) div_σF = 0; (iii) $L$ is the skew-adjoint with respect to the σ inner product. More generally, on $H = L^{2} (M, σ)$ ⁠, we find that

L + L^{†} = - {div}_{σ} F

so that the numerical abscissa ω^37,38 (i.e., logarithmic norm^14,32) of $L$ is given by

ω = sup_{0 \neq u \in H} \frac{R {⟨ u, L u ⟩}_{σ}}{{⟨ u, u ⟩}_{σ}} = sup_{0 \neq u \in H} \frac{{⟨ u, (L + L^{†}) u ⟩}_{σ}}{2 {⟨ u, u ⟩}_{σ}} = sup_{0 \neq u \in H} \frac{{⟨ u, - u {div}_{σ} F ⟩}_{σ}}{2 {⟨ u, u ⟩}_{σ}} = - \frac{1}{2} inf_{x} {div}_{σ} F (x) .

(17)

Using this numerical abscissa ω, we obtain the following $L_{σ}^{2}$ estimation of the Koopman semigroup:

‖ e^{t L} ‖_{L_{σ}^{2}} \leq e^{ω t},

(18)

and moreover, e^ωt is the smallest exponential function that bounds $‖ e^{t L} ‖_{σ}$ ⁠.¹⁴ When $P$ and $Q = I - P$ are the orthogonal projections on $L^{2} (M, σ)$ ⁠, we can observe that the numerical abscissa of $Q L Q$ is bounded by that of $L$ ⁠. In fact,

sup_{0 \neq u \in H} \frac{R {⟨ u, Q L Q u ⟩}_{σ}}{{⟨ u, u ⟩}_{σ}} = sup_{0 \neq u \in H} \frac{R {⟨ Q u, L Q u ⟩}_{σ}}{{⟨ u, u ⟩}_{σ}} = sup_{0 \neq u \in Im Q} \frac{R {⟨ u, L u ⟩}_{σ}}{{⟨ u, u ⟩}_{σ}} \leq sup_{0 \neq u \in H} \frac{R {⟨ u, L u ⟩}_{σ}}{{⟨ u, u ⟩}_{σ}} = ω

[see Eq. (17)] so that

‖ e^{t Q L Q} ‖_{L_{σ}^{2}} \leq e^{ω t} .

(19)

It should be noticed that this bound for the orthogonal semigroup is not necessarily tight. The tightness of the bound depends on the choice of projection $P$ and comes down to whether functions in the image of $Q$ can be chosen localized to regions where div_σF is close to its infimal value.

When estimating the MZ memory integral, we need to deal with the semigroup $e^{t L Q}$ ⁠. It turns out to be extremely difficult to prove strong continuity of such a semigroup in general, due to the unboundedness of $L Q$ ⁠. It is shown in Appendix A that when $P L Q$ is an unbounded operator, as is typical when $P$ is a conditional expectation, the semigroup $e^{t L Q}$ can only be bounded as

‖ e^{t L Q} ‖_{σ} \leq M_{Q} e^{t ω_{Q}}

(20)

for some $M_{Q} > 1$ ⁠, due to the fact that $‖ e^{t L Q} ‖_{σ}$ has an infinite slope at t = 0. More work is needed to obtain satisfactory, computable values for $M_{Q}$ and $ω_{Q}$ ⁠, in the case where $P$ and $Q$ are the infinite-rank projections. It is also shown that when either $P L Q$ or $L P$ is bounded, for example, when $P$ is a finite-rank projection, we can get computable semigroup bounds of the form

‖ e^{t L Q} ‖_{σ} \leq e^{ω_{Q} t} \leq e^{\frac{1}{2} (\sqrt{ω^{2} + ‖ P L Q ‖_{σ}^{2}} + ω) t},

(21a)

‖ e^{t L Q} ‖_{σ} \leq e^{ω_{Q} t} \leq e^{(ω + ‖ L P ‖_{σ}) t},

(21b)

where ω = −inf div_σF if we use the $L_{σ}^{2}$ estimation of $e^{t L}$ ⁠.

B. Memory growth

We begin by seeking to bound the MZ memory integral (B4) as a whole and build our analysis from there. A key assumption of our analysis is that the semigroup $e^{t L Q}$ is strongly continuous, i.e., the map $t \mapsto e^{t L Q} g$ is continuous in the norm topology on the space of observables for each fixed g.¹⁶ As we pointed out in Sec. III A, it is very difficult to rigorously prove a strong continuity of $e^{t L Q}$ and to obtain a computable upper bound for unbounded generators of the form $L Q$ ⁠. We leave this problem for future research, and assume for now that that there exist constants $M_{Q}$ and $ω_{Q}$ such that $‖ e^{t L Q} ‖ \leq M_{Q} e^{t ω_{Q}}$ ⁠. Throughout this section, ∥·∥ denotes a general Banach norm. We begin with the following simple estimate:

Theorem 3.1

(Memory growth). Let

e^{t L}

and

e^{t L Q}

be strongly continuous semigroups with upper bounds

‖ e^{t L} ‖ \leq M e^{t ω}

and

‖ e^{t L Q} ‖ \leq M_{Q} e^{t ω_{Q}}

. Then

∥\int_{0}^{t} P e^{s L} P L e^{(t - s) Q L} Q L u_{0} d s∥ \leq M_{0} (t),

(22)

where

M_{0} (t) = \{\begin{matrix} C_{1} t e^{t ω_{Q}}, & ω = ω_{Q}, \\ \frac{C_{1}}{ω - ω_{Q}} [e^{t ω} - e^{t ω_{Q}}], & ω \neq ω_{Q} \end{matrix}

(23)

and

C_{1} = M M_{Q} ‖ L Q L u_{0} ‖

is a constant. Clearly, lim_t→0M₀(t) = 0.

Proof.

We first rewrite the memory integral in the equivalent form

\int_{0}^{t} P e^{s L} P L e^{(t - s) Q L} Q L u_{0} d s = \int_{0}^{t} P e^{s L} P e^{(t - s) L Q} L Q L u_{0} d s .

Since

e^{t L}

and

e^{t L Q}

are assumed to be strongly continuous semigroups, we have the upper bounds

‖ e^{t L} ‖ \leq M e^{t ω}

⁠,

‖ e^{t L Q} ‖ \leq M_{Q} e^{t ω_{Q}}

⁠. Therefore

\begin{matrix} ∥ \int_{0}^{t} P e^{s L} P e^{(t - s) L Q} L Q L u_{0} d s ∥ & \leq & \int_{0}^{t} ‖ e^{s L} P e^{(t - s) L Q} L Q L u_{0} ‖ d s \\ \leq & M M_{Q} ‖ L Q L u_{0} ‖ \int_{0}^{t} e^{s (ω - ω_{Q})} d s \\ = & \{\begin{matrix} C_{1} t e^{t ω_{Q}}, & ω = ω_{Q}, \\ \frac{C_{1}}{ω - ω_{Q}} [e^{t ω} - e^{t ω_{Q}}], & ω \neq ω_{Q}, \end{matrix} \end{matrix}

where

C_{1} = M M_{Q} ‖ P ‖^{2} ‖ L Q L u_{0} ‖

⁠.

Theorem 3.1 provides an upper bound for the growth of the memory integral based on the assumption that $e^{t L}$ and $e^{t L Q}$ are strongly continuous semigroups. We emphasize that only for simple cases can such upper bounds be computed analytically (we will compute one of the cases later in Sec. IV) because of the fundamental difficulties in computing the upper bound of $e^{t L Q}$ ⁠. However, it will be shown later that, although the specific expression for M₀(t) is unknown, the form of it is already useful as it enables us to derive some verifiable theoretical predictions for general nonlinear systems.

C. Short memory approximation and the t-model

Theorem 3.1 can be employed to obtain upper bounds for well-known approximations of the memory integral. Let us begin with the t-model proposed in Ref. 12. This model relies on the approximation

\int_{0}^{t} P e^{s L} P L e^{(t - s) Q L} Q L u_{0} d s ≃ t e^{t L} P L Q L u_{0} (t - model) .

(24)

Theorem 3.2

(Memory approximation via the t-model¹²). Let

e^{t L}

and

e^{t L Q}

be strongly continuous semigroups with upper bounds

‖ e^{t L} ‖ \leq M e^{t ω}

and

‖ e^{t L Q} ‖ \leq M_{Q} e^{t ω_{Q}}

. Then

∥ \int_{0}^{t} P e^{s L} P L e^{(t - s) L Q} L Q L u_{0} d s - t P e^{t L} L Q L u_{0} ∥ \leq M_{1} (t),

where

M_{1} (t) = \{\begin{matrix} C_{1} (\frac{e^{t ω_{Q}} - e^{t ω}}{ω_{Q} - ω} + \frac{t e^{t ω}}{M_{Q}}), & ω \neq ω_{Q}, \\ C_{1} \frac{M_{Q} + 1}{M_{Q}} t e^{t ω}, & ω = ω_{Q} \end{matrix}

and

C_{1} = M M_{Q} ‖ P ‖^{2} ‖ L Q L u_{0} ‖

.

Proof.

By applying the triangle inequality, we obtain that

\begin{matrix} ∥\int_{0}^{t} P e^{s L} P e^{(t - s) L Q} L Q L u_{0} d s - t P e^{t L} P L Q L u_{0}∥ & \leq & (\int_{0}^{t} ‖ P ‖ ∥e^{s L}∥ ‖ P ‖ ∥e^{(t - s) L Q}∥ d s + t ‖ P ‖ ∥e^{t L}∥ ‖ P ‖) ‖ L Q L u_{0} ‖ \\ \leq & ‖ P ‖^{2} ‖ L Q L u_{0} ‖ (M M_{Q} \int_{0}^{t} e^{s ω} e^{(t - s) ω_{Q}} d s + t M e^{t ω}) \\ = & C_{1} e^{t ω} (\int_{0}^{t} e^{s (ω_{Q} - ω)} d s + \frac{t}{M_{Q}}) \\ = & \{\begin{matrix} C_{1} (\frac{e^{t ω_{Q}} - e^{t ω}}{ω_{Q} - ω} + \frac{t e^{t ω}}{M_{Q}}), & ω \neq ω_{Q}, \\ C_{1} \frac{M_{Q} + 1}{M_{Q}} t e^{t ω}, & ω = ω_{Q}, \end{matrix} \end{matrix}

where

C_{1} = M M_{Q} ‖ P ‖^{2} ‖ L Q L u_{0} ‖

⁠.

Theorem 3.2 provides an upper bound for the error associated with the t-model. The limit

lim_{t \to 0} M_{1} (t) = 0

(25)

guarantees the convergence of the t-model for short integration times. On the other hand, depending on the semigroup constants M, ω,

M_{Q}

⁠, and

ω_{Q}

(which may be estimated numerically), the error of the t-model may remain small for longer integration times (see the numerical results in Sec. IV B 2). Next, we study the short-memory approximation proposed in Ref. 33. The main idea is to replace the integration interval [0, t] in () by a shorter time interval [t − Δt, t], i.e.,

\int_{0}^{t} P e^{s L} P L e^{(t - s) Q L} Q L u_{0} d s ≃ \int_{t - Δ t}^{t} P e^{s L} P L e^{(t - s) Q L} Q L u_{0} d s (short - memory approximation),

where Δt ∈ [0, t] identifies the effective memory length. The following result provides an upper bound to the error associated with the short-memory approximation.

Theorem 3.3

(Short memory approximation³³). Let

e^{t L}

and

e^{t L Q}

be strongly continuous semigroups with upper bounds

‖ e^{t L} ‖ \leq M e^{t ω}

and

‖ e^{t L Q} ‖ \leq M_{Q} e^{t ω_{Q}}

. Then the following error estimate holds true:

∥ \int_{0}^{t} P e^{s L} P L e^{(t - s) Q L} Q L u_{0} d s - \int_{t - Δ t}^{t} P e^{s L} P L e^{(t - s) Q L} Q L u_{0} d s ∥ \leq M_{2} (t - Δ t, t),

where

M_{2} (Δ t, t) = \{\begin{matrix} C_{1} (t - Δ t) e^{t ω_{Q}}, & ω = ω_{Q}, \\ C_{1} e^{Δ t ω_{Q}} \frac{e^{(t - Δ t) ω} - e^{(t - Δ t) ω_{Q}}}{ω - ω_{Q}}, & ω \neq ω_{Q} \end{matrix}

and

C_{1} = M M_{Q} ‖ P ‖^{2} ‖ L Q L u_{0} ‖

.

We omit the proof due to its similarity to that of Theorem 3.1. Note that $lim_{Δ t \to t}$ M₂(Δt, t) = 0 for all finite t > 0.

D. Hierarchical memory approximation

An alternative way to approximate the memory integral (13) was proposed by Stinis in Ref. 35. The key idea is to repeatedly differentiate (13) with respect to time and establish a hierarchy of PDEs which can eventually be truncated or approximated at some level to provide an approximation of the memory. In this section, we derive this hierarchy of memory equations and perform a thorough theoretical analysis to establish the accuracy and convergence of the method. To this end, let us first define

w_{0} (t) = \int_{0}^{t} P e^{s L} P L e^{(t - s) Q L} Q L u_{0} d s

(26)

to be the memory integral (13). By differentiating $w$ ₀(t) with respect to time, we obtain

\frac{d w_{0} (t)}{d t} = P e^{t L} P L Q L u_{0} + w_{1} (t),

where

w_{1} (t) = \int_{0}^{t} P e^{s L} P L e^{(t - s) Q L} {(Q L)}^{2} u_{0} d s .

By iterating this procedure n times, we obtain

\frac{d w_{n - 1} (t)}{d t} = P e^{t L} P L {(Q L)}^{n - 1} u_{0} + w_{n} (t),

(27)

where

w_{n} (t) = \int_{0}^{t} P e^{s L} P L e^{(t - s) Q L} {(Q L)}^{n + 1} u_{0} d s .

(28)

The hierarchy of Eqs. (27) and (28) is equivalent to the following infinite-dimensional system of PDEs:

\{\begin{matrix} \frac{d w_{0} (t)}{d t} & = & P e^{t L} P L Q L u_{0} + w_{1} (t), \\ \frac{d w_{1} (t)}{d t} & = & P e^{t L} P L Q L Q L u_{0} + w_{2} (t), \\ ⋮ \\ \frac{d w_{n - 1} (t)}{d t} & = & P e^{t L} P L {(Q L)}^{n} u_{0} + w_{n} (t), \\ ⋮ \end{matrix}

(29)

evolving from the initial condition $w$ _i(0) = 0, i = 1, 2, … [see Eq. (28)]. With such an initial condition available, we can solve (29) with backward substitution, i.e., from the last equation to the first one, to obtain the following (exact) Dyson series representation of the memory integral (26):

\begin{matrix} w_{0} (t) & = & \int_{0}^{t} P e^{s L} P L Q L u_{0} d s + \int_{0}^{t} \int_{0}^{τ_{1}} P e^{s L} P L Q L Q L u_{0} d s d τ_{1} \\ + \dots + \int_{0}^{t} \int_{0}^{τ_{n - 1}} \dots \int_{0}^{τ_{1}} P e^{s L} P L {(Q L)}^{n} u_{0} d s d τ_{1} \dots d τ_{n - 1} + \dots . \end{matrix}

(30)

So far no approximation was introduced, i.e., the infinite-dimensional system (29) and the corresponding formal solution (30) are exact. To make progress in developing a computational scheme to estimate the memory integral (26), it is necessary to introduce approximations. The simplest of these rely on truncating the hierarchy (29) after n equations, while simultaneously introducing an approximation of the nth order memory integral $w$ _n(t). We denote such an approximation as $w_{n}^{e_{n}} (t)$ ⁠. The truncated system takes the form

\{\begin{matrix} \frac{d w_{0}^{n} (t)}{d t} & = & P e^{t L} P L Q L u_{0} + w_{1}^{n} (t), \\ \frac{d w_{1}^{n} (t)}{d t} & = & P e^{t L} P L Q L Q L u_{0} + w_{2}^{n} (t), \\ ⋮ \\ \frac{d w_{n - 1}^{n} (t)}{d t} & = & P e^{t L} P L {(Q L)}^{n} u_{0} + w_{n}^{e_{n}} (t) . \end{matrix}

(31)

The notation $w_{j}^{n} (t)$ (j = 0, …, n − 1) emphasizes that the solution to (31) is, in general, different from the solution to (29). The initial condition of the system can be set as $w_{i}^{n} (0) = 0$ ⁠, for all i = 0, …, n − 1. By using backward substitution, this yields the following formal solution:

\begin{matrix} w_{0}^{n} (t) & = & \int_{0}^{t} P e^{s L} P L Q L u_{0} d s + \int_{0}^{t} \int_{0}^{τ_{1}} P e^{s L} P L Q L Q L u_{0} d s d τ_{1} \\ + \dots + \int_{0}^{t} \int_{0}^{τ_{n - 1}} \dots \int_{0}^{τ_{1}} P e^{s L} P L {(Q L)}^{n} u_{0} d s d τ_{1} \dots d τ_{n - 1} \\ + \int_{0}^{t} \int_{0}^{τ_{n - 1}} \dots \int_{0}^{τ_{1}} w_{n}^{e_{n}} (s) d s d τ_{1} \dots d τ_{n - 1}, \end{matrix}

(32)

representing an approximation of the memory integral (26). Note that, for a given system, such approximation depends only on the number of equations n in (31) and on the choice of the approximation $w_{n}^{e_{n}} (t)$ ⁠. In the present paper, we consider the following choices:

Approximation by truncation (H-model),

w_{n}^{e_{n}} (t) = 0 .

(33)

Type-I finite memory approximation (FMA),

w_{n}^{e_{n}} (t) = \int_{max (0, t - Δ t_{n})}^{t} P e^{s L} P L e^{(t - s) Q L} {(Q L)}^{n + 1} u_{0} d s .

(34)

Type-II finite memory approximation,

w_{n}^{e_{n}} (t) = \int_{min (t, t_{n})}^{t} P e^{s L} P L e^{(t - s) Q L} {(Q L)}^{n + 1} u_{0} d s .

(35)

H_t-model,

w_{n}^{e_{n}} (t) = t P e^{t L} P L {(Q L)}^{n + 1} u_{0} .

(36)

The quantities t_n and Δt_n appearing in (34) and (35) will be defined in Theorems 3.5 and 3.6, respectively. The first approximation, i.e., (33), is a truncation of the hierarchy obtained by assuming that $w$ _n(t) = 0. Such approximation was originally proposed by Stinis in Ref. 35, and we shall call it the H-model. The Type-I finite memory approximation (FMA) is obtained by applying the short memory approximation to the nth order memory integral $w$ _n(t). The Type-II finite memory approximation (FMA) is a modified version of the Type-I, with a larger memory band. The H_t-model approximation is based on replacing the nth order memory integral $w$ _n(t) with a classical t-model. Note that in this setting, the classical t-model approximation proposed by Chorin and Stinis¹² is equivalent to a zeroth-order H_t-model approximation.

Hereafter, we present a thorough mathematical analysis that aims at estimating the error $‖ w_{0} (t) - w_{0}^{n} (t) ‖$ ⁠, where $w$ ₀(t) is the full memory at time t [see (26) or (30)], while $w_{0}^{n} (t)$ is the solution of the truncated hierarchy (31), with $w_{n}^{e_{n}} (t)$ given by (33), (34), (35), or (36). With such error estimates available, we can infer whether the approximation of the full memory $w$ ₀(t) with $w_{0}^{n} (t)$ is accurate and, more importantly, if the algorithm to approximate the memory integral converges. To the best of our knowledge, this is the first time a rigorous convergence analysis is performed on various approximations of the MZ memory integral. It turns out that the distance $‖ w_{0} (t) - w_{0}^{n} (t) ‖$ can be controlled through the construction of the hierarchy under some constraint on the initial condition.

1. The H-model

Setting $w_{n}^{e_{n}} (t) = 0$ in (31) yields an approximation by truncation, which we will refer to as the H-model (hierarchical model). Such a model was originally proposed by Stinis in Ref. 35. Hereafter we provide error estimates and convergence results for this model. In particular, we derive an upper bound for the error $‖ w_{0} (t) - w_{0}^{n} (t) ‖$ and sufficient conditions for the convergence of the reduced-order dynamical system. Such conditions are problem dependent, i.e., they involve the Liouvillian $L$ ⁠, the initial condition u₀, and the projection $P$ ⁠.

Theorem 3.4

(Accuracy of the H-model). Let

e^{t L}

and

e^{t L Q}

be strongly continuous semigroups with upper bounds

‖ e^{t L} ‖ \leq M e^{t ω}

and

‖ e^{t L Q} ‖ \leq M_{Q} e^{t ω_{Q}}

, and let T > 0 be a fixed integration time. For some fixed n, let

α_{j} = \frac{‖ {(L Q)}^{j + 1} L u_{0} ‖}{‖ {(L Q)}^{j} L u_{0} ‖}, 1 \leq j \leq n .

(37)

Then, for any 1 ≤ p ≤ n and all t ∈ [0, T], we have

‖ w_{0} (t) - w_{0}^{p} (t) ‖ \leq M_{3}^{p} (t) \leq M_{3}^{p} (T),

where

M_{3}^{p} (t) = C_{1} A_{1} A_{2} \frac{t^{p + 1}}{(p + 1)!} \prod_{j = 1}^{p} α_{j}, C_{1} = ‖ L Q L u_{0} ‖ M M_{Q},

and

A_{1} = max_{s \in [0, T]} e^{s (ω - ω_{Q})} = \{\begin{matrix} 1, & ω \leq ω_{Q} \\ e^{T (ω - ω_{Q})}, & ω \geq ω_{Q}, \end{matrix} A_{2} = max_{s \in [0, T]} e^{s ω_{Q}} = \{\begin{matrix} 1, & ω_{Q} \leq 0 \\ e^{T ω_{Q}}, & ω_{Q} \geq 0 . \end{matrix}

(38)

Proof.

We begin with the expression for the difference between the memory term

w

₀ and its approximation

w_{0}^{p}

⁠,

w_{0} (t) - w_{0}^{p} (t) = \int_{0}^{t} \int_{0}^{τ_{p}} \dots \int_{0}^{τ_{2}} \int_{0}^{τ_{1}} P e^{s L} P e^{(τ_{1} - s) L Q} {(L Q)}^{n + 1} L u_{0} d s d τ_{1} \dots d τ_{p} .

(39)

Since

e^{t L}

and

e^{t L Q}

are strongly continuous semigroups, we have

‖ e^{t L} ‖ \leq M e^{ω t}

and

‖ e^{t L Q} ‖ \leq M_{Q} e^{ω_{Q} t}

⁠. By using Cauchy’s formula for repeated integration, we bound the norm of the error () as

\begin{matrix} ‖ w_{0} (t) - w_{0}^{p} (t) ‖ & \leq & \int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} \int_{0}^{σ} ‖ P e^{s L} P e^{(σ - s) L Q} {(L Q)}^{p + 1} L u_{0} ‖ d s d σ \\ \leq & ‖ P ‖^{2} M M_{Q} ‖ {(L Q)}^{p + 1} L u_{0} ‖ \int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} \int_{0}^{σ} e^{s ω} e^{(σ - s) ω_{Q}} d s d σ \\ \leq & C_{1} (\prod_{j = 1}^{p} α_{j}) \underset{f_{p} (t, ω, ω_{Q})}{\underset{︸}{\int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} \int_{0}^{σ} e^{s ω} e^{(σ - s) ω_{Q}} d s d σ}} \\ = & C_{1} (\prod_{j = 1}^{p} α_{j}) f_{p} (t, ω, ω_{Q}), \end{matrix}

(40)

where

C_{1} = ‖ P ‖^{2} ‖ L Q L u_{0} ‖ M M_{Q}

as before. The function

f_{p} (t, ω, ω_{Q})

may be bounded from above as

\begin{matrix} f_{p} (t, ω, ω_{Q}) & \leq & A_{1} A_{2} \int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} \int_{0}^{σ} d s d σ \\ = & A_{1} A_{2} \frac{t^{p + 1}}{(p + 1)!} . \end{matrix}

Hence, we have

‖ w_{0} (t) - w_{0}^{p} (t) ‖ \leq C_{1} A_{1} A_{2} (\prod_{j = 1}^{p} α_{j}) \frac{t^{p + 1}}{(p + 1)!} = M_{3}^{p} (t) .

Theorem 3.4 states that for a given dynamical system (represented by $L$ ⁠) and quantity of interest (represented by $P$ ⁠), the error bound $M_{3}^{p} (t)$ is strongly related to {α_j} which is ultimately determined by the initial condition x₀. It turns out that by bounding {α_j}, we can control $M_{3}^{p} (t)$ and therefore the overall error $‖ w_{0} (t) - w_{0}^{p} (t) ‖$ ⁠. The following corollaries discuss sufficient conditions such that the error $‖ w_{0} (T) - w_{0}^{n} (T) ‖$ decays as we increase the differentiation order n for fixed time T > 0.

Corollary 3.4.1

(Uniform convergence of the H-model). If {α_j} in Theorem 3.4 satisfy

α_{j} < \frac{j + 1}{T}, 1 \leq j \leq n,

(41)

for any fixed time T > 0, then there exists a sequence of constants δ₁ > δ₂ > ⋯ > δ_n such that

‖ w_{0} (T) - w_{0}^{p} (T) ‖ \leq δ_{p}, 1 \leq p \leq n .

Proof.

Evaluating () at any fixed (finite) time, T > 0 yields

\begin{matrix} ‖ w_{0} (T) - w_{0}^{p} (T) ‖ & \leq & C_{2} (\prod_{j = 1}^{p} α_{i}) f_{p} (T, ω, ω_{Q}) \leq C_{2} (\prod_{j = 1}^{p} α_{j}) \frac{T^{p + 1}}{(p + 1)!}, \\ ‖ w_{0} (T) - w_{0}^{p + 1} (T) ‖ & \leq & C_{2} (\prod_{j = 1}^{p + 1} α_{j}) \frac{T^{p + 2}}{(p + 2)!}, \end{matrix}

where C₂ = C₂(T) = C₁A₁A₂. If there exists δ_p ≥ 0 such that

‖ w_{0} (T) - w_{0}^{p} (T) ‖ \leq C_{2} (\prod_{j = 1}^{p} α_{j}) \frac{T^{p + 1}}{(p + 1)!} \leq δ_{p},

then there exists a δ_p+1 such that

‖ w_{0} (T) - w_{0}^{p + 1} (T) ‖ \leq C_{2} (\prod_{j = 1}^{p} α_{j}) \frac{T^{p + 1}}{(p + 1)!} \frac{α_{p + 1} T}{p + 2} \leq δ_{p + 1} < δ_{p}

since α_p+1 < (p + 2)/T. Moreover, the condition α_j < (j + 1)/T holds for all 1 ≤ j ≤ n. Therefore, we conclude that for any fixed time T > 0, there exists a sequence of constants δ₁ > δ₂ > ⋯ > δ_n such that

‖ w_{0} (T) - w_{0}^{p} (T) ‖ \leq δ_{p}

⁠, where 1 ≤ p ≤ n.

Corollary 3.4.1 provides a sufficient condition for the error $‖ w_{0} (t) - w_{0}^{p} (t) ‖$ to decrease monotonically as we increase p in (31). A stronger condition that yields an asymptotically decaying error bound is given by the following Corollary.

Corollary 3.4.2

(Asymptotic convergence of the H-model). If α_j in Theorem 3.4 satisfies

α_{j} < C, 1 \leq j < + \infty

(42)

for some positive constant C, then for any fixed time T > 0 and arbitrary δ > 0, there exists a constant 1 ≤ p < +∞ such that for all n > p,

‖ w_{0} (T) - w_{0}^{n} (T) ‖ \leq δ .

Proof.

By introducing the condition α_j < C in the Proof of Theorem 3.4, we obtain

‖ w_{0} (T) - w_{0}^{p} (T) ‖ \leq C_{2} (\prod_{j = 1}^{p} α_{j}) \frac{T^{p + 1}}{(p + 1)!} \leq C_{2} T \frac{{(C T)}^{p}}{(p + 1)!} for all 1 < p < + \infty .

The limit

lim_{p \to + \infty} C_{2} T \frac{{(C T)}^{p}}{(p + 1)!} = 0

allows us to conclude that there exists a constant 1 < p < +∞ such that for all n > p,

‖ w_{0} (T) - w_{0}^{n} (T) ‖ \leq δ

⁠.

An interesting consequence of Corollary 3.4.2 is the existence of a convergence barrier, i.e., a “hump” in the error plot $‖ w_{0} (T) - w_{0}^{p} (T) ‖$ versus p generated by the H-model. While Corollary 3.4.2 only shows that behavior for an upper bound of the error, not directly the error itself, the feature is often found in the actual errors associated with numerical methods based on these ideas. The following corollary shows that the requirements on {α_j} can be dropped (we still need α_j < +∞) if we consider relatively short integration times T.

Corollary 3.4.3

(Short-time convergence of the H-model). For any integer n for which α_j < ∞ for 1 ≤ j ≤ n and any sequence of constants δ₁ > δ₂ > ⋯ > δ_n > 0, there exists a fixed time T > 0 such that

‖ w_{0} (T) - w_{0}^{p} (T) ‖ \leq δ_{p}

for 1 ≤ p ≤ n.

Proof.

Since α_j < +∞, we can choose

C = max_{1 \leq j \leq n}

α_j. By following the same steps we used in the Proof of Theorem 3.4, we conclude that for

T \leq \frac{1}{C} min_{1 \leq p \leq n} {[\frac{C (p + 1)!}{C_{2}} δ_{p}]}^{\frac{1}{p + 1}},

the errors satisfy

‖ w_{0} (T) - w_{0}^{p} (T) ‖ \leq C_{2} (\prod_{j = 1}^{p} α_{j}) \frac{T^{p + 1}}{(p + 1)!} \leq \frac{C_{2}}{C} \frac{{(C T)}^{p + 1}}{(p + 1)!} \leq δ_{p}

as desired, for all 1 ≤ p ≤ n.

Corollaries 3.4.1 and 3.4.2 provide sufficient conditions for the error $‖ w_{0} (T) - w_{0}^{n} (T) ‖$ generated by the H-model to decay as we increase the truncation order n. However, we still need to answer the important question of whether the H-model actually provides accurate results for a given nonlinear dynamics (⁠ $L$ ⁠), quantity of interest (⁠ $P$ ⁠) and initial state x₀. Corollary 3.4.3 provides a partial answer to this question by showing that, at least in the short time period, condition (41) is always satisfied (assuming that {α_j} are finite). This guarantees the short-time convergence of the H-model for any reasonably smooth nonlinear dynamical system and almost any observable. However, for longer integration times T, convergence of the H-model for arbitrary nonlinear dynamical systems cannot be established in general, which means that we need to proceed on a case-by-case basis by applying Theorem 3.4 or by checking whether the hypotheses of Corollary 3.4.1 or Corollary 3.4.2 are satisfied. On the other hand, the convergence of the H-model can be established for any finite integration time in the case of linear dynamical systems, as we have recently shown in Ref. 44. From a practical viewpoint, the implementation of the H-model requires computing ${(L Q)}^{n} L x_{0}$ to high-order in n. This is not straightforward if the dynamical systems are nonlinear. However, for linear systems, such terms can be easily computed (see Sec. III E).

2. Type-I finite memory approximation (FMA)

The Type-I finite memory approximation is obtained by solving the system (31) with $w_{n}^{e_{n}} (t)$ given by (34). As before, we first derive an upper bound for $‖ w_{0} (t) - w_{0}^{n} (t) ‖$ and then discuss sufficient conditions for convergence. Such conditions basically control the growth of an upper bound on $‖ w_{0} (t) - w_{0}^{n} (t) ‖$ ⁠.

Theorem 3.5

(Accuracy of the Type-I FMA). Let

e^{t L}

and

e^{t L Q}

be strongly continuous semigroups and let T > 0 be a fixed integration time. If

α_{j} = \frac{‖ {(L Q)}^{j + 1} L u_{0} ‖}{‖ {(L Q)}^{j} L u_{0} ‖}, 1 \leq j \leq n,

(43)

then for each 1 ≤ p ≤ n and for Δt_p ≤ t ≤ T,

‖ w_{0} (t) - w_{0}^{p} (t) ‖ \leq M_{4}^{p} (t),

where

M_{4}^{p} (t) = C_{1} A_{1} A_{2} (\prod_{i = 1}^{p} α_{i}) \frac{{(t - Δ t_{p})}^{p + 1}}{(p + 1)!}

and C₁, A₁, A₂ are as in Theorem 3.4.

Proof.

The error at the pth level is of the form

w_{p} (t) - w_{p}^{e_{p}} (t) = \int_{0}^{max (0, t - Δ t_{p})} P e^{s L} P L e^{(t - s) Q L} {(Q L)}^{p + 1} u_{0} d s,

and the error at the zeroth level is

\begin{matrix} w_{0} (t) - w_{0}^{p} (t) & = & \int_{0}^{t} \int_{0}^{τ_{p}} \dots \int_{0}^{τ_{2}} [w_{n} (τ_{1}) - w_{n}^{e_{p}} (τ_{1})] d τ_{1} \dots d τ_{p} \\ = & \int_{0}^{t} \int_{0}^{τ_{p}} \dots \int_{0}^{τ_{2}} \int_{0}^{max (0, τ_{1} - Δ t_{p})} P e^{s L} P e^{(τ_{1} - s) L Q} {(L Q)}^{p + 1} L u_{0} d s d τ_{1} \dots d τ_{p} \\ = & \int_{Δ t_{p}}^{t} \int_{Δ t_{p}}^{τ_{p}} \dots \int_{Δ t_{p}}^{τ_{2}} \int_{0}^{τ_{1} - Δ t_{p}} P e^{s L} P e^{(τ_{1} - s) L Q} {(L Q)}^{p + 1} L u_{0} d s d τ_{1} \dots d τ_{p} \\ = & \int_{Δ t_{p}}^{t} \int_{Δ t_{p}}^{τ_{p}} \dots \int_{0}^{τ_{2} - Δ t_{p}} \int_{0}^{{\tilde{τ}}_{1}} P e^{s L} P e^{({\tilde{τ}}_{1} + Δ t_{p} - s) L Q} {(L Q)}^{p + 1} L u_{0} d s d {\tilde{τ}}_{1} \dots d τ_{p} \\ ⋮ \\ = & \int_{0}^{max (0, t - Δ t_{p})} \int_{0}^{{\tilde{τ}}_{p}} \dots \int_{0}^{{\tilde{τ}}_{2}} \int_{0}^{{\tilde{τ}}_{1}} P e^{s L} P e^{({\tilde{τ}}_{1} + Δ t_{p} - s) L Q} {(L Q)}^{p + 1} L u_{0} d s d {\tilde{τ}}_{1} \dots d {\tilde{τ}}_{p} . \end{matrix}

The norm of this error may be bounded as

\begin{matrix} ‖ w_{0} (t) - w_{0}^{p} (t) ‖ & \leq & \int_{0}^{max (0, t - Δ t_{p})} \int_{0}^{{\tilde{τ}}_{p}} \dots \int_{0}^{{\tilde{τ}}_{2}} \int_{0}^{{\tilde{τ}}_{1}} ∥P e^{s L} P e^{({\tilde{τ}}_{1} + Δ t_{p} - s) L Q} {(L Q)}^{p + 1} L u_{0}∥ d s d {\tilde{τ}}_{1} \dots d {\tilde{τ}}_{p} \\ \leq & C_{1} (\prod_{j = 1}^{p} α_{j}) \int_{0}^{max (0, t - Δ t_{p})} \int_{0}^{{\tilde{τ}}_{p}} \dots \int_{0}^{{\tilde{τ}}_{2}} \int_{0}^{{\tilde{τ}}_{1}} e^{s (ω - ω_{Q})} e^{({\tilde{τ}}_{1} + Δ t_{p}) ω_{Q}} d s d {\tilde{τ}}_{1} \dots d {\tilde{τ}}_{p} \\ \leq & C_{1} (\prod_{j = 1}^{p} α_{j}) f_{p} (t, Δ t_{p}, ω, ω_{Q}), \end{matrix}

where

f_{p} (t, Δ t_{p}, ω, ω_{Q}) = \int_{0}^{max (0, t - Δ t_{p})} \int_{0}^{{\tilde{τ}}_{p}} \dots \int_{0}^{{\tilde{τ}}_{2}} \int_{0}^{{\tilde{τ}}_{1}} e^{s (ω - ω_{Q})} e^{({\tilde{τ}}_{1} + Δ t_{p}) ω_{Q}} d s d {\tilde{τ}}_{1} \dots d {\tilde{τ}}_{p} .

If we bound f_p as

\begin{matrix} f_{p} (t, Δ t_{p}, ω, ω_{Q}) & \leq & A_{1} A_{2} \int_{0}^{max (0, t - Δ t_{p})} \int_{0}^{{\tilde{τ}}_{p}} \dots \int_{0}^{{\tilde{τ}}_{2}} \int_{0}^{{\tilde{τ}}_{1}} d s d {\tilde{τ}}_{1} \dots d {\tilde{τ}}_{p} \\ = & \{\begin{matrix} 0, & 0 \leq t \leq Δ t_{p}, \\ A_{1} A_{2} \frac{{(t - Δ t_{p})}^{p + 1}}{(p + 1)!}, & t \geq Δ t_{p}, \end{matrix} \end{matrix}

where A₁, A₂ are defined in (), then we have that

‖ w_{0} (t) - w_{0}^{p} (t) ‖ \leq C_{1} A_{1} A_{2} (\prod_{j = 1}^{p} α_{j}) \frac{{(t - Δ t_{p})}^{p + 1}}{(p + 1)!} = M_{4}^{p} (t) .

We notice that if the effective memory band at each level decreases as we increase the differentiation order p, then we can control the error $‖ w_{0} (t) - w_{0}^{n} (t) ‖$ ⁠. The following corollary provides a sufficient condition that guarantees this sort of control of the error.

Corollary 3.5.1

(Uniform convergence of the Type-I FMA). If α_j in Therorem 3.5 satisfy

α_{j} < (j + 1) {[\frac{δ j!}{C_{1} A_{1} A_{2} (\prod_{k = 1}^{j - 1} α_{k})}]}^{- \frac{1}{j}}, 1 \leq j \leq n,

(44)

then for any T > 0 and δ > 0, there exists an ordered sequence Δt_n < Δt_n−1 < ⋯ < Δt₁ < T such that

‖ w_{0} (T) - w_{0}^{p} (T) ‖ \leq δ, 1 \leq p \leq n,

which satisfies

Δ t_{p} \leq T - {[\frac{δ (p + 1)!}{C_{1} A_{1} A_{2} (\prod_{j = 1}^{p} α_{j})}]}^{\frac{1}{p + 1}} .

(45)

Proof.

For 1 ≤ p ≤ n, we set

‖ w_{0} (t) - w_{0}^{p} (t) ‖ \leq C_{1} A_{1} A_{2} (\prod_{j = 1}^{p} α_{j}) \frac{{(t - Δ t_{p})}^{p + 1}}{(p + 1)!} \leq δ .

This yields the following requirement on Δt_p:

Δ t_{p} \geq T - {[\frac{δ (p + 1)!}{C_{1} A_{1} A_{2} (\prod_{j = 1}^{p} α_{j})}]}^{\frac{1}{p + 1}} .

(46)

Since hypothesis () holds, it is easy to check that the lower bound on each Δt_p satisfies

T - {[\frac{δ (p + 1)!}{C_{1} A_{1} A_{2} (\prod_{j = 1}^{p} α_{j})}]}^{\frac{1}{p + 1}} < T - {[\frac{δ p!}{C_{1} A_{1} A_{2} (\prod_{j = 1}^{p - 1} α_{j})}]}^{\frac{1}{p}}, Δ t_{p} > Δ t_{p - 1} .

Therefore, by using the equality in () to define a sequence of Δt_n, we find that it is a decreasing time sequence 0 < Δt_n < Δt_n−1 < ⋯ < Δt₁ < T such that

‖ w_{0} (T) - w_{0}^{n} (T) ‖ \leq δ

holds for all t ∈ [0, T] and which satisfies ().

Remark.

The sufficient condition provided in Corollary 3.5.1 guarantees uniform convergence of the Type-I finite memory approximation. If we replace condition () with

α_{j} < C, for all 1 \leq j < + \infty,

where C is a positive constant (independent on T), then we obtain asymptotic convergence. In other words, for each δ > 0, there exists an integer p such that for all n > p, we have

∥w_{0} (t) - w_{0}^{n} (t)∥ < δ

⁠. This result is based on the limit

lim_{p \to + \infty} \frac{δ (p + 1)!}{C_{1} A_{1} A_{2} (\prod_{j = 1}^{p} α_{j})} > lim_{p \to + \infty} \frac{δ (p + 1)!}{C_{1} A_{1} A_{2} C^{p}} = + \infty

which guarantees the existence of an integer p for which the upper bound on Δt_p is smaller or equal to zero. In such a case, the Type I FMA degenerates to the H-model, for which Corollary 3.4.2 holds.

3. Type-II finite memory approximation

The Type-II finite memory approximation is obtained by solving the system (31) with $w_{n}^{e_{n}} (t)$ given in (35). We first derive an upper bound for $‖ w_{0} (t) - w_{0}^{n} (t) ‖$ and then discuss sufficient conditions for convergence.

Theorem 3.6

(Accuracy of the Type-II FMA). Let

e^{t L}

and

e^{t L Q}

be strongly continuous semigroups with upper bounds

‖ e^{t L} ‖ \leq M e^{t ω}

and

‖ e^{t L Q} ‖ \leq M_{Q} e^{t ω_{Q}}

. If

α_{j} = \frac{‖ {(L Q)}^{j + 1} L u_{0} ‖}{‖ {(L Q)}^{j} L u_{0} ‖}, 1 \leq j \leq n,

(47)

then for 1 ≤ p ≤ n,

‖ w_{0} (t) - w_{0}^{p} (t) ‖ \leq M_{5}^{p} (t),

where

M_{5}^{p} (t) = C_{1} (\prod_{j = 1}^{p} α_{j}) f_{p} (ω_{Q}, t) h (ω - ω_{Q}, t_{p}),

f_{p} (ω_{Q}, t) = \int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} e^{σ ω_{Q}} d σ, h (ω - ω_{Q}, t_{p}) = \int_{0}^{t_{p}} e^{s (ω - ω_{Q})} d s,

and

C_{1} = M M_{Q} ‖ P ‖^{2} ‖ {(L Q)}^{p + 1} L u_{0} ‖

.

Proof.

By following the same procedure as in the Proof of the Theorem 3.4, we obtain

w_{0} (t) - w_{0}^{p} (t) = \int_{0}^{t} \int_{0}^{τ_{p}} \dots \int_{0}^{τ_{2}} \int_{0}^{min (τ_{1}, t_{p})} P e^{s L} P e^{(τ_{1} - s) L Q} {(L Q)}^{p + 1} L u_{0} d s d τ_{1} \dots d τ_{p} .

By applying Cauchy’s formula for repeated integration, this expression may be simplified to

w_{0} (t) - w_{0}^{p} (t) = \int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} \int_{0}^{min (σ, t_{p})} P e^{s L} P e^{(σ - s) L Q} {(L Q)}^{p + 1} L u_{0} d s d σ .

Thus,

\begin{matrix} ‖ w_{0} (t) - w_{0}^{p} (t) ‖ & \leq & \int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} \int_{0}^{min (σ, t_{p})} ∥P e^{s L} P e^{(σ - s) L Q} {(L Q)}^{p + 1} L u_{0}∥ d s d σ \\ \leq & M M_{Q} ‖ P ‖^{2} ‖ {(L Q)}^{p + 1} L u_{0} ‖ \int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} \int_{0}^{t_{p}} e^{s ω} e^{(σ - s) ω_{Q}} d s d σ \\ \leq & C_{1} (\prod_{j = 1}^{p} α_{j}) (\int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} e^{σ ω_{Q}} d σ) (\int_{0}^{t_{p}} e^{s (ω - ω_{Q})} d s) \\ = & C_{1} (\prod_{j = 1}^{p} α_{j}) f_{p} (ω_{Q}, t) h (ω - ω_{Q}, t_{p}) = M_{5}^{p} (t), \end{matrix}

where

C_{1} = M M_{Q} ‖ P ‖^{2} ‖ {(L Q)}^{p + 1} L u_{0} ‖

⁠,

f_{p} (ω_{Q}, t) = \int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} e^{σ ω_{Q}} d σ = \{\begin{matrix} \frac{t^{p}}{p!}, & ω_{Q} = 0, \\ \frac{1}{ω_{Q}^{p}} [e^{t ω_{Q}} - \sum_{k = 0}^{p - 1} \frac{{(t ω_{Q})}^{k}}{k!}], & ω_{Q} \neq 0, \end{matrix}

(48)

and

h (ω - ω_{Q}, t_{p}) : = \int_{0}^{t_{p}} e^{s (ω - ω_{Q})} d s = \{\begin{matrix} t_{p}, & ω = ω_{Q}, \\ \frac{e^{t_{p} (ω - ω_{Q})} - 1}{ω - ω_{Q}}, & ω \neq ω_{Q} \end{matrix}

are both strictly increasing functions of t and t_p, respectively.

Corollary 3.6.1

(Uniform convergence of the Type-II FMA). If α_j in Theorem 3.6 satisfy

α_{j} < \frac{j}{T} (ω_{Q} = 0) o r α_{j} < ω_{Q} \frac{e^{T ω_{Q}} - \sum_{k = 0}^{j - 2} \frac{{(T ω_{Q})}^{k}}{k!}}{e^{T ω_{Q}} - \sum_{k = 0}^{j - 1} \frac{{(T ω_{Q})}^{k}}{k!}} (ω_{Q} \neq 0)

(49)

for 1 ≤ j ≤ n, then for any arbitrarily small δ > 0, there exists an ordered sequence 0 < t₀ < t₁ < ⋯ < t_n ≤ T such that

‖ w_{0} (T) - w_{0}^{p} (T) ‖ \leq δ, 1 \leq p \leq n,

which satisfies

t_{j} \geq \{\begin{matrix} \frac{j! δ}{C_{1} (\prod_{i = 1}^{j} α_{i}) T^{j}}, & ω_{Q} = 0, \\ \frac{ω_{Q}^{j} δ}{C_{1} (\prod_{i = 1}^{j} α_{i}) [e^{T ω_{Q}} - \sum_{k = 0}^{j - 1} \frac{{(T ω_{Q})}^{k}}{k!}]}, & ω_{Q} \neq 0, \end{matrix}

when

ω = ω_{Q}

and

t_{j} \geq \{\begin{matrix} \frac{1}{ω} \ln [1 + \frac{j! ω δ}{C_{1} (\prod_{i = 1}^{j} α_{i}) T^{j}}], & ω_{Q} = 0, \\ \frac{1}{ω - ω_{Q}} \ln [1 + \frac{(ω - ω_{Q}) ω_{Q}^{j} δ}{C_{1} (\prod_{i = 1}^{j} α_{i}) [e^{T ω_{Q}} - \sum_{k = 0}^{j - 1} \frac{{(T ω_{Q})}^{k}}{k!}]}], & ω_{Q} \neq 0, \end{matrix}

when

ω \neq ω_{Q}

.

Proof.

We now consider separately the two cases where

ω = ω_{Q}

and where

ω \neq ω_{Q}

⁠. If

ω = ω_{Q}

⁠, then

\begin{matrix} ‖ w_{0} (t) - w_{0}^{p} (t) ‖ & \leq & C_{1} (\prod_{i = 1}^{p} α_{i}) \int_{0}^{t} \int_{0}^{τ_{p}} \dots \int_{0}^{τ_{2}} \int_{0}^{t_{p}} e^{τ_{1} ω_{Q}} e^{s (ω - ω_{Q})} d s d τ_{1} \dots d τ_{p} \\ = & t_{p} C_{1} (\prod_{i = 1}^{p} α_{i}) \int_{0}^{t} \int_{0}^{τ_{p}} \dots \int_{0}^{τ_{2}} e^{τ_{1} ω_{Q}} d τ_{1} \dots d τ_{p} \\ = & t_{p} C_{1} (\prod_{i = 1}^{p} α_{i}) f_{p} (ω_{Q}, t), \end{matrix}

where

f_{p} (ω_{Q}, t)

is defined in (). To ensure that

‖ w_{0} (t) - w_{0}^{p} (t) ‖ \leq δ

for all 0 ≤ t ≤ T, we can take

t_{p} C_{1} (\prod_{i = 1}^{p} α_{i}) f_{p} (ω_{Q}, T) = max_{t \in [0, T]} t_{p} C_{1} (\prod_{i = 1}^{p} α_{i}) f_{p} (ω_{Q}, t) \leq δ

so that

t_{p} \leq \frac{δ}{C_{1} (\prod_{i = 1}^{p} α_{i}) f_{p} (ω_{Q}, T)} = \{\begin{matrix} \frac{p! δ}{C_{1} (\prod_{i = 1}^{p} α_{i}) T^{p}}, & ω_{Q} = 0, \\ \frac{ω_{Q}^{p} δ}{C_{1} (\prod_{i = 1}^{p} α_{i}) [e^{T ω_{Q}} - \sum_{k = 0}^{p - 1} \frac{{(T ω_{Q})}^{k}}{k!}]}, & ω_{Q} \neq 0 . \end{matrix}

On the other hand, if

ω \neq ω_{Q}

⁠, then

\begin{matrix} ‖ w_{0} (t) - w_{0}^{p} (t) ‖ & \leq & C_{1} (\prod_{i = 1}^{p} α_{i}) \int_{0}^{t} \int_{0}^{τ_{p}} \dots \int_{0}^{τ_{2}} \int_{0}^{t_{p}} e^{τ_{1} ω_{Q}} e^{s (ω - ω_{Q})} d s d τ_{1} \dots d τ_{p} \\ = & \frac{e^{t_{p} (ω - ω_{Q})} - 1}{ω - ω_{Q}} C_{1} (\prod_{i = 1}^{p} α_{i}) \int_{0}^{t} \int_{0}^{τ_{p}} \dots \int_{0}^{τ_{2}} e^{τ_{1} ω_{Q}} d τ_{1} \dots d τ_{p} \\ = & \frac{e^{t_{p} (ω - ω_{Q})} - 1}{ω - ω_{Q}} C_{1} (\prod_{i = 1}^{p} α_{i}) f_{p} (ω_{Q}, t) . \end{matrix}

To ensure that

‖ w_{0} (t) - w_{0}^{p} (t) ‖ \leq δ

for all 0 ≤ t ≤ T, we can take

\frac{e^{t_{p} (ω - ω_{Q})} - 1}{ω - ω_{Q}} C_{1} (\prod_{i = 1}^{p} α_{i}) f_{p} (ω_{Q}, T) = max_{t \in [0, T]} \frac{e^{t_{p} (ω - ω_{Q})} - 1}{ω - ω_{Q}} C_{1} (\prod_{i = 1}^{p} α_{i}) f_{p} (ω_{Q}, t) \leq δ .

Let us now consider the two cases

ω > ω_{Q}

and

ω < ω_{Q}

separately. When

ω > ω_{Q}

⁠, we have

e^{t_{p} (ω - ω_{Q})} \leq 1 + \frac{(ω - ω_{Q}) δ}{C_{1} (\prod_{i = 1}^{p} α_{i}) f_{p} (ω_{Q}, T)}

and

t_{p} \leq \{\begin{matrix} \frac{1}{ω} \ln [1 + \frac{p! ω δ}{C_{1} (\prod_{i = 1}^{p} α_{i}) T^{p}}], & ω_{Q} = 0, \\ \frac{1}{ω - ω_{Q}} \ln [1 + \frac{(ω - ω_{Q}) ω_{Q}^{p} δ}{C_{1} (\prod_{i = 1}^{p} α_{i}) [e^{T ω_{Q}} - \sum_{k = 0}^{p - 1} \frac{{(T ω_{Q})}^{k}}{k!}]}], & ω_{Q} \neq 0 . \end{matrix}

On the other hand, when

ω < ω_{Q}

⁠, we have

\frac{1 - e^{- t_{p} (ω_{Q} - ω)}}{ω_{Q} - ω} C_{1} (\prod_{i = 1}^{p} α_{i}) f_{p} (ω_{Q}, T) = max_{t \in [0, T]} \frac{1 - e^{- t_{p} (ω_{Q} - ω)}}{ω_{Q} - ω} C_{1} (\prod_{i = 1}^{p} α_{i}) f_{p} (ω_{Q}, t) \leq δ

so that

e^{- t_{p} (ω_{Q} - ω)} \geq 1 - \frac{(ω_{Q} - ω) δ}{C_{1} (\prod_{i = 1}^{p} α_{i}) f_{p} (ω_{Q}, T)},

i.e.,

- t_{p} (ω_{Q} - ω) \geq \ln [1 - \frac{(ω_{Q} - ω) δ}{C_{1} (\prod_{i = 1}^{p} α_{i}) f_{p} (ω_{Q}, T)}] .

Hence,

t_{p} \leq \{\begin{matrix} \frac{1}{ω} \ln [1 - \frac{p! (- ω) δ}{C_{1} (\prod_{i = 1}^{p} α_{i}) T^{p}}], & ω_{Q} = 0, \\ - \frac{1}{ω_{Q} - ω} \ln [1 - \frac{(ω_{Q} - ω) ω_{Q}^{p} δ}{C_{1} (\prod_{i = 1}^{p} α_{i}) [e^{T ω_{Q}} - \sum_{k = 0}^{p - 1} \frac{{(T ω_{Q})}^{k}}{k!}]}], & ω_{Q} \neq 0 . \end{matrix}

For all the four cases, if

ω_{Q} = 0

⁠, then we have condition α_p < p/T, and the upper bound of the time sequence satisfies

\frac{p! δ}{C_{1} (\prod_{i = 1}^{p} α_{i}) T^{p}} < \frac{(p - 1)! δ}{C_{1} (\prod_{i = 1}^{p - 1} α_{i}) T^{p - 1}}, p \geq 2 .

If

ω_{Q} \neq 0

⁠, then we have the condition

α_{p} < ω_{Q} \frac{e^{T ω_{Q}} - \sum_{k = 0}^{p - 2} \frac{{(T ω_{Q})}^{k}}{k!}}{e^{T ω_{Q}} - \sum_{k = 0}^{p - 1} \frac{{(T ω_{Q})}^{k}}{k!}}

and the upper bound of the time sequence satisfies

\frac{ω_{Q}^{p - 1}}{(e^{T ω_{Q}} - \sum_{k = 0}^{p - 2} \frac{{(T ω_{Q})}^{k}}{k!}) \prod_{i = 1}^{p - 1} α_{i}} < \frac{ω_{Q}^{p}}{(e^{T ω_{Q}} - \sum_{k = 0}^{p - 1} \frac{{(T ω_{Q})}^{k}}{k!}) \prod_{i = 1}^{p} α_{i}}, p \geq 2 .

Therefore, there always exists a increasing time sequence 0 < t₁ < ⋯ < t_n such that

‖ w_{0} (t) - w_{0}^{p} (t) ‖ \leq δ

for all 0 ≤ t ≤ T. And since we have proved that this δ-bound on the error holds for all t_n upper bounded as in the two cases above, there exists such an increasing time sequence 0 < t₁ < ⋯ < t_n with t_n lower-bounded by the same quantities. Indeed, because of the coarseness of the approximations applied in the proof, there may exist such a time sequence with significantly larger t_i.

Remark.

If we replace () with the stronger condition

\{\begin{matrix} α_{j} < \frac{j}{ϵ T}, & ω_{Q} = 0, \\ α_{j} < \frac{1}{ϵ T}, & ω_{Q} \neq 0, \end{matrix} 1 \leq j < \infty,

(50)

where ϵ is some arbitrary constant satisfying ϵ > 1, then we have

lim_{j \to + \infty} t_{j} \geq lim_{j \to + \infty} \{\begin{matrix} \frac{j! δ}{C_{1} (\prod_{i = 1}^{j} α_{i}) T^{j}} \geq \frac{δ}{C_{1}} ϵ^{j} = + \infty, & ω_{Q} = 0, \\ \frac{ω_{Q}^{j} δ}{C_{1} (\prod_{i = 1}^{j} α_{i}) [e^{T ω_{Q}} - \sum_{k = 0}^{j - 1} \frac{{(T ω_{Q})}^{k}}{k!}]} = \frac{ω_{Q}^{j} δ}{C_{1} (\prod_{i = 1}^{j} α_{i}) o (T^{j} ω_{Q}^{j})} = + \infty, & ω_{Q} \neq 0 \end{matrix}

for

ω = ω_{Q}

and

lim_{j \to + \infty} t_{j} \geq lim_{j \to + \infty} \{\begin{matrix} \frac{j}{ω} \ln [\frac{δ}{C_{1}} ω ϵ] = + \infty, & ω_{Q} = 0, \\ \frac{j}{ω - ω_{Q}} \ln [\frac{(ω - ω_{Q}) T^{j}}{C_{1} o (T^{j})}] = + \infty, & ω_{Q} \neq 0 . \end{matrix}

Hence, there exists a j such that the upper bound for t_j is greater than or equal to T. For such a case, the Type II FMA degenerates to the truncation approximation (H-model), for which Corollary 3.4.2 grants us asymptotic convergence.

4. H_t-model

The H_t-model is obtained by solving the system (31) with $w_{n}^{e_{n}} (t)$ approximated using Chorin’s t-model¹² [see Eq. (36)]. Convergence analysis can be performed by using the mathematical methods we employed for the proofs of the H-model. Note that the classical t-model is equivalent to a zeroth-order H_t-model.

Theorem 3.7

(Accuracy of the H_t-model). Let

e^{t L}

and

e^{t L Q}

be strongly continuous semigroups with upper bounds

‖ e^{t L} ‖ \leq M e^{t ω}

and

‖ e^{t L Q} ‖ \leq M_{Q} e^{t ω_{Q}}

, and let T > 0 be a fixed integration time. For some fixed n, let

α_{j} = \frac{‖ {(L Q)}^{j + 1} L u_{0} ‖}{‖ {(L Q)}^{j} L u_{0} ‖}, 1 \leq j \leq n .

(51)

Then, for any 1 ≤ p ≤ n and all t ∈ [0, T], we have

‖ w_{0} (t) - w_{0}^{p} (t) ‖ \leq M_{6}^{p} (t) \leq M_{6}^{p} (T),

where

M_{6}^{p} (t) = C_{4} (\prod_{j}^{p} α_{j}) \frac{t^{p + 1}}{(p + 1)!}, C_{4} = [C_{1} A_{1} A_{2} + \frac{C_{1}}{M_{Q} A_{3}}], A_{3} = max_{s \in [0, T]} s e^{s ω} = \{\begin{matrix} 1, & ω \leq 0, \\ e^{T ω}, & ω > 0 \end{matrix}

and C₁, A₁, A₂ are the same as before.

Proof.

For pth order H_t-model, the difference between the memory term

w

₀ and its approximation

w_{0}^{p}

is

\begin{matrix} w_{0} (t) - w_{0}^{p} (t) & = & \int_{0}^{t} \int_{0}^{τ_{p}} \dots \int_{0}^{τ_{2}} [\int_{0}^{τ_{1}} P e^{s L} P e^{(τ_{1} - s) L Q} {(L Q)}^{p + 1} L u_{0} d s - τ_{1} P e^{τ_{1} L} P {(L Q)}^{p + 1} L u_{0}] \\ \times d τ_{1} \dots d τ_{p} . \end{matrix}

(52)

Using Cauchy’s formula for repeated integration, we can bind the norm of the second term in () as

\begin{matrix} ∥\int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} σ P e^{σ L} P {(L Q)}^{p + 1} L u_{0} d σ∥ & \leq & \int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} ‖ σ P e^{σ L} P {(L Q)}^{p + 1} L u_{0} ‖ d σ \\ \leq & ‖ P ‖^{2} M ‖ {(L Q)}^{p + 1} L u_{0} ‖ \underset{g_{p} (t, ω)}{\underset{︸}{\int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} σ e^{σ ω} d σ}} \\ = & \frac{C_{1}}{M_{Q}} (\prod_{j = 1}^{p} α_{j}) g_{p} (t, ω), \end{matrix}

(53)

where

C_{1} = ‖ P ‖^{2} ‖ L Q L u_{0} ‖ M M_{Q}

as before. The function g_p(t, ω) may be bounded from above as

g_{p} (t, ω) \leq A_{3} \int_{0}^{t} \frac{{(t - σ)}^{p - 1}}{(p - 1)!} σ d σ = A_{3} \frac{t^{p + 1}}{(p + 1)!}, A_{3} = max_{s \in [0, T]} e^{s ω} = \{\begin{matrix} 1, & ω \leq 0 \\ e^{T ω}, & ω > 0 . \end{matrix}

By applying the triangle inequality to () and taking () into account, we obtain

‖ w_{0} (t) - w_{0}^{p} (t) ‖ \leq C_{1} A_{1} A_{2} (\prod_{j = 1}^{p} α_{j}) \frac{t^{p + 1}}{(p + 1)!} + \frac{C_{1}}{M_{Q}} A_{3} (\prod_{j = 1}^{p} α_{j}) \frac{t^{p + 1}}{(p + 1)!} = M_{6}^{p} (t) .

One can see that the upper bounds $M_{6}^{p} (t)$ and $M_{3}^{p} (t)$ (see Theorem 3.4) share the same structure, the only difference being the constant out front. Hence by changing C₂ to C₄, we can prove a series of corollaries similar to 3.4.1–3.4.3. In summary, what holds for the H-model also holds for the H_t-model. For the sake of brevity, we omit the statement and proofs of those corollaries.

E. Linear dynamical systems

The upper bounds we obtained above are not easily computable for general nonlinear systems and infinite-rank projections, e.g., Chorin’s projection (7). However, if the dynamical system is linear, then such upper bounds are explicitly computable and the convergence of the H-model can be established for linear phase space functions in any finite integration time T. To this end, consider the linear system ẋ = Ax with a random initial condition x(0) sampled from the joint probability density function

ρ_{0} (x_{0}) = δ (x_{01} - x_{1} (0)) \prod_{j = 2}^{N} ρ_{0 j} (x_{0 j}) .

(54)

In other words, the initial condition for the quantity of interest u(x) = x₁(t) is set to be deterministic, while all other variables x₂, …, x_N are zero-mean and statistically independent at t = 0. Here we assume for simplicity that ρ_0j (j = 2, …, N) are independent identically distributed (i.i.d.) standard normal distributions. Observe that the Liouville operator associated with the linear system ẋ = Ax is

L = \sum_{i = 1}^{N} \sum_{j = 1}^{N} A_{i j} x_{j} \frac{\partial}{\partial x_{i}},

(55)

where A_ij are the entries of the matrix A. If we choose observable u = x₁(t), then Chorin’s projection operator (9) yields the evolution equation for the conditional expectation $E [x_{1} (t) | x_{1} (0)]$ ⁠, i.e., the conditional mean path (11), which can be explicitly written as

\frac{d}{d t} E [x_{1} | x_{1} (0)] = A_{11} E [x_{1} | x_{1} (0)] + w_{0} (t),

(56)

where $A_{11} = P L x_{1} (0)$ is the first entry of the matrix A and $w$ ₀ represents the memory integral (26). Next, we explicitly compute the upper bounds for the memory growth and the error in the H-model for this system. To this end, we first notice that the domain of the Liouville operator can be restricted to the linear space

V = span {x_{1}, \dots, x_{N}} .

(57)

In fact, V is invariant under $L$ ⁠, $P$ ⁠, and $Q$ ⁠, i.e., $L V \subseteq V$ ⁠, $P V \subseteq V$ ⁠, and $Q V \subseteq V$ ⁠. These operators have the following matrix representations:

L ≃ A^{T} ≃ [\begin{matrix} a_{11} & b^{T} \\ a & M_{11}^{T} \end{matrix}], P ≃ [\begin{matrix} 1 & 0 & \dots & 0 \\ 0 & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & 0 \end{matrix}], Q ≃ [\begin{matrix} 0 & 0 & \dots & 0 \\ 0 & 1 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & 1 \end{matrix}],

where M₁₁ is the minor of the matrix of A obtained by removing the first column and the first row, while

a = {[A_{12} \dots A_{1 N}]}^{T}, b^{T} = [A_{21} \dots A_{N 1}] .

(58)

Therefore,

L Q ≃ [\begin{matrix} 0 & b^{T} \\ 0 & M_{11}^{T} \end{matrix}], L {(Q L)}^{n} x_{1} (0) ≃ [\begin{matrix} b^{T} {(M_{11}^{T})}^{n - 1} a \\ {(M_{11}^{T})}^{n} a \end{matrix}] .

(59)

At this point, we set x₀₁ = x₁(0) and

q (t, x_{01}, {\tilde{x}}_{0}) = \int_{0}^{t} e^{s L} P L e^{(t - s) Q L} Q L x_{01} d s .

Since ${\tilde{x}}_{0} = (x_{2} (0), \dots, x_{N} (0))$ is random, $q (t, x_{01}, {\tilde{x}}_{0})$ is a random variable. By using Jensen’s inequality ${[E (X)]}^{2} \leq E [X^{2}]$ ⁠, we have the following L^∞ estimate:

‖ (P q) (t, x_{01}) ‖_{L^{\infty}} \leq ‖ q (t, x_{01}, \cdot) ‖_{L_{ρ_{0}}^{2}} .

(60)

On the other hand, we have

‖ e^{t L} ‖_{L_{ρ_{0}}^{2} (V)} \leq ‖ e^{t L} ‖_{L_{ρ_{0}}^{2}} \leq e^{t ω}, ω = - \frac{1}{2} inf {div}_{ρ_{0}} F .

(61)

For linear dynamical systems, both $‖ \cdot ‖_{L_{ρ_{0}}^{2} (V)}$ and $‖ \cdot ‖_{L_{ρ_{0}}^{2}}$ upper bounds can be used to estimate the norm of the semigroup $e^{t L}$ ⁠. However, for the semigroup $e^{t L Q}$ ⁠, we can only obtain the explicit form of the $‖ \cdot ‖_{L_{ρ_{0}}^{2} (V)}$ bound, which is given by the following perturbation theorem¹⁶ (see also Appendix A):

‖ e^{t L Q} ‖_{L_{ρ_{0}}^{2} (V)} \leq e^{t ω_{Q}}, where ω_{Q} = ω + \sqrt{A_{11}^{2} + \sum_{i = 2}^{N} A_{1 i}^{2} \frac{{⟨ x_{i}^{2} (0) ⟩}_{ρ_{0}}}{x_{1}^{2} (0)}} \geq ω + ‖ L P ‖_{L_{ρ_{0}}^{2} (V)} .

(62)

Memory growth. It is straightforward at this point to compute the upper bound of the memory growth we obtained in Theorem 3.1. Since $‖ P ‖_{L_{ρ_{0}}^{2}} = ‖ Q ‖_{L_{ρ_{0}}^{2}} = 1$ (⁠ $P$ and $Q$ are the orthogonal projections relative to ρ₀), we have the following result:

| w_{0} (t) | \leq ‖ L Q L x_{1} (0) ‖ \frac{e^{t ω} - e^{t ω_{Q}}}{ω - ω_{Q}} = \sqrt{{(b^{T} a)}^{2} x_{1}^{2} (0) + {∥Λ_{x_{i + 1} (0)} M_{11}^{T} a∥}_{2}^{2}} \frac{e^{t ω} - e^{t ω_{Q}}}{ω - ω_{Q}},

(63)

where $Λ_{x_{i + 1} (0)}$ is a N − 1 × N − 1 diagonal matrix with $Λ_{i i} = {⟨ x_{i + 1} (0) ⟩}_{ρ_{0}}$ and ∥·∥₂ is the vector 2-norm. Accuracy of theH-model. We are interested in computing the upper bound of the approximation error generated by the H-model (see Sec. III D 1–Theorem 3.4). By using the matrix representation of $L$ ⁠, $P$ ⁠, and $Q$ ⁠, the nth order H-model MZ Eq. (56) for the linear system can be explicitly written as

\{\begin{matrix} \frac{d}{d t} E [x_{1} | x_{1} (0)] = A_{11} E [x_{1} | x_{1} (0)] + w_{0}^{n} (t) (MZ equation), \\ \frac{d w_{j}^{n} (t)}{d t} = b^{T} {(M_{11}^{T})}^{j} a^{T} E [x_{1} | x_{1} (0)] + w_{j + 1}^{n} (t), j = 0,1, \dots, n - 1, \\ \frac{d w_{n}^{n} (t)}{d t} = b^{T} {(M_{11}^{T})}^{n} a^{T} E [x_{1} | x_{1} (0)], \end{matrix}

(64)

where M₁₁, a, and b are defined as before [see Eq. (58)]. The upper bound for the memory term approximation error is explicitly obtained as

\begin{matrix} | w_{0} (t) - w_{0}^{n} (t) | & \leq & A_{1} A_{2} ‖ L {(Q L)}^{n} x_{1} (0) ‖ \frac{t^{n + 1}}{(n + 1)!} \\ = & A_{1} A_{2} \sqrt{{[b^{T} {(M_{11}^{T})}^{n} a]}^{2} x_{1}^{2} (0) + {∥Λ_{x_{i + 1} (0)} {(M_{11}^{T})}^{n + 1} a∥}_{2}^{2}} \frac{t^{n + 1}}{(n + 1)!}, \end{matrix}

(65)

where A₁, A₂ are defined in (38), while ω and $ω_{Q}$ are given in (61) and (62), respectively. Note that the error bound (65) is slightly different from the one we obtained in Theorem 3.4. The reason is that here we choose to bound $‖ L {(Q L)}^{n} u_{0} ‖$ ⁠, instead of the quotient $α_{n} = | L {(Q L)}^{n + 1} u_{0} ‖ / ‖ L {(Q L)}^{n} u_{0} ‖$ ⁠. For each fixed integration time T, the upper bound (65) goes to zero as we send n to infinity, i.e.,

lim_{n \to + \infty} | w_{0} (T) - w_{0}^{n} (T) | = 0 .

This means that the H-model converges for all linear dynamical systems with observables in the linear space (57).

F. Memory estimates for finite-rank projections and Hamiltonian systems

The semigroup estimates we obtained in Sec. III A allow us to compute explicitly an a priori estimate of the memory kernel in the Mori-Zwanzig equation if we employ finite-rank projections. In this section, we outline the procedure to obtain such an estimate for Hamiltonian dynamical systems. We begin by recalling that, in general, Hamiltonian systems are necessarily divergence-free, i.e.,

{div}_{ρ_{e q}} (F) = 0 .

(66)

Here, F(x) is the velocity field at the right-hand side of (1), while $ρ_{e q} = e^{- β H} / Z$ is the canonical Gibbs distribution. Equation (66) can be easily obtained by noticing that

\begin{matrix} {div}_{ρ_{e q}} (F) & = e^{β H} \nabla \cdot (e^{- β H} F), \\ = e^{β H} \sum_{i = 1}^{N} (\frac{\partial}{\partial q_{i}} [e^{- β H} \frac{\partial H}{\partial p_{i}}] - \frac{\partial}{\partial p_{i}} [e^{- β H} \frac{\partial H}{\partial q_{i}}]), \\ = 0 . \end{matrix}

By combining the estimate (18) with the divergence-free constraint (66), we find that the Koopman semigroup of a Hamiltonian dynamical system is a contraction in the $L_{ρ_{e q}}^{2}$ norm, i.e.,

{∥e^{t L}∥}_{L_{ρ_{e q}}^{2}} \leq 1 .

(67)

Moreover, the MZ equation (6) with a finite-rank projection $P$ of the form (12) (with σ = ρ_eq) can be reduced to the following Volterra integro-differential equation:

\frac{d}{d t} P u_{i} (t) = \sum_{j = 1}^{M} Ω_{i j} P u_{j} (t) - \sum_{j = 1}^{M} \int_{0}^{t} K_{i j} (t - s) P u_{j} (s) d s, i = 1, \dots, M,

(68)

where

G_{i j} = {⟨ u_{i}, u_{j} ⟩}_{e q},

(69a)

Ω_{i j} = \sum_{k = 1}^{M} {(G^{- 1})}_{j k} {⟨ u_{k}, L u_{i} ⟩}_{e q},

(69b)

K_{i j} (t - s) = - \sum_{k = 1}^{M} {(G^{- 1})}_{j k} {⟨ Q L u_{k}, e^{(t - s) Q L} Q L u_{i} ⟩}_{e q} .

(69c)

Equation (68) is often called the generalized Langevin equation (GLE)^13,31 for the projected quantity of interest u_i(t). To derive (69a)–(69c), we used the fact that $L$ is skew-adjoint and $Q$ is self-adjoint with respect to the $L_{ρ_{e q}}^{2}$ inner product and that $Q^{2} = Q$ ⁠. In classical statistical physics, Eq. (69c) is known as the second fluctuation dissipation theorem. Next, define the temporal-correlation matrix

C_{i j} (t) = {⟨ u_{j} (0), u_{i} (t) ⟩}_{e q} = {⟨ u_{j} (0), P u_{i} (t) ⟩}_{e q} .

(70)

By applying ${⟨ u_{j}, (\cdot) ⟩}_{e q}$ to both sides of Eq. (68), we obtain the following exact evolution equation for C_ij(t):

\frac{d C_{i j}}{d t} = \sum_{k = 1}^{M} Ω_{i k} C_{k j} - \sum_{k = 1}^{M} \int_{0}^{t} K_{i k} (t - s) C_{k j} (s) d s .

(71)

Moreover, if we employ a one-dimensional Mori’s basis, i.e., M = 1, then we obtain the simplified equation

\frac{d C (t)}{d t} = Ω_{1} C (t) - \int_{0}^{t} K (t - s) C (s) d s,

(72)

where $C (t) = {⟨ u_{1} (0), u_{1} (t) ⟩}_{e q}$ ⁠. The main difficulty in solving the GLE (71) [or (72)] lies in computing the memory kernels K_ij(t). Hereafter we prove that such memory kernels can be uniformly bounded by a computable quantity that depends only on the initial condition of the system. For the sake of simplicity, we will focus on the one-dimensional GLE (72).

Theorem 3.8.

The memory kernel K(t) in the one-dimensional GLE () is uniformly bounded by

‖ {\dot{u}}_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2} / ‖ u_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}

, i.e.,

| K (t) | \leq \frac{‖ {\dot{u}}_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}}{‖ u_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}}, \forall t \geq 0 .

(73)

Proof.

From the second-fluctuation dissipation theorem (), the memory kernel K(t) satisfies

| K (t) | = |\frac{{⟨ e^{t Q L} Q L u_{1} (0), Q L u_{1} (0) ⟩}_{e q}}{{⟨ u_{1} (0), u_{1} (0) ⟩}_{e q}}| \leq ‖ e^{t Q L} Q ‖_{L_{ρ_{e q}}^{2}} \frac{‖ L u_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}}{‖ u_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}} = ‖ e^{t Q L} Q ‖_{L_{ρ_{e q}}^{2}} \frac{‖ {\dot{u}}_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}}{‖ u_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}} .

On the other hand, by using the numerical abscissa () and formula (), we see that the semigroup

e^{t Q L Q}

is contractive, i.e.,

‖ e^{t Q L Q} ‖_{L_{ρ_{e q}}^{2}} \leq 1

⁠. Since

Q

is an orthogonal projection with respect to ρ_eq, we have

‖ e^{t Q L} Q ‖_{L_{ρ_{e q}}^{2}} = ‖ Q e^{t Q L Q} ‖_{L_{ρ_{e q}}^{2}} \leq ‖ Q ‖_{L_{ρ_{e q}}^{2}} ‖ e^{t Q L Q} ‖_{L_{ρ_{e q}}^{2}} \leq 1

⁠. This yields

| K (t) | \leq ‖ e^{t Q L} Q ‖_{L_{ρ_{e q}}^{2}} \frac{‖ {\dot{u}}_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}}{‖ u_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}} \leq \frac{‖ \dot{u_{1}} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}}{‖ u_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}} .

Theorem 3.8 provides an a priori and easily computable upper bound for the memory kernel defining the dynamics of any quantity of interest u₁ that is initially in the Gibbs canonical ensemble $ρ_{e q} = e^{- β H} / Z$ ⁠. In Sec. IV, we will calculate the upper bound (73) analytically and compare it with the exact memory kernel we obtain in prototype linear and nonlinear Hamiltonian systems.

Remark.

We emphasized in Sec. III A that the semigroup estimate for

e^{t Q L Q}

is not necessarily tight. In the context of high-dimensional Hamiltonian systems (e.g., molecular dynamics), it is often empirically assumed that the semigroup

e^{t Q L} Q

is dissipative, i.e.,

‖ e^{t Q L Q} ‖ \leq e^{t ω_{0}}

⁠, where ω₀ < 0. In this case, the memory kernel turns out to be uniformly bounded by an exponentially decaying function since

| K (t) | \leq ‖ e^{t Q L} Q ‖_{L_{ρ_{e q}}^{2}} \frac{‖ L u_{1} ‖_{L_{ρ_{e q}}^{2}}^{2}}{‖ u_{1} ‖_{L_{ρ_{e q}}^{2}}^{2}} \leq e^{t ω_{0}} \frac{‖ {\dot{u}}_{1} ‖_{L_{ρ_{e q}}^{2}}^{2}}{‖ u_{1} ‖_{L_{ρ_{e q}}^{2}}^{2}} .

IV. NUMERICAL EXAMPLES

In this section, we provide simple numerical examples of the MZ memory approximation methods we discussed throughout the paper. Specifically, we study Hamiltonian systems (linear and nonlinear) with finite-rank projections (Mori’s projection) and non-Hamiltonian systems with infinite-rank projections (Chorin’s projection). In both cases, we demonstrate the accuracy of the a priori memory estimation method we developed in Secs. III F and III E. We also compute the solution to the MZ equation for non-Hamiltonian systems with the t-model, the H-model, and the H_t-model.

A. Hamiltonian dynamical systems with finite-rank projections

In this section, we consider dimension reduction in linear and nonlinear Hamiltonian dynamical systems with finite-rank projection. In particular, we consider the Mori projection and study the MZ equation for the temporal auto-correlation function of a scalar quantity of interest.

1. Harmonic chains of oscillators

Consider a one-dimensional chain of harmonic oscillators. This is a simple but illustrative example of a linear Hamiltonian dynamical system which has been widely studied in statistical mechanics, mostly in relation with the microscopic theory of Brownian motion.^3,18,17 The Hamiltonian of the system can be written as

H (p, q) = \frac{1}{2 m} \sum_{i = 1}^{N} p_{i}^{2} + \frac{k}{2} \sum_{\begin{array}{c} i, j = 0 \\ i < j \end{array}}^{N + 1} {(q_{i} - q_{j})}^{2},

(74)

where q_i and p_i are, respectively, the displacement and momentum of the ith particle, m is the mass of the particles (assumed constant throughout the network), and k is the elasticity constant that modulates the intensity of the quadratic interactions. We set fixed boundary conditions at the endpoints of the chain, i.e., q₀(t) = q_N+1(t) = 0 and p₀(t) = p_N+1(t) = 0 (particles are numbered from left to right) and m = k = 1. The Hamilton’s equations are

\frac{d q_{i}}{d t} = \frac{\partial H}{\partial p_{i}}, \frac{d p_{i}}{d t} = - \frac{\partial H}{\partial q_{i}},

(75)

which can be written in a matrix-vector form as

[\begin{matrix} \dot{p} \\ \dot{q} \end{matrix}] = [\begin{matrix} 0 & k B - k D \\ I / m & 0 \end{matrix}] [\begin{matrix} p \\ q \end{matrix}],

(76)

where B is the adjacency matrix of the chain and D is the degree matrix (see Ref. 4). Note that (76) is a linear dynamical system. We are interested in the velocity auto-correlation function of a tagged oscillator, say, the one at location j = 1. Such an auto-correlation function is defined as

C_{p_{1}} (t) = \frac{{⟨ p_{1} (0) p_{1} (t) ⟩}_{e q}}{{⟨ p_{1} (0) p_{1} (0) ⟩}_{e q}},

(77)

where the average is with respect to the Gibbs canonical distribution $ρ_{e q} = e^{- β H} / Z$ ⁠. It was shown in Ref. 18 that $C_{p_{1}} (t)$ can be obtained analytically by employing Lee’s continued fraction method. The result is the well-known J₀ − J₄ solution

C_{p_{1}} (t) = J_{0} (2 t) - J_{4} (2 t),

(78)

where J_i(t) is the ith Bessel function of the first kind. On the other hand, the Mori-Zwanzig equation derived by the following Mori’s projection

P (\cdot) = \frac{{⟨ (\cdot), p_{1} (0) ⟩}_{e q}}{{⟨ p_{1} (0), p_{1} (0) ⟩}_{e q}} p_{1} (0)

(79)

yields the following GLE for $C_{p_{1}} (t)$ ⁠:

\frac{d C_{p_{1}} (t)}{d t} = Ω_{p_{1}} C_{p_{1}} (t) - \int_{0}^{t} K (s) C_{p_{1}} (t - s) d s .

(80)

Here,

Ω_{p_{1}} = \frac{{⟨ L p_{1} (0), p_{1} (0) ⟩}_{e q}}{{⟨ p_{1} (0), p_{1} (0) ⟩}_{e q}} = 0

since ${⟨ p_{i} (0), q_{j} (0) ⟩}_{e q} = 0$ ⁠, while K(t) is the MZ memory kernel. For the J₀ − J₄ solution, it is possible to derive the memory kernel K(t) analytically. To this end, we simply insert (78) into (80) and apply the Laplace transform

L [\cdot] (s) = \int_{0}^{\infty} (\cdot) e^{- s t} d t

to obtain

\hat{K} (s) = - s + \frac{1}{Ĉ (s)},

(81)

where $Ĉ (s) = L [C_{p_{1}} (t)]$ and $\hat{K} (s) = L [K (t)]$ ⁠. The inverse Laplace transform of (81) can be computed analytically as

K (t) = \frac{J_{1} (2 t)}{t} + 1 .

(82)

With K(t) available, we can verify the memory estimated we derived in Theorem 3.8. To this end,

| K (t) | \leq \frac{‖ {\dot{p}}_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}}{‖ p_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}} = \frac{‖ q_{2} (0) - 2 q_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}}{‖ p_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}} = 2 .

(83)

Here we used the exact solution of the velocity auto-correlation function and displacement auto-correlation function of the fixed-end harmonic chain given by (see Ref. 18)

{⟨ p_{i} (0), p_{j} (0) ⟩}_{e q} = \frac{k_{B} T}{π} \int_{0}^{π} \sin (i x) \sin (j x) d x, {⟨ q_{i} (0), q_{j} (0) ⟩}_{e q} = \frac{k_{B} T}{π} \int_{0}^{π} \frac{\sin (i x) \sin (j x)}{4 \sin^{2} (x / 2)} d x .

In Fig. 1, we plot the absolute value of the memory kernel K(t) together with the theoretical bound (83). It is seen that the upper bound we obtain in this case is of the same order of magnitude as the memory kernel.

FIG. 1.

View large Download slide

Harmonic chain of oscillators. (a) Velocity auto-correlation function $C_{p_{1}} (t)$ and (b) memory kernel K(t) of the corresponding MZ equation. It is seen that our theoretical estimate (83) (dashed line) correctly bounds the MZ memory kernel. Note that the upper bound we obtain is of the same order of magnitude as the memory kernel.

2. Hald system

In this section, we study the Hamiltonian system studied by Chorin et al. in Refs. 12 and 9. The Hamiltonian function is defined as

H (p, q) = \frac{1}{2} (q_{1}^{2} + p_{1}^{2} + q_{2}^{2} + p_{2}^{2} + q_{1}^{2} q_{2}^{2}),

(84)

while the corresponding Hamilton’s equations of motion are

\{\begin{matrix} {\dot{q}}_{1} & = p_{1}, \\ {\dot{p}}_{1} & = - q_{1} (1 + q_{2}^{2}), \\ {\dot{q}}_{2} & = p_{2}, \\ {\dot{p}}_{2} & = - q_{2} (1 + q_{1}^{2}) . \end{matrix}

(85)

We assume that the initial state is distributed according to canonical Gibbs distribution $ρ_{e q} = e^{- H (p, q)} / Z$ ⁠. The partition function Z is given by

Z = e^{1 / 4} {(2 π)}^{3 / 2} K_{0} (\frac{1}{4}),

(86)

where K₀(t) is the modified Bessel function of the second kind. We aim to study the properties of the autocorrelation function of the first component q₁, which is defined as

C_{q_{1}} (t) = \frac{{⟨ q_{1} (0), q_{1} (t) ⟩}_{e q}}{{⟨ q_{1} (0), q_{1} (0) ⟩}_{e q}} .

Obviously, $C_{q_{1}} (0) = 1$ ⁠. The evolution equation for $C_{q_{1}} (t)$ is obtained by using the MZ formulation with the Mori’s projection

P (\cdot) = \frac{{⟨ (\cdot), q_{1} (0) ⟩}_{e q}}{{⟨ q_{1} (0), q_{1} (0) ⟩}_{e q}} q_{1} (0) .

(87)

This yields the GLE

\frac{d C_{q_{1}} (t)}{d t} = Ω_{q_{1}} C_{q_{1}} (t) - \int_{0}^{t} K (s) C_{q_{1}} (t - s) d s .

(88)

The streaming term $Ω_{q_{1}} C_{q_{1}} (t)$ is again identically zero since

Ω_{q_{1}} = \frac{{⟨ L q_{1} (0), q_{1} (0) ⟩}_{e q}}{{⟨ q_{1} (0), q_{1} (0) ⟩}_{e q}} = 0 .

Theorem 3.8 provides the following computable upper bound for the modulus of K(t):

| K (t) | \leq \frac{‖ {\dot{q}}_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}}{‖ q_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}} = \frac{‖ p_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}}{‖ q_{1} (0) ‖_{L_{ρ_{e q}}^{2}}^{2}} = \frac{e^{1 / 4} K_{0} (1 / 4)}{\sqrt{π} U (1 / 2,0,1 / 2)} \approx 1.39786,

(89)

where U(a, b, y) is the confluent hypergeometric function of the second kind. In Fig. 2(a), we plot the correlation function $C_{q_{1}} (t)$ we obtained by sampling (85) and then averaging over all realizations. The samples of the initial condition are generated from the Gibbs equilibrium distribution by using Markov Chain Monte Carlo (MCMC). In Fig. 2(b), we plot the memory kernel K(t) we computed numerically based on $C_{q_{1}} (t)$ ⁠. To compute such a kernel, we inverted numerically the Laplace transform of (88), i.e.,

K (t) = L^{- 1} [- s + \frac{1}{Ĉ (s)}],

(90)

where $Ĉ (s) = L [C_{q_{1}} (t)]$ ⁠. In practice, we replaced $C_{q_{1}} (t)$ within the time interval [0, 20] with a high-order polynomial interpolant at Gauss-Chebyshev-Lobatto nodes, computed the Laplace transform of such polynomial analytically, and then computed the inverse Laplace transform (90) numerically with the Talbot algorithm.¹

FIG. 2.

View large Download slide

Hald Hamiltonian system (85). (a) Autocorrelation function of the displacement q₁(t) and (b) memory kernel of the governing MZ equation. Here $C_{q_{1}} (t)$ is computed by Markov chain Monte-Carlo (MCMC) while K(t) is determined by inverting numerically the Laplace transform in (90) with the Talbot algorithm. It is seen that the theoretical upper bound (89) (dashed line) is of the same order of magnitude as the memory kernel.

B. Non-Hamiltonian systems with infinite-rank projections

In this section, we study the accuracy of the t-model, the H-model, and the H_t model in predicting scalar quantities of interest in non-Hamiltonian systems. In particular, we consider the MZ formulation with Chorin’s projection operator. For the particular case of linear dynamical systems, we also compute the theoretical upper bounds we obtained in Sec. III E for the memory growth and the error in the H-model and compare such bounds with exact results.

1. Linear dynamical systems

We begin by considering a low-dimensional linear dynamical system ẋ = Ax evolving from a random initial state with density ρ₀(x) to verify the MZ memory estimates we obtained in Sec. III E. For simplicity, we choose A to be negative definite

A = e^{C} B e^{- C}, B = [\begin{matrix} - \frac{1}{8} & 0 & 0 \\ 0 & - \frac{2}{3} & 0 \\ 0 & 0 & - \frac{1}{2} \end{matrix}], C = [\begin{matrix} 0 & 1 & 0 \\ - 1 & 0 & 1 \\ 0 & - 1 & 0 \end{matrix}] .

(91)

In this case, the origin of the phase space is a stable node and it is easy to estimate ${∥e^{t L}∥}_{ρ_{0}}$ ⁠. We set x₁(0) = 1 and x₂(0), x₃(0) independent standard normal random variables. In this setting, the semigroup estimates (61) and (62) are explicit,

\begin{matrix} ‖ e^{t L} ‖ & \leq & e^{t ω}, ω = - \frac{1}{2} Tr (A) = 0.6458, \\ ‖ e^{t L Q} ‖ & \leq & e^{t ω_{Q}}, ω_{Q} = ω + \sqrt{A_{11}^{2} + \sum_{i = 2}^{N} \frac{A_{1 i}^{2}}{x_{1}^{2} (0)}} = 1.1621 . \end{matrix}

Therefore, we obtain the following explicit upper bounds for the memory integral and the error of the H-model [see Eqs. (63) and (65)]:

| w_{0} (t) | \leq 0.1964 (e^{1.1621 t} - e^{0.6458 t}),

(92)

| w_{0} (t) - w_{0}^{n} (t) | \leq e^{1.1624 t} \sqrt{{(b^{T} {(M_{11}^{T})}^{n} a_{1})}^{2} x_{1}^{2} (0) + {∥{(M_{11}^{T})}^{n + 1} a∥}_{2}^{2}} \frac{t^{n + 1}}{(n + 1)!} .

(93)

Next, we compare these error bounds with numerical results obtained by solving numerically the H-model (64). For example, the second-order H-model reads

\{\begin{matrix} \frac{d}{d t} E [x_{1} (t) | x_{1} (0)] = - 0.4560 E [x_{1} (t) | x_{1} (0)] + w_{0}^{2} (t), \\ \frac{d w_{0}^{2} (t)}{d t} = 0.0586 E [x_{1} (t) | x_{1} (0)] + w_{1}^{2} (t), \\ \frac{d w_{1}^{2} (t)}{d t} = - 0.0192 E [x_{1} (t) | x_{1} (0)] . \end{matrix}

(94)

In Fig. 3, we demonstrate the convergence of the H-model to the benchmark solution computed by Monte-Carlo simulation as we increase the H-model differentiation order. In Fig. 4, we plot the bound on the memory growth [Eq. (92)] and the bound in the memory error [Eq. (93)] together with exact results.

FIG. 3.

View large Download slide

Convergence of the H-model for the linear dynamical system with matrix (91). The benchmark solution is computed with Monte-Carlo (MC) simulation. Also, the zero-order H-model represents the Markovian approximation to the MZ equation, i.e., the MZ equation without the memory term.

FIG. 4.

View large Download slide

Linear dynamical system with matrix (91). In (a), we plot the memory term $w$ ₀(t) we obtain from Monte Carlo simulation together with the estimated upper bound (92). In (b) and (c), we plot the H-model approximation error $| w_{0} (T) - w_{0}^{n} (T) |$ together with the upper bound (93) for different differentiation orders n and at different times t.

Remark.

The results we just obtained can be obviously extended to higher-dimensional linear dynamical systems. In Fig. 5, we plot the benchmark conditional mean path we obtained through Monte Carlo simulation together with the solution of the H-model () for the 100-dimensional linear dynamical system defined by the matrix (N = 100)

A = [\begin{matrix} - 1 & 1 & \dots & {(- 1)}^{N} \\ 1 \\ ⋮ & B \\ 1 \end{matrix}],

(95)

where B = e^CΛe^−C and

Λ = [\begin{matrix} - \frac{1}{8} & 0 & \dots & 0 \\ 0 & - \frac{2}{9} & ⋮ \\ ⋮ & ⋱ & 0 \\ 0 & \dots & 0 & - \frac{N - 1}{N + 6} \end{matrix}], C = [\begin{matrix} 0 & 1 & 0 \\ - 1 & 0 & ⋱ \\ ⋱ & ⋱ & 1 \\ 0 & - 1 & 0 \end{matrix}] .

It is seen that the H-model converges as we increase the differentiation order in any finite time interval, in agreement with the theoretical prediction of Sec. III E.

FIG. 5.

View large Download slide

Linear dynamical system with matrix A (95). Convergence of the H-model to the conditional mean path solution $E [x_{1} (t) | x_{1} (0)]$ ⁠. The initial condition is set as x₁(0) = 3, while {x₂(0), …, x₁₀₀(0)} are i.i.d. normals.

2. Nonlinear dynamical systems

The hierarchical memory approximation method we discussed in Sec. III D can be applied to nonlinear dynamical systems in the form (1). As we will see, if we employ the H_t-model, then the nonlinearity introduces a closure problem that needs to be addressed properly.

a. Lorenz-63 system.

Consider the classical Lorenz-63 model

\{\begin{matrix} {\dot{x}}_{1} = σ (x_{2} - x_{1}), \\ {\dot{x}}_{2} = x_{1} (r - x_{3}) - x_{2}, \\ {\dot{x}}_{3} = x_{1} x_{2} - β x_{3}, \end{matrix}

(96)

where σ = 10 and β = 8/3. The phase space Liouville operator for this ordinary differential equation (ODE) is

L = σ (x_{2} - x_{1}) \frac{\partial}{\partial x_{1}} + (x_{1} (r - x_{3}) - x_{2}) \frac{\partial}{\partial x_{2}} + (x_{1} x_{2} - β x_{3}) \frac{\partial}{\partial x_{3}} .

We choose the resolved variables to be $\hat{x} = {x_{1}, x_{2}}$ and aim at formally integrating out $\tilde{x} = x_{3}$ by using the Mori-Zwanzig formalism. To this end, we set $x_{3} (0) \sim N (0,1)$ and consider the zeroth-order H_t-model (t-model)

\{\begin{matrix} \frac{d x_{1 m}}{d t} = σ (x_{1 m} - x_{2 m}), \\ \frac{d x_{2 m}}{d t} = - x_{2 m} + r x_{1 m} - t x_{1 m}^{2} x_{2 m}, \end{matrix}

(97)

where $x_{1 m} (t) = E [x_{1} (t) | x_{1} (0), x_{2} (0)]$ and $x_{2 m} (t) = E [x_{2} (t) | x_{1} (0), x_{2} (0)]$ are the conditional mean paths. To obtain this system, we introduced the following mean field closure approximation:

\begin{matrix} t P e^{t L} P L Q L x_{2} (0) & = & - t E [x_{1} {(t)}^{2} x_{2} (t) | x_{1} (0), x_{2} (0)] \\ ≃ & - t E {[x_{1} (t) | x_{1} (0), x_{2} (0)]}^{2} E [x_{2} (t) | x_{1} (0), x_{2} (0)] \\ = & - t x_{1 m}^{2} x_{2 m} . \end{matrix}

(98)

Higher-order H_t-models can be derived based on (98). As is well known, if r < 1, the fixed point (0, 0, 0) is a global attractor and exponentially stable. In this case, the t-model (zeroth-order H_t-model) yields an accurate prediction of the conditional mean path for long time (see Fig. 6). On the other hand, if we consider the chaotic regime at r = 28, then the t-model and its higher-order extension, i.e., the H_t-model, are accurate only for relatively short time. This is in agreement with our theoretical predictions. In fact, different from linear systems where the hierarchical representation of the memory integral can be proven to be convergent for long time, in nonlinear systems, the memory hierarchy is, in general, provably convergent only in a short time period (Theorem 3.7 and Corollary 3.4.3). This does not mean that the H-model or the H_t-model is not accurate for nonlinear systems. It just means that the accuracy depends on the system, the quantity of interest, and the initial condition.

FIG. 6.

View large Download slide

Accuracy of the H_t model in representing the conditional mean path in the Lorenz-63 system (96). It is seen that if r = 0.5 (first row), then the zeroth-order H_t-model, i.e., the t-model is accurate for long integration times. On the other hand, if we consider the chaotic regime at r = 28 (second row), then we see that the t-model and its high-order extension (H_t-model) are accurate only for relatively short time.

b. Modified Lorenz-96 system.

As an example of a high-dimensional nonlinear dynamical system, we consider the following modified Lorenz-96 system:^22,26

\{\begin{matrix} {\dot{x}}_{1} = - x_{1} + x_{1} x_{2} + F, \\ {\dot{x}}_{2} = - x_{2} + x_{1} x_{3} + F, \\ ⋮ \\ \dot{x_{i}} = - x_{i} + (x_{i + 1} - x_{i - 2}) x_{i - 1} + F, \\ ⋮ \\ {\dot{x}}_{N} = x_{N} - x_{N - 2} x_{N - 1} + F, \end{matrix}

(99)

where F is constant. As is well known, depending on the values of N and F, this system can exhibit a wide range of behaviors.²² Suppose we take the resolved variables to be $\hat{x} = {x_{1}, x_{2}}$ ⁠. Correspondingly, the unresolved ones, i.e., those we aim at integrating through the MZ framework, are $\tilde{x} = {x_{3}, \dots, x_{N}}$ ⁠, which we set to be independent standard normal random variables. By using the mean field approximation (98), we obtain the following zeroth-order H_t-model (t-model) of the modified Lorenz-96 system (99):

\{\begin{matrix} {\dot{x}}_{1 m} = - x_{1 m} + x_{1 m} x_{2 m} + F, \\ {\dot{x}}_{2 m} = - x_{2 m} + F + t (x_{1 m}^{2} x_{2 m} - x_{1 m} F) . \end{matrix}

(100)

In Fig. 7, we study the accuracy of the H_t-model in representing the conditional mean path with F = 5 and N = 100. It is seen that the the H_t-model converges only for short time (in agreement with the theoretical predictions), and it provides results that are more accurate than the classical t-model.

FIG. 7.

View large Download slide

Accuracy of the H_t-model in representing the conditional mean path in the Lorenz-96 system (96). Here we set F = 5 and N = 100. It is seen that the H_t-model converges only for short time and provides results that are more accurate that the classical t-model.

V. SUMMARY

In this paper, we developed a thorough mathematical analysis to deduce conditions for the accuracy and convergence of different approximations of the memory integral in the Mori-Zwanzig (MZ) equation. In particular, we studied the short memory approximation, the t-model, and various hierarchical memory approximation techniques. We also derived useful upper bounds for the MZ memory integral which allowed us to estimate a priori its contribution to the dynamics. Such upper bounds are easily computable for systems with finite-rank projections—such as Mori’s projection. Furthermore, if the system is Hamiltonian, then the MZ memory bounds we derived are tight. This can be attributed to the fact that the Koopman operator of a Hamiltonian dynamical system is always a contraction relative to the $L_{ρ_{e q}}^{2}$ norm weighted by the Gibbs canonical distribution ρ_eq. In the more general case of dissipative systems, however, we do not have any guarantee that the MZ memory bounds we obtained are tight. Indeed, for dissipative systems, the numerical abscissa in Eq. (17) [see also Eqs. (21a) and (21b)] can be positive, yielding an exponential growth of the memory bounds. In the presence of infinite-rank projections—such as Chorin’s projection—we found that it is extremely difficult to obtain explicitly computable upper bounds for the MZ memory or any of its approximations discussed in this paper. The main difficulties are related to the lack of an explicitly computable upper bound for the operator semigroup $e^{t L Q}$ ⁠, which involves the unbounded generator $L Q$ (see the discussion in Sec. III A and Appendix A). To demonstrate the accuracy of the approximation methods and error bounds we developed in this paper, we provided numerical examples involving Hamiltonian dynamical systems (linear and nonlinear) and more general dissipative nonlinear systems evolving from random initial states. The numerical results we obtained are found to be in agreement with our theoretical predictions.

ACKNOWLEDGMENTS

This work was supported by the Air Force Office of Scientific Research Grant No. FA9550-16-586-1-0092.

APPENDIX A: SEMIGROUP BOUNDS VIA FUNCTION DECOMPOSITION

In looking for the numerical abscissa¹⁴ (i.e., the logarithmic norm) of $L Q$ ⁠, we seek to bound

sup_{D (L Q) ∋ x \neq 0} R \frac{{⟨ x, L Q x ⟩}_{σ}}{{⟨ x, x ⟩}_{σ}} .

Notice that if $P L Q = 0$ ⁠, then $L Q = Q L Q$ so that we have the previously proven bound

sup R \frac{{⟨ x, Q L Q x ⟩}_{σ}}{{⟨ x, x ⟩}_{σ}} \leq - \frac{1}{2} inf {div}_{σ} F .

In this section, we consider what happens when $P L Q \neq 0$ ⁠. To that end, let us note that $x \in D (L Q)$ may be decomposed as

x = Q x + α P L Q x + P y,

where $α \in C$ and $P y$ are orthogonal to $P L Q x$ ⁠. In other words, we define $P y$ as

P y = P x - \frac{{⟨ P L Q x, P x ⟩}_{σ}}{{⟨ P L Q x, P L Q x ⟩}_{σ}} P L Q x

and define α as

α = \frac{{⟨ P L Q x, P x ⟩}_{σ}}{{⟨ P L Q x, P L Q x ⟩}_{σ}} .

Then

\begin{matrix} R \frac{{⟨ x, L Q x ⟩}_{σ}}{{⟨ x, x ⟩}_{σ}} & = & R \frac{{⟨ Q x + α P L Q x + P y, L Q x ⟩}_{σ}}{{⟨ x, x ⟩}_{σ}} \\ = & \frac{R {⟨ Q x, L Q x ⟩}_{σ} + R (α) ‖ P L Q x ‖_{σ}^{2}}{‖ Q x ‖_{σ}^{2} + | α |^{2} ‖ P L Q x ‖_{σ}^{2} + ‖ P y ‖_{σ}^{2}} \\ \leq & max [0, \frac{R {⟨ Q x, L Q x ⟩}_{σ} + R (α) ‖ P L Q x ‖_{σ}^{2}}{‖ Q x ‖_{σ}^{2} + | α |^{2} ‖ P L Q x ‖_{σ}^{2}}] . \end{matrix}

Since we assume $P L Q \neq 0$ ⁠, there exists x such that $P L Q x \neq 0$ and then for any α such that

R (α) \geq \frac{R {⟨ Q x, L Q x ⟩}_{σ}}{‖ P L Q x ‖_{σ}^{2}},

we have

R {⟨ Q x, L Q x ⟩}_{σ} + R (α) ‖ P L Q x ‖_{σ}^{2} \geq 0

so that for $P L Q \neq 0$ ⁠,

sup_{D (L Q) ∋ x \neq 0} R \frac{{⟨ x, L Q x ⟩}_{σ}}{{⟨ x, x ⟩}_{σ}} = sup_{D (L Q) ∋ x \neq 0} \frac{R {⟨ Q x, L Q x ⟩}_{σ} + R (α) ‖ P L Q x ‖_{σ}^{2}}{‖ Q x ‖_{σ}^{2} + | α |^{2} ‖ P L Q x ‖_{σ}^{2}} .

Now, fix any $Q x \in D (L) \neq 0$ and consider the expression

\frac{R {⟨ Q x, L Q x ⟩}_{σ} + R (α) ‖ P L Q x ‖_{σ}^{2}}{‖ Q x ‖_{σ}^{2} + | α |^{2} ‖ P L Q x ‖_{σ}^{2}} = \frac{ξ + R (α) β^{2}}{1 + | α |^{2} β^{2}},

where

ξ = ξ (Q x) = R \frac{{⟨ Q x, L Q x ⟩}_{σ}}{‖ Q x ‖_{σ}^{2}}, β = β (Q x) = \frac{‖ P L Q x ‖_{σ}}{‖ Q x ‖_{σ}} .

Then, for this fixed $Q x$ ⁠,

\frac{ξ + R (α) β^{2}}{1 + | α |^{2} β^{2}} \leq max_{a \in R} \frac{ξ + a β^{2}}{1 + a^{2} β^{2}} .

Differentiating with respect to a and setting equal to zero, we find that the latter expression is extremized when

0 = β^{2} (1 + a^{2} β^{2}) - 2 a β^{2} (ξ + a β^{2}),

i.e., when

β^{2} a^{2} + 2 ξ a - 1 = 0, a = \frac{- ξ \pm \sqrt{ξ^{2} + β^{2}}}{β^{2}} .

Since β² > 0, $\frac{ξ + a β^{2}}{1 + a^{2} β^{2}}$ is maximized at $â = \frac{- ξ + \sqrt{ξ^{2} + β^{2}}}{β^{2}}$ ⁠. Then

ξ + â β^{2} = \sqrt{ξ^{2} + β^{2}}, 1 + â^{2} β^{2} = 2 (1 - ξ â) = 2 (\frac{ξ^{2} + β^{2} - ξ \sqrt{ξ^{2} + β^{2}}}{β^{2}})

so that

max_{a \in R} \frac{ξ + a β^{2}}{1 + a^{2} β^{2}} = \frac{ξ + â β^{2}}{1 + â^{2} β^{2}} = \frac{1}{2} \frac{β^{2}}{\sqrt{ξ^{2} + β^{2}} - ξ} = \frac{1}{2} (\sqrt{ξ^{2} + β^{2}} + ξ) .

Therefore,

sup_{D (L Q) ∋ x \neq 0} R \frac{{⟨ x, L Q x ⟩}_{σ}}{{⟨ x, x ⟩}_{σ}} = sup_{D (L) ∋ (Q x) \neq 0} \frac{1}{2} [\sqrt{ξ^{2} (Q x) + β^{2} (Q x)} + ξ (Q x)] .

When $P L Q$ is unbounded, which is the typical case when $P$ is an infinite-rank projection, such as most conditional expectations, there is unlikely to be a finite numerical abscissa for $L Q$ ⁠. In particular, notice that if div_σ(F) is a bounded function (bounded both above and below), then $ξ (Q x)$ is bounded for all $Q x$ while $β (Q x)$ is unbounded, in which case

sup_{D (L Q) ∋ x \neq 0} R \frac{{⟨ x, L Q x ⟩}_{σ}}{{⟨ x, x ⟩}_{σ}} = \infty .

It follows Refs. 37 and 30 that in these cases, $‖ e^{t L Q} ‖_{σ}$ has an infinite slope at t = 0, and therefore there is no finite ω such that $‖ e^{t L Q} ‖_{σ} \leq e^{ω t}$ for all t ≥ 0 (see Ref. 14). Assuming still that $L Q$ generates a strongly continuous semigroup, we must look to bound the semigroup as $‖ e^{t L Q} ‖_{σ} \leq M e^{ω t}$ ⁠, where

ω > ω_{0} = lim_{t \to \infty} \frac{\ln ‖ e^{t L Q} ‖_{σ}}{t}

and

M \geq M (ω) = sup {‖ e^{t L Q} ‖_{σ} e^{- ω t} : t \geq 0} .

On the other hand, if $P L Q$ is a bounded operator, for example, when $P$ is a finite-rank projection [e.g., Mori’s projection (12)], there exists a finite value for the numerical abscissa. Indeed, in this case, since ζ is bounded by ω = −inf div_σ(F), the numerical abscissa of $L Q$ may be bounded as

ω_{L Q} = sup_{D (L Q) ∋ x \neq 0} R \frac{{⟨ x, L Q x ⟩}_{σ}}{{⟨ x, x ⟩}_{σ}} \leq \frac{1}{2} [\sqrt{ω^{2} + ‖ P L Q ‖_{σ}^{2}} + ω] .

(A1)

Alternatively, in the case of finite rank $P$ ⁠, the operator $L Q$ may be thought of as a bounded perturbation of $L$ ⁠, i.e., $L Q = L - L Q$ ⁠, and the numerical abscissa of $L Q$ can be bounded using the bounded perturbation theorem Ref. 16, III.1.3, obtaining

ω_{L Q} \leq ω + ‖ L P ‖_{σ} .

(A2)

Either of these bounds for $ω_{L Q}$ can be used to bound the semigroup norm

\begin{matrix} ‖ e^{t L Q} ‖_{σ} & \leq & e^{ω_{L Q} t} \leq e^{\frac{1}{2} (\sqrt{ω^{2} + ‖ P L Q ‖_{σ}^{2}} + ω) t}, \\ ‖ e^{t L Q} ‖_{σ} & \leq & e^{ω_{L Q} t} \leq e^{(ω + ‖ L P ‖_{σ}) t} . \end{matrix}

Which of these two estimates gives the tighter bound will generally depend on the values of $‖ P L Q ‖_{σ}$ and $‖ L P ‖_{σ}$ ⁠. It may be noted, however, that, when σ is invariant, ω = 0 and $L$ is skew-adjoint so that

{‖ P L Q ‖}_{σ} = {‖ Q L^{†} P ‖}_{σ} = {‖ Q L P ‖}_{σ} \leq {‖ L P ‖}_{σ}

and therefore the bound in (A1) is half that of (A2).

APPENDIX B: THE MORI-ZWANZIG FORMULATION IN PDF SPACE

It was shown in Ref. 15 that the Banach dual of (4) defines an evolution in the probability density function space. Specifically, the joint probability density function of the state vector u(t) that solves Eq. (1) is pushed forward by the Frobenius-Perron operator $F (t, 0)$ [Banach dual of the Koopman operator (3)],

p (x, t) = F (t, s) p (x, s), F (t, s) = e^{(t - s) M},

(B1)

where

M (x) p (x, t) = - \nabla \cdot (F (x) p (x, t)) .

(B2)

By introducing a projection $P$ in the space of probability density functions and its complement $Q = I - P$ ⁠, it is easy to show that the projected density $P p$ satisfies the MZ equation⁴²

\frac{\partial P p (t)}{\partial t} = P M P p (t) + P e^{t Q M} Q p (0) + \int_{0}^{t} P M e^{(t - s) Q M} Q M P p (s) d s .

(B3)

In Appendixes B 1 and B 2, we perform an analysis of different types of approximations of the MZ memory integral

\int_{0}^{t} P M e^{(t - s) Q M} Q M P p (s) d s .

(B4)

The main objective of such analysis is to establish rigorous error bounds for widely used approximation methods and also propose new provably convergent approximation schemes.

1. Analysis of the memory integral

In this section, we develop the analysis of the memory integral arising in the PDF formulation of the MZ equation. The starting point is the definition (B4). As before, we begin with the following estimate of upper bound estimation of the integral:

Theorem B.1

(Memory growth). Let

e^{t M Q}

and

e^{t M}

be strongly continuous semigroups with upper bounds

‖ e^{t M} ‖ \leq M e^{t ω}

and

‖ e^{t M Q} ‖ \leq M_{Q} e^{t ω_{Q}}

, and let T > 0 be a fixed integration time. Then for any 0 ≤ t ≤ T, we have

∥ \int_{0}^{t} P M e^{(t - s) Q M} Q M P p (s) d s ∥ \leq N_{0} (t),

where

N_{0} (t) = \{\begin{matrix} t C_{4}, & ω_{Q} = 0, \\ \frac{C_{4}}{ω_{Q}} (e^{t ω_{Q}} - 1), & ω_{Q} \neq 0 \end{matrix}

and

C_{4} = max_{0 \leq s \leq T} ‖ P ‖ ‖ M P M p (s) ‖

. Moreover, N(t) satisfies

lim_{t \to 0}

N(t) = 0.

Proof.

Consider

\begin{matrix} ∥ \int_{0}^{t} P M e^{(t - s) Q M} Q M P p (s) d s ∥ & = & ∥ \int_{0}^{t} P e^{(t - s) M Q} M Q M P p (s) d s ∥ \\ \leq & C_{4} M_{Q} \int_{0}^{t} e^{(t - s) ω_{Q}} d s \\ = & \{\begin{matrix} t C_{4}, & ω_{Q} = 0, \\ \frac{C_{4}}{ω_{Q}} (e^{t ω_{Q}} - 1), & ω_{Q} \neq 0, \end{matrix} \end{matrix}

where

C_{4} = max_{0 \leq s \leq T} ‖ P ‖ ‖ M Q M P p (s) ‖

⁠.

Theorem B.2

(Memory approximation via the t-model). Let

e^{t M Q}

and

e^{t M}

be strongly continuous semigroups with bounds

‖ e^{t M} ‖ \leq M e^{t ω}

and

‖ e^{t M Q} ‖ \leq M_{Q} e^{t ω_{Q}}

, and let T > 0 be a fixed integration time. If the function

k (s, t) = P M e^{(t - s) Q M} Q M P p (s)

(integrand of the memory term) is at least twice differentiable respect to s for all t ≥ 0, then

∥ \int_{0}^{t} P M e^{(t - s) Q M} Q M P p (s) d s - t P M Q M P p (t) ∥ \leq N_{1} (t),

where N₁(t) is defined as

N_{1} (t) = \{\begin{matrix} t (M_{Q} + 1) C_{4}, & ω_{Q} = 0, \\ \frac{C_{2} M_{Q}}{ω_{Q}} (e^{t ω_{Q}} - 1) + t C_{4}, & ω_{Q} \neq 0 \end{matrix}

and C₄ is as in Theorem B.1.

Proof.

\begin{matrix} ∥ \int_{0}^{t} P M e^{(t - s) Q M} Q M P p (s) d s - t P M Q M P p (t) ∥ & \leq & C_{4} M_{Q} \int_{0}^{t} e^{(t - s) ω_{Q}} d s + C_{4} t \\ = & \{\begin{matrix} t (M_{Q} + 1) C_{4}, & ω_{Q} = 0, \\ \frac{C_{4} M_{Q}}{ω_{Q}} (e^{t ω_{Q}} - 1) + t C_{4}, & ω_{Q} \neq 0 . \end{matrix} \end{matrix}

2. Hierarchical memory approximation in PDF space

The hierarchical memory approximation methods we discussed in Sec. III D can also be developed in the PDF space. To this end, let us first define

v_{0} (t) = \int_{0}^{t} P M e^{(t - s) Q M} Q M P p (s) d s .

(B5)

By repeatedly differentiating $v$ ₀(t) with respect to time [assuming $v$ ₀(t) smooth enough], we obtain the hierarchy of equations

\frac{\partial}{\partial t} v_{i - 1} (t) = P M {(Q M)}^{i} P p (t) + v_{i} (t), i = 1, \dots, n,

where

v_{i} (t) = \int_{0}^{t} P M e^{(t - s) Q M} {(Q M)}^{i + 1} P p (s) d s .

By following closely the discussion in Sec. III D, we introduce the hierarchy of memory equations

\{\begin{matrix} \frac{d v_{0}^{n} (t)}{d t} = P M Q M P p (t) + v_{1}^{n} (t), \\ \frac{d v_{1}^{n} (t)}{d t} = P M {(Q M)}^{2} P p (t) + v_{2}^{n} (t), \\ ⋮ \\ \frac{d v_{n - 1}^{n} (t)}{d t} = P M {(Q M)}^{n} P p (t) + v_{n}^{e_{n}} (t) \end{matrix}

(B6)

and approximate the last term in such hierarchy in different ways. Specifically, we consider

\begin{matrix} v_{n}^{e_{n}} (t) & = & \int_{t}^{t} P M e^{(t - s) Q M} {(Q M)}^{n + 1} P p (s) d s = 0 (H - model), \\ v_{n}^{e_{n}} (t) & = & \int_{max (0, t - Δ t)}^{t} P M e^{(t - s) Q M} {(Q M)}^{n + 1} P p (s) d s (Type I Finite Memory Approximation), \\ v_{n}^{e_{n}} (t) & = & \int_{min (t, t_{n})}^{t} P M e^{(t - s) Q M} {(Q M)}^{n + 1} P p (s) d s (Type II Finite Memory Approximation), \\ v_{n}^{e n} (t) & = & t P M {(Q M)}^{n + 1} P p (t) (H_{t} - model) . \end{matrix}

Hereafter we establish the accuracy of the approximation schemes resulting from the substitution of each $v_{n}^{e_{n}} (t)$ above into (B6).

Theorem B.3

(Accuracy of the H-model). Let

e^{t M}

and

e^{t M Q}

be strongly continuous semigroups, T > 0 a fixed integration time, and

β_{i} = \frac{sup_{s \in [0, T]} ‖ {(M Q)}^{i + 1} M P p (s) ‖}{sup_{s \in [0, T]} ‖ {(M Q)}^{i} M P p (s) ‖}, 1 \leq i \leq n .

(B7)

Then for 1 ≤ q ≤ n, we have

‖ v_{0} (t) - v_{0}^{q} (t) ‖ \leq N_{2}^{q} (t),

where

N_{2}^{q} (t) = A_{2} C_{4} (\prod_{i = 1}^{q} β_{i}) \frac{t^{q + 1}}{(q + 1)!},

A_{2} = max_{s \in [0, T]} e^{s ω_{Q}}

, and C₄ is as in Theorem B.1.

Proof.

The error at the nth level can be bounded as

\begin{matrix} ‖ v_{0} (t) - v_{0}^{q} (t) ‖ & \leq & \int_{0}^{t} \int_{0}^{τ_{q}} \dots \int_{0}^{τ_{2}} ‖ P M e^{(τ_{1} - s) Q M} {(Q M)}^{p + 1} P p (s) ‖ d s d τ_{1} \dots d τ_{q} \\ \leq & A_{2} C_{4} (\prod_{i = 1}^{q} β_{i}) \frac{t^{q + 1}}{(q + 1)!}, \end{matrix}

where

A_{2} = max_{s \in [0, T]} e^{s ω_{Q}} = \{\begin{matrix} 1, & ω_{Q} \leq 0 \\ e^{T ω_{Q}}, & ω_{Q} \geq 0 . \end{matrix}

(B8)

Let

β_{i} = \frac{sup_{s \in [0, T]} ‖ {(M Q)}^{i + 1} M P p (s) ‖}{sup_{s \in [0, T]} ‖ {(M Q)}^{i} M P p (s) ‖},

(B9)

under the assumption that these quantities are finite. Then we have

‖ v_{0} (t) - v_{0}^{q} (t) ‖ \leq A_{2} C_{4} (\prod_{i = 1}^{q} β_{i}) \frac{t^{q + 1}}{(q + 1)!} .

Corollary B.3.1

(Uniform convergence of the H-model). If β_i in Theorem B.3 satisfy

β_{i} < \frac{i + 1}{T}, 1 \leq i \leq n

for any fixed time T, then there exits a sequence δ₁ > δ₂ > ⋯ > δ_n such that

‖ v_{0} (T) - v_{0}^{q} (T) ‖ \leq δ_{q},

where 1 ≤ q ≤ n.

Corollary B.3.2

(Asymptotic convergence of the H-model). If β_i in Theorem B.3 satisfy

β_{i} < C, 1 \leq i < + \infty

for some constant C, then for any fixed time T and arbitrary δ > 0, there exits an integer q such that for all n > q,

‖ v_{0} (T) - v_{0}^{n} (T) ‖ \leq δ .

The Proofs of the Corollaries B.3.1 and B.3.2 closely follow the Proofs of Corollaries 3.4.1 and 3.4.2. Therefore we omit details here.

Theorem B.4

(Accuracy of Type-I FMA). Let

e^{t M}

and

e^{t M Q}

be strongly continuous semigroups, T > 0 be a fixed integration time, and let

β_{i} = \frac{sup_{s \in [0, T]} ‖ {(M Q)}^{i + 1} M P p (s) ‖}{sup_{s \in [0, T]} ‖ {(M Q)}^{i} M P p (s) ‖}, 1 \leq i \leq n .

(B10)

Then for 1 ≤ q ≤ n,

‖ v_{0} (t) - v_{0}^{q} (t) ‖ \leq N_{3}^{q} (t),

where

N_{3}^{q} (t) = A_{2} C_{4} (\prod_{i = 1}^{q} β_{i}) \frac{{(t - Δ t_{q})}^{q + 1}}{(q + 1)!}

and C₄ is as in Theorem B.1.

Proof.

The proof is very similar with the Proof of Theorem 3.5. We begin with the estimate of

v_{0} (t) - v_{0}^{q} (t)

⁠,

\begin{matrix} \int_{0}^{max (0, t - Δ t_{q})} \int_{0}^{{\tilde{τ}}_{q}} \dots \int_{0}^{{\tilde{τ}}_{2}} \int_{0}^{{\tilde{τ}}_{1}} P e^{({\tilde{τ}}_{1} + Δ t_{q} - s) M Q} {(M Q)}^{q + 1} M Q p (s) d s d {\tilde{τ}}_{1} \dots d {\tilde{τ}}_{q}, \\ = \{\begin{matrix} 0, & 0 \leq t \leq Δ t_{q}, \\ \int_{0}^{t - Δ t_{q}} \int_{0}^{σ} \frac{{(t - Δ t_{q} - σ)}^{q - 1}}{(q - 1)!} P e^{(σ + Δ t_{q} - s) M Q} {(M Q)}^{q + 1} M Q p (s) d s d σ, & t \geq Δ t_{q} . \end{matrix} \end{matrix}

This can be bounded by following the technique in the Proof of Theorem B.3. This yields

‖ v_{0} (t) - v_{0}^{q} (t) ‖ \leq \{\begin{matrix} 0, & 0 \leq t \leq Δ t_{q}, \\ A_{2} C_{4} (\prod_{i = 1}^{q} β_{i}) \frac{{(t - Δ t_{q})}^{q + 1}}{(q + 1)!}, & t \geq Δ t_{q} . \end{matrix}

(B11)

Corollary B.4.1

(Uniform convergence of Type-I FMA). If β_i in Theorem B.4 satisfy

β_{i} < (i + 1) {[\frac{δ i!}{C_{4} A_{2} (\prod_{k = 1}^{i} β_{k})}]}^{- \frac{1}{i}},

(B12)

then for any δ > 0, there exists an ordered time sequence Δt_n < Δt_n−1 < ⋯ < Δt₁ < T such that

‖ w_{0} (T) - w_{0}^{q} (T) ‖ \leq δ, 1 \leq q \leq n

which satisfies

Δ t_{q} \leq T - {[\frac{δ (q + 1)!}{C_{4} A_{2} (\prod_{i = 1}^{q} β_{i})}]}^{\frac{1}{q + 1}} .

The proof is very similar with the Proof of Corollary 3.5.1, and therefore we omit it.

Theorem B.5

(Accuracy of Type-II FMA). Let

e^{t M}

and

e^{t M Q}

be strongly continuous semigroups and T > 0 be a fixed integration time. Set

β_{i} = \frac{sup_{s \in [0, T]} ‖ {(M Q)}^{i + 1} M P p (s) ‖}{sup_{s \in [0, T]} ‖ {(M Q)}^{i} M P p (s) ‖}, 1 \leq i \leq n .

(B13)

Then for 1 ≤ q ≤ n,

‖ v_{0} (t) - v_{0}^{q} (t) ‖ \leq N_{4}^{q} (t),

where

N_{4}^{q} (t) = \frac{C_{4}}{ω_{Q}} [1 - e^{- t_{q} ω_{Q}}] (\prod_{i = 1}^{q} β_{i}) f_{q} (ω_{Q}, t),

f_{q} (ω_{Q}, t)

is defined in (), and C₄ is as in Theorem B.1.

Proof.

The proof is very similar with the Proof of Theorem 3.6. Hereafter we provide the proof for the case when

ω_{Q} > 0

⁠. Other cases can be easily obtained by using the same method. First of all, we have error estimation

\begin{matrix} ‖ v_{0} (t) - v_{0}^{q} (t) ‖ & \leq & \int_{0}^{t} \int_{0}^{τ_{q}} \dots \int_{0}^{τ_{2}} \int_{0}^{t_{q}} ‖ P M e^{(τ_{1} - s) Q M} {(Q M)}^{q + 1} P p (s) ‖ d s d τ_{1} \dots d τ_{q} \\ \leq & C_{4} (\prod_{i = 1}^{q} β_{i}) \int_{0}^{t} \int_{0}^{τ_{q}} \dots \int_{0}^{τ_{2}} \int_{0}^{t_{q}} e^{(τ_{1} - s) ω_{Q}} d s d τ_{1} \dots d τ_{q} \\ = & \frac{C_{4}}{ω_{Q}} [1 - e^{- t_{q} ω_{Q}}] (\prod_{i = 1}^{q} β_{i}) f_{q} (ω_{Q}, t), \end{matrix}

where

f_{q} (ω_{Q}, t)

is as in ().

Corollary B.5.1

(Uniform convergence of Type-II FMA). If β_i in Theorem B.5 satisfy

\{\begin{matrix} β_{i} < \frac{i}{T}, & ω_{Q} = 0, \\ β_{i} < ω_{Q} \frac{e^{T ω_{Q}} - \sum_{k = 0}^{i - 2} \frac{{(T ω_{Q})}^{k}}{k!}}{e^{T ω_{Q}} - \sum_{k = 0}^{i - 1} \frac{{(T ω_{Q})}^{k}}{k!}}, & ω_{Q} \neq 0 \end{matrix}

for all t ∈ [0, T], then for arbitrarily small δ > 0, there exists an ordered time sequence 0 < t₀ < t₁ < … t_n < T such that

‖ v_{0} (T) - v_{0}^{q} (T) ‖ \leq δ, 1 \leq q \leq n,

which satisfies

t_{q} \geq \frac{1}{ω_{Q}} \ln [1 - \frac{δ ω_{Q}^{q + 1}}{C_{4} (\prod_{i = 1}^{q} β_{i}) [e^{T ω_{Q}} - \sum_{k = 0}^{q - 1} \frac{{(T ω_{Q})}^{k}}{k!}]}] .

Proof.

To ensure that

‖ v_{0} (t) - v_{0}^{q} (t) ‖ \leq δ

for all 0 ≤ t ≤ T, we can take (for

ω_{Q} > 0

⁠)

\frac{C_{4}}{ω_{Q}} [1 - e^{- t_{q} ω_{Q}}] (\prod_{i = 1}^{q} β_{i}) f_{q} (ω_{Q}, T) = max_{t \in [0, T]} \frac{C_{2}}{ω_{Q}} [1 - e^{- t_{q} ω_{Q}}] (\prod_{i = 1}^{q} β_{i}) f_{q} (ω_{Q}, t) \leq δ .

Therefore

\begin{matrix} e^{- t_{q} ω_{Q}} & \geq & 1 - \frac{δ ω_{Q}}{C_{4} (\prod_{i = 1}^{q} β_{i}) f_{q} (ω_{Q}, T)}, i.e., \\ t_{q} & \leq & \frac{1}{ω_{Q}} \ln [1 - \frac{δ ω_{Q}^{q + 1}}{C_{4} (\prod_{i = 1}^{q} β_{i}) [e^{T ω_{Q}} - \sum_{k = 0}^{q - 1} \frac{{(T ω_{Q})}^{k}}{k!}]}] . \end{matrix}

Since for

ω_{Q} > 0

⁠, we have have condition

β_{i} (t) < ω_{Q} \frac{e^{T ω_{Q}} - \sum_{k = 0}^{i - 2} \frac{{(T ω_{Q})}^{k}}{k!}}{e^{T ω_{Q}} - \sum_{k = 0}^{i - 1} \frac{{(T ω_{Q})}^{k}}{k!}} .

Thus, there exists an ordered time sequence 0 < t₁ < ⋯ < t_n such that

‖ v_{0} (T) - v_{0}^{n} (T) ‖ \leq δ

⁠. As in Theorem 3.6, this δ-bound on the error holds for all t_n (with the upper bound as above), which implies the existence of such an increasing time sequence 0 < t₁ < ⋯ < t_n with t_n bounded from below by the same quantities.

REFERENCES

1.

J.

Abate

and

W.

Whitt

, “

A unified framework for numerically inverting Laplace transforms

,”

INFORMS J. Comput.

18

(

4

),

408

–

421

(

2006

).

https://doi.org/10.1287/ijoc.1050.0137

Google Scholar

Crossref

2.

B. J.

Alder

and

T. E.

Wainwright

, “

Decay of the velocity autocorrelation function

,”

Phys. Rev. A

1

(

1

),

18

(

1970

).

https://doi.org/10.1103/physreva.1.18

Google Scholar

Crossref

3.

R. J.

Baxter

,

Exactly Solved Models in Statistical Mechanics

(

Elsevier

,

2016

).

Google Scholar

4.

N.

Biggs

,

Algebraic Graph Theory

(

Cambridge University Press

,

1993

).

Google Scholar

5.

C.

Brennan

and

D.

Venturi

, “

Data-driven closures for stochastic dynamical systems

,”

J. Comput. Phys.

372

,

281

–

298

(

2018

).

https://doi.org/10.1016/j.jcp.2018.06.038

Google Scholar

Crossref

6.

H.-P.

Breuer

,

B.

Kappler

, and

F.

Petruccione

, “

The time-convolutionless projection operator technique in the quantum theory of dissipation and decoherence

,”

Ann. Phys.

291

(

1

),

36

–

70

(

2001

).

https://doi.org/10.1006/aphy.2001.6152

Google Scholar

Crossref

7.

A.

Chertock

,

D.

Gottlieb

, and

A.

Solomonoff

, “

Modified optimal prediction and its application to a particle-method problem

,”

J. Sci. Comput.

37

(

2

),

189

–

201

(

2008

).

https://doi.org/10.1007/s10915-008-9242-4

Google Scholar

Crossref

8.

H.

Cho

,

D.

Venturi

, and

G. E.

Karniadakis

, “

Statistical analysis and simulation of random shocks in Burgers equation

,”

Proc. R. Soc. A

470

(

2171

),

20140080

(

2014

).

https://doi.org/10.1098/rspa.2014.0080

Google Scholar

Crossref

9.

A.

Chorin

,

O.

Hald

, and

R.

Kupferman

, “

Optimal prediction with memory

,”

Physica D

166

(

3-4

),

239

–

257

(

2002

).

https://doi.org/10.1016/s0167-2789(02)00446-3

Google Scholar

Crossref

10.

A. J.

Chorin

,

O. H.

Hald

, and

R.

Kupferman

, “

Optimal prediction and the Mori-Zwanzig representation of irreversible processes

,”

Proc. Natl. Acad. Sci. U. S. A.

97

(

7

),

2968

–

2973

(

2000

).

https://doi.org/10.1073/pnas.97.7.2968

Google Scholar

Crossref

PubMed

11.

A. J.

Chorin

,

R.

Kupferman

, and

D.

Levy

, “

Optimal prediction for Hamiltonian partial differential equations

,”

J. Comput. Phys.

162

(

1

),

267

–

297

(

2000

).

https://doi.org/10.1006/jcph.2000.6536

Google Scholar

Crossref

12.

A. J.

Chorin

and

P.

Stinis

, “

Problem reduction, renormalization and memory

,”

Commun. Appl. Math. Comput. Sci.

1

(

1

),

1

–

27

(

2006

).

https://doi.org/10.2140/camcos.2006.1.1

Google Scholar

Crossref

13.

E.

Darve

,

J.

Solomon

, and

A.

Kia

, “

Computing generalized Langevin equations and generalized Fokker-Planck equations

,”

Proc. Natl. Acad. Sci. U. S. A.

106

(

27

),

10884

–

10889

(

2009

).

https://doi.org/10.1073/pnas.0902633106

Google Scholar

Crossref

PubMed

14.

E. B.

Davies

, “

Semigroup growth bounds

,”

J. Oper. Theory

53

(

2

),

225

–

249

(

2005

).

Google Scholar

15.

J.

Dominy

and

D.

Venturi

, “

Duality and conditional expectations in the Nakajima-Mori-Zwanzig formulation

,”

J. Math. Phys.

58

,

082701

(

2017

).

https://doi.org/10.1063/1.4997015

Google Scholar

Crossref

16.

K.-J.

Engel

and

R.

Nagel

,

One-Parameter Semigroups for Linear Evolution Equations

(

Springer

,

1999

), Vol. 194.

Google Scholar

17.

P.

Español

, “

Dissipative particle dynamics for a harmonic chain: A first-principles derivation

,”

Phys. Rev. E

53

(

2

),

1572

(

1996

).

https://doi.org/10.1103/physreve.53.1572

Google Scholar

Crossref

18.

J.

Florencio

and

H. M.

Lee

, “

Exact time evolution of a classical harmonic-oscillator chain

,”

Phys. Rev. A

31

(

5

),

3231

(

1985

).

https://doi.org/10.1103/physreva.31.3231

Google Scholar

Crossref

19.

R. F.

Fox

, “

Functional-calculus approach to stochastic differential equations

,”

Phys. Rev. A

33

(

1

),

467

–

476

(

1986

).

https://doi.org/10.1103/physreva.33.467

Google Scholar

Crossref

20.

S.

Gudder

, “

A Radon-Nikodým theorem for *-algebras

,”

Pac. J. Math.

80

(

1

),

141

–

149

(

1979

).

https://doi.org/10.2140/pjm.1979.80.141

Google Scholar

Crossref

21.

G. D.

Harp

and

B. J.

Berne

, “

Time-correlation functions, memory functions, and molecular dynamics

,”

Phys. Rev. A

2

(

3

),

975

(

1970

).

https://doi.org/10.1103/physreva.2.975

Google Scholar

Crossref

22.

A.

Karimi

and

M. R.

Paul

, “

Extensive chaos in the Lorenz-96 model

,”

Chaos

20

(

4

),

043105

(

2010

).

https://doi.org/10.1063/1.3496397

Google Scholar

Crossref

PubMed

23.

J.

Kim

and

I.

Sawada

, “

Dynamics of a harmonic oscillator on the Bethe lattice

,”

Phys. Rev. E

61

(

3

),

R2172

(

2000

).

https://doi.org/10.1103/physreve.61.r2172

Google Scholar

Crossref

24.

B. O.

Koopman

, “

Hamiltonian systems and transformation in Hilbert spaces

,”

Proc. Natl. Acad. Sci. U. S. A.

17

(

5

),

315

–

318

(

1931

).

https://doi.org/10.1073/pnas.17.5.315

Google Scholar

Crossref

PubMed

25.

H.

Lei

,

N. A.

Baker

, and

X.

Li

, “

Data-driven parameterization of the generalized Langevin equation

,”

Proc. Natl. Acad. Sci. U. S. A.

113

(

50

),

14183

–

14188

(

2016

).

https://doi.org/10.1073/pnas.1609587113

Google Scholar

Crossref

PubMed

26.

E. N.

Lorenz

, “

Predictability—A problem partly solved

,” in

ECMWF Seminar on Predictability

(

ECMWF

,

1996

), Vol. 1, pp.

1

–

18

.

Google Scholar

Crossref

27.

H.

Mori

, “

A continued-fraction representation of the time-correlation functions

,”

Prog. Theor. Phys.

34

(

3

),

399

–

416

(

1965

).

https://doi.org/10.1143/ptp.34.399

Google Scholar

Crossref

28.

H.

Mori

, “

Transport, collective motion, and Brownian motion

,”

Prog. Theor. Phys.

33

(

3

),

423

–

455

(

1965

).

https://doi.org/10.1143/ptp.33.423

Google Scholar

Crossref

29.

Noise in Nonlinear Dynamical Systems

, Theory of Continuous Fokker-Planck Systems, edited by

F.

Moss

and

E. V. P.

McClintock

(

Cambridge University Press

,

1995

), Vol. 1.

Google Scholar

30.

A.

Pazy

,

Semigroups of Linear Operators and Applications to Partial Differential Equations

(

Springer

,

1992

).

Google Scholar

31.

I.

Snook

,

The Langevin and Generalised Langevin Approach to the Dynamics of Atomic, Polymeric and Colloidal Systems

, 1st ed. (

Elsevier

,

2007

).

Google Scholar

32.

G.

Söderlind

, “

The logarithmic norm. History and modern theory

,”

BIT Numer. Math.

46

(

3

),

631

–

652

(

2006

).

https://doi.org/10.1007/s10543-006-0069-9

Google Scholar

Crossref

33.

P.

Stinis

, “

Stochastic optimal prediction for the Kuramoto–Sivashinsky equation

,”

Multiscale Model. Simul.

2

(

4

),

580

–

612

(

2004

).

https://doi.org/10.1137/030600424

Google Scholar

Crossref

34.

P.

Stinis

, “

A comparative study of two stochastic model reduction methods

,”

Physica D

213

,

197

–

213

(

2006

).

https://doi.org/10.1016/j.physd.2005.11.010

Google Scholar

Crossref

35.

P.

Stinis

, “

Higher order Mori-Zwanzig models for the Euler equations

,”

Multiscale Model. Simul.

6

(

3

),

741

–

760

(

2007

).

https://doi.org/10.1137/06066504x

Google Scholar

Crossref

36.

P.

Stinis

, “

Renormalized Mori–Zwanzig-reduced models for systems without scale separation

,”

Proc. R. Soc. A

471

(

2176

),

20140446

(

2015

).

https://doi.org/10.1098/rspa.2014.0446

Google Scholar

Crossref

37.

L. N. P.

Trefethen

, “

Pseudospectra of linear operators

,”

SIAM Rev.

39

(

3

),

383

–

406

(

1997

).

https://doi.org/10.1137/s0036144595295284

Google Scholar

Crossref

38.

L. N.

Trefethen

and

M.

Embree

,

Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators

(

Princeton University Press

,

2005

).

Google Scholar

Crossref

39.

U.

Umegaki

, “

Conditional expectation in an operator algebra. I

,”

Tohoku Math. J.

6

(

2

),

177

–

181

(

1954

).

https://doi.org/10.2748/tmj/1178245177

Google Scholar

Crossref

40.

D.

Venturi

, “

The numerical approximation of nonlinear functionals and functional differential equations

,”

Phys. Rep.

732

,

1

–

102

(

2018

).

https://doi.org/10.1016/j.physrep.2017.12.003

Google Scholar

Crossref

41.

D.

Venturi

,

H.

Cho

, and

G. E.

Karniadakis

, “

The Mori-Zwanzig approach to uncertainty quantification

,” in

Handbook of Uncertainty Quantification

, edited by

R.

Ghanem

,

D.

Higdon

, and

H.

Owhadi

(

Springer

,

2016

).

Google Scholar

42.

D.

Venturi

and

G. E.

Karniadakis

, “

Convolutionless Nakajima-Zwanzig equations for stochastic analysis in nonlinear dynamical systems

,”

Proc. R. Soc. A

470

(

2166

),

20130754

(

2014

).

https://doi.org/10.1098/rspa.2013.0754

Google Scholar

Crossref

43.

D.

Venturi

,

T. P.

Sapsis

,

H.

Cho

, and

G. E.

Karniadakis

, “

A computable evolution equation for the joint response-excitation probability density function of stochastic dynamical systems

,”

Proc. R. Soc. A

468

(

2139

),

759

–

783

(

2012

).

https://doi.org/10.1098/rspa.2011.0186

Google Scholar

Crossref

44.

Y.

Zhu

and

D.

Venturi

, “

Faber approximation of the Mori-Zwanzig equation

,”

J. Comput. Phys.

372

,

694

–

718

(

2018

).

https://doi.org/10.1016/j.jcp.2018.06.047

Google Scholar

Crossref

45.

R.

Zwanzig

, “

Memory effects in irreversible thermodynamics

,”

Phys. Rev.

124

,

983

–

992

(

1961

).

https://doi.org/10.1103/physrev.124.983

Google Scholar

Crossref

46.

R.

Zwanzig

,

Nonequilibrium Statistical Mechanics

(

Oxford University Press

,

2001

).

Google Scholar

2018

Author(s)

On the estimation of the Mori-Zwanzig memory integral

I. INTRODUCTION

II. THE MORI-ZWANZIG FORMULATION

A. Projection operators

1. Infinite-rank projections

2. Finite-rank projections

III. ANALYSIS OF THE MEMORY INTEGRAL

A. Semigroup estimates

B. Memory growth

C. Short memory approximation and the t-model

D. Hierarchical memory approximation

1. The H-model

2. Type-I finite memory approximation (FMA)

3. Type-II finite memory approximation

4. H_t-model

E. Linear dynamical systems

F. Memory estimates for finite-rank projections and Hamiltonian systems

IV. NUMERICAL EXAMPLES

A. Hamiltonian dynamical systems with finite-rank projections

1. Harmonic chains of oscillators

2. Hald system

B. Non-Hamiltonian systems with infinite-rank projections

1. Linear dynamical systems

2. Nonlinear dynamical systems

a. Lorenz-63 system.

b. Modified Lorenz-96 system.

V. SUMMARY

ACKNOWLEDGMENTS

APPENDIX A: SEMIGROUP BOUNDS VIA FUNCTION DECOMPOSITION

APPENDIX B: THE MORI-ZWANZIG FORMULATION IN PDF SPACE

1. Analysis of the memory integral

2. Hierarchical memory approximation in PDF space

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

On the estimation of the Mori-Zwanzig memory integral

I. INTRODUCTION

II. THE MORI-ZWANZIG FORMULATION

A. Projection operators

1. Infinite-rank projections

2. Finite-rank projections

III. ANALYSIS OF THE MEMORY INTEGRAL

A. Semigroup estimates

B. Memory growth

C. Short memory approximation and the t-model

D. Hierarchical memory approximation

1. The H-model

2. Type-I finite memory approximation (FMA)

3. Type-II finite memory approximation

4. Ht-model

E. Linear dynamical systems

F. Memory estimates for finite-rank projections and Hamiltonian systems

IV. NUMERICAL EXAMPLES

A. Hamiltonian dynamical systems with finite-rank projections

1. Harmonic chains of oscillators

2. Hald system

B. Non-Hamiltonian systems with infinite-rank projections

1. Linear dynamical systems

2. Nonlinear dynamical systems

a. Lorenz-63 system.

b. Modified Lorenz-96 system.

V. SUMMARY

ACKNOWLEDGMENTS

APPENDIX A: SEMIGROUP BOUNDS VIA FUNCTION DECOMPOSITION

APPENDIX B: THE MORI-ZWANZIG FORMULATION IN PDF SPACE

1. Analysis of the memory integral

2. Hierarchical memory approximation in PDF space

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

4. H_t-model