Electronic Journal of Differential Equations, Vol. 2021 (2021), No. 55, pp. 1–17. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu TURING INSTABILITY ANALYSIS OF A SINGULAR CROSS-DIFFUSION PROBLEM GONZALO GALIANO, VÍCTOR GONZÁLEZ-TABERNERO Abstract. The population model by Busenberg and Travis is a paradigmatic model in ecology and tumor modeling because its ability to capture interest- ing phenomena such as segregation of populations. Its singular mathematical structure enforces the consideration of regularized problems to deduce prop- erties as fundamental as the existence of solutions. In this article we perform a weakly nonlinear stability analysis of a general class of regularized problems to study the convergence of the instability modes in the limit of the regu- larization parameter. We demonstrate with some specific examples that the pattern formation observed in the regularized problems, with unbounded wave numbers, is not present in the limit problem because of the amplitude decay of the oscillations. We also check the results of the stability analysis with direct finite element simulations of the problem. 1. Introduction Busenberg and Travis [4] introduced a class of singular cross-diffusion problems under the assumption that the spatial relocation of each species is caused by a diffusion flow which depends on the densities of all the involved species. In the case of two species, if u1 and u2 denote their densities, the flow, in its simplest form, may be assumed to be determined by the total population u1 + u2, and thus the conservation laws for both species lead to the system ∂tu1 − div ( u1(∇u1 +∇u2) ) = f1(u1, u2), (1.1) ∂tu2 − div ( u2(∇u1 +∇u2) ) = f2(u1, u2). (1.2) Here, the functions f1 and f2 capture some ecological features of the populations, such as growth, competition, etc., and the equations are complemented with non- negative initial data and non-flow boundary conditions. System (1.1)-(1.2) is called a cross-diffusion system because the flow of each species depend upon the densities of the other species. We call it singular because the resulting diffusion matrix is singular. Indeed, when rewriting (1.1)-(1.2) in matrix form, for u = (u1, u2), we obtain the equation ∂tu− div(A(u)∇ bu) = (u), with A(u) = ( u1 u1 u2 u2 ) , 2010 Mathematics Subject Classification. 35K55, 35B36, 92D25. Key words and phrases. Cross-diffusion; Turing instability; weakly nonlinear equation; finite element method. c©2021. This work is licensed under a CC BY 4.0 license. Submitted December 8, 2020. Published June 21, 2021. 1 2 G. GALIANO, V. GONZÁLEZ-TABERNERO EJDE-2021/55 where the divergence is applied by rows. The full and singular structure of A introduces serious difficulties in the mathematical analysis of the problem, as we shall comment later. In his seminal paper [14], Turing introduced a mechanism explaining how spa- tially uniform equilibriums may evolve, small perturbations mediating, into stable equilibriums with non-trivial spatial structure. He considered a system of the type ∂tu1 −∆u1 = f1(u1, u2), (1.3) ∂tu2 − σ∆u2 = f2(u1, u2), (1.4) with σ > 0, and proved that when σ is small or large enough then the stable equilibrium of the dynamical system ∂tv1 = f1(v1, v2), (1.5) ∂tv2 = f2(v1, v2) (1.6) is not stable for the diffusion system (1.3)-(1.4), being replaced by a non-uniform equilibrium with spatial structure or patterns. This mechanism is known as Turing instability or Turing bifurcation. In this article we study Turing instability for the cross-diffusion singular sys- tem (1.1)-(1.2). We already know that some cross-diffusion systems, such as the paradigmatic SKT model introduced by Shigesada, Kawasaki and Teramoto [13], exhibit Turing instability when cross-diffusion coefficients are large in comparison with self-diffusion coefficients, see e.g. [10, 11]. However, the singularity of the dif- fusion matrix of the system (1.1)-(1.2), not present in the SKT model, introduces important mathematical difficulties to the analysis of this system. Regarding the existence of solutions of (1.1)-(1.2), it has been proved only in some special situations: for a bounded spatial domain Ω ⊂ R (Bertsch et al. [2]) and for Ω = Rn (Bertsch et al. [3]). In their proofs, the following observation is crucial: adding the two equations of (1.1)-(1.2) we deduce that if a solution of this system does exist then the total population, u = u1 +u2, satisfies the porous medium type equation ∂tu− div(u∇u) = f(u), (1.7) for which the theory of existence and uniqueness of solutions is well established. In particular, if the initial data of the total population is bounded away from zero and if f is regular enough with f(0) ≥ 0, it is known that the solution of (1.7) remains positive and smooth for all time. This allows to introduce the change of unknowns w1 = u1/u into the original problem (1.1)-(1.2) to deduce the equivalent formulation ∂tu− div(u∇u) = F1(u,w1), (1.8) ∂tw1 −∇u · ∇w1 = F2(u,w1), (1.9) for certain well-behaved functions F1 and F2. This change reveals the parabolic- hyperbolic structure of the problem, which is handled by Bertsch et al. [2, 3] by the vanishing viscosity method, adding the regularization term −δ∆w1 to the left hand side of (1.9), and justifying the passing to the limit δ → 0 in the corresponding sequence of solutions. In [9] we followed a different approach to prove the existence of solutions of the original system (1.1)-(1.2) for a bounded domain Ω ⊂ R. We introduced a parabolic regularization of the original system by adding a cross-diffusion perturbation term EJDE-2021/55 TURING INSTABILITY ANALYSIS OF A CROSS-DIFFUSION PROBLEM 3 while keeping the porous medium type equation satisfied by u. More concretely, we considered the system ∂tu1 − div ( u1(∇u1 +∇u2) ) − δ 2 ∆(u1(u1 + u2)) = f1(u1, u2), (1.10) ∂tu2 − div ( u2(∇u1 +∇u2) ) − δ 2 ∆(u2(u1 + u2)) = f2(u1, u2), (1.11) and then used previous results for cross-diffusion systems [7, 6, 12] to establish the existence of solutions of the approximated problems. Then, estimates of the bounded variation (BV) similar to those in [2] allowed us to prove the convergence of the sequence (u (δ) 1 , u (δ) 2 ) to a solution of the original problem. Let us mention that the system (1.1)-(1.2) is a limit case of a general type of problems with diffusion matrix given by A(u) = ( a11u1 a12u1 a21u2 a22u2 ) , for which, if aii > 0, for i = 1, 2, and a11a22 > a12a21, then the existence of solutions is ensured for any spatial dimension of Ω, see [8]. In addition, it has been recently shown that this kind of systems, when set in the whole space Ω = Rn, may be rigorously deduced as mean field limits [5]. Concerning Turing instability, since the diffusion matrix corresponding to the system (1.1)-(1.2), A(u), is singular, the linearization of this system at an equilib- rium of the dynamical system (1.5)-(1.6) does not provide any information on the behavior of the equilibrium in the spatial dependent case. Thus, our approach to Turing instability for the system (1.1)-(1.2) relies on the study of this property for regularized problems like (1.10)-(1.11) and their limit behavior. For the regularized problems, we prove that linear instability is always present for δ small enough. Interestingly, the linear analysis also establishes that the main instability wave number is unbounded as δ → 0. For understanding this convergence of a increasingly oscillating sequence of func- tions to a BV function (the solution of (1.1)-(1.2) ensured in [2, 9]), we perform a weakly nonlinear analysis (WNA) which allows to gain insight into the behavior of the amplitude of the main instability mode as δ → 0. As expected, we find that the amplitude of the main instability mode vanishes in the limit δ → 0 and is, therefore, coherent with the BV convergence. In addition, this result also suggests that the uniform equilibrium is stable for the original problem. We, finally, illustrate our analytical results with some numerical examples. 2. Main results For simplicity, we study Turing instability for the one-dimensional spatial setting which has also the advantage of a well established existence theory for the case of a bounded domain [2, 9]. Redefining the functions f1 and f2, we can fix without loss of generality Ω = (0, π) and then rewrite problem (1.1)-(1.2) together with the auxiliary conditions as ∂tu1 − ∂x ( u1(∂xu1 + ∂xu2) ) = f1(u1, u2) in QT , (2.1) ∂tu2 − ∂x ( u2(∂xu1 + ∂xu2) ) = f2(u1, u2) in QT , (2.2) u1(∂xu1 + ∂xu2) = u2(∂xu1 + ∂xu2) = 0 on ΓT , (2.3) u1(0, ·) = u10, u2(0, ·) = u20 in Ω, (2.4) 4 G. GALIANO, V. GONZÁLEZ-TABERNERO EJDE-2021/55 where QT = (0, T )×Ω and the initial data u10, u20 are non-negative functions. We assume a competitive Lotka-Volterra form for the reaction term, this is, fi(u1, u2) = ui(αi − βi1u1 − βi2u2), for i = 1, 2, and for some non-negative parameters αi, βij , for i, j = 1, 2. To deal with several types of regularized problems we introduce, for positive δ and b, the uniformly parabolic cross-diffusion system ∂tu1 − ∂x ( dδ11(u)∂xu1 + dδ12(u)∂xu2 ) = f b1(u) in QT , (2.5) ∂tu2 − ∂x ( dδ21(u)∂xu1 + dδ22(u)∂xu2 ) = f b2(u) in QT , (2.6) dδ11(u)∂xu1 + dδ12(u)∂xu2 = dδ21(u)∂xu1 + dδ22(u)∂xu2 = 0 on ΓT , (2.7) u1(0, ·) = u10, u2(0, ·) = u20 in Ω, (2.8) where the diffusion matrix Dδ(u) = (dδij(u)) and the Lotka-Volterra function b(u) = (fb1 (u), fb2 (u)) satisfy the following assumptions. Assumption (H1): (1) Dδ(u) is linear in u and affine in δ, so that it allows the decomposition Dδ(u) = D0(u) + δD1(u) = Dδ1u1 +Dδ2u2, (2.9) for some matrices Dδi for i = 1, 2, being the coefficients of Dδ(u) given by dδij(u) = d10 ij u1 + d11 ij u1δ + d20 ij u2 + d21 ij u2δ, for some non-negative constants dmnij , for i, j,m = 1, 2 and n = 0, 1. (2) We assume that dδii(u) > 0 for i = 1, 2, and that det(Dδ(u)) is an increasing function with respect to δ satisfying det(Dδ(u)) > 0 if δ > 0 and u ∈ R2 +. (3) For i, j = 1, 2, f bi (u1, u2) = ui(α b i − βbi1u1 − βbi2u2) for some non-negative αbi , β b ij such that αbi → αi and βbij → βij as b → 0. Moreover, using the notation α0 i = αi and β0 ij = βij , we assume, for b ≥ 0, βb22α b 1 − βb12α b 2 > 0, βb11α b 2 − βb21α b 1 > 0, det(Bb) > 0, tr(Bb) ≥ 0, where Bb = ( βbij ) . (2.10) Observe that (2.10) ensures the existence of a stable coexistence equilibrium for the dynamical system (1.5)-(1.6), given by u∗ = (u∗1, u ∗ 2) = ( βb22α b 1 − βb12α b 2 βb11β b 22 − βb12β b 21 , βb11α b 2 − βb21α b 1 βb11β b 22 − βb12β b 21 ) . (2.11) There are two examples of Dδ(u) in which we are specially interested. The first is chosen because of its simplicity for the calculations, Dδ(u) = ( (1 + δ)u1 u1 u2 (1 + δ)u2 ) , (2.12) with det(Dδ(u)) = δ(2 + δ)u1u2. According to [8], the second hypothesis of (H1) guarantees the well-posedness of the problem (2.5)-(2.8) corresponding to this dif- fusion matrix. The second example corresponds to the approximation used in [9] for proving the existence of BV solutions of the original problem (2.1)-(2.4): Dδ(u) = ( (1 + δ)u1 + δ 2u2 (1 + δ 2 )u1 (1 + δ 2 )u2 δ 2u1 + (1 + δ)u2 ) , (2.13) for which det(Dδ(u)) = 1 2δ(1 + δ)(u1 + u2)2. EJDE-2021/55 TURING INSTABILITY ANALYSIS OF A CROSS-DIFFUSION PROBLEM 5 Our first result provides conditions under which linear instability arises. The following notation is used, K = Db(u∗) = ( −βb11u ∗ 1 −βb12u ∗ 1 −βb21u ∗ 2 −βb22u ∗ 2 ) . (2.14) Theorem 2.1 (Linear instability). Assume (H1), with b ≥ 0. Let u∗ be the coex- istence equilibrium defined by (2.11). If tr(K−1Dδ(u∗)) > 0 ∀δ ≥ 0 (2.15) then there exists δc > 0 such that if δ < δc then u∗ is a linearly unstable equilibrium for problem (2.5)-(2.8). In such situation, the wave number of the main instability mode tends to infinity as δ → 0. Condition (2.15) is equivalent to dδ11(u∗)βb22u ∗ 2 + dδ22(u∗)βb11u ∗ 1 < dδ12(u∗)βb21u ∗ 2 + dδ21(u∗)βb12u ∗ 1 (2.16) and introduces a further restriction on the competence coefficients. Roughly speak- ing, for Bb to fulfill (2.10) and (2.16), its elements must be such that intra- population joint competence is larger than inter-population joint competence (con- dition (2.10)) and one of the inter-population competence coefficients is large in comparison with the others (condition (2.16)). A numeric example we shall work with along the article is Bb = ( 1 b 2 2 1 ) , with b ∈ (0, 1 2 ). (2.17) Assuming the forms of Dδ(u∗) given in Experiments 1 and 2, see (2.12) and (2.13), we have that the conditions (2.10) and (2.16) on Bb are satisfied if δ < b/4 (Ex- periment 1) or δ < bu∗1u ∗ 2/(u ∗ 1 +u∗2)2 (Experiment 2), which are always true for the case of interest: δ close to zero. Our second result allows to estimate not only the instability wave numbers pro- vided by the linear analysis but also the amplitude corresponding to these modes. The approximation of the steady state solution is obtained by using a weakly non- linear analysis (WNA) based on the method of multiple scales. Theorem 2.2. Assume the hypothesis of Theorem 2.1 and let ε2 = (δc − δ)/δc be a small number. Then, there exist sets of data problem such that the stationary WNA approximation of the solution u of problem (2.5)-(2.8) is given by v(x) = u∗ + ερ √ A∞ cos(kcx) + ε2A∞ ( v20 + v22 cos(2kcx) ) +O(ε3), (2.18) where kc ∈ Z is the critical wave number corresponding to δc, A∞ is a positive constant and ρ,v20 and v22 are constant vectors. Our third result focuses on the limit behavior of the critical parameters and the amplitude when δ → 0, this is, when the solutions of the regularized problems converge to the solution of the original singular problem. For the sake of simplicity, we limit our study to the following example: ∂tu1 − ∂x ( u1(∂xu1 + ∂xu2) ) = u1(1− u1), (2.19) ∂tu2 − ∂x ( u2(∂xu1 + ∂xu2) ) = u2(4− (2u1 + u2)), (2.20) 6 G. GALIANO, V. GONZÁLEZ-TABERNERO EJDE-2021/55 whose solutions we approximate by the two-parameter family of solutions of ∂tu1 − ∂x ( u1((1 + δ)∂xu1 + ∂xu2) ) = u1(1− (u1 + b 2 u2)), (2.21) ∂tu2 − ∂x ( u2(∂xu1 + (1 + δ)∂xu2) ) = u2(4− (2u1 + u2)). (2.22) On one hand, Theorem 2.1 ensures the existence of δc > 0 such that, for any b ≥ 0, the equilibrium u∗ = 1 1−b (1− 2b, 2) of (2.21)-(2.22) becomes unstable for δ < δc, with an associated critical wave number such that kc →∞ as δ → 0. On the other hand, for δ < b/4 and b → 0, the sequence of solutions of (2.21)- (2.22) converges to a solution of (2.19)-(2.20) in the space BV (0, T, L∞(Ω)) ∪ L∞(0, T ;BV (Ω)). Therefore, for the approximation (2.18) provided by the weakly nonlinear analysis to remain valid for all δ > 0, the corresponding amplitude A∞ must vanish in the limit δ → 0, making in this way compatible the increase of oscillations with its BV regularity. Theorem 2.3. Set α = (1, 4), and let Dδ(u) and Bb be given by (2.12) and (2.17), respectively, for b < 1/2 and 0 < δ < b/4. Then, there exists δc(b) > 0 such that if δ < δc(b) then u∗ = 1 1−b (1 − 2b, 2) is linearly unstable for problem (2.5)-(2.8). In addition, lim b→0 δc(b) = 0, lim b→0 kc(b) =∞, and the amplitude provided by the weakly nonlinear analysis satisfies lim b→0 A∞(b) = 0. In particular, the weakly nonlinear approximation v given by (2.18) satisfies v→ u∗ uniformly in Ω as b→ 0. 3. Numerical experiments To analyze the quality of the approximation provided by the WNA, and the properties stated in Theorems 2.1–2.3, we compare it to a numerical approximation of the evolution problem computed with the finite element method (FEM). For the FEM approximation, we used the open source software deal.II [1] to implement a time semi-implicit scheme with a spatial linear-wise finite element discretization. For the time discretization, we take in the experiments a uniform time partition of time step τ = 0.01. For the spatial discretization, we take a uniform partition of the interval Ω = (0, π) with spatial step depending on the predicted wave number of the pattern, see Table 1. Let, initially, t = t0 = 0 and set (u0 1, u 0 2) = (u10, u20). For n ≥ 1, the discrete problem is: Find un1 , u n 2 ∈ Sh such that 1 τ ( un1 − un−1 1 , χ)h + ( dδ11(un)∂xu n 1 + dδ12(un)∂xu n 2 , ∂xχ )h = ( f b1(un1 , u n 2 ), χ)h, (3.1) 1 τ ( un2 − un−1 2 , χ)h + ( dδ21(un)∂xu n 1 + dδ22(un)∂xu n 2 , ∂xχ )h = ( f b2(un1 , u n 2 ), χ)h, (3.2) for every χ ∈ Sh, the finite element space of piece-wise Q1-elements. Here, (·, ·)h stands for a discrete semi-inner product on C(Ω). EJDE-2021/55 TURING INSTABILITY ANALYSIS OF A CROSS-DIFFUSION PROBLEM 7 Since (3.1)-(3.2) is a nonlinear algebraic problem, we use a fixed point argument to approximate its solution, (un1 , u n 2 ), at each time slice t = tn, from the previous approximation (un−1 1 , un−1 2 ). Let un,01 = un−1 1 and un,02 = un−1 2 . Then, for k ≥ 1 the linear problem to solve is: Find (un,k1 , un,k2 ) such that for for all χ ∈ Sh 1 τ ( un,k1 − un−1 1 , χ)h + ( dδ11(un,k−1)∂xu n,k 1 + dδ12(un,k−1)∂xu n,k 2 , ∂xχ )h = ( un,k1 (αb1 − βb11u n,k−1 1 − βb12u n,k−1 2 ), χ)h, 1 τ ( un,k2 − un−1 2 , χ)h + ( dδ21(un,k−1)∂xu n,k 1 + dδ22(un,k−1)∂xu n,k 2 , ∂xχ )h = ( un,k2 (αb2 − βb21u n,k−1 1 − βb22u n,k−1 2 )χ)h. We use the fixed point stopping criteria max ( ‖un,k1 − un,k−1 1 ‖2, ‖un,k2 − un,k−1 2 ‖2 ) < tolFP , for values of tolFP chosen empirically, and set (un1 , u n 2 ) = (un,k1 , un,k2 ). Finally, we integrate in time until a numerical stationary solution, (uS1 , u S 2 ), is achieved. This is determined by max ( ‖un,11 − un−1 1 ‖2, ‖un,12 − un−1 2 ‖2 ) < tolS , where tolS is chosen empirically too. In the following experiments we always fix tolFP = 1.e− 07 and tolS = 1.e− 12. Experiment 1. We investigate the behavior of the instabilities arising in the solu- tions of the approximated problems (2.5)-(2.8) when δ → 0. Our main objective is checking if the predictions of the weakly nonlinear analysis stated in Theorem 2.3 are captured by the FEM approximation too. Thus, we use the diffusion matrix Dδ(u) and the competence parameters Bb given by (2.12) and (2.17), respectively. We run three simulations according to the choice of b, see Table 1, and fix δ = 0.95δc(b) in all of them, so that u∗ is unstable and pattern formation follows. Table 1. Data set for the Experiment 1. Wave numbers and times are rounded. Execution time measured for a standard laptop with i7 processor. Simulation 1 Simulation 2 Simulation 3 b 3.85e-02 9.91e-03 4.42e-03 δ(b) 4.53e-05 2.94e-06 5.83e-07 kc(b) 10 20 30 A∞(b) 1.21e-02 3.1e-03 1.4e-03 Number of nodes 128 256 512 Time steps to stationary 3.e+04 1.9e+05 4.4e+05 Execution time (hours) 1.67 19.26 84.77 In Figure 1 we show the typical onset and transmission of disturbances found in all the experiments. In this figure and in the following we plot only the first com- ponent of the solution, being the behavior of the second component similar. After a fast decay of the initial data towards the unstable equilibrium, a perturbation 8 G. GALIANO, V. GONZÁLEZ-TABERNERO EJDE-2021/55 with the wave number predicted by the linear analysis grows from one side of the boundary to the rest of the domain until reaching the steady state, see Figure 2. In the latter figure, we may check the good accordance between the FEM and the WNA approximations which, in numeric figures, have a relative difference of the order 10−5. In Figure 3 we show three interesting behaviors of solutions when δ → 0. In the left panel, the shrinking amplitude of the stationary patterns while the wave number increases. The equilibrium has been subtracted from the solution to center the pattern in y = 0. The center panel shows the time evolution of the amplitude (log scale) as given by the exact solution of the Stuart-Landau equation (4.17). We readily see that the stabilization time is a decreasing function of δ. This fact together with the increment of the wave number when δ → 0 results in very high execution times, see Table 1. Finally, the third panel shows how the variation of the numerical stationary solution∫ Ω |∂xu1(T, x)|dx is an increasing function of δ and tends to zero as δ → 0, in agreement with the regularity of solutions stated by the theoretical results. Experiment 2. We repeated Experiment 1 replacing the diffusion matrix Dδ(u) by that defined in (2.13). In Table 2 we show the relative differences in Lp, given by RDp(φ1, φ2) = ‖φ1 − φ2‖Lp ‖φ1‖Lp , (3.3) in terms of the critical bifurcation parameter, δc, the stationary solution of the FEM approximation, u(T, ·), the WNA approximation, v, and the pattern ampli- tude, A∞, corresponding to both approximations of the original diffusion matrix. We see that although the critical bifurcation parameter is clearly affected by the ap- proximation scheme, the FEM and WNA approximations provided by both schemes are in a very good agreement, as well as the amplitudes of the instability patterns, suggesting that in the limit δ → 0 both sequences of approximations converge to the same limit. Table 2. Comparison between the results obtained with the ap- proximated diffusion matrices corresponding to Experiments 1 and 2, given by (2.12) and (2.13) respectively. RDp denotes the relative difference in Lp, see (3.3). Simulation 1 Simulation 2 Simulation 3 RD∞(δ (E1) c , δ (E2) c ) 0.136 0.117 0.113 RD2(u(E1)(T, ·),u(E2)(T, ·)) 3.74e-06 8.68e-07 5.41e-07 RD2(v(E1),v(E2)) 3.46e-06 2.13e-07 4.19e-08 RD∞(A (E1) ∞ , A (E2) ∞ ) 2.90e-03 6.82e-04 2.99e-4 EJDE-2021/55 TURING INSTABILITY ANALYSIS OF A CROSS-DIFFUSION PROBLEM 9 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.9894 0.9897 0.9900 0.9903 0.9906 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Figure 1. Typical evolution of disturbances in Simulation 2. 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.9576 0.9588 0.9600 0.9612 0.9623 WNA FEM 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.9894 0.9897 0.9900 0.9903 0.9906 WNA FEM 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.9953 0.9954 0.9956 0.9957 0.9958 WNA FEM Figure 2. Experiment 1. WNA and FEM approximations corre- sponding to Simulations 1 to 3 (left to right). Notice the different scales in the ordinates axis showing the decreasing amplitude of the oscillations. 0.0 0.5 1.0 1.5 2.0 2.5 3.0 -0.0024 -0.0012 0.0000 0.0012 0.0024 Patterns around equilibrium = 4.53e-05 = 2.94e-06 = 5.84e-07 0 200 400 600 800 1000 t 8 7 6 5 4 3 lo g( A) Evolution of the amplitude = 4.53e-05 = 2.94e-06 = 5.84e-07 0.0000 0.0005 0.0010 0.0015 0.0020 0.00 0.02 0.04 0.06 0.08 0.10 Variation of the patterns Figure 3. Experiment 1. Behavior of the patterns as δ → 0. 4. Proofs We use the decomposition of the nonlinear problem (2.5)-(2.8) in terms of its linear and nonlinear parts. Let v = u − u∗, where u is a solution of (2.5)-(2.8). Then, v satisfies ∂tv = Lδv +N δ(v), (4.1) where we split the reaction-diffusion terms into their linear parts Lδv = Dδ(u∗)∂xxv +Kv, with K given by (2.14), and their nonlinear parts N δ = ∂x ( Dδ(v)∂xv ) +˜b(v), (4.2) being f̃ bi (v) = −βbiiv2 i − βbijvivj , for i, j = 1, 2 and i 6= j. 10 G. GALIANO, V. GONZÁLEZ-TABERNERO EJDE-2021/55 Proof of Theorem 2.1. We study the linearization of (4.1), this is, the equation ∂tw = Lδw, (4.3) satisfying Neumann homogeneous boundary conditions and with initial data w0 = u0 − u∗. This linear problem is well-posed due to the second assumption of (H1). The type of boundary conditions lead to seek for solutions of the form w = eλt cos(kx)w, with k = 1, 2, . . ., where w is a constant vector. Replacing w in (4.3) we obtain the matrix eigenvalue problem Akw = λw, with Ak = K − k2Dδ(u∗). Since, by hypotheses, tr(Ak) = tr(K) − k2 tr(Dδ(u∗)) < 0 for all k = 0, 1, . . . ; an eigenvalue with positive real part (instability) may exist only if det(Ak) is negative for some wave number k. We introduce the notation h(k2) = det(Ak): h(k2) = det(Dδ(u∗))k4 + qδ(u ∗)k2 + det(K), where qδ(u ∗) = dδ11(u∗)βb22u ∗ 2 +dδ22(u∗)βb11u ∗ 1− (dδ12(u∗)βb21u ∗ 2 +dδ21(u∗)βb12u ∗ 1). The minimum of the convex parabola h is attained at k2 m(δ) = − qδ(u ∗) 2 det(Dδ(u∗)) , requiring qδ(u ∗) < 0, which is true in view of (2.15). A necessary condition for linear instability is h(k2 m(δ)) < 0, where h(k2 m(δ)) = det(K)− qδ(u ∗)2 4 det(Dδ(u∗)) . In this expression, det(K) is a positive constant and qδ(u ∗)2 > 0 for all δ ≥ 0. Thus, since qδ(u ∗) and Dδ(u∗) are monotone with respect to δ and det(Dδ(u∗)) → 0 as δ → 0, we deduce the existence of an unique δ̄c > 0 such that h(k2 m(δ̄c)) = 0. Therefore, for δ < δ̄c we have h(k2(δ)) < 0 if k2(δ) ∈ (k2 −(δ), k2 +(δ)), where k2 ±(δ) = −qδ(u∗)± √ (qδ(u∗))2 − 4 det(Dδ(u∗)) det(K) 2 det(Dδ(u∗)) . Because of the boundary conditions, the onset of instabilities only occurs when one of the extreme values of the interval (k−(δ), k+(δ)) is an integer number. Since k+(δ) → ∞ as δ → 0, this will certainly hold for δ small enough. We define the critical bifurcation parameter, δc, as such number, and the critical wave number, kc ∈ Z, as the corresponding root of h(k2). Finally, the last assertion of the theorem is a consequence of the infinite limit of k+(δ) as δ → 0. � Proof of Theorem 2.2. We retake the whole nonlinear equation (4.1) for v = u−u∗. The idea of the weakly nonlinear analysis is to look for an approximation of v for a value of δ near the critical bifurcation parameter δc. This approximation is defined as an expansion in terms of the small parameter ε2 = (δc − δ)/δc, for δ < δc. We consider the expansions δ = δc − εδ1 − ε2δ2 − ε3δ3 +O(ε4), t = εt1 + ε2t2 + ε3t3 +O(ε4), v = εv1 + ε2v2 + ε3v3 +O(ε4), and then introduce these expressions in equation (4.1) and collect the resulting equations in terms of powers of ε. Since this procedure is standard, we give the EJDE-2021/55 TURING INSTABILITY ANALYSIS OF A CROSS-DIFFUSION PROBLEM 11 results and omit intermediate calculations for the sake of brevity. We obtain Order ε: Lδcv1 = 0. (4.4) Order ε2: Lδcv2 = ∂t1v1 + δ1D 1(u∗)∂xxv1 − 1 2 ( QK(v1,v1) + ∂xxQDδc (v1,v1) ) − SDδc (v1) =: F. (4.5) Order ε3: Lδcv3 = (∂t1v2 + ∂t2v1) +D1(u∗)∂xx(δ1v2 + δ2v1) −QK(v1,v2)− ∂xxQDδc (v1,v2) + 1 2 δ1∂xxR1(v1)−R2(v1,v2) + δ1R3(v1) =: G. (4.6) Here, D1(u∗) is given by (2.9), QK(x,y) = − ( 2βb11x1y1 + βb12(x1y2 + x2y1) 2βb22x2y2 + βb21(x1y2 + x2y1) ) , QDδ(x,y) = ( dδ111x1y1 + dδ212x2y2 dδ121x1y1 + dδ222x2y2 ) , SDδ(v) = ∂x ( dδ211v2∂xv1 + dδ112v1∂xv2 dδ221v2∂xv1 + dδ122v1∂xv2 ) , R1(v1) = ( d11 11(v11)2 + d21 12(v12)2 d11 21(v11)2 + d21 22(v12)2 ) , R3(v1) = ∂x ( d21 11v12∂xv11 + d11 12v11∂xv12 d21 21v12∂xv11 + d11 22v11∂xv12 ) , R2(v1,v2) = ∂x ( dδc211 (v12∂xv21 + v22∂xv11) + dδc112 (v11∂xv22 + v21∂xv12) dδc221 (v12∂xv21 + v22∂xv11) + dδc122 (v11∂xv22 + v21∂xv12) ) , where we introduced the notation dδmij = dm0 ij + δdm1 ij , for m = 1, 2 so that dδij(v) = dδ1ij v1 + dδ2ij v2. Observe that dδ1ij , d δ2 ij are the elements of the matrices Dδ1, Dδ2 introduced in the first assumption of (H1). Observe also that (4.2) may be written as N δv = 1 2 ( QK(v,v) + ∂xxQDδ(v,v) ) + SDδ(v). We now compute the solutions corresponding to each order in the expansion. Order ε: The solution of (4.4) is given by v1(t1, t2, x) = A(t1, t2)ρ cos(kcx), with ρ ∈ ker(K − k2 cD δc(u∗)), where A is the amplitude of the pattern, unknown at the moment. Observe that ker(Aδckc) is a one-dimensional subspace, implying that the vector ρ is defined up to a multiplicative constant. We shall fix this constant later. Order ε2: We start by expressing F in terms of A and ρ. We have ∂t1v1 = ∂t1A cos(kcx)ρ δ1D 1(u∗)∂xxv1 = −δ1Ak2 c cos(kcx)D1(u∗)ρ 12 G. GALIANO, V. GONZÁLEZ-TABERNERO EJDE-2021/55 On noting that QU (v1,v1) = A2QU (ρ,ρ) cos2(kcx), for U = K, Dδc , we find 1 2 ( QK(v1,v1) + ∂xxQDδc (v1,v1) ) = 1 4 A2 ∑ j=0,2 Mj(ρ,ρ) cos(jkcx), with Mj = QK − j2k2 cQDδc . Using standard trigonometric identities, we obtain SDδc (v1) = −k2 cA 2ρ1ρ2 cos(2kcx)d, where d = (dδc211 + dδc112 , d δc2 21 + dδc122 ). Gathering the above expressions, we obtain F = [ ∂t1Aρ− δ1Ak2 cD 1(u∗)ρ ] cos(kcx)− 1 4 A2 ∑ j=0,2 Mj(ρ,ρ) cos(jkcx) + k2 cA 2ρ1ρ2d cos(2kcx). By Fredholm’s alternative, (4.5) admits a solution if and only if 〈F,ψ〉L2 = 0, where 〈·, ·〉L2 denotes the scalar product in L2(0, π), and ψ ∈ ker((Lδc)∗) is of the form ψ = η cos(kcx), with η ∈ ker((K − k2 cD δc(u∗))∗). (4.7) Observe that η, for similar reasons than ρ, is defined up to a multiplicative constant. We fix η at the end of this proof, and also show that 〈ρ,η〉 6= 0. The compatibility condition implies ∂t1A(t1, t2) = δ1k 2 c 〈D1(u∗)ρ,η〉 〈ρ,η〉 A(t1, t2). Since the solution to this equation is an exponential function, we do not obtain from it any useful indication on the asymptotic behavior of the pattern amplitude. Therefore, to suppress the secular terms appearing in F, we impose t1 ≡ 0 and δ1 ≡ 0. (4.8) In particular, this implies A ≡ A(t2). Assuming these restrictions, the Fredholm’s alternative is satisfied, and moti- vated by the functional form of F, we seek for a solution of (4.5) of the form v2(t2, x) = A2(t2) ∑ j=0,2 v2j cos(jkcx), where v2j are constant vectors. The linear operator Lδc may be decomposed as Lδcv2 = A2 ∑ j=0,2 Ljv2j cos(jkcx), with Lj = K − j2k2 cD δc(u∗). Then, Lδcv2 = F if the vectors v2j are the solutions of the linear systems L0v20 = −1 4 M0(ρ,ρ), L2v22 = k2 cρ1ρ2d− 1 4 M2(ρ,ρ). Order ε3: We have to solve Lδcv3 = G, where, taking into account (4.8), G = ∂t2v1 + δ2D 1(u∗)∂xxv1 −QK(v1,v2)− ∂xxQDδc (v1,v2)−R2(v1,v2). Replacing the solutions obtained for the orders ε and ε2, i.e. v1 = A(t2)ρ cos(kcx) and v2 = A(t2)2(v20 + v22 cos(2kcx)) in G yields G = ( ρ∂t2A−Ak2 cδ2D 1(u∗)ρ−A3 ( M1(ρ,v20) + 1 2 M1(ρ,v22) + k2 cR1 )) cos(kcx) −A3 (1 2 M3(ρ,v22) + k2 cR2 ) cos(3kcx), EJDE-2021/55 TURING INSTABILITY ANALYSIS OF A CROSS-DIFFUSION PROBLEM 13 where R (i) 1 = dδc2i1 [ ρ1 (1 2 v (2) 22 − v (2) 20 ) − ρ2v (1) 22 ] + dδc1i2 [ ρ2 (1 2 v (1) 22 − v (1) 20 ) − ρ1v (2) 22 ] , R (i) 2 = −3 ( dδc2i1 [ ρ2v (1) 22 + 1 2 ρ1v (2) 22 ] + dδc1i2 [ ρ1v (2) 22 + 1 2 ρ2v (1) 22 ]) . The solvability condition for problem (4.6) is 〈G,ψ〉L2 = 0, with ψ = η cos(kcx) given by (4.7). This condition leads to the differential equation 〈ρ,η〉∂t2A = 〈G1,η〉A+ 〈G2,η〉A3, where G1 = k2 cδ2D 1(u∗)ρ, G2 =M1(ρ,v20) + 1 2 M1(ρ,v22) + k2 cR1, (4.9) Thus, we deduce the cubic Stuart-Landau equation for the amplitude ∂t2A = σA− `A3, (4.10) with σ = 〈G1,η〉 〈ρ,η〉 , ` = −〈G2,η〉 〈ρ,η〉 . (4.11) We, finally, fix the vectors ρ ∈ ker(K − k2 cD δc(u∗)) , and η ∈ ker((K − k2 cD δc(u∗))∗). Since all the elements of both matrices are negative, we may set ρ = (1,M)t and η = (1,M∗)t for some M,M∗ < 0, implying 〈ρ,η〉 > 0. Thus, the asymptotic behavior of the solution to (4.10) is fully determined by the signs of the numerators in (4.11). When σ and ` are positive, the amplitude stabilizes to a positive value, this is, A(t2) → A∞ := √ σ/` as t2 → ∞. Therefore, in this case, the corresponding solution v = εv1 + ε2v2 +O(ε3), is given by v = ερ √ σ ` cos(kcx) + ε2 σ ` ( v20 + v22 cos(2kcx) ) +O(ε3). An example of this situation is studied in Theorem 2.3. � Proof of Theorem 2.3. Our objective is to compute the coefficients of the Stuart- Landau equation (4.10). Specifically, we are interested in the ratio σ ` = −〈G1,η〉 〈G2,η〉 . Determination of 〈G1,η〉. For the given data, we obtain qδ(u ∗) = −u∗1u∗2 ( b 2 − 2δ ) , which is negative if b > 4δ. The corresponding roots of h(k2 m) are positive and, therefore, we take δc = δ−, so that for any δ < δc we have h(k2 m) < 0. The corresponding critical wave number is the minimum of h(k2), given by k2 c = b− 4δc 4δc(2 + δc) . The vectors ρ = (1,M) and η = (1,M∗) are elements of ker(Aδckc) and ker(Aδckc) ∗, respectively. Thus, M = −1 + k2 c (1 + δc) b 2 + k2 c M∗ = −u ∗ 1(1 + k2 c (1 + δc)) u∗2(2 + k2 c ) . 14 G. GALIANO, V. GONZÁLEZ-TABERNERO EJDE-2021/55 From (4.9), we obtain G1 = k2 cδ2(u∗1, u ∗ 2M), and then 〈G1,η〉 = k2 cδ2(u∗1+u∗2MM∗). � Lemma 4.1. Let εM = 1 +M . Then we have: lim b→0 δc(b) = lim b→0 δ′c(b) = lim b→0 δck 2 c = 0, lim b→0 k2 c =∞, (4.12) lim b→0 M = −1, lim b→0 M∗ = −1 2 , (4.13) lim b→0 k2 c εM = −1, lim b→0 k2 c (1 + 2M∗) = 9. (4.14) Taking into account that δ2 ≈ δc, a first consequence of Lemma 4.1 is lim b→0 〈G1,η〉 = 0. (4.15) Proof of Lemma 4.1. For proving (4.12), we use L’Hôpital’s rule to obtain lim b→0 δc(b) = 1 4 lim b→0 ( − 3− −6(4− 3b)− 3b2 2 √ (4− 3b)2 − b3 ) = 0. Let φ(b, δ) = 4 det(Dδ(u∗) det(K) − q2 δ (u∗) = −4bδ2 + (8 − 6b)δ − b2/4. By the definition of δc, we have φ(b, δc(b)) = 0 for all b ∈ (0, 1). Thus 0 = d db φ(b, δc(b)) = ∂bφ(b, δc(b)) + ∂δφ(b, δc(b))δ ′ c(b). (4.16) Since ∂bφ(0, 0) = 0 and ∂δφ(0, 0) = 8, we deduce that δ′c(0) = 0. We then have lim b→0 k2 c = lim b→0 b− 4δc 4δc(2 + δc) = lim b→0 1− 4δ′c 8δ′c(1 + δc)) =∞, lim b→0 δck 2 c = lim b→0 b− 4δc 4(2 + δc) = 0. The limits in (4.13) follow easily from the definitions of M and M∗. Finally, to prove (4.14), we use the definition of M to obtain k2 c εM = k2 c b 2 + k2 c − 1− k2 c (1 + δc) b 2 + k2 c = b 2 − 1− k2 cδc b 2k2c + 1 , from where the first limit follows. The second limit is computed in a similar way. We write k2 c (1 + 2M∗) = k2 c ( 1 + u∗2 u∗1 M∗ +M∗ ( 2− u∗2 u∗1 )) . On the one hand, we have k2 c ( 1 + u∗2 u∗1 M∗ ) = k2 c 2 + k2 c − (1 + k2 c (1 + δc)) 2 + k2 c = k2 c 1− k2 cδc 2 + k2 c → 1 as b→ 0. On the other hand, using the definition of u∗ and k2 c , we obtain k2 c ( 2− u∗2 u∗1 ) = − 1 1− 2b 4k2 cb, with 4k2 cb = 1 2 + δc (b2 δc − 4b ) . Using a concatenation of L’Hôpital’s rule, limb→0 b 2/δc(b) = 2 limb→0 1/δ′′c (0). Determination of 〈G2,η〉. The following lemma gives the expression of this scalar product. Since the calculation is straightforward, we omit the details. EJDE-2021/55 TURING INSTABILITY ANALYSIS OF A CROSS-DIFFUSION PROBLEM 15 Lemma 4.2. We have −〈G2,η〉 = S1v (1) 20 + S2v (2) 20 + T1v (1) 22 + T2v (2) 22 , where v20 = 1 4u∗1u ∗ 2(1− b) ( −u∗2(2 + bM) + bu∗1M(M + 2) 2u∗2(2 + bM)− 2u∗1M(M + 2) ) , v22 = 1 det(L2) ( w2k 4 c + w1k 2 c + w0 ) , with w2 = 4 ( −(1 + δc)u ∗ 2(εM + δc) + u∗1M(εM +Mδc) u∗2(εM + δc)− (1 + δc)u ∗ 1M(εM +Mδc) ) , w1 = ( −4(1 + δc)u ∗ 2( 1 2 + bM 4 )− u∗2(εM + δc) + 4u∗1(M + M2 2 ) + b 2u ∗ 1M(εM +Mδc) 4u∗2( 1 2 + bM 4 ) + 2u∗2(εM + δc)− 4(1 + δc)u ∗ 1(M + M2 2 )− u∗1M(εM +Mδc) ) , w0 = ( −u∗2( 1 2 + bM 4 ) + b 2u ∗ 1(M + M2 2 ) 2u∗2( 1 2 + bM 4 )− u∗1(M + M2 2 ) ) , being det(L2) = 9u∗1u ∗ 2(2 + δc)k 4 cδc, Ti = T (i) 0 + T (i) 1 k2 c , for i = 1, 2, and S1 = 2 + Mb 2 + 2MM∗ + k2 c (εM + δc), S2 = b 2 + 2M∗εM + k2 cM ∗(εM +Mδc), T (1) 0 = 1 2 (2 + Mb 2 + 2MM∗), T (1) 1 = 1 2 (1−M + 2MM∗ + δc), T (2) 0 = b 4 +M∗εM , T (2) 1 = 1 + M∗ 2 (M − 1 +Mδc). Lemma 4.3. There exists a constant C < 0 such that k2 cδc〈G2,η〉 → C as b→ 0. Consequently, 〈G2,η〉 → −∞ as b→ 0. This result together with (4.15) implies that the solution of the Stuart-Landau equation (4.10), given by A2(t) = σ ` 1 1 + ( A−2 0 σ ` − 1 ) e−2σt , A0 = A(0), (4.17) satisfies A→ 0 uniformly in (0,∞) as b→ 0, which completes the proof. � Proof of Lemma 4.3. We set ζk2 cδc(T1v (1) 22 + T2v (2) 22 ) = I4k 4 c + I2k 2 c + I0 + I−2k −2 c , with ζ = 9u∗1u ∗ 2(2 + δc), I4 = w (1) 2 T (1) 1 + w (2) 2 T (2) 1 , I−2 = w (1) 0 T (1) 0 + w (2) 0 T (2) 0 , I2 = w (1) 2 T (1) 0 + w (1) 1 T (1) 1 + w (2) 2 T (2) 0 + w (2) 1 T (2) 1 , I0 = w (1) 1 T (1) 0 + w (1) 0 T (1) 1 + w (2) 1 T (2) 0 + w (2) 0 T (2) 1 . 16 G. GALIANO, V. GONZÁLEZ-TABERNERO EJDE-2021/55 As b → 0, we have w0 → (−1, 3/2) and (T (1) 0 , T (2) 0 ) → (3/2, 0), implying that I−2k −2 c → 0 as b→ 0. For I4, we have 1 4 I4 = ( (−u∗2 + u∗1M)εM + (−u∗2 + u∗1M 2)δc )(M∗ − 1 2 εM + 1−M∗M 2 δc ) +O(εMδc) +O(δ2 c ). Since k4 cδ 2 c → 0, k4 cδcεM → 0 and k4 c ε 2 M → 1, we deduce lim b→0 I4k 4 c = 4 lim b→0 (−u∗2 + u∗1M) M∗ − 1 2 = 9. For I2, we have w (1) 1 T (1) 1 + w (2) 1 T (2) 1 = ( − 2u∗2 + 4u∗1M ( 1 + M 2 ) − u∗2bM )(M∗ − 1 2 εM + 1−M∗M 2 δc ) + ( − u∗2T (1) 1 + (2u∗2 − u∗1M)T (2) 1 ) εM +O(δc) and w (1) 2 T (1) 0 + w (2) 2 T (2) 0 = (−u∗2 + u∗1M)εM (1−M∗) +O(δc). Since k2 cδc → 0 and k2 c εM → −1, we have lim b→0 I2k 2 c = − lim b→0 [( − 2u∗2 + 4u∗1M ( 1 + M 2 ))M∗ − 1 2 + (−u∗2T (1) 1 + (2u∗2 − u∗1M)T (2) 1 ) + (−u∗2 + u∗1M)(1−M∗) ] = 5 2 . Finally, for I0, we have w (1) 0 T (1) 1 + w (2) 0 T (2) 1 = ( − u∗2 2 +O(b) ) T (1) 1 + ( u∗2 − u∗1(M + M2 2 ) +O(b) ) T (2) 1 , w (1) 1 T (1) 0 + w (2) 1 T (2) 0 = ( − 2u∗2 + 4u∗1M ( 1 + M 2 )) (T (1) 0 − T (2) 0 ) +O(b) +O(εM ) +O(δc) = O(b) +O(εM ) +O(δc). Therefore, lim b→0 I0 = lim b→0 ( − u∗2 2 T (1) 1 + ( u∗2 − u∗1(M + M2 2 ) ) T (2) 1 ) = 9 4 implying limb→0 k 2 cδc(T 1v (1) 22 + T 2v (2) 22 ) = 55/288. � Acknowledgments. G. Galiano was supported by the Spanish MEC Project no. MTM2017-87162-P. References [1] W. Bangerth, T. Heister, L. Heltai, G. Kanschat, M. Kronbichler, M. Maier, B. Turcksin; The deal.II Library, Version 8.3, Arch. Numer. Software, 4(100) (2016), 1-11. [2] M. Bertsch, R. Dal Passo, M. Mimura; A free boundary problem arising in a simplified tumour growth model of contact inhibition, Interfaces and Free Bound., 12 (2010), 235–250. [3] M. Bertsch, D. Hilhorst, H. Izuhara, M. Mimura; A nonlinear parabolic-hyperbolic system for contact inhibition of cell-growth, Differ. Equ. Appl., 4(1) (2012), 137–157. [4] S. N. Busenberg, C. C. Travis; Epidemic models with spatial spread due to population mi- gration, J. Math. Biol., 16 (1983), 181–198. [5] L. Chen, E. S. Daus, A. Jüngel; Rigorous mean-field limit and cross-diffusion, Z. Angew. Math. Phys., (2019), 70–122. EJDE-2021/55 TURING INSTABILITY ANALYSIS OF A CROSS-DIFFUSION PROBLEM 17 [6] L. Chen, A. Jüngel; Analysis of a multidimensional parabolic population model with strong cross-diffusion, SIAM J. Math. Anal., 36 (2004), 301–322. [7] G. Galiano, M. L. Garzón, A. Jüngel; Semi-discretization in time and numerical convergence of solutions of a nonlinear cross-diffusion population model, Numer. Math., 93 (2003), 655– 673. [8] G. Galiano, V. Selgas; On a cross-diffusion segregation problem arising from a model of interacting particles, Nonlinear Anal. Real World Appl., 18 (2014), 34–49. [9] G. Galiano, V. Selgas; Deterministic particle method approximation of a contact inhibition cross-diffusion problem, Appl. Numer. Math., 95 (2015), 229–237. [10] G. Gambino, M. C. Lombardo, M. Sammartino; Turing instability and traveling fronts for a nonlinear reaction–diffusion system with cross-diffusion, Math. Comput. Simul., 82 (2012), 1112–1132. [11] G. Gambino, M. C. Lombardo, M. Sammartino; Pattern formation driven by cross-diffusion in a 2D domain, Nonlinear Anal. Real World Appl., 14 (2013), 1755–1779. [12] A. Jüngel; The boundedness-by-entropy method for cross-diffusion systems, Nonlinearity, 28 (2015), 1963. [13] N. Shigesada, K. Kawasaki, E. Teramoto; Spatial segregation of interacting species, J. The- oret. Biol., 79 (1979), 83–99. [14] A. M. Turing; The chemical basis of morphogenesis; Philosophical Trans. Royal Soc. London. Series B, Biol. Sciences, 237 (1952), 37–72. Gonzalo Galiano Department of Mathematics, University of Oviedo, Oviedo, Spain Email address: galiano@uniovi.es V́ıctor González-Tabernero Department of Mathematics, University of Santiago de Compostela, Santiago de Com- postela, Spain Email address: victor.gonzalez.tabernero@rai.usc.es 1. Introduction 2. Main results 3. Numerical experiments Experiment 1 Experiment 2 4. Proofs Acknowledgments References