1 Introduction
The interaction between species has been widely studied with reaction–diffusion models. Crossdiffusion systems are quasilinear parabolic equations in which the gradient of one variable induces a flux of another variable. They arise in multicomponent systems from physics, chemistry, and biology. The correlation between diffusion and crossdiffusion terms may cause an unstable steadystate, called Turing instability or diffusiondriven instability, leads to formation of patterns [15, 16]. In this paper we consider a wellknown crossdiffusion system from population dynamics, the ShigesadaKawasakiTeramoto (SKT) with LotkaVolterra kinetics [30]
. Reactiondiffusion systems as SKT equation have to be computed for many parameters to predict the patterns. Therefore the numerical simulations are computationally costly. Reducedorder models (ROMs) have emerged as a powerful approach to reduce the cost of evaluating large systems of partial differential equations (PDEs) in multiquery scenario for different parameters. The proper orthogonal decomposition (POD) has been widely used and is a computationally efficient reducedorder modeling technique in largescale numerical simulations of nonlinear PDEs. Applying POD Galerkin projection, dominant modes of the PDEs are extracted from the snapshots of the fullorder solutions and the reducedorder solutions are computed in a linear reduced space. During the offline stage, a set of reduced basis is extracted from a collection of highfidelity solutions of the fullorder model (FOM). In the online stage, the reducedorder solutions are computed in loworder reduced space, spanned by a set of basis functions.
The SKT equation has linear and quadratic nonlinear terms, both in the crossdiffusion and in the LoktaVolterra parts. Consequently, the semidiscretization of the SKT equation by second order finite difference results in a system of linearquadratic ordinary differential equations (ODEs). For time discretization, we use secondorder 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 CrankNicolson 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 lowdimensional 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 hyperreduction techniques like discrete empirical interpolation (DEIM)
[12]. Online computation of the ROMs are further accelerated by matricizations of tensors [6, 8, 23]. Applying tensorial POD (TPOD) to the SKT equation recovers an efficient offlineonline decomposition. Here we make use of the sparse matrix technique MULTIPROD [24] to speed up the tensor calculations.For smooth systems where the systems’ energetics can be characterized by using few modes, the global POD (GPOD) method in the whole time interval provides a very efficient way to generate reducedorder 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 crossdiffusion 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 subintervals in the principal decomposition framework (PPOD). We show for one and two dimensional SKT equations, the patterns can be more efficiently and accurately computed with PPOD than with the GPOD. Moreover, crossdiffusion 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, GPOD, PPOD and TPOD are described in Section 4. Numerical experiments for one and twodimensional SKT equations are presented in Section 5. Finally, we provide brief conclusions and directions for future work.
2 ShigesadaKawasakiTeramoto equation
The interaction between species has been widely studied with reactiondiffusion models. The most prominent example is the LotkaVolterra 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 crossdiffusion 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 crossdiffusion induced instability. When cross and selfdiffusion are absent, for linear diffusion in a convex domain, the only stable equilibrium solutions are spatially homogenous. In reactiondiffusion with crossdiffusion, the destabilization of a constant steadystate, is followed by the transition to a nonhomogeneous steadystate, i.e., formation of patterns. The linear stability analysis shows that the crossdiffusion is the key mechanism for the formation of spatial patterns through Turing instability [32]. In this paper, we consider the strongly coupled reactiondiffusion system with nonlinear selfand crossdiffusion terms. One of the most popular model in population ecology with pattern formation is the SKT crossdiffusion system [30] with the LotkaVolterra reaction terms
(1)  
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
where
is the unit outward normal vector to the boundary
. The homogeneous Neumann (zeroflux) boundary conditions impose the weakest constraint on formation of selforganizing patterns [15, 16]The parameters and are selfdiffusion and linear diffusion coefficients, respectively, while the parameters are the crossdiffusion coefficients. The parameters are assumed to be nonnegative. The parameters denote the intrinsic growth rates, the intraspecific competition coefficients, and are interspecific 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]. Crossdiffusion destabilizes the uniform equilibrium leading to traveling fonts [15] in the onedimensional SKT equation (1) and formation of patterns in the twodimensional SKT equation (1) [16]. 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
(2) 
The critical value of bifurcation parameter is calculated using the Turing instability analysis for one and twodimensional SKT equations (1) in [15, 16].
The entropy structure is crucial to understand various theoretical properties of crossdiffusion 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]
(3) 
with the diffusion matrix
where and are the LotkaVolterra reaction terms, . A characteristic feature of the crossdiffusion 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
(4) 
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 [4], semiimplicit finite difference method [14], semiimplicit spectral method [16], explicit Euler and finite difference method [15]. In [27, 26], the SKT system (1) is transformed to a semilinear PDE through replacing the nonlinear self and crossdiffusion terms by linear reactiondiffusion terms. The resulting semilinear equations with LotkaVolterra 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
(5)  
where are semidiscrete approximations to the exact solutions and at spatial grid nodes, and the powers together with the multiplication operator are driven elementwise. The number of spatial grid nodes differ for one and twodimensional regions. The components of the semidiscrete solution vectors () in the case of one and twodimensional regions are given respectively by
with on onedimensional regions and on twodimensional 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 matrix
given bycorresponds to the centered finite differences approximation to the Laplace operator under homogeneous Neumann boundary conditions with grid nodes. Then, we have on the onedimensional regions
whereas on the twodimensional 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 linearquadratic ODE system
where we set
In compact form, we can also write as
(6) 
where is the state vector, represents the matrix of linear terms, and () includes the quadratic terms given by
(7)  
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 semidiscrete linearquadratic ODE system (6) in time by linearly implicit Kahan’s method [21, 22]
(8) 
where are the symmetric bilinear forms obtained by the polarization of the quadratic vector fields [10]:
and is the time step size, and is the full discrete solution vector at time . Kahan’s method is timereversal, 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 RungeKutta method of the form [11]
When fully implicit time integrators are used, at each time step nonlinear equations have to be solved iteratively by Newton’s method or by fixedpoint iteration. The linearly implicit Kahan’s methods is much faster than the fully implicit time integrators like the implicit Euler and midpoint 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, GlobalPOD (GPOD). The solutions of the SKT equation (1) converge quickly to steadystate after a short transient phase. We have constructed partitionedPOD (PPOD) in two subintervals respecting different behavior of the average densities of the FOM solutions in the transient and steadystate phase. The computation of the quadratic terms scale in the reduced order model with . We have constructed an efficient tensorialPOD (TPOD) exploiting the quadratic form of the semidiscrete SKT equation (6).
4.1 GlobalPOD
The POD basis vectors are computed using the method of snapshots. The POD basis for the semidiscrete 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 aswhere 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 linearquadratic ODE as the FOM (6)
(9) 
where , , , and
The number of POD modes for each component () is determined usually by the relative information content (RIC) defined by
(10) 
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 nonoverlapping 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 GPOD subspace into a few of locally optimal subspaces to obtain accurate partitioned ROMs with smaller sizes in each individual subinterval. The PID was first studied in [18] and then is applied successfully to nonlinear convective fluid problems: Burgers equation, Boussinesq equation and NavierStokes equation [9, 29, 1, 2] and twodimensional turbulence flow [3]. Adaptive partitioning and clustering techniques can be used to construct the subintervals of different size, but also the time domain may be decomposed into equidistant subintervals.
The crossdiffusion 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 subintervals 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
(11) 
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 subintervals 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 timeinterval into two subintervals 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 subintervals, 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. TPOD 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
(12) 
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 [5].
Recently tensorial algorithms are developed by exploiting the particular structure of Kronecker product [8, 7]. The reduced matrix can be given in MATLAB notation as follows
(13) 
which utilizes the structure of without constructing explicitly. In [8, 7] the CUR matrix approximation [25] of is used to increase computational efficiency. Instead, here we make use of the ”MULTIPROD” [24] to increase computational efficiency. The MULTIPROD ^{1}^{1}1https://www.mathworks.com/matlabcentral/fileexchange/8773multiplematrixmultiplicationswitharrayexpansionenabled handles multiple multiplications of the multidimensional 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
(14) 
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 oneand two dimensional SKT system (1). We compare the FOM and ROM solutions by GPOD and PPOD, 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 steadystates, 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
(15) 
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 timedependent 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 [17] which only needs to perform SVD of small matrices, to efficiently generate a reduced basis with large snapshot matrices. In the ROMs, TPOD framework is applied in both the GPOD and PPOD approaches.
5.1 Onedimensional SKT equation
Our first example is the onedimensional SKT equation (1) with the parameters are taken as in [15]
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 steadystates are reached at .
In Figure 1, left, the patterns at the steadystate 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 [15]. 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 subintervals and .
In Figure 2, singular values are plotted in the whole time interval, in the intervals of the transient and of the steadystate 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 steadystate 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 PPOD requires fewer POD modes in the interval of the steadystates comparing with the ones required in the interval of the transition phase. The time averaged relative errors obtained by the PPOD are smaller than the errors obtained by the GPOD approach. The ROM solutions in Figure 3 computed using the RIC tolerance are very close to the FOM solutions.
GPOD  PPOD  

#Modes  Error  #Modes  #Modes  Error  
5(5)  5.35e03(7.20e03)  6(5)  2(2)  1.84e03(2.35e03)  
9(7)  9.23e04(1.15e03)  9(8)  2(2)  7.74e04(9.16e04)  
12(10)  4.14e04(5.01e04)  12(10)  2(2)  2.15e04(2.62e04)  
15(14)  1.34e04(1.64e04)  15(13)  3(4)  9.27e05(1.16e04) 
5.2 Twodimensional SKT equation
Our second example is the twodimensional SKT equation (1) with the parameters are taken as in [16]
Spatial domain is set to . We take in both space directions the same number of partition , and time step size is set to . The steadystates are reached at .
In Figure 4, the densities start almost unchanged around . By the PID tolerance , the time interval is split into subintervals and . Decay of the singular vales in Figure 5 is similar to the onedimensional 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 PPOD approach are almost the same as the FOM solutions.
The errors by GPOD ad PPOD in Table 2 show the same behavior as the errors in Table 1 for the onedimensional case.
GPOD  PPOD  

#Modes  Error  #Modes  #Modes  Error  
4(3)  1.97e03(2.18e03)  5(4)  1(1)  1.26e03(1.40e03)  
6(5)  9.17e04(1.02e03)  7(6)  2(2)  3.53e04(3.89e04)  
9(8)  2.14e04(2.30e04)  10(8)  2(2)  2.20e04(2.46e04)  
12(11)  6.73e05(7.29e05)  11(10)  3(3)  5.76e05(6.42e05) 
5.3 Entropy decrease
The entropy (4) of the SKT equation is defined with the LoktaVolterra kinetics terms , . The entropies are computed with the same diffusion coefficients for oneand two dimensional SKT equation (1) setting the . Initial conditions are taken from [31] are given by
for oneand 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 GPOD approach.
6 Conclusions
Exploiting the different behavior of transient and steadystate 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 semidiscrete linearquadratic ODE system using finitedifferences, which enables separation of the offlineonline 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.
References
 [1] 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.
 [2] S. E. Ahmed, Sk. M. Rahman, O. San, A. Rasheed, and I. M. Navon. Memory embedded nonintrusive reduced order modeling of nonergodic flows. Physics of Fluids, 31(12):126602, 2019.
 [3] Shady E. Ahmed and Omer San. Breaking the Kolmogorov barrier in model reduction of fluid flows. Fluids, 5(1):26, 2020.
 [4] B. Andreianov, M. Bendahmane, and R. RuizBaier. Analysis of a finite volume method for a crossdiffusion model in population dynamics. Mathematical Models and Methods in Applied Sciences, 21(02):307–344, 2011.
 [5] P. Benner and T. Breiten. Twosided projection methods for nonlinear model order reduction. SIAM Journal on Scientific Computing, 37(2):B239–B260, 2015.
 [6] P. Benner and L. Feng. Model order reduction for coupled problems (survey). Appl. Comput. Math., 14(1):3–22, 2015.
 [7] P. Benner and P. Goyal. Interpolationbased model order reduction for polynomial parametric systems. arXiv, 2019.
 [8] P. Benner, P. Goyal, and S. Gugercin. quasioptimal model order reduction for quadraticbilinear control systems. SIAM Journal on Matrix Analysis and Applications, 39(2):983–1032, 2018.
 [9] J. Borggaard, A. Hay, and D. Pelletier. Intervalbased reducedorder models for unsteady fluid flow. International Journal of Numerical Analysis and Modeling, 4:353–367, 2007.
 [10] 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.
 [11] 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.

[12]
Saifon Chaturantabut and Danny C. Sorensen.
A state space error estimate for PODDEIM nonlinear model reduction.
SIAM Journal on Numerical Analysis, 50(1):46–63, 2012.  [13] X. Chen and A. Jüngel. When do crossdiffusion systems have an entropy structure? arXiv, 2019.
 [14] G. Galiano, M. L. Garzón, and A. Jüngel. Analysis and numerical solution of a nonlinear crossdiffusion system arising in population dynamics. RACSAM, 95(2):281–295, 2001.
 [15] G. Gambino, M.C. Lombardo, and M. Sammartino. Turing instability and traveling fronts for a nonlinear reactiondiffusion system with crossdiffusion. Mathematics and Computers in Simulation, 82(6):1112 – 1132, 2012.
 [16] G. Gambino, M.C. Lombardo, and M. Sammartino. Pattern formation driven by crossdiffusion in a 2D domain. Nonlinear Analysis: Real World Applications, 14(3):1755 – 1779, 2013.
 [17] 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.
 [18] W.L. IJzerman and E. van Groesen. Lowdimensional model for vortex merging in the twodimensional temporal mixing layer. European Journal of Mechanics  B/Fluids, 20(6):821 – 840, 2001.
 [19] A. Jüngel. The boundednessbyentropy method for crossdiffusion systems. Nonlinearity, 28(6):1963–2001, 2015.
 [20] A. Jüngel. Crossdiffusion systems with entropy structure. In Conference on Differential Equations and Their Applications, pages 181–190. Bratislava, 2017.
 [21] 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.
 [22] William Kahan and RenChang Li. Unconventional schemes for a class of ordinary differential equationswith applications to the Kortewegde Vries equation. Journal of Computational Physics, 134(2):316 – 331, 1997.
 [23] Boris Kramer and Karen E. Willcox. Nonlinear model order reduction via lifting transformations and proper orthogonal decomposition. AIAA Journal, 57(6):2297–2307, 2019.
 [24] P. D. Leva. MULTIPROD TOOLBOX, multiple matrix multiplications, with array expansion enabled. Technical report, University of Rome Foro Italico, Rome, 2008.
 [25] 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.
 [26] H. Murakawa. A linear finite volume method for nonlinear crossdiffusion systems. Numerische Mathematik, 136(1):1–26, 2017.
 [27] Murakawa, H. A linear scheme to approximate nonlinear crossdiffusion systems. ESAIM: M2AN, 45(6):1141–1161, 2011.
 [28] T. Reis and T. Stykel. Stability analysis and model order reduction of coupled systems. Math. Comput. Model. Dyn. Syst., 13(5):413–436, 2007.
 [29] O. San and J. Borggaard. Principal interval decomposition framework for POD reducedorder modeling of convective Boussinesq flows. International Journal for Numerical Methods in Fluids, 78(1):37–62, 2015.
 [30] N. Shigesada, K. Kawasaki, and E. Teramoto. Spatial segregation of interacting species. Journal of Theoretical Biology, 79(1):83–99, 1979.
 [31] Zheng Sun, José A. Carrillo, and ChiWang Shu. An entropy stable highorder discontinuous galerkin method for crossdiffusion gradient flow systems. Kinetic & Related Models, 12(19375093_2019_4_885):885, 2019.
 [32] 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.
Comments
There are no comments yet.