# Optimal Combination of EEFs for the Model Reduction of Nonlinear Partial Differential Equations

Authors: Jun Shuai; Xuli Han
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

Proper orthogonal decomposition is a popular approach for determining the principal spatial structures from the measured data. Generally, model reduction using empirical eigenfunctions (EEFs) can generate a relatively low-dimensional model among all linear expansions. However, the neglectful modes representing only a tiny amount of energy will be crucial in the modeling for certain type of nonlinear partial differential equations (PDEs). In this paper, an optimal combination of EEFs is proposed for model reduction of nonlinear partial differential equations (PDEs),  obtained by the basis function transformation from the initial EEFs. The transformation matrix is derived from straightforward optimization techniques. The present new EEFs can keep the dynamical information of neglectful modes and generate a lower-dimensional and more precise dynamical system for the PDEs. The numerical example shows its effectiveness and feasibility for model reduction of the nonlinear PDEs.


---

## Body

## 1. Introduction

Many problems in science and engineering are reduced to a set of nonlinear partial differential equations (PDEs) [1] through a process of mathematical modeling. In this paper, we focus on the modeling of the dynamics for the following kind of nonlinear PDEs:
(1)∂u(x,t)∂t=Qu(x,t)+F(u(x,t),∂u(x,t)∂x,…,x,t).Equation (1) will satisfy many types of the boundary conditions. It can use the following description to denote boundary conditions containing the Dirichlet boundary condition and Neumann boundary condition of the nonlinear PDEs (1):
(2)B1u(α,t)+C1∂u(α,t)∂x=D1,B2u(β,t)+C2∂u(β,t)∂x=D2.
The initial conditions are also given:
(3)u(x,0)=u0(x),
where u(x,t) denotes the spatiotemporal state variable, x∈Ω=[α,β]⊂R is the spatial coordinate, and t∈[0,∞) is the time variable. Q is a linear operator involving linear spatial derivatives on the state variable. F(u(x,t),∂u(x,t)/∂x,…,x,t) is a nonlinear function which contains spatial derivatives of the state variable. B1, C1, D1, B2, C2, and D2 are constants and u0(x) is a function of x.Generally, two kinds of methods are mainly used to model the dynamics of nonlinear PDEs: conventional discretization approaches and advanced methods based on spatial basis function expansions [2–4]. However, conventional discretization approaches such as finite difference method (FDM) [5] and finite element method (FEM) [6] often lead to high-order dynamical models of nonlinear PDEs. The advanced methods based on spatial basis function expansions such as spectral methods [2, 3] can derive a low-order ODE system to model the dynamics of nonlinear PDEs.Empirical eigenfunctions (EEFs) [7, 8] which are identified by proper orthogonal decomposition are one of the spatial basis functions used for the time/space separation.With the help of the traditional identification technique [9], EEFs are widely used for PDE model identification. It is well known that the amount of energy of a system represented by the EEFs is often taken as an indication of the quality of a reduced model using those first several EEFs. However, the linear approximation for the nonlinear system cannot guarantee the assumption that minor components do not contain important information [10]. Modes representing only a tiny amount of energy can also be crucial in the generation of certain types of dynamics. For instance, the EEF-based models can have difficulties in reproducing behavior dominated by irregular transitions between different dynamical states, which have been pointed out in the past studies [11–13]. For the more precise approximation of the PDEs, it is necessary to keep the influence of the neglectful modes representing only a tiny amount of energy on retained dominated modes. Adding the dynamical information by the combination of the initial EEFs will give both qualitatively and quantitatively better results than the common EEFs-based models [8].In this paper, an optimal combination of EEFs is proposed for model reduction of nonlinear partial differential equations (PDEs), which is obtained by the basis function transformation from the initial EEFs. High-order EEFs are initially obtained by the proper orthogonal decomposition of the measured data. Thus, a new set of basis functions is constructed by the basis function transformation of the initial EEFs. The transformation matrix is derived from straightforward optimization techniques. The present new EEFs can keep the dynamical information of neglectful modes and generate a lower dimensional and more precise dynamical system for the PDEs. The numerical example shows its effectiveness and feasibility for model reduction of the nonlinear PDEs.


## 2. EEFs for Model Reduction

The proper orthogonal decomposition (POD) [14], also known as the Karhunen-Loeve (K-L) decomposition, principal component analysis (PCA), and empirical orthogonal function analysis (EOF) can facilitate the modal projections of partial differential equations into reduced-order models, via a Galerkin projection [15]. The POD is a mean of extracting spatial structure from a set of time series available on a domain. The most attracting property of the POD is that it can minimize the average squared distance between the original measured signal and its reduced representation. POD allows the identification of a useful set of empirical eigenfunctions and the dimension of the subspace necessary to achieve a satisfactory approximation of nonlinear PDEs.For simplicity, assume the measured spatial-temporal variable{u(xi,tj)}i=1,t=1n,K (called snapshots) is uniformly sampled in the space and time coordinates. Define the inner product, norm, and ensemble average as (f(x),g(x))=∫Ωf(x)g(x)dx, ∥f(x)∥=(f(x),f(x))1/2, and 〈f(x,t)〉=1/K∑t=1Kf(x,t). Motivated by the Fourier series, the spatiotemporal variable u(x,t) can be expanded onto an infinite number of orthonormal spatial basis functions ϕi(x) with temporal coefficients ui(t):
(4)u(x,t)=∑i=1∞ui(t)ϕi(x).
The spatial basis functions are orthonormal; the following are derived:
(5)(ϕi(x),ϕj(x))=∫Ωϕi(x)ϕj(x)dx={0,i≠j,1,i=j.
The temporal coefficients can be calculated from
(6)ui(t)=(ϕi(x),u(x,t)),i=1,…,∞.If a large enoughM is given (i.e., M→∞), u(x,t) can be approximated by the following M-order approximation, which is truncated to be a finite dimension:
(7)u(x,t)≈uM(x,t)=∑i=1Mui(t)ϕi(x),
where uM(x,t) denotes the M-order approximation.The main problem of using proper orthogonal decomposition for time/space separation is to obtain the most characteristic spatial structure{ϕi(x)}i=1N from the spatiotemporal snapshots {u(xi,tj)}i=1,t=1n,K. Finding the typical {ϕi(x)}i=1N can be performed by using the following eigenvalue problems:
(8)∫ΩR(x,μ)ϕi(μ)dμ=λiϕi(x),(ϕi,ϕi)=1,
where R(x,μ)=〈u(x,t)u(μ,t)〉  denotes the spatial two-point correlation function.Since the data are always discrete in space, one must solve numerically the integral equation (8). Discretizing the integral equation gives an n×n matrix eigenvalue problem. Thus, at most eigenfunctions at n sampled spatial locations can be obtained. Then one can interpolate the eigenfunctions to locations where the data are not available.Assume that the maximum number of nonzero eigenvalues isM≤n. Let λ1>λ2>⋯>λM be the eigenvalues with the corresponding eigenfunctions {ϕ1(x),ϕ2(x),…,ϕM(x)}. The eigenfunction that corresponds to the first eigenvalue is considered to be the most “energetic.” The total “energy” is defined as being the sum of the eigenvalues. To each eigenfunction, assign an “energy” percentage based on the associated eigenvalue:
(9)Ei=λi∑j=1Mλj.Usually, the sufficient number of eigenfunctions that capture 99% of the system “energy” is used to determine the value ofN. Experiences show that only the first few basis functions can represent the dominant dynamics of many PDE systems.As to the model reduction for PDE (1), the spatiotemporal variables u(x,t) can be approximated by expanding onto the first N spatial-dependent EEFs {ϕ1(x),ϕ2(x),…,ϕN(x)} with the corresponding temporal coefficients ui(t):
(10)u(x,t)≈uN(x,t)=∑i=1Nui(t)ϕi(x).
For PDE (1), the following can be derived:
(11)∑i=1Nu˙i(t)ϕi(x)=Q(∑i=1Nui(t)ϕi(x))+F(∑i=1Nui(t)ϕi(x),∑i=1Nui(t)∂ϕi(x)∂x,…,x,t).
Via a Galerkin projection, the following N-order dynamical systems are then obtained:
(12)u˙i(t)=λiui(t)+fi(u(t),t),
where λi,i=1,2…,N denote the eigenvalues, fi(u(t),t) denote the nonlinear terms, and u(t)=[u1(t),u2(t),…,uN(t)]T. The series of dynamical systems can be rewritten in a general form as follows:
(13)u˙(t)=Lu(t)+f(u(t),t),
where L=diag(λ1,λ2,…,λN);  f(u(t),t)=[f1(u(t),t),f2(u(t),t),…,fN(u(t),t)]T  and fi(u(t),t) = (ϕi(x),F(∑i=1Nui(t)ϕi(x),∑i=1Nui(t)(∂ϕi(x)/∂x),…,x,t)).


## 3. Line Combination of EEFs

Considering the influence of neglectful EEFs with the tiny amount energy, a new kind of EEFs can be derived by line combinations of the initial EEFs. As to obtain the lower-order dynamical systems of nonlinear PDEs, a smaller set of new basis functions can be derived. Each new basis function is the line combination of theM EEFs from POD:
(14)φj(x)=a1jϕ1(x)+a2jϕ2(x)+⋯+aMjϕM(x),=a1jϕ1(x)+aaaaa2jϕ2(x)+⋯j=1,2,…,N,
where N<M. Equation (14) can also be rewritten as follows:
(15)[φ1(x),φ2(x),…,φN(x)]=[ϕ1(x),ϕ2(x),…,ϕM(x)]A,
where the transformation matrix A=[a11a12⋯a1Na21a22⋯a2N⋮⋮⋱⋮aM1aM2⋯aMN].Using spatial basis functions expansion based on the set of new basis functions (14), the approximation of spatiotemporal variables u(x,t) can be expanded as follows:
(16)u(x,t)≈∑i=1Nu-i(t)φi(x).
Combining (7), (15), and (16), the time coefficients have
(17)u-(t)=A+u(t),
where u(t)=[u1(t),u2(t),…,uM(t)]T;  u-(t)=[u-1(t),u-2(t),…,u-N(t)]T, where A+ denotes the generalized inverse of transformation matrix A.For PDE (1), the following is derived:
(18)∑i=1Nu-˙i(t)φi(x)=Q(∑i=1Nu-i(t)φi(x))+F(∑i=1Nu-i(t)φi(x),∑i=1Nu-i(t)∂φi(x)∂x,…,x,t).Using the Galerkin method, a lower-dimensional dynamical system of nonlinear PDE (1) in a general form can be derived as follows:
(19)u-˙(t)=L-u-(t)+f-(u-(t),t),
where u-(t)=[u-1(t),u-2(t),…,u-N(t)]T;
(20)L-={L-ij},i,j=1,2,…,N;L-ij=(φi,L(φj));f-(u-(t),t)=[f-1(u-(t),t),f-2(u-(t),t),…,f-N(u-(t),t)]T;f-i(u-(t),t)=(∑i=1Nu-i(t)∂φi(x)∂xφi(x),F(∑i=1Nu-i(t)φi(x),mmmmmmm∑i=1Nu-i(t)∂φi(x)∂x,…,x,t)).


## 4. Optimization for the Transformation Matrix

An error function that measures energy error between the ODE system (13) and the lower-dimensional dynamical system (19) is introduced by
(21)Error=∫0∞(∑i=1M(ui(t))2-∑i=1N(u-i(t))2)2dt.
Combining the vector expressions of u(t),  u-(t) with (17), the error (21) can be strictly derived as follows:
(22)Error=∫0∞(∑i=1M(ui(t))2-∑i=1N(u-i(t))2)2dt=∫0∞(u(t)Tu(t)-u-(t)Tu-(t))2dt=∫0∞(u(t)Tu(t)-u(t)T((A+)TA+)u(t))2dt=∫0∞(u(t)T(I-(A+)TA+)u(t))2dt.The combination matrixA is determined by minimizing the error function (22). To evaluate the error function, a temporal interval [0,Tmax] is considered for (22) and approximated by a finite sum, where the max integral time Tmax is selected as free parameter. The temporal integral interval [0,Tmax] generally contains the processes from initial state to steady state of the nonlinear PDEs. The restricted condition for the minimizing problem (22) is that the combination matrix should be column orthogonal. Thus, the error functions for combination matrix can be approximated as follows.Since(23)Error(A)≈∫0Tmax(u(t)T(I-(A+)TA+)u(t))2dt≈(∑i(ti-ti-1)(u(ti)T(I-(A+)TA+)u(ti)))2,
we consider approximate error
(24)E(A)=(∑i(ti-ti-1)(u(ti)T(I-(A+)TA+)u(ti)))2,mmmmmmmmmmsubjecttoA(:,i)TA(:,j)=kiδij.The optimization ofE(A) (24) mainly contains several steps. Firstly, the fourth-order Runge-Kutta method is used to calculate the ODE system (13); thus ti, u(ti) of Error (24) are known.Secondly, the particle swarm optimization (PSO) algorithm [16] is used to optimize the restricted error functions (24). PSO is an efficient, robust, and simple optimization algorithm for solving many optimization problems. In nature, PSO is a stochastic optimization approach which maintains a swarm of candidate solutions, referred to as particles. The method is inspired by the movement of particles and their interactions with their neighbors in the group. The principle and details of PSO algorithm are given in [17].For the optimization of (24) by PSO directly, the obtained combination matrix A~ is not orthogonal. To keep the minimization of the error function (24), a satisfied combination matrix A can be derived by singular value decomposition for A~. Let the singular value decomposition of A~ be A~=USVT, where U,  S are orthogonal and diagonal matrices, respectively. Thus combination matrix A=US is orthogonal and it can be derived that (A~+)TA~+=US(STS)-1(STS)-1STUT=(A+)TA+, which means that the minimum of the error function (24) is retained.


## 5. Numerical Example

In this section, the optimal combination of EEFs is applied for model reduction of a simple PDE model. Note that this example serves as an illustration of the proposed methodology rather than a demonstration of the reduction of a very large and complex system.In order to evaluate the performance of optimal combination of EEFs for model reduction of nonlinear PDEs, the Burgers equation [18, 19] is studied. The Burgers equation is a one-dimensional spatial model of a variety of three-dimensional physical phenomena, greatly simplifying the problems while retaining many of the complex behavior characteristics. This equation contains nonlinear convection and diffusion terms and retains many of the interesting features of the Navier-Stokes equation. The governing equation may be written as follows:
(25)∂u(x,t)∂t+u∂u(x,t)∂x=ε∂2u(x,t)∂x2,
where the viscosity ε is 0.2 and (25) subjects to the following boundary conditions:
(26)u(0,t)=u(2π,t)=0,
and an initial condition:
(27)u(x,0)=sin(x)+sin(2x)+sin(3x).Thirty-two sensors uniformly distributed in the space are used for measurement. A noise-free dataset of 300 data is collected from (25). The sampling interval Δt is 0.01 s and the simulation time is 3 s.


### 5.1. Comparisons with the EEFs

The EEFs are obtained from the singular value decomposition for the sampling data. The energy values of the dynamical systems based on the EEFs are shown in Table1.Table 1
Energy value of the dynamical systems based on EEFs.
EEFs
1 order
2 order
3 order
4 order
5 order
6 order
7 order
8 order
Energy value
0.0026
0.0074
0.0220
0.0340
0.1176
0.1245
0.1279
0.1281Using spatial basis functions expansion based on the new basis functions and the Galerkin methods, low-dimensional dynamical systems can be derived. Meanwhile, the energy functions of the obtained dynamical systems are calculated. The results in Table2 show that the low-dimensional dynamical system based on the first several optimal combinations of EEFs can capture the dynamics of the nonlinear PDE (25) perfectly.Table 2
Energy value of the dynamical systems based on new basis functions.
New basis functions
1 order
2 order
3 order
Energy value
0.1277
0.1280
0.1281


### 5.2. Low-Dimensional Approximation Based on 3 Optimal Combined EEFs

Using the line combinations of the EEFs and optimization for the energy error functions of the line combination coefficients, the optimal combination of EEFs for model reduction of nonlinear PDEs is obtained, where the first 3 initial EEFs and 3 optimal combined EEFs are shown in Figures1 and 2.Figure 1
The first 3 initial EEFs.Figure 2
The first 3 optimal combined EEFs.A new set of data is collected for testing to demonstrate the spatial-temporal performance of the two kinds of EEFs. For practical industrial control, the low-dimensional model based on optimal combined EEFs is very useful because it cannot only capture the dominant dynamics of the nonlinear PDEs but simplify the control design as well. The spatiotemporal errors of the predicted output of two models combined with two kinds of EEFs are shown in Figures3 and 4, respectively. As given in Figures 3 and 4, the predicted distributed error of the approximate model based on the optimal combination of EEFs is shown to be significantly smaller than that of the same order approximate model based on initial EEFs.Figure 3
Predicted error based on 3 initial EEFs.Figure 4
Predicted error based on 3 optimal combined EEFs.


## 5.1. Comparisons with the EEFs

The EEFs are obtained from the singular value decomposition for the sampling data. The energy values of the dynamical systems based on the EEFs are shown in Table1.Table 1
Energy value of the dynamical systems based on EEFs.
EEFs
1 order
2 order
3 order
4 order
5 order
6 order
7 order
8 order
Energy value
0.0026
0.0074
0.0220
0.0340
0.1176
0.1245
0.1279
0.1281Using spatial basis functions expansion based on the new basis functions and the Galerkin methods, low-dimensional dynamical systems can be derived. Meanwhile, the energy functions of the obtained dynamical systems are calculated. The results in Table2 show that the low-dimensional dynamical system based on the first several optimal combinations of EEFs can capture the dynamics of the nonlinear PDE (25) perfectly.Table 2
Energy value of the dynamical systems based on new basis functions.
New basis functions
1 order
2 order
3 order
Energy value
0.1277
0.1280
0.1281


## 5.2. Low-Dimensional Approximation Based on 3 Optimal Combined EEFs

Using the line combinations of the EEFs and optimization for the energy error functions of the line combination coefficients, the optimal combination of EEFs for model reduction of nonlinear PDEs is obtained, where the first 3 initial EEFs and 3 optimal combined EEFs are shown in Figures1 and 2.Figure 1
The first 3 initial EEFs.Figure 2
The first 3 optimal combined EEFs.A new set of data is collected for testing to demonstrate the spatial-temporal performance of the two kinds of EEFs. For practical industrial control, the low-dimensional model based on optimal combined EEFs is very useful because it cannot only capture the dominant dynamics of the nonlinear PDEs but simplify the control design as well. The spatiotemporal errors of the predicted output of two models combined with two kinds of EEFs are shown in Figures3 and 4, respectively. As given in Figures 3 and 4, the predicted distributed error of the approximate model based on the optimal combination of EEFs is shown to be significantly smaller than that of the same order approximate model based on initial EEFs.Figure 3
Predicted error based on 3 initial EEFs.Figure 4
Predicted error based on 3 optimal combined EEFs.


## 6. Conclusions

In this note, an optimal combination of EEFs is proposed for model reduction of nonlinear partial differential equations (PDEs), which is obtained by the basis function transformation from the initial EEFs. The transformation matrix is determined by the optimization for the energy error functions of the dynamical systems; the present new EEFs can keep the dynamical information of neglectful modes and generate a lower-dimensional and more precise dynamical system for the PDEs. The numerical example shows its effectiveness and feasibility for model reduction of the nonlinear PDEs.



---

*Source: 347248-2013-07-14.xml*