Mori-Zwanzig reduced models for uncertainty quantification

Introduction

The problem of quantifying the uncertainty of the solution of systems of partial or ordinary differential equations has become in recent years a rather active area of research. The realization that more often than not, for problems of practical interest, one is not able to determine the parameters, initial conditions, boundary conditions etc. to within high enough accuracy, has led to a flourishing literature of methods for quantifying the impact that this uncertainty imposes on the solution of the problems under investigation (see e.g. ). However, despite the increase in computational power and the development of various techniques for uncertainty quantification there is still a wealth of problems where reliable uncertainty quantification is beyond reach. One way to address this problem is to look for reduced models for a subset of the variables needed for a complete description of the uncertainty. The effect of all types of uncertainty is intimately connected with the inherent instabilities that may be present in the underlying system which we subject to the uncertainty. These considerations remain equally, if not more, important when we attempt to construct reduced models for uncertainty quantification.

In the current work, we are concerned with the construction of reduced models for systems of differential equations that arise from polynomial chaos expansions of solutions of a PDE or ODE system. In particular, we focus on the case that the given PDE or ODE system contains uncertain parameters or initial conditions and we want to construct a reduced model for the evolution of a subset of the polynomial chaos expansions that are needed for a complete description of the uncertainty caused by the uncertain parameters. There are different methods to construct reduced models for PDE or ODE systems (see e.g. and references therein). We choose to use the Mori-Zwanzig formalism in order to construct the reduced model .

The main issue with all model reduction approaches is the computation of the memory caused by the process of eliminating variables from the given system (referred to as the full system from this point on) . The memory terms are, in general, integral terms which account for the history of the variables that are not resolved. In principle the integrands appearing in the memory terms can be computed through the solution of the orthogonal dynamics equation . We present some examples where this procedure can be implemented and the resulting reduced model can be estimated. Those examples highlight the definite improvement in accuracy of a reduced model when it includes a memory term. However, it is also easy to come up with examples where the solution of the orthogonal dynamics equation becomes prohibitively expensive.

For such cases we present a Markovian reformulation of the MZ formalism which allows the calculation of the memory terms through the solution of ordinary differential equations instead of the computation of convolution integrals as they appear in the original formulation. We present an algorithm which allows the estimation of the necessary parameters on the fly. This means that one starts evolving the full system and use it to estimate the reduced model parameters. Once this is achieved, the simulation continues by evolving only the reduced model with the necessary parameters set equal to their estimated values from the first part of the algorithm. Of course, such an approximation of the memory term cannot work under all circumstances. We present results for a nontrivial problem where it does yield a reduced model with improved behavior compared to a model that ignores the memory terms altogether.

We should note that this alternative approach to computing the memory term fits in the renormalization framework advocated recently by one of the authors in order to construct reduced models for singular PDEs. In particular, the idea is that one embeds the MZ reduced model in a larger family of reduced models which share the same functional form but may have additional parameters for enhanced flexibility. These extra parameters are determined so that the reduced model reproduces some dynamic features of the full system. After this is done, one can switch to the reduced model for the rest of the simulation. In the current work, the extra parameters are the lengths of the memory appearing in the MZ reduced model.

Section 2 presents a brief introduction to the MZ formalism for the construction of reduced models of systems of ODEs. In Section 3 we develop the Markovian reformulation of the MZ formalism and show how one can estimate adaptively the parameters appearing in the reduced model. Section 4 presents numerical results both for the original MZ formalism (Sections 4.1-4.2) and its Markovian reformulation (Section 4.3). Finally, in Section 5 we discuss directions for future work.

Mori-Zwanzig formalism

We begin with a brief presentation of the Mori-Zwanzig formalism . Suppose we are given the system

where $u=(\{u_{k}\}),\;k\in H\cup G$ with initial condition $u(0)=u_{0}.$ The unknown variables (modes) are divided into two groups, one group is indexed in H and the order indexed in G. Our goal is to construct a reduced model for the modes in the set $H.$ The system of ordinary differential equations we are given can be transformed into a system of linear partial differential equations

where $L=\sum_{k\in H\cup G}R_{i}(u_{0})\frac{\partial}{\partial u_{0i}}.$ The solution of (2) is given by $u_{k}(u_{0},t)=\varphi_{k}(u_{0},t)$. Using semigroup notation we can rewrite (2) as

Equation (3) is the Mori-Zwanzig identity. Note that this relation is exact and is an alternative way of writing the original PDE. It is the starting point of our approximations. Of course, we have one such equation for each of the resolved variables $u_{k},k\in H$. The first term in (3) is usually called Markovian since it depends only on the values of the variables at the current instant, the second is called “noise” and the third “memory”.

since $PQ=0$. Also for the initial condition

by the same argument. Thus, the solution of (5) is at all times orthogonal to the range of $P.$ We call (5) the orthogonal dynamics equation. Since the solutions of the orthogonal dynamics equation remain orthogonal to the range of $P$, we can project the Mori-Zwanzig equation (3) and find

We will not present here more details about how to start from Eq. (6) and construct reduced models of different orders for a general system of ODEs. Such constructions have been documented thoroughly elsewhere (see e.g. ). However, we will provide such details for the specific numerical examples in Sections 4.1-4.2.

Markovian reformulation of the MZ formalism

While the MZ model given by Eq. (6) is exact, its construction can be involved and most importantly, very costly. The main source of computational expense is the memory term. Technically, the cost associated with the memory term comes from two sources: i) the presence of the orthogonal dynamics equation solution operator $e^{sQL}$ and ii) the need to find an expression in terms of the resolved variables and time for $PLe^{sQL}QLu_{0k}$ which appears in the memory integrand. The presence of $e^{sQL}$ is problematic because the orthogonal dynamics equation is, for the general case, a PDE in as many dimensions as the original system of ODEs. Also, finding an expression for $PLe^{sQL}QLu_{0k}$ is problematic because, in general, it is not possible to separate the dependence of the expression on time and on the resolved variables. Both are formidable tasks and we will show with several examples how they can increase the cost of constructing the reduced model. For some cases (see e.g. Section 4.1 and 4.2) both tasks can be tackled through the use of a finite-rank projection for the operator $P.$ However, we will show with a simple example (see Section 4.3) that the use of a finite-rank projection may be too costly itself. For such cases, we need an alternative approach to the construction of the memory term. In this section we describe a reformulation of the problem of computing the memory term which can alleviate some of these issues. Also, we present numerical results from the application of this approach in Section 4.3.

We focus on the case when the memory has a finite extent only. The case of infinite memory is simpler and is a special case of the formulation presented below. Also, the current reformulation allows us to comment on what happens in the case when the memory is very short.

Let $w_{0k}(t)=P\int_{0}^{t}e^{(t-s)L}PLe^{sQL}QLu_{0k}ds=P\int_{0}^{t}e^{sL}PLe^{(t-s)QL}QLu_{0k}ds,$ by the change of variables $t^{\prime}=t-s.$ Note, that $w_{0k}$ depends both on $t$ and the resolved part of the initial conditions $\hat{u}_{0}.$ We have suppressed the $\hat{u}_{0}$ dependence for simplicity of notation. If the memory extends only for $t_{0}$ units in the past (with $t_{0}\leq t,$) then

To allow for more flexibility, let us assume that the integrand in the formula for $w_{1k}(t)$ contributes only for $t_{1}$ units with $t_{1}\leq t_{0}.$ Then

We can proceed and write an equation for the evolution of $w_{1k}(t)$ which reads

Similarly, if this integral extends only for $t_{2}$ units in the past with $t_{2}\leq t_{1},$ then

This hierarchy of equations continues indefinitely. Also, we can assume for more flexibility that at every level of the hierarchy we allow the interval of integration for the integral term to extend to fewer or the same units of time than the integral in the previous level. If we keep, say, $n$ terms in this hierarchy, the equation for $w_{(n-1)k}(t)$ will read

Note that the last term in (9) involves the unknown evolution operator for the orthogonal dynamics equation. This situation is the well-known closure problem. We can stop the hierarchy at the $n$th term by assuming that $w_{nk}(t)=0.$

In addition to the closure problem, the unknown evolution operator for the orthogonal dynamics equation appears in the equations for the evolution of the quantities $w_{0k}(t),$ $\ldots,$ $w_{(n-1)k}(t)$ through the various terms $Pe^{(t-t_{0})L}PLe^{t_{0}QL}QLu_{0k},$ $\ldots,$ $Pe^{(t-t_{0})L}PLe^{t_{0}QL}(QL)^{n-1}QLu_{0k}$ respectively.

We describe now a way to express these terms involving the unknown orthogonal dynamics operator through known quantities so that we obtain a closed system for the evolution of $w_{0k}(t),\ldots,w_{(n-1)k}(t).$

Since we want to treat the case where $t_{0}$ is not necessarily small, we divide the interval $[t-t_{0},t]$ in $n_{0}$ subintervals. Define

where $n_{0}\Delta t_{0}=t_{0}$ and $w_{0k}(t)=\sum_{i=1}^{n_{0}}w_{0k}^{(i)}(t).$ Similarly, we can define the quantities $w_{1k}^{(1)}(t),\ldots,w_{1k}^{(n_{1})}(t)$

where $n_{1}\Delta t_{1}=t_{1}$ and $w_{1k}(t)=\sum_{i=1}^{n_{1}}w_{1k}^{(i)}(t).$ In a similar fashion we can define corresponding quantities for all the memory terms up to $w_{(n-1)k}(t)=\sum_{i=1}^{n_{n-1}}w_{(n-1)k}^{(i)}(t).$

In order to proceed we need to make an approximation for the integrals over the subintervals.

2 Trapezoidal rule approximation

By dropping the $O((\Delta t_{0})^{2}),\ldots,O((\Delta t_{n-1})^{2})$ terms we obtain a system of $n_{0}+n_{1}+\ldots+n_{n-1}$ differential equations for the evolution of the quantities $w_{0k}^{(1)}(t),\ldots,w_{(n-1)k}^{(n_{n-1})}.$ This system allows us to determine the memory term $w_{0k}(t).$ Since the approximation we have used for the integral leads to an error $O(\Delta t)^{2},$ the ODE solver should also be $O(\Delta t)^{2}.$ We have used the modified Euler method to solve numerically the equations for the reduced model.

Note that the implementation of the above scheme requires the knowledge of the expressions for $Pe^{tL}PLQLu_{0k},\ldots,Pe^{tL}PL(QL)^{n-1}QLu_{0k}.$ Since the computation of these expressions for large $n$ can be rather involved for nonlinear systems (see Section 4.3), we expect that the above scheme will be used with a small to moderate value of $n.$ Finally, we mention that the above construction can be carried out for integration rules of higher order e.g. Simpson’s rule.

3 Estimation of the memory length

The construction presented above relies on an accurate determination of the memory lengths $t_{0},t_{1},\ldots,t_{n-1}.$ We present in this section a way to estimate these quantities on the fly. This means that we start evolving the full system, use it to estimate $t_{0},t_{1},\ldots,t_{n-1}$ and then switch to the reduced model with the estimated values for $t_{0},t_{1},\ldots,t_{n-1}.$

For simplicity of presentation we assume that we evolve only $w_{0k}(t).$ If we use the trapezoidal rule to discretize $w_{0k}(t)$ and eliminate the term $Pe^{(t-t_{0})L}PLe^{t_{0}QL}QLu_{0k}$ from (7), the reduced model reads

for $k\in H.$ We can solve (13) formally and substitute in (12) to get

where $\lambda_{0}=2/t_{0}.$ Recall that, for the resolved variables, we have from the full system

We would like to estimate the memory decay parameter $t_{0}$ so that the reduced equation (14) for $u_{k}$ reproduces the behavior of $u_{k}$ as predicted by the full system (15). We can do that by requiring that the evolution of some integral quantity of the solution is the same when predicted by the reduced and full systems.

We begin by discretizing the integral term in (14). Suppose that we are evolving the full system with a step size $\delta t,$ where $t=n_{t}\delta t$ (note that $n_{t}$ increases as $t$ increases). If we discretize the integral with the trapezoidal rule we find

where $f_{k}(j\delta t,\hat{u}_{0})=2Pe^{j\delta tL}PLQLu_{0k}$ for $j=0,\ldots,n_{t}.$ The quantities $f_{k}(j\delta t,\hat{u}_{0})$ can be computed from the full system.

There is freedom in the choice of the integral quantity whose evolution the reduced model should be able to reproduce. For example, we can use $\sum_{k\in H}|Pu_{k}(t)|^{2}$ the squared $l_{2}$ norm of the resolved variables. If we use this integral quantity, then from (16) and (15) we find that the unknown parameter $t_{0}$ must satisfy

Let $y=\exp[-\lambda_{0}\delta t].$ Then,

With this identification, equation (17) becomes a polynomial equation for $y$ with $y\in.$ It is not difficult to solve equation (17) with an iterative method, for example Newton’s method. For the numerical results we present in Section 4.3, Newton’s method converged to double precision accuracy within 4-5 iterations. After an estimate $\hat{y}$ has been obtained, we can find the estimate $\hat{t}_{0}$ of $t_{0}$ (recall $\lambda_{0}=2/t_{0}$) from

For each time instant $t$ we can obtain through equations (17) and (19), an estimate $\hat{t}_{0}(t)$ for $t_{0}.$ Thus, the most important issue that we have to address is that of deciding which is the best estimate of $t_{0}.$ In other words, at what time $t_{f}$ should we stop estimating the value of $t_{0}$ so that we can use the estimated value $\hat{t}_{0}(t_{f})$ to evolve the reduced model from then on.

We define $\epsilon(t)=\underset{l\in[1,n_{t}]}{\max}|\hat{y}^{l}(t+\delta t)-\hat{y}^{l}(t)|.$ The quantity $\epsilon(t)$ monitors the convergence of not only the value of the estimate $\hat{y}$ as a function of the time $t$, but of the whole function $e^{-\lambda_{0}(t-s)}.$ Ideally, $\epsilon(t)$ converges to zero with increasing $t.$ That will be the case if the approximation of the memory term only through $Pe^{tL}PLQLu_{0kr}$ is enough (see (12)-(13)). However, this will not always be the case. If keeping $Pe^{tL}PLQLu_{0kr}$ is not enough, then $\epsilon(t)$ will decrease with increasing $t$ up to some time $t_{min}$ when it will reach a nonzero minimum. After that time, it starts increasing. This signals that keeping only $Pe^{tL}PLQLu_{0kr}$ is not enough to describe accurately the memory.

In order to proceed we have two options: (i) construct a higher order model and (ii) identify $t_{f}=t_{min}$ and thus $\hat{t}_{0}(t_{f})=\hat{t}_{0}(t_{min}).$ Results for higher order models will be presented elsewhere (see also discussion in Section 5). In the numerical experiments we present in the next section we have chosen $\hat{t}_{0}(t_{f})=\hat{t}_{0}(t_{min}).$ Note that the procedure just outlined allows the automation of the algorithm. This means that there is no adjustable reduced model parameter that needs to be specified at the onset of the algorithm.

We are now in a position to state the adaptive Mori-Zwanzig algorithm which constructs a reduced model with the necessary memory term parameter $t_{0}$ estimated on the fly.

Evolve the full system and compute, at every step, the estimate $\hat{t}_{0}(t).$ Use estimates of $t_{0}$ from successive steps to calculate $\epsilon(t)=\underset{l\in[1,n_{t}]}{\max}|\hat{y}^{l}(t+\delta t)-\hat{y}^{l}(t)|.$

When $\epsilon(t)$ reaches a minimum (possibly non zero) value at some instant $t_{min}$, pick $\hat{t}_{0}(t_{min})$ as the final estimate of $t_{0}.$

For the remaining simulation time ($t>t_{min}$), switch from the full system to the reduced model. The reduced model is evolved with the necessary parameter $t_{0}$ set to its estimated value $\hat{t}_{0}(t_{min}).$

This procedure can be extended to the computation of optimal estimates for $t_{1},t_{2},\ldots,$ i.e. when we evolve, in addition to $w_{0k}(t),$ the quantities $w_{1k}(t),w_{2k}(t),\ldots.$ Results for such higher order models will be presented elsewhere.

Numerical Examples

Consider the following linear ordinary equation with an uncertain coefficient

where $\kappa\sim U$. This equation has the solution $u=u^{\circ}exp(-\kappa t)$. To represent the dependence of the solution of (20) on $\kappa,$ we can expand it in a general polynomial chaos (gPC) expansion , say using Legendre polynomials. Let $u(t,\cdot)\approx\sum_{i=0}^{M}u_{i}(t)\phi_{i}(\xi)$, where $\xi\sim U$ and $\{\phi_{i}\}$ are normalized Legendre polynomials which are orthonormal with respect to the uniform distribution of $\xi$, i.e.,

We can write $\kappa$ as $\kappa=\frac{1}{2}\xi+\frac{1}{2}=\sum_{i=0}^{1}{k_{i}}\phi_{i}(\xi)$. We substitute this expansion in (20) and obtain (through Galerkin projection) the (truncated) system up to order $M$

where $e_{ijk}=\int_{-1}^{1}\phi_{i}(\xi)\phi_{j}(\xi)\phi_{k}(\xi)\frac{1}{2}d\xi$ and $u_{00}=u^{\circ}$, $u_{0r}=0$ for $r=1,\dots,M$ (for details, see e.g ).

where $h^{\nu}(\hat{u}_{0})$ are tensor product Hermite polynomials up to some order $p$, $\nu$ is the multi-index $\nu=(\nu_{0},\dots,\nu_{\Lambda})$ with $|\nu|=\sum_{i=0}^{\Lambda}\nu_{i}$ and $I$ is the index set up to order $p$, i.e., I=\{\mu\big{|}|\mu|\leq p\}. The order $p$ for the basis functions was set to 5 for a total of 21 basis functions. In formula (22) the inner product is defined as

For each $j\leq\Lambda$, the component $F_{j}(u_{0},t)$ denotes the solution of the orthogonal dynamics

Eq. (27) is a Volterra integral equation for the function $a^{\nu}_{j}(t)$, which can be rewritten as follows:

The functions $f^{\nu}_{j}(t)$, $g^{\mu\nu}(t)$ can be found by averaging over a collection of experiments or simulations, with initial conditions drawn from the initial distribution. In this example, we use a sparse grid quadrature rule for the multi-dimensional integrals .

Finally, we perform one more projection to eliminate the noise term (see Section 2) and the memory term becomes

After calculating $a_{i}^{\mu}$ and $\gamma^{\mu\nu}$ we obtain the following reduced system,

here $A$ and $\Gamma$ are the matrix form of $a^{\mu}_{i}$ and $\gamma^{\mu\nu}$, $\hat{u}_{0}$ is the initial condition of resolved variables.

Fig. 1 shows the evolution of the memory kernel $(Le^{tQL}QLu_{1},h^{01})$ which is indicative of the behavior of the memory kernels. The basis function $h^{01}$ is the product of the zero order Hermite polynomial in the variable $u_{0}$ and the first order Hermite polynomial in the variable $u_{1}.$ We see that the memory kernel is rather slowly decaying which means that the resulting reduced order model will have a long memory. Fig. 2 shows the solution for the resolved variables as predicted by the full system and two different reduced order models, the Markovian model which results from dropping the memory term in (29) and the non-Markovian reduced model given by (29). It is obvious from Fig. 2 that the Markovian model loses accuracy quickly. On the other hand, the non-Markovian model retains its accuracy for the length of the simulation interval. This difference in behavior is quantified in Fig. 3 where we see that for both resolved variables the relative error of the Markovian model becomes greater than $50\%$ by the end of the simulation interval. On the other hand, the error of the non-Markovian model remains less than $1\%$ for the whole simulation interval.

2 Nonlinearly damped and randomly forced particle

Consider the following equation describing a particle moving in a double well potential and driven by a force term (see )

where $f=\sin(t+t_{0})\xi$ and $\xi\sim U$. We use order $M=6$ polynomials in $\xi$ to approximate the full system solution up to time $10$. We want to construct a reduced model for the first 2 coefficients of the polynomial expansion ($\Lambda=1$). As before, we let $u(t,\xi)\approx\sum_{i=0}^{M}u_{i}(t)\phi_{i}(\xi)$ and we obtain through Galerkin projection the system

As can be seen from Fig. 5, the difference between the (memoryless) Markovian and non-Markovian reduced models is even more pronounced than in the case of the linear ODE. The inclusion of the memory term is indeed crucial for maintaining the accuracy of the reduced model for long times. For the case of the resolved variable $u_{1},$ the relative error spikes at a couple of points even for the otherwise very accurate non-Markovian reduced model. As can be seen from Fig. 4, this is because the exact value of $u_{1}$ becomes zero at these points so that the relative error becomes very large even for an accurate approximation. However, the significant improvement in accuracy with the inclusion of the memory term is evident in Fig. 5 which plots the error in a logarithmic scale.

3 Viscous 1D Burgers with uncertain initial conditions

In this section we show how the above MZ formulation can be used for uncertainty quantification for the one-dimensional Burgers equation with uncertain initial condition. As is explained at the end of this section, the calculation of the MZ memory term cannot proceed as for the last two examples. The reason is that it is prohibitively expensive due to the number of basis functions needed. Thus, we will apply the alternative construction that was presented in Section 3.

where $\nu>0.$ Equation (32) should be supplemented with an initial condition $u(x,0)=u_{0}(x)$ and boundary conditions. We solve (32) in the interval $[0,2\pi]$ with periodic boundary conditions. This allows us to expand the solution in Fourier series

where $F=[-\frac{N}{2},\frac{N}{2}-1].$ The equation of motion for the Fourier mode $u_{k}$ becomes

We assume that the initial condition $u_{0}(x)$ is uncertain (random) and can be expanded as $u_{0}(x,\xi)=(\alpha_{0}+\alpha_{1}\xi)v_{0}(x)$ where $\xi$ is uniformly distributed in $$and$v_{0}(x)$a given function. In the numerical experiments we have taken$\alpha_{0}=\alpha_{1}=1$and$v_{0}(x)=\sin x.$Thus, the initial condition varies “uniformly” between the functions 0 and$2\sin x.$

To proceed we expand the solution $u_{k}(t,\xi)$ for $k\in F$ in a polynomial chaos expansion using the standard Legendre polynomials which are orthogonal in the interval $.$ In particular, we have that

where $\phi_{i}(\xi)$ is the standard Legendre polynomial of order $i.$ For each wavenumber $k$ we expand the solution $u_{k}(t,\xi)$ of (36) in Legendre polynomials and keep the first $M$ polynomials

Similarly, the initial condition can be written as $u_{0}(x,\xi)=\sin x\sum_{i=0}^{1}\alpha_{i}\phi_{i}(\xi)$ since $\phi_{0}(\xi)=1$ and $\phi_{1}(\xi)=\xi.$ Substitution of (34) in (33) and use of the orthogonality property of the Legendre polynomials gives

where the expectation $E[\cdot]$ is taken with respect to the uniform density on $.$ The Legendre polynomial triple product integral defines a tensor which has the following sparsity pattern: $E[\phi_{l}(\xi)\phi_{m}(\xi)\phi_{r}(\xi)]=0,$ if $l+m<r$ or $l+r<m$ or $m+r<l$ or $l+m+r=\text{odd}$ . Due to this sparsity pattern, for a given value of $M$ only about $1/4$ of the $M^{3}$ tensor entries are different from zero.

Before we proceed we have to comment on the cost of applying the MZ formalism to construct a reduced model. We have set the viscosity coefficient to $\nu=0.03.$ The solution of the full system was computed with $N=196$ Fourier modes ($F=$) and the first 7 Legendre polynomials ($M=7$). The first 7 Legendre polynomials were enough to obtain converged statistics for the full system. We want to construct reduced models for the evolution of the coefficients of the first 2 Legendre polynomials i.e., $u_{k0},u_{k1}$ for $k\in F.$ If we want to apply the MZ formalism in the way we did for the previous two examples (employing a finite-rank projection etc.) we would need to construct a basis in $2\times 98$ dimensions (exploiting the fact that the solution of the Burgers equation is real-valued). Any attempt to use basis functions up to a high order is infeasible for such a high-dimensional situation. We have attempted to use only low order basis functions but they are not enough to guarantee accuracy of the reduced model. Thus, we turn to the reformulated reduced model that was presented in Section 3.

To conform with the Mori-Zwanzig formalism we set

where $u=\{u_{kr}\}$ for $k\in F$ and $r=0,\ldots,M-1.$ Thus, we have

The Markovian term has the same functional form as the RHS of the full system but is restricted to a sum over only the first $\Lambda$ Legendre expansion coefficients for each Fourier mode.

Finally, to implement any method to solve equation (17) for the estimation of $t_{0}$ we need to specify the RHS of the equation (17). This requires the evaluation of the expression $Pe^{tL}QLu_{0kr}.$ For the case of the viscous Burgers equation, we find

Note that since we restrict attention to initial conditions for which the unresolved variables are zero and the projection sets the unresolved variables to zero, the quantity $Pe^{tL}QLu_{0kr}$ can be computed through the evolution of the full system (36).

The full system was solved with the modified Euler method with $\delta t=0.001.$ The reduced model uses $N=196$ Fourier modes but only the first two Legendre polynomials, so $\Lambda=2.$ It was solved using the modified Euler method with $\delta t=0.001.$ The parameter $t_{0}$ needed for the evolution of the memory term was found to be 0.3783 through the procedure described in Section 3.3.1.

Figure 6 shows the evolution of the mean energy of the solution

as computed from the full system (with $M=7$ Legendre polynomials), the MZ reduced model with $\Lambda=2$ without memory (keeping only the Markovian term) and the MZ reduced model with $\Lambda=2$ with memory. Figure 7 shows the evolution of the standard deviation of the energy of the solution. The variance of the energy is given by

The reduced model performs equally well with or without memory. Of course, the reduced model with memory is slower than the reduced model without memory. However, the reduced model with memory is still about 4 times faster than the the full system.

Figure 8 shows the evolution of the mean squared $l_{2}$ norm of the gradient of the solution

as computed from the full system (with $M=7$ Legendre polynomials), the MZ reduced model with $\Lambda=1$ without memory (keeping only the Markovian term) and the MZ reduced model with $\Lambda=1$ with memory. Figure 9 shows the evolution of the standard deviation. The variance is given by

The large values of the standard deviation of the mean squared $l_{2}$ norm of the gradient are justified by the uncertainty in the initial condition. Recall that we have chosen an initial condition which can vary “uniformly” between the functions 0 and $2\sin x.$ As a result, the standard deviation is large because it has to account for a wide range of possible initial conditions.

It is evident from the figures that the inclusion of the memory term improves the performance of the reduced model. Also, it is evident that there is room for improvement of the reduced model with memory. In particular, more terms are needed in the reformulated MZ model to approximate better the memory.

Recall that the solution of Burgers equation is a contraction . Eventually, the complete description of the uncertainty caused by the uncertainty in the initial condition requires only a few polynomial chaos expansion coefficients. This happens at a time scale that is dictated by the magnitude of the viscosity coefficient. That is why for long times the reduced model with and without memory have comparable behavior to that of the full system. However, for short times, the inclusion of the memory term does make a difference because information from the higher chaos expansion coefficients is needed. The higher chaos expansion coefficients will have a more prolonged contribution for systems that possess unstable modes. In such cases, the inclusion of the memory term becomes imperative for short as well long times. Results for such cases will be presented elsewhere.

Discussion and future work

We have examined the application of the Mori-Zwanzig formalism to the problem of constructing reduced models for uncertainty quantification. In particular, we have constructed reduced models for subsets of the polynomial chaos expansion coefficients needed to describe fully the uncertainty. We have examined cases of parametric or initial condition uncertainty. The main conclusion from the current work is that while the MZ formalism can be applied for the construction of reduced models, the task of constructing an efficient (or even feasible) reduced model can be involved. For cases where the straightforward application of the MZ formalism is not possible, we have offered an alternative construction. The implementation of this alternative construction is reminiscent of renormalization constructions used to describe the evolution of complex solutions of PDEs .

The current work opens several directions for future work. First, we should investigate whether there is a more economical way of choosing the basis functions for cases when the basis functions have many arguments (as was the case for the Burgers example). This is important because the calculation of the memory kernels through the finite-rank projection is well defined and the solution of the corresponding Volterra equations can be performed with high accuracy. A related question is whether there is sparsity in the coefficients of the basis functions. It is plausible that even though in principle the number of basis functions to reach a specific order may be very large, many of them may not contribute to the representation. A related approach would be the use of machine learning algorithms to obtain a more efficient representation of the memory term. Finally, a related issue to be investigated is how to ensure the stability of the reduced model when the finite-rank projection is employed. For example, for the nonlinearly damped and forced particle case, we had to assign smaller variances for the higher coefficients to stabilize the reduced model. This procedure needs to be investigated and, if possible, automated.

A second interesting research direction has to do with the representation of the memory term when the finite-rank projection is not possible due to a prohibitively large number of basis functions. We have explored here an expansion of the memory term that involves, in essence, a Taylor expansion of the orthogonal dynamics operator. Such an expansion seems more plausible when the timescale of the orthogonal dynamics is slower than that of the resolved variables. However, there is an alternative way of performing the expansion of the memory term that is more suited to the case when the orthogonal dynamics is faster than the resolved variables. Such an expansion leads to a Taylor expansion of the whole memory term, not just the orthogonal dynamics operator. If the memory kernel becomes insignificant after a time interval $t_{0},$ then one can use the full system up to time $t_{0},$ estimate the Taylor expansion of the whole memory term around time $t_{0}$ and then switch to the reduced model with the memory given by the Taylor expansion. We will investigate this alternative memory representation and report the results elsewhere.

Acknowledgements

The authors would like to thank D. Barajas-Solano, H. Lei and A. Tartakovsky for useful discussions and comments. This research at Pacific Northwest National Laboratory (PNNL) was partially supported by the U.S. Department of Energy (DOE) Office of Advanced Scientific Computing Research (ASCR) Collaboratory on Mathematics for Mesoscopic Modeling of Materials (CM4), under Award Number DE-SC0009280 and partially by the U.S. DOE ASCR project “Uncertainty Quantification For Complex Systems Described by Stochastic Partial Differential Equations”. PNNL is operated by Battelle for the DOE under Contract DE-AC05-76RL01830.

References