The interaction between species has been widely studied with reaction–diffusion models. Cross-diffusion systems are quasilinear parabolic equations in which the gradient of one variable induces a flux of another variable. They arise in multi-component systems from physics, chemistry, and biology. The correlation between diffusion and cross-diffusion terms may cause an unstable steady-state, called Turing instability or diffusion-driven instability, leads to formation of patterns [15, 16]. In this paper we consider a well-known cross-diffusion system from population dynamics, the Shigesada-Kawasaki-Teramoto (SKT) with Lotka-Volterra kinetics 
. Reaction-diffusion systems as SKT equation have to be computed for many parameters to predict the patterns. Therefore the numerical simulations are computationally costly. Reduced-order models (ROMs) have emerged as a powerful approach to reduce the cost of evaluating large systems of partial differential equations (PDEs) in multi-query scenario for different parameters. The proper orthogonal decomposition (POD) has been widely used and is a computationally efficient reduced-order modeling technique in large-scale numerical simulations of nonlinear PDEs. Applying POD Galerkin projection, dominant modes of the PDEs are extracted from the snapshots of the full-order solutions and the reduced-order solutions are computed in a linear reduced space. During the offline stage, a set of reduced basis is extracted from a collection of high-fidelity solutions of the full-order model (FOM). In the online stage, the reduced-order solutions are computed in low-order reduced space, spanned by a set of basis functions.
The SKT equation has linear and quadratic nonlinear terms, both in the cross-diffusion and in the Lokta-Volterra parts. Consequently, the semi-discretization of the SKT equation by second order finite difference results in a system of linear-quadratic ordinary differential equations (ODEs). For time discretization, we use second-order linearly implicit Kahan’s method [22, 11]
, which is designed for ODEs with quadratic nonlinear terms, as the SKT equation. In contrast to the fully implicit schemes, such as the Crank-Nicolson scheme, Kahan’s method requires only one step Newton iteration at each time step. When nonlinear PDEs like the SKT equation have polynomial structure, projecting the FOM onto the reduced space yields low-dimensional matrix operators that preserve the polynomial structure of the FOMs, such that the offline and online phases are separated. This enables construction of computationally efficient ROMs without using hyper-reduction techniques like discrete empirical interpolation (DEIM). Online computation of the ROMs are further accelerated by matricizations of tensors [6, 8, 23]. Applying tensorial POD (T-POD) to the SKT equation recovers an efficient offline-online decomposition. Here we make use of the sparse matrix technique MULTIPROD  to speed up the tensor calculations.
For smooth systems where the systems’ energetics can be characterized by using few modes, the global POD (G-POD) method in the whole time interval provides a very efficient way to generate reduced-order systems. However, its applicability to complex, nonlinear PDEs is often limited by the errors associated with the finite truncation in POD modes and unsteadiness of the problem. An alternative (and complementary) approach is the principal interval decomposition (PID) [29, 2, 3], which optimizes the length of time windows over which to perform the POD procedure in such systems. The cross-diffusion systems with pattern formation as the SKT equation have a rapidly changing short transient phase and long stationary phase. This two phases provide as natural decomposition of the whole time domain into two sub-intervals in the principal decomposition framework (P-POD). We show for one- and two dimensional SKT equations, the patterns can be more efficiently and accurately computed with P-POD than with the G-POD. Moreover, cross-diffusion systems have an entropy structure. We show that dissipation of reduced entropy is well preserved for the SKT equation.
The organization of this paper is as follows. In Section 2, we briefly describe the SKT equation. The fully discrete model in space and time is derived in Section 3. The three reduced order methods, G-POD, P-POD and T-POD are described in Section 4. Numerical experiments for one- and two-dimensional SKT equations are presented in Section 5. Finally, we provide brief conclusions and directions for future work.
2 Shigesada-Kawasaki-Teramoto equation
The interaction between species has been widely studied with reaction-diffusion models. The most prominent example is the Lotka-Volterra competition diffusion system which has been extensively investigated in population ecology. When the diffusion of one of the species depends not only on the density of these species but also on the density of the other species, then cross-diffusion occurs, which may rise to formation of patterns. The species with high densities diffuse faster than predicted by the usual linear diffusion towards lower density areas, that leads to the coexistence of two spatial segregated competing species, known as cross-diffusion induced instability. When cross- and self-diffusion are absent, for linear diffusion in a convex domain, the only stable equilibrium solutions are spatially homogenous. In reaction-diffusion with cross-diffusion, the destabilization of a constant steady-state, is followed by the transition to a non-homogeneous steady-state, i.e., formation of patterns. The linear stability analysis shows that the cross-diffusion is the key mechanism for the formation of spatial patterns through Turing instability . In this paper, we consider the strongly coupled reaction-diffusion system with nonlinear self-and cross-diffusion terms. One of the most popular model in population ecology with pattern formation is the SKT cross-diffusion system  with the Lotka-Volterra reaction terms
in a convex bounded domain with a smooth boundary on the time interval with . In (1), and with denote population densities of two competing species and is the Laplace operator. The initial and boundary conditions are
is the unit outward normal vector to the boundary. The homogeneous Neumann (zero-flux) boundary conditions impose the weakest constraint on formation of self-organizing patterns [15, 16]
The parameters and are self-diffusion and linear diffusion coefficients, respectively, while the parameters are the cross-diffusion coefficients. The parameters are assumed to be non-negative. The parameters denote the intrinsic growth rates, the intra-specific competition coefficients, and are inter-specific competition rates. The parameter represents the relative strength of reaction terms.
Pattern formation in the SKT system (1) was investigated using linear and weakly nonlinear stability analysis in [15, 16]. Cross-diffusion destabilizes the uniform equilibrium leading to traveling fonts  in the one-dimensional SKT equation (1) and formation of patterns in the two-dimensional SKT equation (1) . In both papers, it was shown, for parameter values , patterns start to emerge, from an initial condition with a random periodic perturbation of the equilibrium
The entropy structure is crucial to understand various theoretical properties of cross-diffusion systems, such as existence, regularity and long time asymptotic weak solutions of the STK equation (1). The SKT equation (1) can be written alternatively as [19, 20, 13]
with the diffusion matrix
where and are the Lotka-Volterra reaction terms, . A characteristic feature of the cross-diffusion is that the diffusion matrix is generally neither symmetric nor positive definite which complicates the mathematical analysis. However, using a transformation of variables (called entropy variable), the transformed diffusion matrix becomes positive definite and sometimes even symmetric. Hence the existence of global solutions can be established. The entropy for the SKT equation (3) is given without the reaction terms [19, 20, 13] as
when two constants and exist satisfying . The entropy decreases in time .
3 Full order model
The SKT system (1) has been solved by various numerical methods: fully implicit finite volume method , semi-implicit finite difference method , semi-implicit spectral method , explicit Euler and finite difference method . In [27, 26], the SKT system (1) is transformed to a semi-linear PDE through replacing the nonlinear self and cross-diffusion terms by linear reaction-diffusion terms. The resulting semi-linear equations with Lotka-Volterra reaction terms and linear diffusion terms are solved by explicit Euler method with finite difference or finite volume discretization in space.
Here we discretize the SKT equation (1) in space by finite differences, which leads to a system of ODEs of the form
where are semi-discrete approximations to the exact solutions and at spatial grid nodes, and the powers together with the multiplication operator are driven element-wise. The number of spatial grid nodes differ for one- and two-dimensional regions. The components of the semi-discrete solution vectors () in the case of one- and two-dimensional regions are given respectively by
with on one-dimensional regions and on two-dimensional regions, where and are the number of partition in and -directions, respectively.
In the ODE system (5), the matrix represents the finite difference matrix related to the second order centered finite differences approximation to the Laplace operator under homogeneous Neumann boundary conditions. More clearly, let denotes the
-dimensional identity matrix, and let the matrixgiven by
corresponds to the centered finite differences approximation to the Laplace operator under homogeneous Neumann boundary conditions with grid nodes. Then, we have on the one-dimensional regions
whereas on the two-dimensional regions we get that
where and are the mesh sizes in and -directions, respectively, and denotes the Kronecker product.
Collecting linear and quadratic parts, the system (5) can be written as the following linear-quadratic ODE system
where we set
In compact form, we can also write as
where is the state vector, represents the matrix of linear terms, and () includes the quadratic terms given by
In (7), stands for the matricized tensor so that it satisfies the identity for any vector .
The ODE system (6) can be solved in time using explicit or implicit methods. Explicit methods require small time steps and can produce unstable solutions. The implicit integrators require at each time step solution of nonlinear equations iteratively. We solve the semi-discrete linear-quadratic ODE system (6) in time by linearly implicit Kahan’s method [21, 22]
where are the symmetric bilinear forms obtained by the polarization of the quadratic vector fields :
and is the time step size, and is the full discrete solution vector at time . Kahan’s method is time-reversal, symmetric and linearly implicit, i.e., can be computed by solving a single linear system
where is the Jacobian matrix of evaluated at .
Kahan’s method can also be written as a second order convergent Runge-Kutta method of the form 
When fully implicit time integrators are used, at each time step nonlinear equations have to be solved iteratively by Newton’s method or by fixed-point iteration. The linearly implicit Kahan’s methods is much faster than the fully implicit time integrators like the implicit Euler and mid-point rule.
4 Reduced order model
In this section, we introduce three different types of ROMs. The standard way of constructing ROMs is the use of POD with Galerkin projection on the whole time interval, Global-POD (G-POD). The solutions of the SKT equation (1) converge quickly to steady-state after a short transient phase. We have constructed partitioned-POD (P-POD) in two sub-intervals respecting different behavior of the average densities of the FOM solutions in the transient and steady-state phase. The computation of the quadratic terms scale in the reduced order model with . We have constructed an efficient tensorial-POD (T-POD) exploiting the quadratic form of the semi-discrete SKT equation (6).
The POD basis vectors are computed using the method of snapshots. The POD basis for the semi-discrete SKT system (5) can be obtained by stacking all and in one vector and to determine the common subspace
by taking the singular value decomposition (SVD) of that data. But this may produce unstable ROMs such that the resulting ROMs do not preserve the coupling structure of the PDE[6, 28]. In order to maintain the coupling structure in ROMs of the coupled SKT equation, we compute the snapshot matrices and the POD basis vectors separately for the state components and . Consider the discrete state vectors and of (5). The snapshot matrix corresponding to the state is defined as
where each column vector is the full discrete solution vector at the time instance , . Assuming , we then expand the SVD of the snapshot matrix
where the columns of and are the left and right singular vectors of , respectively, and is the diagonal matrix whose diagonal elements are the singular values .
The -POD basis matrix minimizes the least squares error of the snapshot reconstruction
where denotes the Euclidean -norm and denotes the Frobenius norm. The optimal solution of basis matrix to this problem is given by the left singular vectors of corresponding to the largest singular values. The POD state approximation is then , where is the reduced state vector. Throughout the paper, we omit the subscript for easy notation and we denote by simply the -POD basis matrix corresponding to the state . Moreover, we may choose different number of POD modes and related to each state component and , respectively.
Once the POD basis matrices are found, the ROM for the SKT system is obtained as a linear-quadratic ODE as the FOM (6)
where , , , and
The number of POD modes for each component () is determined usually by the relative information content (RIC) defined by
with a tolerance .
4.2 Partitioned POD
POD depends on a global approximation of the snapshot data, which can result in overall deformation of the obtained modes for systems with fast variations in state and the constructed POD cannot capture any dominant structure at all. The PID approach is developed [18, 9, 29, 1, 2, 3] as an alternative approach which preserves the optimality of POD respecting local characteristics of the solutions. In PID, the time domain is divided into non-overlapping intervals, each characterizing a specific stage in the system’s dynamics and evolution. The same POD algorithm is applied within each subinterval to generate a set of basis functions that best fit the respective partition locally. In short, the PID approach can be viewed as decomposing the G-POD subspace into a few of locally optimal subspaces to obtain accurate partitioned ROMs with smaller sizes in each individual sub-interval. The PID was first studied in  and then is applied successfully to nonlinear convective fluid problems: Burgers equation, Boussinesq equation and Navier-Stokes equation [9, 29, 1, 2] and two-dimensional turbulence flow . Adaptive partitioning and clustering techniques can be used to construct the sub-intervals of different size, but also the time domain may be decomposed into equidistant sub-intervals.
The cross-diffusion with pattern formation like SKT equation (1) is characterized by short transient phase and long stationary phase. This leads to a natural decomposition of the whole time domain into two sub-intervals in the PID approach, one for the transient phase and the other for the stationary phase. The transition from transient phase to stationary phase can be determined by using average densities
When the difference of the both average densities at two consecutive time instances are lower than a prespecified tolerance , the stationary phase occur. This transition point is then used as the interface of two sub-intervals in the PID approach. More clearly, let denotes the transition point, i.e., the index is the minimum integer such that
for a given tolerance . Then, we decompose the whole time-interval into two sub-intervals and with the common interface . According to the PID approach, we set different POD basis matrices and through the snapshot matrices composed of the full solution vectors on either intervals and , respectively. Finally, we should enforce the interface constraint that the full discrete solution vectors from the first interval agree with the initial vectors on the second interval . Using the POD basis matrices defined on two sub-intervals, this requires the following identity
where and are reduced solution vectors at the interface time on the intervals and , respectively. Using the orthonormality of the POD basis matrices, the initial reduced solution vector on the interval can be recovered from the reduced solution vector at the interface on the interval as
4.3 Tensorial POD
Although the dimension of the ROM (9) is small compared to the dimension of the FOM (5), the computation of the quadratic nonlinearities still scale with the dimension of FOM. T-POD approach can handle this computational inefficiency utilizing the matricized tensor together with the properties of Kronecker product.
The explicit form of the reduced quadratic terms in the SKT equation (9) are given by
It is clear that all the terms above are of the form
The terms with tensor in (12) are computed using the properties of Kronecker product, which depends on the computation of the reduced tensor so that we get
where the small matrix can be precomputed in the offline stage. Although is computed offline, the explicit computation of may be inefficient since it depends on the full dimension . In order to avoid from this computational burden, is computed in an efficient way using -mode matricizations of tensors .
which utilizes the structure of without constructing explicitly. In [8, 7] the CUR matrix approximation  of is used to increase computational efficiency. Instead, here we make use of the ”MULTIPROD”  to increase computational efficiency. The MULTIPROD 111https://www.mathworks.com/matlabcentral/fileexchange/8773-multiple-matrix-multiplications-with-array-expansion-enabled handles multiple multiplications of the multi-dimensional arrays via virtual array expansion. It is a fast and memory efficient generalization for arrays of the MATLAB matrix multiplication operator. For any given two vectors and , the Kronecker product satisfies
where vec denotes the vectorization of a matrix. Using the identity in (14), the matrix can be constructed as
Reshaping the matrix as and computing MULTIPROD of and in and dimensions, we obtain
5 Numerical results
In this section we present results of the numerical experiments for the one-and two dimensional SKT system (1). We compare the FOM and ROM solutions by G-POD and P-POD, computational gain by TPOD over POD, and show the decreasing structure of the entropy.
Initial conditions are taken as random periodic perturbation around the equilibrium given in (2). We stop the computation of the FOMs when the solutions are sufficiently close to the steady-states, e.g., when the termination condition
is satisfied for a prescribed tolerance , where denotes the usual -norm over the domain , and calculated by the trapezoidal quadrature. We take in all simulations .
The accuracy of the ROM solutions are measured using the time averaged relative -errors defined by
All the simulations are performed on a machine with Intel: CoreTM i7 2.5 GHz 64 bit CPU, 8 GB RAM, Windows 10, using 64 bit MatLab R2014. For the time-dependent problems with many time steps, such as the SKT system, the snapshot matrix is large, leading to an expensive SVD. We use randomized SVD (rSVD) algorithm  which only needs to perform SVD of small matrices, to efficiently generate a reduced basis with large snapshot matrices. In the ROMs, T-POD framework is applied in both the G-POD and P-POD approaches.
5.1 One-dimensional SKT equation
where, is taken larger than the critical value of the bifurcation parameter is , so that pattern formation can occur. Spatial interval is set to with the mesh size , and time step size is taken as . The steady-states are reached at .
In Figure 1, left, the patterns at the steady-state are shown, that are formed starting from an initial datum which is a random periodic perturbation of the equilibrium (2). The FOM solutions are very close to those in . In Figure 1, right, the densities start a plateau around . Accordingly, using the PID tolerance , we obtain the transition point , and the time interval is split into sub-intervals and .
In Figure 2, singular values are plotted in the whole time interval, in the intervals of the transient and of the steady-state phases. The singular values decay at the same rate, slowly in the whole time interval and in the first interval , whereas the decay is faster in the second interval of the steady-state phase.
The number of POD modes and the time averaged relative -errors (15) between FOM and ROM approximations for different RIC tolerances in (10) are listed in Table 1. For the same RIC tolerances, the P-POD requires fewer POD modes in the interval of the steady-states comparing with the ones required in the interval of the transition phase. The time averaged relative -errors obtained by the P-POD are smaller than the errors obtained by the G-POD approach. The ROM solutions in Figure 3 computed using the RIC tolerance are very close to the FOM solutions.
5.2 Two-dimensional SKT equation
Spatial domain is set to . We take in both space directions the same number of partition , and time step size is set to . The steady-states are reached at .
In Figure 4, the densities start almost unchanged around . By the PID tolerance , the time interval is split into sub-intervals and . Decay of the singular vales in Figure 5 is similar to the one-dimensional SKT equation in the previous example. The FOM and ROM solutions in Figure 6 computed with the RIC tolerance agree well, and that the ones by the P-POD approach are almost the same as the FOM solutions.
5.3 Entropy decrease
The entropy (4) of the SKT equation is defined with the Lokta-Volterra kinetics terms , . The entropies are computed with the same diffusion coefficients for one-and two dimensional SKT equation (1) setting the . Initial conditions are taken from  are given by
for one-and two dimensional problems, respectively. Since the SKT equation is solved without the reaction terms, the transient phase is absent . Therefore, the ROMs are computed only by the G-POD approach.
Exploiting the different behavior of transient and steady-state solutions of the SKT equation, reduced solutions are obtained in a computationally efficient way. The quadratic nonlinear terms of SKT equation are reflected in the semi-discrete linear-quadratic ODE system using finite-differences, which enables separation of the offline-online computation. The ROM solutions depend affinely on the parameters in both of the linear and quadratic parts. This allows the prediction of patterns for different parameter values without interpolation. We plan to investigate the bifurcation behavior of ROM using parametrized ROM techniques.
-  M. Ahmed and O. San. Stabilized principal interval decomposition method for model reduction of nonlinear convective systems with moving shocks. Computational and Applied Mathematics, 37(5):6870–6902, 2018.
-  S. E. Ahmed, Sk. M. Rahman, O. San, A. Rasheed, and I. M. Navon. Memory embedded non-intrusive reduced order modeling of non-ergodic flows. Physics of Fluids, 31(12):126602, 2019.
-  Shady E. Ahmed and Omer San. Breaking the Kolmogorov barrier in model reduction of fluid flows. Fluids, 5(1):26, 2020.
-  B. Andreianov, M. Bendahmane, and R. Ruiz-Baier. Analysis of a finite volume method for a cross-diffusion model in population dynamics. Mathematical Models and Methods in Applied Sciences, 21(02):307–344, 2011.
-  P. Benner and T. Breiten. Two-sided projection methods for nonlinear model order reduction. SIAM Journal on Scientific Computing, 37(2):B239–B260, 2015.
-  P. Benner and L. Feng. Model order reduction for coupled problems (survey). Appl. Comput. Math., 14(1):3–22, 2015.
-  P. Benner and P. Goyal. Interpolation-based model order reduction for polynomial parametric systems. arXiv, 2019.
-  P. Benner, P. Goyal, and S. Gugercin. -quasi-optimal model order reduction for quadratic-bilinear control systems. SIAM Journal on Matrix Analysis and Applications, 39(2):983–1032, 2018.
-  J. Borggaard, A. Hay, and D. Pelletier. Interval-based reduced-order models for unsteady fluid flow. International Journal of Numerical Analysis and Modeling, 4:353–367, 2007.
-  E. Celledoni, R. I. McLachlan, D. I. McLaren, B. Owren, and G. R. W. Quispel. Discretization of polynomial vector fields by polarization. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 471(2184), 2015.
-  E. Celledoni, R. I McLachlan, B. Owren, and G. R. W. Quispel. Geometric properties of Kahan’s method. Journal of Physics A: Mathematical and Theoretical, 46(2):025201, 2013.
Saifon Chaturantabut and Danny C. Sorensen.
A state space error estimate for POD-DEIM nonlinear model reduction.SIAM Journal on Numerical Analysis, 50(1):46–63, 2012.
-  X. Chen and A. Jüngel. When do cross-diffusion systems have an entropy structure? arXiv, 2019.
-  G. Galiano, M. L. Garzón, and A. Jüngel. Analysis and numerical solution of a nonlinear cross-diffusion system arising in population dynamics. RACSAM, 95(2):281–295, 2001.
-  G. Gambino, M.C. Lombardo, and M. Sammartino. Turing instability and traveling fronts for a nonlinear reaction-diffusion system with cross-diffusion. Mathematics and Computers in Simulation, 82(6):1112 – 1132, 2012.
-  G. Gambino, M.C. Lombardo, and M. Sammartino. Pattern formation driven by cross-diffusion in a 2D domain. Nonlinear Analysis: Real World Applications, 14(3):1755 – 1779, 2013.
-  N. Halko, P. G. Martinsson, and J. A. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review, 53(2):217–288, 2011.
-  W.L. IJzerman and E. van Groesen. Low-dimensional model for vortex merging in the two-dimensional temporal mixing layer. European Journal of Mechanics - B/Fluids, 20(6):821 – 840, 2001.
-  A. Jüngel. The boundedness-by-entropy method for cross-diffusion systems. Nonlinearity, 28(6):1963–2001, 2015.
-  A. Jüngel. Cross-diffusion systems with entropy structure. In Conference on Differential Equations and Their Applications, pages 181–190. Bratislava, 2017.
-  W. Kahan. Unconventional numerical methods for trajectory calculations. Technical report, Computer Science Division and Department of Mathematics, University of California, Berkeley, 1993. Unpublished lecture notes.
-  William Kahan and Ren-Chang Li. Unconventional schemes for a class of ordinary differential equations-with applications to the Korteweg-de Vries equation. Journal of Computational Physics, 134(2):316 – 331, 1997.
-  Boris Kramer and Karen E. Willcox. Nonlinear model order reduction via lifting transformations and proper orthogonal decomposition. AIAA Journal, 57(6):2297–2307, 2019.
-  P. D. Leva. MULTIPROD TOOLBOX, multiple matrix multiplications, with array expansion enabled. Technical report, University of Rome Foro Italico, Rome, 2008.
-  M. W. Mahoney and P. Drineas. CUR matrix decompositions for improved data analysis. Proceedings of the National Academy of Sciences, 106(3):697–702, 2009.
-  H. Murakawa. A linear finite volume method for nonlinear cross-diffusion systems. Numerische Mathematik, 136(1):1–26, 2017.
-  Murakawa, H. A linear scheme to approximate nonlinear cross-diffusion systems. ESAIM: M2AN, 45(6):1141–1161, 2011.
-  T. Reis and T. Stykel. Stability analysis and model order reduction of coupled systems. Math. Comput. Model. Dyn. Syst., 13(5):413–436, 2007.
-  O. San and J. Borggaard. Principal interval decomposition framework for POD reduced-order modeling of convective Boussinesq flows. International Journal for Numerical Methods in Fluids, 78(1):37–62, 2015.
-  N. Shigesada, K. Kawasaki, and E. Teramoto. Spatial segregation of interacting species. Journal of Theoretical Biology, 79(1):83–99, 1979.
-  Zheng Sun, José A. Carrillo, and Chi-Wang Shu. An entropy stable high-order discontinuous galerkin method for cross-diffusion gradient flow systems. Kinetic & Related Models, 12(1937-5093_2019_4_885):885, 2019.
-  A. M. Turing. The chemical basis of morphogenesis. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 237(641):37–72, 1952.