# Adaptive Optimalm-Stage Runge-Kutta Methods for Solving Reaction-Diffusion-Chemotaxis Systems

Authors: Jui-Ling Yu
Journal: Journal of Applied Mathematics
 (2011)

---

## Abstract

We present a class of numerical methods for the reaction-diffusion-chemotaxis system which is significant for biological and chemistry pattern formation problems. To solve reaction-diffusion-chemotaxis systems, efficient and reliable numerical algorithms are essential for pattern generations. Along with the implementation
of the method of lines, implicit or semi-implicit schemes are typical time stepping solvers to reduce the effect on time step constrains due to the stability condition. However, these two schemes are usually difficult to employ. In this paper, we propose an adaptive optimal time stepping strategy for the explicitm-stage Runge-Kutta method to solve reaction-diffusion-chemotaxis systems. Instead of relying on empirical approaches to control the time step size, variable
time step sizes are given explicitly. Yet, theorems about stability and convergence of the algorithm are provided in analyzing robustness and efficiency. Numerical experiment results on a testing problem and a real application problem are shown.


---

## Body

## 1. Introduction and Background

The reaction-diffusion-chemotaxis system is the most common model used to describe chemical and biological processes and the most successful model used to describe the generation of pattern formations. Numerical simulation of such system is necessary because complicated chemotaxis terms and highly nonlinear reaction terms increase the difficulty on finding analytical solutions for realistic systems. Even though, some of the early chemotaxis models had analytical solution or one-dimensional simulation with reduced models [1, 2], the development of numerical methods for solving chemotaxis systems does not seem sufficient comparing with other type of PDEs systems. Since interesting patterns may emerge and evolve before they reach steady state [3–8], some particular heed needs to be paid to the numerical methods to ensure that the numerical simulated patterns literally obey primary model equations [5, 8].Chemotaxis is a directed motion, shown in motile organism, which moves toward or away from certain chemical or environment. It has been observed in microorganism, plants, animals, even in blood vessels of cancer tumors. For bacteria, chemotaxis was first discovered in the lab by Wilhelm Pfeffer and his associates inE. coli. in 1884 [9]. Besides the random motion of bacteria, they also show directed random movements in response to some higher concentration of oxygen. Keller and Segel [1] construct the first differential chemotaxis models. They extend the reaction-diffusion model with a convection term to simulate chemotaxis. Based on their work, many models are constructed, especially on the development of bacteria colonies [3–8, 10].The general form of the reaction-diffusion-chemotaxis is given by [11].  (1.1)ut=DuΔu+∇(uχ(p)∇p)+f(p,u).Here,u represents the microorganism (bacteria) density; p stands for chemoattractant concentration and is given or obtained from another equation. f is the reaction term which models the interactions between the chemicals. The diffusion term, DuΔu, imitates the random motion of particles with positive constant diffusivity Du. ∇·(uχ(p)∇p) describes the directed random movement of particles p that are modulated by the concentration gradient of the chemoattractant. The chemotaxis parameter function, χ, indicates the strength of the chemotaxis, which is given by the Keller-Segal model.(1.2)χ(x)=C0(C1+x)s,
where C0 and C1 are positive constants and s≥1 is an integer power. Observe that limx→∞χ(x)=0 which indicates that particles tend to stay where they are when the local concentration of chemoattractant is large.The basic idea for illustrating patterns using reaction-diffusion-chemotaxis systems is to demonstrate that eigenvalues of local linearization that are responsible for the growth (decay) can have positive (negative) real parts. Therefore, it can lead to growth at some places and decay at others, resulting in spatially inhomogeneous patterns [12, 13]. Based on this idea, we consider the method of lines (MOL) approach in this paper because it is not only the first step for construction of many numerical methods, but also provides good explanations about the generation of patterns through analyzing the positive (growth modes) and negative (decay modes) of eigenvalues. Other numerical methods used on the chemotaxis model include the fractional step method [14], the Alternating Direction Implicit (ADI) method [15, 16], and the optimal two-stage scheme [5, 8] listed here for readers's references.The MOL approach reduces the partial differential equations to a system of ordinary differential equations (ODEs) by approximating spatial derivatives with numerous spatial discretization methods, for example, the finite difference, finite element, and spectral methods. The solution can be obtained after following some suitable time integrator. According to various spatial discretization features, the computation work can be facilitated by developing some more sophisticated robust and efficient adaptive time integrators [5, 8].In this paper, we use standard finite differences of center difference for the Laplacian operator and upwind differences for the chemotaxis term. The reason for adapting upwind differences is the stability concern. When the chemotaxis dominates, the system behaves more like a hyperbolic system. Then the upwind differences are more stable methods [17]. It is worth to mention that the following time-stepping scheme can also be applied to other discretization methods when the domain containing the eigenvalues is known or can be estimated. The general form of the resulting ODE system after spatial discretization reads(1.3)dudt=M(t)u+b(t),
where M(t) is a time-dependent matrix and b(t) is a reaction vector which also depends on time.Typically, the time integrator for (1.3) can be classified as implicit or explicit methods. Generally, the numerical scheme for solving implicit method is more complicated than explicit method since it often involves solving systems of algebraic equations or needs to deal with the inverse of matrices. Although implicit methods allow larger step size than explicit ones for those spatial discretization matrices only having eigenvalues with negative real parts, the computational cost may still be high due to the complicated implementation for implicit methods. Most importantly, the problem regarding time constrain still exists if the spatial matrices contain very positive real parts of eigenvalues; even for the implicit method the use of very small time steps is yet unpreventable in order to meet the stability condition. This situation is likely to occur for pattern formation problems for which the spatial discretization matrices comprise both positive and negative real parts for their eigenvalues [5, 8]. Besides, for pattern formation problems, those positive real eigenvalues are the unpredictable ones which evolve along with time. Bear these factors in mind, the explicit methods are actually better ones for the reaction-diffusion-chemotaxis systems.In the present paper, we start with reviewing the explicitm-stage Runge-Kutta methods and then construct their adaptive optimal time steps to enhance the computation efficiency for the chemotaxis system. The optimal time step formulas are presentable and easy to calculate. Only some simple estimations about the bounds for the largest and smallest eigenvalues of the discretization spatial matrices are required at each time step. The numerical simulations on the test and realistic problems are provided. Comparison will be made in the testing problem to demonstrate that our method is not only stable but also more efficient than other similar methods, referring to CPU time and the number of iterations.


## 2. Optimal Time Steps for Explicit Runge-Kutta Methods

### 2.1. Time Independence ODEs Systems

Since interesting patterns are caused by the interaction between growth and decay modes, to analyze system (1.3) and relate it to the pattern formation model (1.1), one can do the local linearization property of the right-hand side of (1.3). The positive and negative real eigenvalues of the linearization matrix then, respectively, correspond to the growth and decay modes. Because of this reason, we will consider the corresponding ODEs induced from the chemotaxis system and introduce our idea by starting with the constant coefficients case, that is,(2.1)dudt=Mu+b
then generalize it to the time-dependent coefficients case. The general formulation of the m-stage Runge-Kutta method is(2.2)k1=f(tn,un),kq=f(tn+αqΔt,un+Δtn∑s=1q-1βqsks),q=2,…,m,αq=∑s=1q-1βqs,un+1=un+Δtn∑s=1qγsks,q=2,…,m,
where m is the stage number and αq=∑j=1q-1βqj,q=2,…,m. The subscripts n and n+1 are the iteration level. For simplicity, we first look at the uniform time step size Δt=tn-tn-1 for all n. Then un≈u(xn) is the numerical solution at the nth step. Equations for the unknowns γs and βqs can be obtained by matching terms with Taylor's expansion.We choose single-step explicit methods, which basically are Runge-Kutta methods, because they are easy to implement and especially efficient for ODEs systems derived from the PDEs. The method is set to be explicit since the range of eigenvalues with negative real parts is measurable for chemotaxis systems. In [5, 8], we have developed adaptive two-stage time stepping scheme with optimum time step sizes for solving chemotaxis system and verified that it is reliable and more economical when comparing with other similar methods. We conjecture that the similar strategy used in [5, 8] could be used to derive optimal adaptive three- and four-stage methods. In the two-stage case, the size of the optimal step is determined by the root of the first-order equation hence it is absolute optimal within one time step. If our conjecture is correct, the sizes of optimal steps for higher-stage methods should be the roots of some higher-order equations. But, are these adaptive step sizes absolute optimal? And are they easy to be found as in the two-stage case?To answer these questions, we begin our study with the three-stage Runge-Kutta method and setb=0 in (2.1). As in [5, 8], we first fix all the intermediate coefficients at t=tn for the three-stage Runge-Kutta method. The general time dependence case is proposed in the next section. Denoting M=M(tn) and applying (2.2) to the three-stage Runge-Kutta method with all the intermediate values being fixed at tn, we arrive at(2.3)k1=Mun,k2=M(1+α2ΔtM)un,k3=M(1+α3ΔtM+α2β32Δt2M2)un,(2.4)un+1=(1+δ1ΔtM+δ2Δt2M2+δ3Δt3M3)un,
where δi,i=1,2,3, are combinations of the coefficients βqs and γs. Since u(tn+1)=eMΔtu(tn) and assuming un is exact up to the order Δt3, we have δ1=1, δ2=(1/2!), δ3=(1/3!). Similar arguments also hold for RK methods with other stages. Thus, if we treat it as a standard iterative method, (2.4) is a one-step method [18]:(2.5)un+1=[1+ΔtM+12!(ΔtM)2+⋯+1m!(ΔtM)m]un.=pm(Δt)un.In the case ofb≠0, (2.5) becomes(2.6)un+1=pm(Δt)un+(1+Δt2!M+⋯+Δtmm!Mm-1)b.Denotev as the solution by solving the system with the perturbated initial condition; we have(2.7)vn+1=pm(Δt)vn+(1+Δt2!M+⋯+Δtmm!Mm-1)b.From (2.6) and (2.7), the accumulation error can be written as (2.8)|un+1-vn+1|=|pm(Δt)|×|un-vn|.Note that(2.9)pm(λ)=1+Δtλ+12!(Δtλ)2+⋯+1m!(Δtλ)m
is called the amplification factor, where λ is an eigenvalue of M. It is obtained when a scalar test problem, u′=λu, is applied to the single-step method.Suppose thatλ=a+bi=Reλ+Imλi. For pattern formation problems, what we are interested in is the long-term behavior of the solution vectors. It is important to ensure that observed patterns are accurate in the sense that they are not coming from the accumulation error of numerical schemes. We therefore make the following definition of stability.Definition 2.1.
An explicit single-step method is said to be stable if|pm(λ)|<1 when a<0 and |pm(λ)|≤1+aΔt+O(Δt) when a≥0.HereO(Δt) means higher order of Δt. Our definition is designed to guarantee appropriate numerical solutions for both decay and growth cases. Based on this definition, we get restrictions on the step size. Typically, time steps are restricted either when the scheme is explicit and a<0 or when the scheme is implicit and a>0 [5, 8]. Since the amplification factors are functions of complex numbers, it is sometimes difficult to analyze. The following theorem shows that the amplification factor can be approximated by simpler functions of a plus some error terms instead of functions of complex eigenvalues λ. Observe that the order of error terms indicated in Theorem 2.2 is the same as the truncation error of the m-stage Runge-Kutta method. Therefore, we can work with this simpler function without losing the order of accuracy.


#### 2.1.1. Stability Analysis and the Optimization Method

The following theorem states that them-stage Runge-Kutta method satisfies the stability condition automatically for a>0 but needs constraints on Δt when a<0.Theorem 2.2.
|pm(λ)|2=|pm(a+bi)|2=(1+Δta+(1/2)Δt2a2+⋯+(1/m!)Δtmam)2+O(Δtm+1), where  m=1,…,4.Proof.
The proof form=1 is straightforward, and the proof for m=2 is in [5, 8]. For m=3,
(2.10)|p3(λ)|2=(1+Δta+Δt22(a2-b2)+Δt36(a3-3ab2))2+(Δtb+Δt2ab+Δt36(3a2b-b3))2=(1+Δta+12(Δta)2+16(Δta)3)2-2(1+Δta)(12Δt2b2+12Δt3ab2)+O(Δt4)+(Δtb+Δt2ab+Δt36(3a2b-b3))2=(1+Δta+12(Δta)2+16(Δta)3)2-Δt2b2-2Δt3ab2+O(Δt4)+Δt2b2+2Δt3ab2+O(Δt4)=(1+Δta+12(Δta)2+16(Δta)3)2+O(Δt4).
For m=4,
(2.11)|p4(λ)|2=(1+Δta+Δt22(a2-b2)+Δt36(a3-3ab2)+Δt424(a4-6a2b2+b4))2+(Δtb+Δt2ab+Δt36(3a2b-b3)+Δt424(4a3b-4ab3))2=(1+Δta+12(Δta)2+16(Δta)3+124(Δta)4)2-2(1+Δta+12(Δta)2+16(Δta)3+124(Δta)4)*(12Δt2b2+12Δt3ab2+14Δt4a2b2-124Δt4b4)+(12Δt2b2+12Δt3ab2+14Δt4a2b2-124Δt4b4)2+(Δtb+Δt2ab+Δt36(3a2b-b3)+Δt424(4a3b-4ab3))2=(1+Δta+12(Δta)2+16(Δta)3+124(Δta)4)2+Δt4(-2a2b2+112b4)-2Δt3ab2-Δt2b2+O(Δt5)+14Δt4b4+O(Δt5)+Δt4(-13b4+2a2b2)+2Δt3ab2+Δt2b2+O(Δt5)=(1+Δta+12(Δta)2+16(Δta)3+124(Δta)4)2+O(Δt5).
The geometric interpretation is that, for spatial matricesM(tn), the eigenvalues with negative real parts are assumed to be contained in the left-plane ellipse. Let us denote 1+Δta+(1/2)Δt2a2+⋯+(1/m!)Δtmam as pm(a,Δt). To examine whether numerically generated patterns truly obey the original models, one needs to verify if the negative modes (a<0) are diminished and positive modes (a>0) are growing accordingly. That is, when the negative modes occur, pm(a,Δt) is less than one, and is greater than one when the spatial matrix M(tn) happens to have positive modes. From pm(a,Δt), we find that when a>0, pm(a,Δt) is definitely greater than unity, hence the stability condition is automatically satisfied; when a<0, pm(a,Δt) is then possible to be less than unity. It is necessary to put the limitation on the time steps to ensure the stability condition pm(a,Δt)<1 is preserved. Besides, the negative modes should be diminished in an effective way so that the computational cost can be reduced. Under these consideration, we derive the optimal time step which not only diminishes decay modes efficiently but also guarantees the stability conditions. Before moving on the theorems about showing how we get the optimal time steps for m-stage Runge-Kutta method, let us recall the Gerschgorin theorem [19]. It says that for matrices with real coefficients, the distribution for the real parts of eigenvalues can be estimated simply by the discs centered at the diagonal elements. So we can assume that the range of negative real parts of eigenvalues of the discretized spatial matrices is known. Here, only negative real eigenvalues are needed as they are the only ones responsible for the stability and restriction of the step sizes for explicit schemes. We are now ready to give relevant theorems for constructing optimal step sizes by effectively damping the negative modes.Theorem 2.3.
Fora<0 and  m being even,  pm(a,Δt)>0 and is concave up for all  Δt>0.Proof.
Form=2, it is a second-degree polynomial with a positive leading coefficient. The proof is straightforward. Note that p2(a,Δt)>0.
Assumepm-2(a,Δt)>0 for all even m. For m being even and m>2, we have
(2.12)∂2pm(a,Δt)∂Δt2=a2pm-2(a,Δt)>0.
It follows from mathematical induction, ∂2/∂Δt2pm(a,Δt)>0 for all even m. Hence, pm(a,Δt) is concave up w.r.t. Δt if m is even.
Moreover,pm(a,Δt)>0 for all even m, because pm(a,Δt)>exp(aΔt)>0 for even m.Corollary 2.4.
Fora<0 and  m being an odd number,  pm(a,Δt)<1 and  pm(a,Δt) is strictly decreasing for all  Δt>0.Proof.
Form=1, the proof is trivial. For m>1,pm(a,0)=1 and
(2.13)∂pm(a,Δt)∂Δt=apm-1(a,Δt)<0
since a<0 and m-1 is an even number.
Therefore, we concludepm(a,Δt)<1 for all odd m and Δt>0. Moreover, pm(a,Δt) is strictly decreasing with respect to Δt.Theorem 2.5 (Case form being odd).
Assume that-amax≤a≤-amin, where amax>amin>0 are two positive numbers. Then the solution of the following min-max problem occurs at Δt=Δtopt in which Δt satisfies fm(-amax,Δt)=0, where m is an odd number. The min-max problem is stated as
(2.14)min⁡Δt>0max⁡-amax⁡≤a≤-amin⁡pm(a,Δt).Proof.
Note thatpm(a1,Δt)>pm(a2,Δt) if a1>a2 and a1,a2<0. For any fixed Δt, the maximum of pm(a,Δt) is obtained at pm(-amin,Δt). So the min-max occurs when pm(-amax,Δt)=0. Since pm(a,Δt) is continuous, pm(a,0)=1>0, and pm(a,Δt)<0 as Δt→∞, the solution exists by Intermediate Value Theorem. Furthermore, from Corollary 2.4, pm(a,Δt) is monotone decreasing w.r.t. Δt, so the solution is unique.Theorem 2.6 (Case form being even).
Assume that-amax≤a≤-amin, where amax>amin>0 are two positive numbers. Then the solution of the following min-max problem occurs at Δt=Δtopt in which Δt satisfies pm(-amax,Δt)=pm(-amin,Δt), where m is an even number (m=2,4). The min-max problem is stated as
(2.15)min⁡Δt>0max⁡-amax⁡≤a≤-amin⁡pm(a,Δt).Proof.
For each fixeda, pm(a,Δt) is a quadratic-like polynomial function in Δt. Then for any fixed Δt, the maximum of pm(a,Δt) is obtained at pm(-amin,Δt),Δt<Δtopt, and at pm(-amax,Δt), Δt>Δtopt. So the minimax occurs when pm(-amax,Δt)=pm(-amin,Δt), that is,
(2.16)∑i=0m(-amax⁡Δt)mm!=∑i=0m(-amin⁡Δt)mm!.
This leads to the solution.Let us denotegm(Δt)=pm(-amax,Δt)-pm(-amin,Δt). The following theorem states that the optimal time step is unique for m being even. For m being odd, it is shown in Theorem 2.5.Theorem 2.7.
Ifm is even, gm(Δt) has only one positive real solution.Proof.
(2.17)gm(Δt)=Δt(amax⁡-amin⁡)[∑i=1m(-Δt)m-1(amax⁡m-amin⁡m)m!(amax⁡-amin⁡)]=Δt(amax⁡-amin⁡)qm(Δt).
It follows from (2.17) and the fact that amax≠amin, to find the zeros of gm(Δt) is enough to solve qm(Δt)=0. Since qm(0)=1 and qm(Δt)<0 as Δt tends to infinity, positive real solutions hence exist for gm(Δt) by the Intermediate Value Theorem. Moreover, qm(Δt) is strictly decreasing (Figures 1 and 2). Therefore, the solution is unique.Figure 1
Graphs forp3(a,Δt) where a=-6: ·, a=-10: *, a=-12: ∘.Figure 2
Graphs forp4(a,Δt) where a=-6: ·, a=-10: *, a=-12: ∘.


### 2.2. Time-Dependent ODEs Systems

In the previous section, we derived simpler forms of amplification factors as well as theorems for setting optimal time steps for constant ODEs systems. The optimal time step is obtained under the assumption that the intermediate values of Runge-Kutta methods are fixed at current time. Under the stability consideration, we continue our study on searching optimal time steps and relevant schemes when the ODEs system and Runge-Kutta methods are time dependent. The general form of the time dependent system is(2.18)dudt=M(t)u+b(t).Derivation of Adaptive Optimal Time Step for the Three-Stage Runge-Kutta Method
We start with the three-stage Runge-kutta method. As before, we first set the reaction vector to be zero.(2.19)k1=Mnun,k2=Mn+α2(un+β21Δtk1),k3=Mn+α3(un+Δt(β31k1+β32k2)),un+1=un+Δt(γ1k1+γ2k2+γ3k3).
Here,Mn=M(tn). Expand Mn+αi at tn by Taylor expansion, then Mn+αi=Mn+αiΔtnMn′+((αiΔtn)2/2!)Mn′′+O(Δt3), where i=2,3.
Therefore,(2.20)k1=Mnun,k2=[Mn+(α2Mn′+β21Mn2)Δtn+(α222Mn′′+α2β21MnMn′)Δtn2]un+O(Δtn3),k3=[Mn+α3Mn′Δtn+(β31+β32)Mn2Δtn+α322Mn′′Δtn2+(β31α3+β32(α2+α3))Mn′MnΔtn2+β32β21Mn3Δtn2α322]un+O(Δtn3).
Since by consistency,(2.21)∑i=13γi=1,αq=∑j=1q-1βqj,q=2,3,γ2α2+γ3α3=12,γ2α22+γ3α32=13,γ3α2β32=16.
Followed by (2.20) and the consistency condition, (2.19) can be simplified as
(2.22)un+1=un+[MnΔtn+12Mn′Δtn2+16Mn′′Δtn3+12Mn2Δtn2+12Mn′MnΔtn3+16Mn3Δtn3]un+O(Δtn4).
Now, we assume(2.23)∑i=13τi=1,τ2μ+τ3=12,τ2μ2+τ3=13.
Then, we can rewrite (2.22) as
(2.24)un+1=un+[(τ1+τ2+τ3)MnΔtn+(τ2μ+τ3)Mn′Δtn2+12(τ2μ2+τ3)Mn′′Δtn3+12(τ1+τ2+τ3)2Mn2Δtn2+(τ1+τ2+τ3)(τ2μ+τ3)Mn′MnΔtn3+16(τ1+τ2+τ3)3Mn3Δtn3]un+O(Δtn4).
This way, we arrive the following desirable form through rearrangement:(2.25)un+1=un+[τ1Mn+τ2(Mn+μΔtnMn′+(μΔtn)22Mn′′)+τ3(Mn+ΔtnMn′+Δtn22Mn′′)+12[(τ1+τ2+τ3)Mn+(τ2+τ3μ)Mn′Δtn+O(Δtn2)]2Δtn2+16(τ1+τ2+τ3)3Mn3Δtn3Mn+μΔtnMn′+(μΔtn)22Mn′′]un+O(Δtn4).
Let(2.26)Jn=τ1Mn+τ2Mn+μ+τ3Mn+1,
where
(2.27)τ2=-16μ(μ-1),τ3=3μ-26(μ-1).
Then the previous equation implies(2.28)un+1=[1+JnΔtn+12(JnΔtn)2+13!(JnΔtn)3]un+O(Δtn4).
The above result illustrates that the time-dependent three-stage Runge-Kutta method can be treated as a single-step method without losing the order of accuracy. In addition, it shares the same form with (2.5) if Jn is replaced by M. Next, we derive a similar argument for the four-stage Runge-Kutta method.Derivation of Adaptive Optimal Time Step for the Four-Stage Runge-Kutta Method
Applying the four-stage Runge-Kutta method to (2.18), we get
(2.29)k1=Mnun,k2=Mn+α2(un+β21Δtk1),k3=Mn+α3(un+Δt(β31k1+β32k2)),k4=Mn+α4(un+Δt(β41k1+β42k2+β43k3)),un+1=un+Δt(r1k1+r2k2+r3k3+r4k4).
SinceMn+αi=Mn+αiΔtnMn′+((αiΔtn)2/2!)Mn′′+((αiΔtn)3/3!)Mn′′′+O(Δt4),
(2.30)k2=[Mn+(α2Mn′+β21Mn2)Δtn+(α222Mn′′+α2β21MnMn′)Δtn2+(α233!Mn′′′+α22β212Mn′′Mn)Δtn3]un+O(Δtn4),k3=[Mn+α3Mn′Δtn+(β31+β32)Mn2Δtn+α322Mn′′Δtn2+(β31α3+β32(α2+α3))Mn′MnΔtn2+β32β21Mn3Δtn2+α333!Mn′′′Δtn3+12(β31α32+β32(α22+α32))Mn′′MnΔtn3+β32α2α3(Mn′)2Δtn3+β32β21(α2+α3)Mn′Mn2Δtn3α322]un+O(Δtn4),k4=[Mn+α4Mn′Δtn+(β41+β42+β43)Mn2Δtn+12α42Mn′′Δtn2+(β41α4+β42(α2+α4)+β43(α3+α4))Mn′MnΔtn2+(β42β21+β43(β31+β32))Mn3Δtn2+13!α43Mn′′′Δtn3+12(β41α42+β42(α22+α42)+β43(α32+α42))Mn′′MnΔtn3+(β42β21(α2+α4)+β43β31(α3+α4)+β43β32(α2+α3+α4))Mn′Mn2Δtn3+(β42α2α4+β43α3α4)(Mn′)2Δtn3+β43β32β21Mn4Δtn312]un+O(Δtn4).
By consistency of the method,(2.31)∑i=14γi=1,αq=∑j=1q-1βqj,q=2,…,4,γ2α2+γ3α3+γ4α4=12,γ2α22+γ3α32+γ4α42=13,γ2α23+γ3α33+γ4α43=14,γ3α2β32+γ4(α2β42+β43α3)=16,γ3β32α2α3+γ4β42α2α4+γ4β43α3α4=18,γ3α22β32+γ4α22β42+γ4α32β43=112,γ4α2β43β32=124.
Thus,(2.32)un+1=un+[MnΔtn+12Mn′Δtn2+16Mn′′Δtn3+124Mn′′′Δtn4+12Mn2Δtn2+12Mn′MnΔtn3+16Mn′′MnΔtn4+18(Mn′)2Δtn4+16Mn3Δtn3+14Mn′Mn2Δtn4+124Mn4Δtn4]un+O(Δtn5).
We further assume(2.33)∑i=14τi=1,τ2μ1+τ3μ2+τ4=12,τ2μ12+τ3μ22+τ4=13,τ2μ13+τ3μ23+τ4=14.
The previous equation hence can be written as(2.34)un+1=un+[124Mn4Δtn4(τ1+τ2+τ3+τ4)MnΔtn+(τ2μ1+τ3μ2+τ4)Mn′Δtn2+12(τ2μ12+τ3μ22+τ4)Mn′′Δtn3+16(τ2μ13+τ3μ23+τ4)Mn′′′Δtn4+12(τ1+τ2+τ3+τ4)2Mn2Δtn2+(τ1+τ2+τ3+τ4)*(τ2μ1+τ3μ2+τ4)Mn′MnΔtn3+12(τ1+τ2+τ3+τ4)(τ2μ12+τ3μ22+τ4)Mn′′MnΔtn4+12(τ2μ1+τ3μ2+τ4)2(Mn′)2Δtn4+16(τ1+τ2+τ3+τ4)3Mn3Δtn3+12(τ1+τ2+τ3+τ4)(τ2μ1+τ3μ2+τ4)Mn′Mn2Δtn4+124Mn4Δtn4]un+O(Δtn5).
That is,(2.35)un+1=[1+JnΔtn+12(JnΔtn)2+13!(JnΔtn)3+14!(JnΔtn)4]un+O(Δtn5),
where
(2.36)Jn=τ1Mn+τ2Mn+μ1+τ3Mn+μ2+τ4Mn+1,τ2=2μ2-112(-μ2+μ1)(μ1-1)μ1,τ3=--1+2μ112μ2(μ1μ2+μ2-μ22-μ1),τ4=6μ1μ2-4μ2+3-4μ112(μ2-1)(μ1-1).
We now turn to examine the more general situation in which the reaction term is not zero, that is,b≠0. In this case, we have(2.37)un+1=(∑i=0m(ΔtnJn)ii!)un+F(bn,bn+μ1,…,bn+1)+O(Δtni+1).
Similar to our discussion in the previous section, the accumulation error for (2.37) is
(2.38)|un+1-vn+1|=|∑i=0m(ΔtnJn)ii!|×|un-vn|+O(Δtni+1).
We remark again that the amplification factor∑i=0m((ΔtnJn)i/i!) in (2.38) at the nth local step has the same form with (2.5) except M being replaced by Jn. Therefore, we can apply the strategy demonstrated in Theorems 2.5 and 2.6 to look for the optimal time steps for Runge-Kutta schemes. The present amax and amin appeared in the optimal step size formulas then, respectively, represent the greatest and smallest real parts of the eigenvalues of the matrix Jn. Notice that by the similar argument outlined in the previous proofs, the stability is preserved with this step size.An Explicit Adaptive Optimal Time Stepping Scheme for Solving (2.18)
(1)
Set up initial conditionu0 and initial step size Δt0.(2)
SetJn in accordance with the stage number m, and apply m-stage Runge-Kutta schemes.(3)
Use the Gerschgorin theorem to estimate the interval,[-amaxn,-aminn] for negative real parts of the eigenvalues of Jn.(4)
PlaceΔtn+1,
(i)
m=2:Δtn+1=2/(amaxn+aminn),(ii)
m=3:Δtn+1=2.512745327/amaxn,(iii)
m=4:(2.39)Δtn+1=13((+30amax⁡amin⁡5+29amin⁡j⁡6)172amax⁡6+192amax⁡5amin⁡+276amax⁡4amin⁡2+448amax⁡3amin⁡3+276amin⁡4amax⁡2+192amax⁡amin⁡5+172amin⁡6+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amax⁡3+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amin⁡amax⁡2+36(29amax⁡6+30amax⁡5+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amax⁡amin⁡2+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amin⁡3(29agg))1/3/(amax⁡3+amin⁡amax⁡2+amax⁡amin⁡2+amin⁡3)-43(5amax⁡4+10amax⁡3amin⁡+6amax⁡2amin⁡2+10amax⁡amin⁡3+5amin⁡4)/(∑222(amax⁡3+amin⁡amax⁡2+amax⁡amin⁡2+amin⁡3)×((+30amax⁡amin⁡5+29amin⁡j⁡6)172amax⁡6+192amax⁡5amin⁡+276amax⁡4amin⁡2+448amax⁡3amin⁡3+276amin⁡4amax⁡2+192amax⁡amin⁡5+172amin⁡6+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2ar3+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amin⁡amax⁡2+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amax⁡amin⁡2+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amin⁡3(29agg))1/3(amin⁡3)3)+43(amax⁡2+amax⁡amin⁡+amin⁡2)/(amax⁡3+amin⁡amax⁡2+amax⁡amin⁡2+amin⁡3),
and then let tn+1=tn+Δtn+1. When an amaxn is equal to or close to zero, simply choose a reasonable step size for Δtn+1.
Note that the time step form=4 can be solved by mathematical software such as MAPLE from the following equation:
(2.40)∑i=141i!(amax⁡i-amin⁡i)Δtn+1i-1=0.


#### 2.2.1. The “Best” Optimal Step Size and the Time-Dependent Homogeneous Systems

We point out that for those eigenvalues with negative real parts, one may also assume they are enclosed in an ellipse with one of the vertex lying at the origin, that is,aminn=0. In this way, only amaxn are needed in forming optimal time step formulas and time steps Δtn+1 in step 4 are switched as(i)
m=2:Δtn+1=2/amaxn,(ii)
m=3:Δtn+1=2.512745327/amaxn,(iii)
m=4:Δtn+1=2.785293563/amaxn,whereamaxn is defined as before in step 2. We remark that when Δtn+1 is specified as above values, pm(a,Δt) is equal to one. The stability condition does not seem to be satisfied. However, this hidden problem is avoidable by requiring Gerschgorin theorem or other equivalent theorems in step 3 as an estimation for the distribution of real part of eigenvalues. This is because the theorem provides an overestimation for the distribution of negative real part of eigenvalues. As a result, pm(a,Δt) is less than one. Hence, it is safe to use the above time step.In particular, when the reaction vector is negligible, we can further write them-stage Runge-kutta method as some equivalent m-stage formulas while the order of accuracy remains invariant. For example, m=3, we have(2.41)un(0)=un+Δtn3Jnun,un(1)=un+Δtn2Jnun(0),un+1=un+ΔtJnun(1).Form=4,
(2.42)un(0)=un+Δtn4Jnun,un(1)=un+Δtn3Jnun(0),un(2)=un+Δtn2Jnun(1),un+1=un+ΔtnJnun(2).Of course,Jn is defined according to their stage number.


## 2.1. Time Independence ODEs Systems

Since interesting patterns are caused by the interaction between growth and decay modes, to analyze system (1.3) and relate it to the pattern formation model (1.1), one can do the local linearization property of the right-hand side of (1.3). The positive and negative real eigenvalues of the linearization matrix then, respectively, correspond to the growth and decay modes. Because of this reason, we will consider the corresponding ODEs induced from the chemotaxis system and introduce our idea by starting with the constant coefficients case, that is,(2.1)dudt=Mu+b
then generalize it to the time-dependent coefficients case. The general formulation of the m-stage Runge-Kutta method is(2.2)k1=f(tn,un),kq=f(tn+αqΔt,un+Δtn∑s=1q-1βqsks),q=2,…,m,αq=∑s=1q-1βqs,un+1=un+Δtn∑s=1qγsks,q=2,…,m,
where m is the stage number and αq=∑j=1q-1βqj,q=2,…,m. The subscripts n and n+1 are the iteration level. For simplicity, we first look at the uniform time step size Δt=tn-tn-1 for all n. Then un≈u(xn) is the numerical solution at the nth step. Equations for the unknowns γs and βqs can be obtained by matching terms with Taylor's expansion.We choose single-step explicit methods, which basically are Runge-Kutta methods, because they are easy to implement and especially efficient for ODEs systems derived from the PDEs. The method is set to be explicit since the range of eigenvalues with negative real parts is measurable for chemotaxis systems. In [5, 8], we have developed adaptive two-stage time stepping scheme with optimum time step sizes for solving chemotaxis system and verified that it is reliable and more economical when comparing with other similar methods. We conjecture that the similar strategy used in [5, 8] could be used to derive optimal adaptive three- and four-stage methods. In the two-stage case, the size of the optimal step is determined by the root of the first-order equation hence it is absolute optimal within one time step. If our conjecture is correct, the sizes of optimal steps for higher-stage methods should be the roots of some higher-order equations. But, are these adaptive step sizes absolute optimal? And are they easy to be found as in the two-stage case?To answer these questions, we begin our study with the three-stage Runge-Kutta method and setb=0 in (2.1). As in [5, 8], we first fix all the intermediate coefficients at t=tn for the three-stage Runge-Kutta method. The general time dependence case is proposed in the next section. Denoting M=M(tn) and applying (2.2) to the three-stage Runge-Kutta method with all the intermediate values being fixed at tn, we arrive at(2.3)k1=Mun,k2=M(1+α2ΔtM)un,k3=M(1+α3ΔtM+α2β32Δt2M2)un,(2.4)un+1=(1+δ1ΔtM+δ2Δt2M2+δ3Δt3M3)un,
where δi,i=1,2,3, are combinations of the coefficients βqs and γs. Since u(tn+1)=eMΔtu(tn) and assuming un is exact up to the order Δt3, we have δ1=1, δ2=(1/2!), δ3=(1/3!). Similar arguments also hold for RK methods with other stages. Thus, if we treat it as a standard iterative method, (2.4) is a one-step method [18]:(2.5)un+1=[1+ΔtM+12!(ΔtM)2+⋯+1m!(ΔtM)m]un.=pm(Δt)un.In the case ofb≠0, (2.5) becomes(2.6)un+1=pm(Δt)un+(1+Δt2!M+⋯+Δtmm!Mm-1)b.Denotev as the solution by solving the system with the perturbated initial condition; we have(2.7)vn+1=pm(Δt)vn+(1+Δt2!M+⋯+Δtmm!Mm-1)b.From (2.6) and (2.7), the accumulation error can be written as (2.8)|un+1-vn+1|=|pm(Δt)|×|un-vn|.Note that(2.9)pm(λ)=1+Δtλ+12!(Δtλ)2+⋯+1m!(Δtλ)m
is called the amplification factor, where λ is an eigenvalue of M. It is obtained when a scalar test problem, u′=λu, is applied to the single-step method.Suppose thatλ=a+bi=Reλ+Imλi. For pattern formation problems, what we are interested in is the long-term behavior of the solution vectors. It is important to ensure that observed patterns are accurate in the sense that they are not coming from the accumulation error of numerical schemes. We therefore make the following definition of stability.Definition 2.1.
An explicit single-step method is said to be stable if|pm(λ)|<1 when a<0 and |pm(λ)|≤1+aΔt+O(Δt) when a≥0.HereO(Δt) means higher order of Δt. Our definition is designed to guarantee appropriate numerical solutions for both decay and growth cases. Based on this definition, we get restrictions on the step size. Typically, time steps are restricted either when the scheme is explicit and a<0 or when the scheme is implicit and a>0 [5, 8]. Since the amplification factors are functions of complex numbers, it is sometimes difficult to analyze. The following theorem shows that the amplification factor can be approximated by simpler functions of a plus some error terms instead of functions of complex eigenvalues λ. Observe that the order of error terms indicated in Theorem 2.2 is the same as the truncation error of the m-stage Runge-Kutta method. Therefore, we can work with this simpler function without losing the order of accuracy.


### 2.1.1. Stability Analysis and the Optimization Method

The following theorem states that them-stage Runge-Kutta method satisfies the stability condition automatically for a>0 but needs constraints on Δt when a<0.Theorem 2.2.
|pm(λ)|2=|pm(a+bi)|2=(1+Δta+(1/2)Δt2a2+⋯+(1/m!)Δtmam)2+O(Δtm+1), where  m=1,…,4.Proof.
The proof form=1 is straightforward, and the proof for m=2 is in [5, 8]. For m=3,
(2.10)|p3(λ)|2=(1+Δta+Δt22(a2-b2)+Δt36(a3-3ab2))2+(Δtb+Δt2ab+Δt36(3a2b-b3))2=(1+Δta+12(Δta)2+16(Δta)3)2-2(1+Δta)(12Δt2b2+12Δt3ab2)+O(Δt4)+(Δtb+Δt2ab+Δt36(3a2b-b3))2=(1+Δta+12(Δta)2+16(Δta)3)2-Δt2b2-2Δt3ab2+O(Δt4)+Δt2b2+2Δt3ab2+O(Δt4)=(1+Δta+12(Δta)2+16(Δta)3)2+O(Δt4).
For m=4,
(2.11)|p4(λ)|2=(1+Δta+Δt22(a2-b2)+Δt36(a3-3ab2)+Δt424(a4-6a2b2+b4))2+(Δtb+Δt2ab+Δt36(3a2b-b3)+Δt424(4a3b-4ab3))2=(1+Δta+12(Δta)2+16(Δta)3+124(Δta)4)2-2(1+Δta+12(Δta)2+16(Δta)3+124(Δta)4)*(12Δt2b2+12Δt3ab2+14Δt4a2b2-124Δt4b4)+(12Δt2b2+12Δt3ab2+14Δt4a2b2-124Δt4b4)2+(Δtb+Δt2ab+Δt36(3a2b-b3)+Δt424(4a3b-4ab3))2=(1+Δta+12(Δta)2+16(Δta)3+124(Δta)4)2+Δt4(-2a2b2+112b4)-2Δt3ab2-Δt2b2+O(Δt5)+14Δt4b4+O(Δt5)+Δt4(-13b4+2a2b2)+2Δt3ab2+Δt2b2+O(Δt5)=(1+Δta+12(Δta)2+16(Δta)3+124(Δta)4)2+O(Δt5).
The geometric interpretation is that, for spatial matricesM(tn), the eigenvalues with negative real parts are assumed to be contained in the left-plane ellipse. Let us denote 1+Δta+(1/2)Δt2a2+⋯+(1/m!)Δtmam as pm(a,Δt). To examine whether numerically generated patterns truly obey the original models, one needs to verify if the negative modes (a<0) are diminished and positive modes (a>0) are growing accordingly. That is, when the negative modes occur, pm(a,Δt) is less than one, and is greater than one when the spatial matrix M(tn) happens to have positive modes. From pm(a,Δt), we find that when a>0, pm(a,Δt) is definitely greater than unity, hence the stability condition is automatically satisfied; when a<0, pm(a,Δt) is then possible to be less than unity. It is necessary to put the limitation on the time steps to ensure the stability condition pm(a,Δt)<1 is preserved. Besides, the negative modes should be diminished in an effective way so that the computational cost can be reduced. Under these consideration, we derive the optimal time step which not only diminishes decay modes efficiently but also guarantees the stability conditions. Before moving on the theorems about showing how we get the optimal time steps for m-stage Runge-Kutta method, let us recall the Gerschgorin theorem [19]. It says that for matrices with real coefficients, the distribution for the real parts of eigenvalues can be estimated simply by the discs centered at the diagonal elements. So we can assume that the range of negative real parts of eigenvalues of the discretized spatial matrices is known. Here, only negative real eigenvalues are needed as they are the only ones responsible for the stability and restriction of the step sizes for explicit schemes. We are now ready to give relevant theorems for constructing optimal step sizes by effectively damping the negative modes.Theorem 2.3.
Fora<0 and  m being even,  pm(a,Δt)>0 and is concave up for all  Δt>0.Proof.
Form=2, it is a second-degree polynomial with a positive leading coefficient. The proof is straightforward. Note that p2(a,Δt)>0.
Assumepm-2(a,Δt)>0 for all even m. For m being even and m>2, we have
(2.12)∂2pm(a,Δt)∂Δt2=a2pm-2(a,Δt)>0.
It follows from mathematical induction, ∂2/∂Δt2pm(a,Δt)>0 for all even m. Hence, pm(a,Δt) is concave up w.r.t. Δt if m is even.
Moreover,pm(a,Δt)>0 for all even m, because pm(a,Δt)>exp(aΔt)>0 for even m.Corollary 2.4.
Fora<0 and  m being an odd number,  pm(a,Δt)<1 and  pm(a,Δt) is strictly decreasing for all  Δt>0.Proof.
Form=1, the proof is trivial. For m>1,pm(a,0)=1 and
(2.13)∂pm(a,Δt)∂Δt=apm-1(a,Δt)<0
since a<0 and m-1 is an even number.
Therefore, we concludepm(a,Δt)<1 for all odd m and Δt>0. Moreover, pm(a,Δt) is strictly decreasing with respect to Δt.Theorem 2.5 (Case form being odd).
Assume that-amax≤a≤-amin, where amax>amin>0 are two positive numbers. Then the solution of the following min-max problem occurs at Δt=Δtopt in which Δt satisfies fm(-amax,Δt)=0, where m is an odd number. The min-max problem is stated as
(2.14)min⁡Δt>0max⁡-amax⁡≤a≤-amin⁡pm(a,Δt).Proof.
Note thatpm(a1,Δt)>pm(a2,Δt) if a1>a2 and a1,a2<0. For any fixed Δt, the maximum of pm(a,Δt) is obtained at pm(-amin,Δt). So the min-max occurs when pm(-amax,Δt)=0. Since pm(a,Δt) is continuous, pm(a,0)=1>0, and pm(a,Δt)<0 as Δt→∞, the solution exists by Intermediate Value Theorem. Furthermore, from Corollary 2.4, pm(a,Δt) is monotone decreasing w.r.t. Δt, so the solution is unique.Theorem 2.6 (Case form being even).
Assume that-amax≤a≤-amin, where amax>amin>0 are two positive numbers. Then the solution of the following min-max problem occurs at Δt=Δtopt in which Δt satisfies pm(-amax,Δt)=pm(-amin,Δt), where m is an even number (m=2,4). The min-max problem is stated as
(2.15)min⁡Δt>0max⁡-amax⁡≤a≤-amin⁡pm(a,Δt).Proof.
For each fixeda, pm(a,Δt) is a quadratic-like polynomial function in Δt. Then for any fixed Δt, the maximum of pm(a,Δt) is obtained at pm(-amin,Δt),Δt<Δtopt, and at pm(-amax,Δt), Δt>Δtopt. So the minimax occurs when pm(-amax,Δt)=pm(-amin,Δt), that is,
(2.16)∑i=0m(-amax⁡Δt)mm!=∑i=0m(-amin⁡Δt)mm!.
This leads to the solution.Let us denotegm(Δt)=pm(-amax,Δt)-pm(-amin,Δt). The following theorem states that the optimal time step is unique for m being even. For m being odd, it is shown in Theorem 2.5.Theorem 2.7.
Ifm is even, gm(Δt) has only one positive real solution.Proof.
(2.17)gm(Δt)=Δt(amax⁡-amin⁡)[∑i=1m(-Δt)m-1(amax⁡m-amin⁡m)m!(amax⁡-amin⁡)]=Δt(amax⁡-amin⁡)qm(Δt).
It follows from (2.17) and the fact that amax≠amin, to find the zeros of gm(Δt) is enough to solve qm(Δt)=0. Since qm(0)=1 and qm(Δt)<0 as Δt tends to infinity, positive real solutions hence exist for gm(Δt) by the Intermediate Value Theorem. Moreover, qm(Δt) is strictly decreasing (Figures 1 and 2). Therefore, the solution is unique.Figure 1
Graphs forp3(a,Δt) where a=-6: ·, a=-10: *, a=-12: ∘.Figure 2
Graphs forp4(a,Δt) where a=-6: ·, a=-10: *, a=-12: ∘.


## 2.1.1. Stability Analysis and the Optimization Method

The following theorem states that them-stage Runge-Kutta method satisfies the stability condition automatically for a>0 but needs constraints on Δt when a<0.Theorem 2.2.
|pm(λ)|2=|pm(a+bi)|2=(1+Δta+(1/2)Δt2a2+⋯+(1/m!)Δtmam)2+O(Δtm+1), where  m=1,…,4.Proof.
The proof form=1 is straightforward, and the proof for m=2 is in [5, 8]. For m=3,
(2.10)|p3(λ)|2=(1+Δta+Δt22(a2-b2)+Δt36(a3-3ab2))2+(Δtb+Δt2ab+Δt36(3a2b-b3))2=(1+Δta+12(Δta)2+16(Δta)3)2-2(1+Δta)(12Δt2b2+12Δt3ab2)+O(Δt4)+(Δtb+Δt2ab+Δt36(3a2b-b3))2=(1+Δta+12(Δta)2+16(Δta)3)2-Δt2b2-2Δt3ab2+O(Δt4)+Δt2b2+2Δt3ab2+O(Δt4)=(1+Δta+12(Δta)2+16(Δta)3)2+O(Δt4).
For m=4,
(2.11)|p4(λ)|2=(1+Δta+Δt22(a2-b2)+Δt36(a3-3ab2)+Δt424(a4-6a2b2+b4))2+(Δtb+Δt2ab+Δt36(3a2b-b3)+Δt424(4a3b-4ab3))2=(1+Δta+12(Δta)2+16(Δta)3+124(Δta)4)2-2(1+Δta+12(Δta)2+16(Δta)3+124(Δta)4)*(12Δt2b2+12Δt3ab2+14Δt4a2b2-124Δt4b4)+(12Δt2b2+12Δt3ab2+14Δt4a2b2-124Δt4b4)2+(Δtb+Δt2ab+Δt36(3a2b-b3)+Δt424(4a3b-4ab3))2=(1+Δta+12(Δta)2+16(Δta)3+124(Δta)4)2+Δt4(-2a2b2+112b4)-2Δt3ab2-Δt2b2+O(Δt5)+14Δt4b4+O(Δt5)+Δt4(-13b4+2a2b2)+2Δt3ab2+Δt2b2+O(Δt5)=(1+Δta+12(Δta)2+16(Δta)3+124(Δta)4)2+O(Δt5).
The geometric interpretation is that, for spatial matricesM(tn), the eigenvalues with negative real parts are assumed to be contained in the left-plane ellipse. Let us denote 1+Δta+(1/2)Δt2a2+⋯+(1/m!)Δtmam as pm(a,Δt). To examine whether numerically generated patterns truly obey the original models, one needs to verify if the negative modes (a<0) are diminished and positive modes (a>0) are growing accordingly. That is, when the negative modes occur, pm(a,Δt) is less than one, and is greater than one when the spatial matrix M(tn) happens to have positive modes. From pm(a,Δt), we find that when a>0, pm(a,Δt) is definitely greater than unity, hence the stability condition is automatically satisfied; when a<0, pm(a,Δt) is then possible to be less than unity. It is necessary to put the limitation on the time steps to ensure the stability condition pm(a,Δt)<1 is preserved. Besides, the negative modes should be diminished in an effective way so that the computational cost can be reduced. Under these consideration, we derive the optimal time step which not only diminishes decay modes efficiently but also guarantees the stability conditions. Before moving on the theorems about showing how we get the optimal time steps for m-stage Runge-Kutta method, let us recall the Gerschgorin theorem [19]. It says that for matrices with real coefficients, the distribution for the real parts of eigenvalues can be estimated simply by the discs centered at the diagonal elements. So we can assume that the range of negative real parts of eigenvalues of the discretized spatial matrices is known. Here, only negative real eigenvalues are needed as they are the only ones responsible for the stability and restriction of the step sizes for explicit schemes. We are now ready to give relevant theorems for constructing optimal step sizes by effectively damping the negative modes.Theorem 2.3.
Fora<0 and  m being even,  pm(a,Δt)>0 and is concave up for all  Δt>0.Proof.
Form=2, it is a second-degree polynomial with a positive leading coefficient. The proof is straightforward. Note that p2(a,Δt)>0.
Assumepm-2(a,Δt)>0 for all even m. For m being even and m>2, we have
(2.12)∂2pm(a,Δt)∂Δt2=a2pm-2(a,Δt)>0.
It follows from mathematical induction, ∂2/∂Δt2pm(a,Δt)>0 for all even m. Hence, pm(a,Δt) is concave up w.r.t. Δt if m is even.
Moreover,pm(a,Δt)>0 for all even m, because pm(a,Δt)>exp(aΔt)>0 for even m.Corollary 2.4.
Fora<0 and  m being an odd number,  pm(a,Δt)<1 and  pm(a,Δt) is strictly decreasing for all  Δt>0.Proof.
Form=1, the proof is trivial. For m>1,pm(a,0)=1 and
(2.13)∂pm(a,Δt)∂Δt=apm-1(a,Δt)<0
since a<0 and m-1 is an even number.
Therefore, we concludepm(a,Δt)<1 for all odd m and Δt>0. Moreover, pm(a,Δt) is strictly decreasing with respect to Δt.Theorem 2.5 (Case form being odd).
Assume that-amax≤a≤-amin, where amax>amin>0 are two positive numbers. Then the solution of the following min-max problem occurs at Δt=Δtopt in which Δt satisfies fm(-amax,Δt)=0, where m is an odd number. The min-max problem is stated as
(2.14)min⁡Δt>0max⁡-amax⁡≤a≤-amin⁡pm(a,Δt).Proof.
Note thatpm(a1,Δt)>pm(a2,Δt) if a1>a2 and a1,a2<0. For any fixed Δt, the maximum of pm(a,Δt) is obtained at pm(-amin,Δt). So the min-max occurs when pm(-amax,Δt)=0. Since pm(a,Δt) is continuous, pm(a,0)=1>0, and pm(a,Δt)<0 as Δt→∞, the solution exists by Intermediate Value Theorem. Furthermore, from Corollary 2.4, pm(a,Δt) is monotone decreasing w.r.t. Δt, so the solution is unique.Theorem 2.6 (Case form being even).
Assume that-amax≤a≤-amin, where amax>amin>0 are two positive numbers. Then the solution of the following min-max problem occurs at Δt=Δtopt in which Δt satisfies pm(-amax,Δt)=pm(-amin,Δt), where m is an even number (m=2,4). The min-max problem is stated as
(2.15)min⁡Δt>0max⁡-amax⁡≤a≤-amin⁡pm(a,Δt).Proof.
For each fixeda, pm(a,Δt) is a quadratic-like polynomial function in Δt. Then for any fixed Δt, the maximum of pm(a,Δt) is obtained at pm(-amin,Δt),Δt<Δtopt, and at pm(-amax,Δt), Δt>Δtopt. So the minimax occurs when pm(-amax,Δt)=pm(-amin,Δt), that is,
(2.16)∑i=0m(-amax⁡Δt)mm!=∑i=0m(-amin⁡Δt)mm!.
This leads to the solution.Let us denotegm(Δt)=pm(-amax,Δt)-pm(-amin,Δt). The following theorem states that the optimal time step is unique for m being even. For m being odd, it is shown in Theorem 2.5.Theorem 2.7.
Ifm is even, gm(Δt) has only one positive real solution.Proof.
(2.17)gm(Δt)=Δt(amax⁡-amin⁡)[∑i=1m(-Δt)m-1(amax⁡m-amin⁡m)m!(amax⁡-amin⁡)]=Δt(amax⁡-amin⁡)qm(Δt).
It follows from (2.17) and the fact that amax≠amin, to find the zeros of gm(Δt) is enough to solve qm(Δt)=0. Since qm(0)=1 and qm(Δt)<0 as Δt tends to infinity, positive real solutions hence exist for gm(Δt) by the Intermediate Value Theorem. Moreover, qm(Δt) is strictly decreasing (Figures 1 and 2). Therefore, the solution is unique.Figure 1
Graphs forp3(a,Δt) where a=-6: ·, a=-10: *, a=-12: ∘.Figure 2
Graphs forp4(a,Δt) where a=-6: ·, a=-10: *, a=-12: ∘.


## 2.2. Time-Dependent ODEs Systems

In the previous section, we derived simpler forms of amplification factors as well as theorems for setting optimal time steps for constant ODEs systems. The optimal time step is obtained under the assumption that the intermediate values of Runge-Kutta methods are fixed at current time. Under the stability consideration, we continue our study on searching optimal time steps and relevant schemes when the ODEs system and Runge-Kutta methods are time dependent. The general form of the time dependent system is(2.18)dudt=M(t)u+b(t).Derivation of Adaptive Optimal Time Step for the Three-Stage Runge-Kutta Method
We start with the three-stage Runge-kutta method. As before, we first set the reaction vector to be zero.(2.19)k1=Mnun,k2=Mn+α2(un+β21Δtk1),k3=Mn+α3(un+Δt(β31k1+β32k2)),un+1=un+Δt(γ1k1+γ2k2+γ3k3).
Here,Mn=M(tn). Expand Mn+αi at tn by Taylor expansion, then Mn+αi=Mn+αiΔtnMn′+((αiΔtn)2/2!)Mn′′+O(Δt3), where i=2,3.
Therefore,(2.20)k1=Mnun,k2=[Mn+(α2Mn′+β21Mn2)Δtn+(α222Mn′′+α2β21MnMn′)Δtn2]un+O(Δtn3),k3=[Mn+α3Mn′Δtn+(β31+β32)Mn2Δtn+α322Mn′′Δtn2+(β31α3+β32(α2+α3))Mn′MnΔtn2+β32β21Mn3Δtn2α322]un+O(Δtn3).
Since by consistency,(2.21)∑i=13γi=1,αq=∑j=1q-1βqj,q=2,3,γ2α2+γ3α3=12,γ2α22+γ3α32=13,γ3α2β32=16.
Followed by (2.20) and the consistency condition, (2.19) can be simplified as
(2.22)un+1=un+[MnΔtn+12Mn′Δtn2+16Mn′′Δtn3+12Mn2Δtn2+12Mn′MnΔtn3+16Mn3Δtn3]un+O(Δtn4).
Now, we assume(2.23)∑i=13τi=1,τ2μ+τ3=12,τ2μ2+τ3=13.
Then, we can rewrite (2.22) as
(2.24)un+1=un+[(τ1+τ2+τ3)MnΔtn+(τ2μ+τ3)Mn′Δtn2+12(τ2μ2+τ3)Mn′′Δtn3+12(τ1+τ2+τ3)2Mn2Δtn2+(τ1+τ2+τ3)(τ2μ+τ3)Mn′MnΔtn3+16(τ1+τ2+τ3)3Mn3Δtn3]un+O(Δtn4).
This way, we arrive the following desirable form through rearrangement:(2.25)un+1=un+[τ1Mn+τ2(Mn+μΔtnMn′+(μΔtn)22Mn′′)+τ3(Mn+ΔtnMn′+Δtn22Mn′′)+12[(τ1+τ2+τ3)Mn+(τ2+τ3μ)Mn′Δtn+O(Δtn2)]2Δtn2+16(τ1+τ2+τ3)3Mn3Δtn3Mn+μΔtnMn′+(μΔtn)22Mn′′]un+O(Δtn4).
Let(2.26)Jn=τ1Mn+τ2Mn+μ+τ3Mn+1,
where
(2.27)τ2=-16μ(μ-1),τ3=3μ-26(μ-1).
Then the previous equation implies(2.28)un+1=[1+JnΔtn+12(JnΔtn)2+13!(JnΔtn)3]un+O(Δtn4).
The above result illustrates that the time-dependent three-stage Runge-Kutta method can be treated as a single-step method without losing the order of accuracy. In addition, it shares the same form with (2.5) if Jn is replaced by M. Next, we derive a similar argument for the four-stage Runge-Kutta method.Derivation of Adaptive Optimal Time Step for the Four-Stage Runge-Kutta Method
Applying the four-stage Runge-Kutta method to (2.18), we get
(2.29)k1=Mnun,k2=Mn+α2(un+β21Δtk1),k3=Mn+α3(un+Δt(β31k1+β32k2)),k4=Mn+α4(un+Δt(β41k1+β42k2+β43k3)),un+1=un+Δt(r1k1+r2k2+r3k3+r4k4).
SinceMn+αi=Mn+αiΔtnMn′+((αiΔtn)2/2!)Mn′′+((αiΔtn)3/3!)Mn′′′+O(Δt4),
(2.30)k2=[Mn+(α2Mn′+β21Mn2)Δtn+(α222Mn′′+α2β21MnMn′)Δtn2+(α233!Mn′′′+α22β212Mn′′Mn)Δtn3]un+O(Δtn4),k3=[Mn+α3Mn′Δtn+(β31+β32)Mn2Δtn+α322Mn′′Δtn2+(β31α3+β32(α2+α3))Mn′MnΔtn2+β32β21Mn3Δtn2+α333!Mn′′′Δtn3+12(β31α32+β32(α22+α32))Mn′′MnΔtn3+β32α2α3(Mn′)2Δtn3+β32β21(α2+α3)Mn′Mn2Δtn3α322]un+O(Δtn4),k4=[Mn+α4Mn′Δtn+(β41+β42+β43)Mn2Δtn+12α42Mn′′Δtn2+(β41α4+β42(α2+α4)+β43(α3+α4))Mn′MnΔtn2+(β42β21+β43(β31+β32))Mn3Δtn2+13!α43Mn′′′Δtn3+12(β41α42+β42(α22+α42)+β43(α32+α42))Mn′′MnΔtn3+(β42β21(α2+α4)+β43β31(α3+α4)+β43β32(α2+α3+α4))Mn′Mn2Δtn3+(β42α2α4+β43α3α4)(Mn′)2Δtn3+β43β32β21Mn4Δtn312]un+O(Δtn4).
By consistency of the method,(2.31)∑i=14γi=1,αq=∑j=1q-1βqj,q=2,…,4,γ2α2+γ3α3+γ4α4=12,γ2α22+γ3α32+γ4α42=13,γ2α23+γ3α33+γ4α43=14,γ3α2β32+γ4(α2β42+β43α3)=16,γ3β32α2α3+γ4β42α2α4+γ4β43α3α4=18,γ3α22β32+γ4α22β42+γ4α32β43=112,γ4α2β43β32=124.
Thus,(2.32)un+1=un+[MnΔtn+12Mn′Δtn2+16Mn′′Δtn3+124Mn′′′Δtn4+12Mn2Δtn2+12Mn′MnΔtn3+16Mn′′MnΔtn4+18(Mn′)2Δtn4+16Mn3Δtn3+14Mn′Mn2Δtn4+124Mn4Δtn4]un+O(Δtn5).
We further assume(2.33)∑i=14τi=1,τ2μ1+τ3μ2+τ4=12,τ2μ12+τ3μ22+τ4=13,τ2μ13+τ3μ23+τ4=14.
The previous equation hence can be written as(2.34)un+1=un+[124Mn4Δtn4(τ1+τ2+τ3+τ4)MnΔtn+(τ2μ1+τ3μ2+τ4)Mn′Δtn2+12(τ2μ12+τ3μ22+τ4)Mn′′Δtn3+16(τ2μ13+τ3μ23+τ4)Mn′′′Δtn4+12(τ1+τ2+τ3+τ4)2Mn2Δtn2+(τ1+τ2+τ3+τ4)*(τ2μ1+τ3μ2+τ4)Mn′MnΔtn3+12(τ1+τ2+τ3+τ4)(τ2μ12+τ3μ22+τ4)Mn′′MnΔtn4+12(τ2μ1+τ3μ2+τ4)2(Mn′)2Δtn4+16(τ1+τ2+τ3+τ4)3Mn3Δtn3+12(τ1+τ2+τ3+τ4)(τ2μ1+τ3μ2+τ4)Mn′Mn2Δtn4+124Mn4Δtn4]un+O(Δtn5).
That is,(2.35)un+1=[1+JnΔtn+12(JnΔtn)2+13!(JnΔtn)3+14!(JnΔtn)4]un+O(Δtn5),
where
(2.36)Jn=τ1Mn+τ2Mn+μ1+τ3Mn+μ2+τ4Mn+1,τ2=2μ2-112(-μ2+μ1)(μ1-1)μ1,τ3=--1+2μ112μ2(μ1μ2+μ2-μ22-μ1),τ4=6μ1μ2-4μ2+3-4μ112(μ2-1)(μ1-1).
We now turn to examine the more general situation in which the reaction term is not zero, that is,b≠0. In this case, we have(2.37)un+1=(∑i=0m(ΔtnJn)ii!)un+F(bn,bn+μ1,…,bn+1)+O(Δtni+1).
Similar to our discussion in the previous section, the accumulation error for (2.37) is
(2.38)|un+1-vn+1|=|∑i=0m(ΔtnJn)ii!|×|un-vn|+O(Δtni+1).
We remark again that the amplification factor∑i=0m((ΔtnJn)i/i!) in (2.38) at the nth local step has the same form with (2.5) except M being replaced by Jn. Therefore, we can apply the strategy demonstrated in Theorems 2.5 and 2.6 to look for the optimal time steps for Runge-Kutta schemes. The present amax and amin appeared in the optimal step size formulas then, respectively, represent the greatest and smallest real parts of the eigenvalues of the matrix Jn. Notice that by the similar argument outlined in the previous proofs, the stability is preserved with this step size.An Explicit Adaptive Optimal Time Stepping Scheme for Solving (2.18)
(1)
Set up initial conditionu0 and initial step size Δt0.(2)
SetJn in accordance with the stage number m, and apply m-stage Runge-Kutta schemes.(3)
Use the Gerschgorin theorem to estimate the interval,[-amaxn,-aminn] for negative real parts of the eigenvalues of Jn.(4)
PlaceΔtn+1,
(i)
m=2:Δtn+1=2/(amaxn+aminn),(ii)
m=3:Δtn+1=2.512745327/amaxn,(iii)
m=4:(2.39)Δtn+1=13((+30amax⁡amin⁡5+29amin⁡j⁡6)172amax⁡6+192amax⁡5amin⁡+276amax⁡4amin⁡2+448amax⁡3amin⁡3+276amin⁡4amax⁡2+192amax⁡amin⁡5+172amin⁡6+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amax⁡3+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amin⁡amax⁡2+36(29amax⁡6+30amax⁡5+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amax⁡amin⁡2+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amin⁡3(29agg))1/3/(amax⁡3+amin⁡amax⁡2+amax⁡amin⁡2+amin⁡3)-43(5amax⁡4+10amax⁡3amin⁡+6amax⁡2amin⁡2+10amax⁡amin⁡3+5amin⁡4)/(∑222(amax⁡3+amin⁡amax⁡2+amax⁡amin⁡2+amin⁡3)×((+30amax⁡amin⁡5+29amin⁡j⁡6)172amax⁡6+192amax⁡5amin⁡+276amax⁡4amin⁡2+448amax⁡3amin⁡3+276amin⁡4amax⁡2+192amax⁡amin⁡5+172amin⁡6+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2ar3+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amin⁡amax⁡2+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amax⁡amin⁡2+36(29amax⁡6+30amax⁡5amin⁡+51amax⁡4amin⁡2+68amax⁡3amin⁡3+51amin⁡4amax⁡2+30amax⁡amin⁡5+29amin⁡6)1/2amin⁡3(29agg))1/3(amin⁡3)3)+43(amax⁡2+amax⁡amin⁡+amin⁡2)/(amax⁡3+amin⁡amax⁡2+amax⁡amin⁡2+amin⁡3),
and then let tn+1=tn+Δtn+1. When an amaxn is equal to or close to zero, simply choose a reasonable step size for Δtn+1.
Note that the time step form=4 can be solved by mathematical software such as MAPLE from the following equation:
(2.40)∑i=141i!(amax⁡i-amin⁡i)Δtn+1i-1=0.


### 2.2.1. The “Best” Optimal Step Size and the Time-Dependent Homogeneous Systems

We point out that for those eigenvalues with negative real parts, one may also assume they are enclosed in an ellipse with one of the vertex lying at the origin, that is,aminn=0. In this way, only amaxn are needed in forming optimal time step formulas and time steps Δtn+1 in step 4 are switched as(i)
m=2:Δtn+1=2/amaxn,(ii)
m=3:Δtn+1=2.512745327/amaxn,(iii)
m=4:Δtn+1=2.785293563/amaxn,whereamaxn is defined as before in step 2. We remark that when Δtn+1 is specified as above values, pm(a,Δt) is equal to one. The stability condition does not seem to be satisfied. However, this hidden problem is avoidable by requiring Gerschgorin theorem or other equivalent theorems in step 3 as an estimation for the distribution of real part of eigenvalues. This is because the theorem provides an overestimation for the distribution of negative real part of eigenvalues. As a result, pm(a,Δt) is less than one. Hence, it is safe to use the above time step.In particular, when the reaction vector is negligible, we can further write them-stage Runge-kutta method as some equivalent m-stage formulas while the order of accuracy remains invariant. For example, m=3, we have(2.41)un(0)=un+Δtn3Jnun,un(1)=un+Δtn2Jnun(0),un+1=un+ΔtJnun(1).Form=4,
(2.42)un(0)=un+Δtn4Jnun,un(1)=un+Δtn3Jnun(0),un(2)=un+Δtn2Jnun(1),un+1=un+ΔtnJnun(2).Of course,Jn is defined according to their stage number.


## 2.2.1. The “Best” Optimal Step Size and the Time-Dependent Homogeneous Systems

We point out that for those eigenvalues with negative real parts, one may also assume they are enclosed in an ellipse with one of the vertex lying at the origin, that is,aminn=0. In this way, only amaxn are needed in forming optimal time step formulas and time steps Δtn+1 in step 4 are switched as(i)
m=2:Δtn+1=2/amaxn,(ii)
m=3:Δtn+1=2.512745327/amaxn,(iii)
m=4:Δtn+1=2.785293563/amaxn,whereamaxn is defined as before in step 2. We remark that when Δtn+1 is specified as above values, pm(a,Δt) is equal to one. The stability condition does not seem to be satisfied. However, this hidden problem is avoidable by requiring Gerschgorin theorem or other equivalent theorems in step 3 as an estimation for the distribution of real part of eigenvalues. This is because the theorem provides an overestimation for the distribution of negative real part of eigenvalues. As a result, pm(a,Δt) is less than one. Hence, it is safe to use the above time step.In particular, when the reaction vector is negligible, we can further write them-stage Runge-kutta method as some equivalent m-stage formulas while the order of accuracy remains invariant. For example, m=3, we have(2.41)un(0)=un+Δtn3Jnun,un(1)=un+Δtn2Jnun(0),un+1=un+ΔtJnun(1).Form=4,
(2.42)un(0)=un+Δtn4Jnun,un(1)=un+Δtn3Jnun(0),un(2)=un+Δtn2Jnun(1),un+1=un+ΔtnJnun(2).Of course,Jn is defined according to their stage number.


## 3. Numerical Examples

### 3.1. Overview

In this section, we demonstrate the implementation and the efficiency of our optimal adaptive time stepping method by setting a one-dimension testing problem and simulating some problems from a series of real biological experiments. The 1D testing problem with the known analytic solution is used to show the implementation process of our scheme. It is also designed to show the rate of convergence as well as the computational cost. Those results are illustrated through comparison with standardm-stage Runge-Kutta methods and the Runge-Kutta-Fehlberg (RKF) method. To test our numerical methods on real biological experiments, we apply our optimal adaptive two-stage Runge-Kutta method to a two-dimensional reaction-diffusion-chemotaxis model for a series of biological experiments on bacterial chemotaxis reported in [9]. Some results will be shown here.


### 3.2. Testing Problems

Consider(3.1)ut=uxx-(txu)x+f,x∈(0,2π),t>0,
where u=u(x,t) and f=(tx-1)e-tcos(x-t)+te-tsin(x-t). The three terms on the right-hand side represent diffusion, chemotaxis, and reaction, respectively. Assume the periodic boundary condition: u(o,t)=u(2π,t),t≥0 and the initial condition: u(x,0)=sinx. The analytic solution of (3.1) is(3.2)u(x,t)=e-tsin(x-t).To apply MOL, we simply use the difference scheme for the spatial discretization. The standard second-order difference is applied to the diffusion term. The upwind scheme is employed in the chemotaxis for some stability consideration [17]. After the discretization, we then apply m-stage Runge-Kutta method with the adaptive optimal time stepping scheme as the time integration strategy. By denoting h=2π/N, xi=ih, i=0,1,…,N, and the numerical solution as ui, the resulting ODE system is  (3.3)dudt=M(t)u+f(t),
where (3.4)u=(u0,u1,…,uN)T,f(t)=(f(x0,t),f(x1,t),…,f(xN,t))T,M(tn)=(-2h2-x0tn-1h1h21h2-tn-101h2+x0tn-1h-2h2-x1tn-1h1h2001h2+x1tn-1h-2h2-x2tn-1h1h200⋱⋱⋱⋱⋱⋱⋱⋱⋱01h201h2+xN-1tn-1h-2h2-xNtn-1h).
Generally, the actual eigenvalues of the spatial matrix M(t) are not available. One may use the power method and the inverse power method [20] to estimate the largest and smallest magnitude of eigenvalues. But these methods are usually time consuming. Moreover, only negative real parts of the eigenvalues are needed for our scheme. Therefore, we suggest using Gerschgorin theorem to estimate the distribution of eigenvalues. As we should show immediately, it is a simple and inexpensive way to find the range for real parts of the eigenvalues. It is therefore adequate for implementing our scheme. In fact, if a formula for the distribution of eigenvalues or actual eigenvalues can be found ahead of time, then it is not necessary to use the theorem to estimate the range at each time step. A simple version of the Gerschgorin theorem is stated as follows.Theorem 3.1 (Gerschgorin).
The union of all discs(3.5)Ki:={μ∈|μ-aii|≤∑k=1,k≠in|aik|}
contains all eigenvalues of the  n×n matrix  A=[aik].The theorem states that ifA is a real matrix, the range of the real parts of the eigenvalues can be easily detected by checking the intersection intervals of these discs with the real axis.From Gerschgorin Theorem, we obtain(3.6)-amax⁡≤Re(λ)<0,
where  (3.7)amax⁡=4Δx2+tn(4πΔx+1).Thus, form-stages, where m=2,…,4, the optimal time steps at the current step are(3.8)m=2:Δtopt=2(4/Δx2)+tn((4π/Δx)+1),m=3:Δtopt=2.512745327(4/Δx2)+tn((4π/Δx)+1),m=4:Δtopt=2.785293563(4/Δx2)+tn((4π/Δx)+1).Figure3 shows the relative error between the numerical and the analytic solutions. Observe that the relative error for the standard two-stage Runge-Kutta method increases dramatically around t=5.2, but the optimal time-stepping m-stage Runge-Kutta methods still converge to the analytic solution.Figure 3
Relative error between numerical simulations and the analytic solution.Tables1, 2, and 3 demonstrate the number of iteration for our scheme (Optimal m-stage RK) with that of two other methods: standard m-stage Runge-Kutta method with fixed step size (Std m-stage RK) and the Runge-Kutta-Fehlberg methods (RKF). These tables basically indicate that while the optimal time stepping scheme always ensure convergence, the standard m-stage Runge-Kutta methods or other adaptive methods like RKF either diverge or converge slower since they require a larger number of iterations. Similar results can also be illustrated by calculating CPU time following from Tables 4, 5, and 6. The computations were done on Dell Latitude D620 2600.Table 1
Comparison for number of iterations among optimal two-stage RK, RKF 23, and standard two-stage RK. The optimal two-stage RK method is counted twice for the function evaluation per iteration, while the RKF 23 is counted four times per iteration. “Div” is the abbreviation for Divergent.
Method/Time12345Optimal two-stage RK6381278191825583200RKF 2311752765431060758250Std two-stage RK (Δt=0.0025)80016002400DivDivTable 2
Comparison for number of iterations among the optimal three-stage RK, RKF 34, and standard three-stage RK. The optimal three-stage RK method is counted three times for the function evaluation per iteration, while the RKF 34 is five times per iteration.
Method/Time12345Optimal three-stage RK6781524253537175067RKF 34392563257935965511570Std three-stage RK (Δt=0.003)100220013000DivDivTable 3
Comparison for number of iterations among the optimal four-stage RK, RKF 45, and standard four-stage RK. The optimal four-stage RK method is counted four times for the function evaluation per iteration while the RKF 45 is six times per iteration.
Method/Time12345Optimal four-stage RK8161832305244726096RKF 45348055987152891010866Std four-stage RK (Δt=0.004)10002000DivDivDivTable 4
Comparison of CPU time among the optimal two-stage RK, RKF 23, and standard two-stage RK withΔt=0.001.
CPU/Time1020304050Optimal two-stage RK0.59401.84303.79606.46909.7820RKF 231.39103.90607.453012.062017.6090Std two-stage RK0.98403.5310DivDivDivTable 5
Comparison of CPU time among the optimal three-stage RK, RKF 34, and standard three-stage RK withΔt=0.0005 and Δt=0.001, respectively.
CPU/Time1020304050Optimal three-stage RK0.89102.82805.79709.797015.0000RKF 341.71904.51508.563013.688020.0780Std three-stage RK3.17106.641011.063013.125019.2340Std three-stage RK1.68803.310DivDivDivTable 6
Comparison of CPU time among the optimal four-stage RK, RKF 45, and standard four-stage RK withΔt=0.001.
CPU/Time1020304050Optimal four-stage RK1.00003.20306.578011.156016.9530RKF 451.87505.00009.500015.328023.6410Std four-stage RK2.26604.6250DivDivDivTables4–6 illustrate the computational efficiency for our method with CPU time.


### 3.3. Further Discussion on the Efficiency of the Optimal RK Method

We continue our discussion about the computational efficiency for optimal RK methods with time steps mentioned in Sections2.2. We are interested in measuring the efficiency difference between RKF and our method. To make one's point, we take the optimal four-stage RK method as an example. Let us start with the “best” optimal step size as in Section 2.2.1. Then the current time step for the optimal four-stage RK is(3.9)Δt=2.785293563amax⁡.To get the “worst" case, setamax=amin. Then  (3.10)Δt=1.596071638amax⁡.Now, let us assignΔt to be the average of (3.9) and (3.10). We obtain(3.11)Δt=2.1906826005amax⁡.Table7 shows that even with the “medium” optimal time step (3.11), our numerical scheme still performs better than the RKF 45. In addition, the computational efficiency of RKF 45 is approximately equal to the optimal four-stage RK when the time step is chosen as (3.11).Table 7
Comparison of CPU time between the optimal four-stage RK with the average time step of the “best” and the “worst” cases, and the RKF 45.
CPU/Time1020304050Optimal four-stage RK1.00003.20306.578011.156016.9530Four-stage RKΔt= (3.11)1.26604.04708.360014.219021.5160RKF 451.87505.00009.500015.328023.6410


### 3.4. Computer Simulations of Chemotaxis in Bacteria

Our numerical experiment is based on the biological experiments reported in (2.28) and the corresponding mathematical model given in (2.28). Interesting patterns are developed when Escherichia coli or Salmonella typhimurium is incoulated on intermediates of the tricarboxylic acid (TCA) cycle. Bacteria form symmetrical rings of spots or strips if they are inoculated in semisolid agar on the intermediates of TCA cycle. When placed in a liquid medium and fed on intermediates of the TCA cycle, they secrete aspartate as a chemoattractant. Cells then arrange themselves into high density aggregates. In this case, aggregates form randomly as time evolves and fade in a time interval less than the bacteria generation. Here, attractant chemicals are excreted by the bacterial themselves instead of being added into the Petri dish artificially.A nondimensionalized mathematical models were proposed in (2.28) for the liquid medium experiments in (2.28). The model is a classical chemotaxis system, which reads as(3.12)∂u∂t=∇⋅(Du∇u-αu(1+v)2∇v),∂v∂t=Δv+ωu2μ+u2,
where u,v, and ω represent cell density, concentration of chemotactic substance, and succinate, respectively. Simulations were done with zero flux boundary condition for both u and v. Initial conditions are set to be zero for v and the random perturbation at u0=1 for the initial condition of u. The noise amplitude is of order 10-1. They were created by the C++ random number generator. The initial, boundary conditions for u and v, and the noise amplitude are chosen to be the same as described in (2.28) except for the random numbers which may not be identical due to the use of different softwares.Numerical simulations are showed in Figure4. Simulations did capture the characteristics as described in (2.28) qualitatively. The aggregates first are divided into small regions and then cojoin with others.Numerical simulations for bacterial aggregates they exposed to TCA cycle: from left to right, top to bottom. The aggregates are separated into small regions and then cojoin other aggregates. The aggregates become fewer but larger.
(a)(b)(c)(d)


## 3.1. Overview

In this section, we demonstrate the implementation and the efficiency of our optimal adaptive time stepping method by setting a one-dimension testing problem and simulating some problems from a series of real biological experiments. The 1D testing problem with the known analytic solution is used to show the implementation process of our scheme. It is also designed to show the rate of convergence as well as the computational cost. Those results are illustrated through comparison with standardm-stage Runge-Kutta methods and the Runge-Kutta-Fehlberg (RKF) method. To test our numerical methods on real biological experiments, we apply our optimal adaptive two-stage Runge-Kutta method to a two-dimensional reaction-diffusion-chemotaxis model for a series of biological experiments on bacterial chemotaxis reported in [9]. Some results will be shown here.


## 3.2. Testing Problems

Consider(3.1)ut=uxx-(txu)x+f,x∈(0,2π),t>0,
where u=u(x,t) and f=(tx-1)e-tcos(x-t)+te-tsin(x-t). The three terms on the right-hand side represent diffusion, chemotaxis, and reaction, respectively. Assume the periodic boundary condition: u(o,t)=u(2π,t),t≥0 and the initial condition: u(x,0)=sinx. The analytic solution of (3.1) is(3.2)u(x,t)=e-tsin(x-t).To apply MOL, we simply use the difference scheme for the spatial discretization. The standard second-order difference is applied to the diffusion term. The upwind scheme is employed in the chemotaxis for some stability consideration [17]. After the discretization, we then apply m-stage Runge-Kutta method with the adaptive optimal time stepping scheme as the time integration strategy. By denoting h=2π/N, xi=ih, i=0,1,…,N, and the numerical solution as ui, the resulting ODE system is  (3.3)dudt=M(t)u+f(t),
where (3.4)u=(u0,u1,…,uN)T,f(t)=(f(x0,t),f(x1,t),…,f(xN,t))T,M(tn)=(-2h2-x0tn-1h1h21h2-tn-101h2+x0tn-1h-2h2-x1tn-1h1h2001h2+x1tn-1h-2h2-x2tn-1h1h200⋱⋱⋱⋱⋱⋱⋱⋱⋱01h201h2+xN-1tn-1h-2h2-xNtn-1h).
Generally, the actual eigenvalues of the spatial matrix M(t) are not available. One may use the power method and the inverse power method [20] to estimate the largest and smallest magnitude of eigenvalues. But these methods are usually time consuming. Moreover, only negative real parts of the eigenvalues are needed for our scheme. Therefore, we suggest using Gerschgorin theorem to estimate the distribution of eigenvalues. As we should show immediately, it is a simple and inexpensive way to find the range for real parts of the eigenvalues. It is therefore adequate for implementing our scheme. In fact, if a formula for the distribution of eigenvalues or actual eigenvalues can be found ahead of time, then it is not necessary to use the theorem to estimate the range at each time step. A simple version of the Gerschgorin theorem is stated as follows.Theorem 3.1 (Gerschgorin).
The union of all discs(3.5)Ki:={μ∈|μ-aii|≤∑k=1,k≠in|aik|}
contains all eigenvalues of the  n×n matrix  A=[aik].The theorem states that ifA is a real matrix, the range of the real parts of the eigenvalues can be easily detected by checking the intersection intervals of these discs with the real axis.From Gerschgorin Theorem, we obtain(3.6)-amax⁡≤Re(λ)<0,
where  (3.7)amax⁡=4Δx2+tn(4πΔx+1).Thus, form-stages, where m=2,…,4, the optimal time steps at the current step are(3.8)m=2:Δtopt=2(4/Δx2)+tn((4π/Δx)+1),m=3:Δtopt=2.512745327(4/Δx2)+tn((4π/Δx)+1),m=4:Δtopt=2.785293563(4/Δx2)+tn((4π/Δx)+1).Figure3 shows the relative error between the numerical and the analytic solutions. Observe that the relative error for the standard two-stage Runge-Kutta method increases dramatically around t=5.2, but the optimal time-stepping m-stage Runge-Kutta methods still converge to the analytic solution.Figure 3
Relative error between numerical simulations and the analytic solution.Tables1, 2, and 3 demonstrate the number of iteration for our scheme (Optimal m-stage RK) with that of two other methods: standard m-stage Runge-Kutta method with fixed step size (Std m-stage RK) and the Runge-Kutta-Fehlberg methods (RKF). These tables basically indicate that while the optimal time stepping scheme always ensure convergence, the standard m-stage Runge-Kutta methods or other adaptive methods like RKF either diverge or converge slower since they require a larger number of iterations. Similar results can also be illustrated by calculating CPU time following from Tables 4, 5, and 6. The computations were done on Dell Latitude D620 2600.Table 1
Comparison for number of iterations among optimal two-stage RK, RKF 23, and standard two-stage RK. The optimal two-stage RK method is counted twice for the function evaluation per iteration, while the RKF 23 is counted four times per iteration. “Div” is the abbreviation for Divergent.
Method/Time12345Optimal two-stage RK6381278191825583200RKF 2311752765431060758250Std two-stage RK (Δt=0.0025)80016002400DivDivTable 2
Comparison for number of iterations among the optimal three-stage RK, RKF 34, and standard three-stage RK. The optimal three-stage RK method is counted three times for the function evaluation per iteration, while the RKF 34 is five times per iteration.
Method/Time12345Optimal three-stage RK6781524253537175067RKF 34392563257935965511570Std three-stage RK (Δt=0.003)100220013000DivDivTable 3
Comparison for number of iterations among the optimal four-stage RK, RKF 45, and standard four-stage RK. The optimal four-stage RK method is counted four times for the function evaluation per iteration while the RKF 45 is six times per iteration.
Method/Time12345Optimal four-stage RK8161832305244726096RKF 45348055987152891010866Std four-stage RK (Δt=0.004)10002000DivDivDivTable 4
Comparison of CPU time among the optimal two-stage RK, RKF 23, and standard two-stage RK withΔt=0.001.
CPU/Time1020304050Optimal two-stage RK0.59401.84303.79606.46909.7820RKF 231.39103.90607.453012.062017.6090Std two-stage RK0.98403.5310DivDivDivTable 5
Comparison of CPU time among the optimal three-stage RK, RKF 34, and standard three-stage RK withΔt=0.0005 and Δt=0.001, respectively.
CPU/Time1020304050Optimal three-stage RK0.89102.82805.79709.797015.0000RKF 341.71904.51508.563013.688020.0780Std three-stage RK3.17106.641011.063013.125019.2340Std three-stage RK1.68803.310DivDivDivTable 6
Comparison of CPU time among the optimal four-stage RK, RKF 45, and standard four-stage RK withΔt=0.001.
CPU/Time1020304050Optimal four-stage RK1.00003.20306.578011.156016.9530RKF 451.87505.00009.500015.328023.6410Std four-stage RK2.26604.6250DivDivDivTables4–6 illustrate the computational efficiency for our method with CPU time.


## 3.3. Further Discussion on the Efficiency of the Optimal RK Method

We continue our discussion about the computational efficiency for optimal RK methods with time steps mentioned in Sections2.2. We are interested in measuring the efficiency difference between RKF and our method. To make one's point, we take the optimal four-stage RK method as an example. Let us start with the “best” optimal step size as in Section 2.2.1. Then the current time step for the optimal four-stage RK is(3.9)Δt=2.785293563amax⁡.To get the “worst" case, setamax=amin. Then  (3.10)Δt=1.596071638amax⁡.Now, let us assignΔt to be the average of (3.9) and (3.10). We obtain(3.11)Δt=2.1906826005amax⁡.Table7 shows that even with the “medium” optimal time step (3.11), our numerical scheme still performs better than the RKF 45. In addition, the computational efficiency of RKF 45 is approximately equal to the optimal four-stage RK when the time step is chosen as (3.11).Table 7
Comparison of CPU time between the optimal four-stage RK with the average time step of the “best” and the “worst” cases, and the RKF 45.
CPU/Time1020304050Optimal four-stage RK1.00003.20306.578011.156016.9530Four-stage RKΔt= (3.11)1.26604.04708.360014.219021.5160RKF 451.87505.00009.500015.328023.6410


## 3.4. Computer Simulations of Chemotaxis in Bacteria

Our numerical experiment is based on the biological experiments reported in (2.28) and the corresponding mathematical model given in (2.28). Interesting patterns are developed when Escherichia coli or Salmonella typhimurium is incoulated on intermediates of the tricarboxylic acid (TCA) cycle. Bacteria form symmetrical rings of spots or strips if they are inoculated in semisolid agar on the intermediates of TCA cycle. When placed in a liquid medium and fed on intermediates of the TCA cycle, they secrete aspartate as a chemoattractant. Cells then arrange themselves into high density aggregates. In this case, aggregates form randomly as time evolves and fade in a time interval less than the bacteria generation. Here, attractant chemicals are excreted by the bacterial themselves instead of being added into the Petri dish artificially.A nondimensionalized mathematical models were proposed in (2.28) for the liquid medium experiments in (2.28). The model is a classical chemotaxis system, which reads as(3.12)∂u∂t=∇⋅(Du∇u-αu(1+v)2∇v),∂v∂t=Δv+ωu2μ+u2,
where u,v, and ω represent cell density, concentration of chemotactic substance, and succinate, respectively. Simulations were done with zero flux boundary condition for both u and v. Initial conditions are set to be zero for v and the random perturbation at u0=1 for the initial condition of u. The noise amplitude is of order 10-1. They were created by the C++ random number generator. The initial, boundary conditions for u and v, and the noise amplitude are chosen to be the same as described in (2.28) except for the random numbers which may not be identical due to the use of different softwares.Numerical simulations are showed in Figure4. Simulations did capture the characteristics as described in (2.28) qualitatively. The aggregates first are divided into small regions and then cojoin with others.Numerical simulations for bacterial aggregates they exposed to TCA cycle: from left to right, top to bottom. The aggregates are separated into small regions and then cojoin other aggregates. The aggregates become fewer but larger.
(a)(b)(c)(d)


## 4. Conclusions and Discussions

We investigate an optimal adaptive time stepping strategy for fully explicitm-stage Runge-Kutta methods to solve general reaction-diffusion-chemotaxis systems. It is designed specially for handling spatially nonhomogeneous pattern formations. The method is based on m-stage Runge-Kutta methods and simple estimation on the optimal step size at each iteration. With some optimization techniques, we efficiently discard the decay modes and only allow the growth modes to grow. In this way, we guarantee that the numerically generated patterns follow the primary model equations. The scheme is proven to be extremely robust and accurate by theoretical analysis, numerical experiments on the testing problem, and real pattern simulations on bacterial chemotaxis. The numerical experiments also show that our method is superior to RKF in solving reaction-diffusion-chemotaxis systems.



---

*Source: 389207-2011-06-28.xml*