Adaptive Optimal 𝑚-Stage Runge-Kutta Methods for Solving Reaction-Diffusion-Chemotaxis Systems

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 explicit 𝑚-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.

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 in E. 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)

Here, 𝑢 represents the microorganism (bacteria) density; 𝑝 stands for chemoattractant concentration and is given or obtained from another equation. 𝑓 is the reaction term which models the interactions between the chemicals. The diffusion term, 𝐷_𝑢Δ𝑢, imitates the random motion of particles with positive constant diffusivity 𝐷_𝑢. ∇⋅(𝑢𝜒(𝑝)∇𝑝) describes the directed random movement of particles 𝑝 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) where 𝐶₀ and 𝐶₁ are positive constants and 𝑠≥1 is an integer power. Observe that lim_𝑥→∞𝜒(𝑥)=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) where 𝑀(𝑡) is a time-dependent matrix and 𝐛(𝐭) 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 explicit 𝑚-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) then generalize it to the time-dependent coefficients case. The general formulation of the 𝑚-stage Runge-Kutta method is𝑘₁=𝐟𝑡_𝑛,𝐮_𝑛,𝑘_𝑞=𝐟𝑡_𝑛+𝛼_𝑞Δ𝑡,𝐮_𝑛+Δ𝑡_𝑛^𝑞−1_𝑠=1𝛽_𝑞𝑠𝑘_𝑠𝛼,𝑞=2,…,𝑚,_𝑞=^𝑞−1_𝑠=1𝛽_𝑞𝑠,𝐮_𝑛+1=𝐮_𝑛+Δ𝑡_𝑛^𝑞_𝑠=1𝛾_𝑠𝑘_𝑠,𝑞=2,…,𝑚, (2.2) where 𝑚 is the stage number and 𝛼_𝑞=∑^𝑞−1_𝑗=1𝛽_𝑞𝑗,𝑞=2,…,𝑚. The subscripts 𝑛 and 𝑛+1 are the iteration level. For simplicity, we first look at the uniform time step size Δ𝑡=𝑡_𝑛−𝑡_𝑛−1 for all 𝑛. Then 𝐮_𝑛≈𝐮(𝐱_𝑛) is the numerical solution at the 𝑛th step. Equations for the unknowns 𝛾_𝑠 and 𝛽_𝑞𝑠 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 set 𝐛=0 in (2.1). As in [5, 8], we first fix all the intermediate coefficients at 𝑡=𝑡_𝑛 for the three-stage Runge-Kutta method. The general time dependence case is proposed in the next section. Denoting 𝑀=𝑀(𝑡_𝑛) and applying (2.2) to the three-stage Runge-Kutta method with all the intermediate values being fixed at 𝑡_𝑛, we arrive at𝑘₁=𝑀𝐮_𝑛,𝑘₂=𝑀1+𝛼₂Δ𝑡𝑀𝐮_𝑛,𝑘₃=𝑀1+𝛼₃Δ𝑡𝑀+𝛼₂𝛽₃₂Δ𝑡²𝑀²𝐮_𝑛,𝐮_𝑛+1=1+𝛿₁Δ𝑡𝑀+𝛿₂Δ𝑡²𝑀²+𝛿₃Δ𝑡³𝑀³𝐮_𝑛, (2.3)(2.4) where 𝛿_𝑖,𝑖=1,2,3, are combinations of the coefficients 𝛽_𝑞𝑠 and 𝛾_𝑠. Since 𝐮(𝑡_𝑛+1)=𝑒^𝑀Δ𝑡𝐮(𝑡_𝑛) and assuming 𝐮_𝑛 is exact up to the order Δ𝑡³, we have 𝛿₁=1, 𝛿₂=(1/2!), 𝛿₃=(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]:𝐮_𝑛+1=11+Δ𝑡𝑀+()2!Δ𝑡𝑀²1+⋯+()𝑚!Δ𝑡𝑀^𝑚𝐮_𝑛.=𝑝_𝑚(Δ𝑡)𝐮_𝑛. (2.5)

In the case of 𝐛≠0, (2.5) becomes𝐮_𝑛+1=𝑝_𝑚(Δ𝑡)𝐮_𝑛+1+Δ𝑡2!𝑀+⋯+Δ𝑡^𝑚𝑀𝑚!^𝑚−1𝐛. (2.6)

Denote 𝐯 as the solution by solving the system with the perturbated initial condition; we have𝐯_𝑛+1=𝑝_𝑚(Δ𝑡)𝐯_𝑛+1+Δ𝑡2!𝑀+⋯+Δ𝑡^𝑚𝑀𝑚!^𝑚−1𝐛. (2.7)

From (2.6) and (2.7), the accumulation error can be written as 𝐮_𝑛+1−𝐯_𝑛+1=𝑝_𝑚(Δ𝑡)×𝐮_𝑛−𝐯_𝑛. (2.8)

Note that 𝑝_𝑚(𝜆)1=1+Δ𝑡𝜆+()2!Δ𝑡𝜆²1+⋯+()𝑚!Δ𝑡𝜆^𝑚 (2.9) is called the amplification factor, where 𝜆 is an eigenvalue of 𝑀. It is obtained when a scalar test problem, 𝑢′=𝜆𝑢, is applied to the single-step method.

Suppose that 𝜆=𝑎+𝑏𝑖=Re𝜆+Im𝜆𝑖. 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 |𝑝_𝑚(𝜆)|<1 when 𝑎<0 and |𝑝_𝑚(𝜆)|≤1+𝑎Δ𝑡+𝑂(Δ𝑡) when 𝑎≥0.

Here 𝑂(Δ𝑡) means higher order of Δ𝑡. 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 𝑎<0 or when the scheme is implicit and 𝑎>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 𝑎 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 𝑚-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 the 𝑚-stage Runge-Kutta method satisfies the stability condition automatically for 𝑎>0 but needs constraints on Δ𝑡 when 𝑎<0.

Theorem 2.2. |𝑝_𝑚(𝜆)|²=|𝑝_𝑚(𝑎+𝑏𝑖)|²=(1+Δ𝑡𝑎+(1/2)Δ𝑡²𝑎²+⋯+(1/𝑚!)Δ𝑡^𝑚𝑎^𝑚)²+𝑂(Δ𝑡^𝑚+1), where 𝑚=1,…,4.

Proof. The proof for 𝑚=1 is straightforward, and the proof for 𝑚=2 is in [5, 8]. For 𝑚=3, 𝑝₃(𝜆)²=1+Δ𝑡𝑎+Δ𝑡²2𝑎²−𝑏²+Δ𝑡³6𝑎³−3𝑎𝑏²²+Δ𝑡𝑏+Δ𝑡²𝑎𝑏+Δ𝑡³63𝑎²𝑏−𝑏³²=11+Δ𝑡𝑎+2(Δ𝑡𝑎)²+16(Δ𝑡𝑎)³²()1−21+Δ𝑡𝑎2Δ𝑡²𝑏²+12Δ𝑡³𝑎𝑏²+𝑂Δ𝑡⁴+Δ𝑡𝑏+Δ𝑡²𝑎𝑏+Δ𝑡³63𝑎²𝑏−𝑏³²=11+Δ𝑡𝑎+2(Δ𝑡𝑎)²+16(Δ𝑡𝑎)³²−Δ𝑡²𝑏²−2Δ𝑡³𝑎𝑏²+𝑂Δ𝑡⁴+Δ𝑡²𝑏²+2Δ𝑡³𝑎𝑏²+𝑂Δ𝑡⁴=11+Δ𝑡𝑎+2(Δ𝑡𝑎)²+16(Δ𝑡𝑎)³²+𝑂Δ𝑡⁴. (2.10) For 𝑚=4, 𝑝₄(𝜆)²=1+Δ𝑡𝑎+Δ𝑡²2𝑎²−𝑏²+Δ𝑡³6𝑎³−3𝑎𝑏²+Δ𝑡⁴𝑎24⁴−6𝑎²𝑏²+𝑏⁴²+Δ𝑡𝑏+Δ𝑡²𝑎𝑏+Δ𝑡³63𝑎²𝑏−𝑏³+Δ𝑡⁴244𝑎³𝑏−4𝑎𝑏³²=11+Δ𝑡𝑎+2(Δ𝑡𝑎)²+16(Δ𝑡𝑎)³+1()24Δ𝑡𝑎⁴²1−21+Δ𝑡𝑎+2(Δ𝑡𝑎)²+16(Δ𝑡𝑎)³+1()24Δ𝑡𝑎⁴∗12Δ𝑡²𝑏²+12Δ𝑡³𝑎𝑏²+14Δ𝑡⁴𝑎²𝑏²−124Δ𝑡⁴𝑏⁴+12Δ𝑡²𝑏²+12Δ𝑡³𝑎𝑏²+14Δ𝑡⁴𝑎²𝑏²−124Δ𝑡⁴𝑏⁴²+Δ𝑡𝑏+Δ𝑡²𝑎𝑏+Δ𝑡³63𝑎²𝑏−𝑏³+Δ𝑡⁴244𝑎³𝑏−4𝑎𝑏³²=11+Δ𝑡𝑎+2(Δ𝑡𝑎)²+16(Δ𝑡𝑎)³+1()24Δ𝑡𝑎⁴²+Δ𝑡⁴−2𝑎²𝑏²+1𝑏12⁴−2Δ𝑡³𝑎𝑏²−Δ𝑡²𝑏²+𝑂Δ𝑡⁵+14Δ𝑡⁴𝑏⁴+𝑂Δ𝑡⁵+Δ𝑡⁴−13𝑏⁴+2𝑎²𝑏²+2Δ𝑡³𝑎𝑏²+Δ𝑡²𝑏²+𝑂Δ𝑡⁵=11+Δ𝑡𝑎+2(Δ𝑡𝑎)²+16(Δ𝑡𝑎)³+1()24Δ𝑡𝑎⁴²+𝑂Δ𝑡⁵. (2.11)
The geometric interpretation is that, for spatial matrices 𝑀(𝑡_𝑛), the eigenvalues with negative real parts are assumed to be contained in the left-plane ellipse. Let us denote 1+Δ𝑡𝑎+(1/2)Δ𝑡²𝑎²+⋯+(1/𝑚!)Δ𝑡^𝑚𝑎^𝑚 as 𝑝_𝑚(𝑎,Δ𝑡). To examine whether numerically generated patterns truly obey the original models, one needs to verify if the negative modes (𝑎<0) are diminished and positive modes (𝑎>0) are growing accordingly. That is, when the negative modes occur, 𝑝_𝑚(𝑎,Δ𝑡) is less than one, and is greater than one when the spatial matrix 𝑀(𝑡_𝑛) happens to have positive modes. From 𝑝_𝑚(𝑎,Δ𝑡), we find that when 𝑎>0, 𝑝_𝑚(𝑎,Δ𝑡) is definitely greater than unity, hence the stability condition is automatically satisfied; when 𝑎<0, 𝑝_𝑚(𝑎,Δ𝑡) is then possible to be less than unity. It is necessary to put the limitation on the time steps to ensure the stability condition 𝑝_𝑚(𝑎,Δ𝑡)<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 𝑚-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. For 𝑎<0 and 𝑚 being even, 𝑝_𝑚(𝑎,Δ𝑡)>0 and is concave up for all Δ𝑡>0.

Proof. For 𝑚=2, it is a second-degree polynomial with a positive leading coefficient. The proof is straightforward. Note that 𝑝₂(𝑎,Δ𝑡)>0.
Assume 𝑝_𝑚−2(𝑎,Δ𝑡)>0 for all even 𝑚. For 𝑚 being even and 𝑚>2, we have 𝜕²𝑝_𝑚(𝑎,Δ𝑡)𝜕Δ𝑡²=𝑎²𝑝_𝑚−2(𝑎,Δ𝑡)>0. (2.12) It follows from mathematical induction, 𝜕²/𝜕Δ𝑡²𝑝_𝑚(𝑎,Δ𝑡)>0 for all even 𝑚. Hence, 𝑝_𝑚(𝑎,Δ𝑡) is concave up w.r.t. Δ𝑡 if 𝑚 is even.
Moreover, 𝑝_𝑚(𝑎,Δ𝑡)>0 for all even 𝑚, because 𝑝_𝑚(𝑎,Δ𝑡)>exp(𝑎Δ𝑡)>0 for even 𝑚.

Corollary 2.4. For 𝑎<0 and 𝑚 being an odd number, 𝑝_𝑚(𝑎,Δ𝑡)<1 and 𝑝_𝑚(𝑎,Δ𝑡) is strictly decreasing for all Δ𝑡>0.

Proof. For 𝑚=1, the proof is trivial. For 𝑚>1,𝑝_𝑚(𝑎,0)=1 and 𝜕𝑝_𝑚(𝑎,Δ𝑡)𝜕Δ𝑡=𝑎𝑝_𝑚−1(𝑎,Δ𝑡)<0 (2.13) since 𝑎<0 and 𝑚−1 is an even number.
Therefore, we conclude 𝑝_𝑚(𝑎,Δ𝑡)<1 for all odd 𝑚 and Δ𝑡>0. Moreover, 𝑝_𝑚(𝑎,Δ𝑡) is strictly decreasing with respect to Δ𝑡.

Theorem 2.5 (Case for 𝑚 being odd). Assume that −𝑎_max≤𝑎≤−𝑎_min, where 𝑎_max>𝑎_min>0 are two positive numbers. Then the solution of the following min-max problem occurs at Δ𝑡=Δ𝑡_opt in which Δ𝑡 satisfies 𝑓_𝑚(−𝑎_max,Δ𝑡)=0, where m is an odd number. The min-max problem is stated as min_Δ𝑡>0max_{−𝑎max≤𝑎≤−amin}𝑝_𝑚(𝑎,Δ𝑡). (2.14)

Proof. Note that 𝑝_𝑚(𝑎₁,Δ𝑡)>𝑝_𝑚(𝑎₂,Δ𝑡) if 𝑎₁>𝑎₂ and 𝑎₁,𝑎₂<0. For any fixed Δ𝑡, the maximum of 𝑝_𝑚(𝑎,Δ𝑡) is obtained at 𝑝_𝑚(−𝑎_min,Δ𝑡). So the min-max occurs when 𝑝_𝑚(−𝑎_max,Δ𝑡)=0. Since 𝑝_𝑚(𝑎,Δ𝑡) is continuous, 𝑝_𝑚(𝑎,0)=1>0, and 𝑝_𝑚(𝑎,Δ𝑡)<0 as Δ𝑡→∞, the solution exists by Intermediate Value Theorem. Furthermore, from Corollary 2.4, 𝑝_𝑚(𝑎,Δ𝑡) is monotone decreasing w.r.t. Δ𝑡, so the solution is unique.

Theorem 2.6 (Case for 𝑚 being even). Assume that −𝑎_max≤𝑎≤−𝑎_min, where 𝑎_max>𝑎_min>0 are two positive numbers. Then the solution of the following min-max problem occurs at Δ𝑡=Δ𝑡_opt in which Δt satisfies 𝑝_𝑚(−𝑎_max,Δ𝑡)=𝑝_𝑚(−𝑎_min,Δ𝑡), where m is an even number (𝑚=2,4). The min-max problem is stated as min_Δ𝑡>0max_{−𝑎max≤𝑎≤−𝑎min}𝑝_𝑚(𝑎,Δ𝑡). (2.15)

Proof. For each fixed 𝑎, 𝑝_𝑚(𝑎,Δ𝑡) is a quadratic-like polynomial function in Δ𝑡. Then for any fixed Δ𝑡, the maximum of 𝑝_𝑚(𝑎,Δ𝑡) is obtained at 𝑝_𝑚(−𝑎_min,Δ𝑡),Δ𝑡<Δ𝑡_opt, and at 𝑝_𝑚(−𝑎_max,Δ𝑡), Δ𝑡>Δ𝑡_opt. So the minimax occurs when 𝑝_𝑚(−𝑎_max,Δ𝑡)=𝑝_𝑚(−𝑎_min,Δ𝑡), that is, ^𝑚_𝑖=0−𝑎_maxΔ𝑡^𝑚=𝑚!^𝑚_𝑖=0−𝑎_minΔ𝑡^𝑚.𝑚! (2.16) This leads to the solution.

Let us denote 𝑔_𝑚(Δ𝑡)=𝑝_𝑚(−𝑎_max,Δ𝑡)−𝑝_𝑚(−𝑎_min,Δ𝑡). The following theorem states that the optimal time step is unique for 𝑚 being even. For 𝑚 being odd, it is shown in Theorem 2.5.

Theorem 2.7. If m is even, 𝑔_𝑚(Δ𝑡) has only one positive real solution.

Proof. 𝑔_𝑚(Δ𝑡)=Δ𝑡𝑎_max−𝑎_min^𝑚_𝑖=1(−Δ𝑡)^𝑚−1𝑎^𝑚_max−𝑎^𝑚_min𝑚!𝑎_max−𝑎_min=Δ𝑡𝑎_max−𝑎_min𝑞_𝑚(Δ𝑡). (2.17)
It follows from (2.17) and the fact that 𝑎_max≠𝑎_min, to find the zeros of 𝑔_𝑚(Δ𝑡) is enough to solve 𝑞_𝑚(Δ𝑡)=0. Since 𝑞_𝑚(0)=1 and 𝑞_𝑚(Δ𝑡)<0 as Δ𝑡 tends to infinity, positive real solutions hence exist for 𝑔_𝑚(Δ𝑡) by the Intermediate Value Theorem. Moreover, 𝑞_𝑚(Δ𝑡) is strictly decreasing (Figures 1 and 2). Therefore, the solution is unique.

Figure 1 

Graphs for 𝑝3(𝑎,Δ𝑡) where 𝑎=−6: ·, 𝑎=−10: *, 𝑎=−12: ∘.

Figure 2 

Graphs for 𝑝4(𝑎,Δ𝑡) where 𝑎=−6: ·, 𝑎=−10: *, 𝑎=−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)

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𝑢_𝑛+𝛽₂₁Δ𝑡𝑘₁,𝑘₃=𝑀_𝑛+𝛼3𝑢_𝑛+Δ𝑡𝛽₃₁𝑘₁+𝛽₃₂𝑘₂,𝐮_𝑛+1=𝐮_𝑛+Δ𝑡𝛾₁𝑘₁+𝛾₂𝑘₂+𝛾₃𝑘₃. (2.19)
Here, 𝑀_𝑛=𝑀(𝑡_𝑛). Expand 𝑀_{𝑛+𝛼𝑖} at 𝑡_𝑛 by Taylor expansion, then 𝑀_{𝑛+𝛼𝑖}=𝑀_𝑛+𝛼_𝑖Δ𝑡_𝑛𝑀^′_𝑛+((𝛼_𝑖Δ𝑡_𝑛)²/2!)𝑀_𝑛^′′+𝑂(Δ𝑡³), where 𝑖=2,3.
Therefore, 𝑘₁=𝑀_𝑛𝐮_𝑛,𝑘₂=𝑀_𝑛+𝛼₂𝑀^′_𝑛+𝛽₂₁𝑀²_𝑛Δ𝑡_𝑛+𝛼²₂2𝑀_𝑛^′′+𝛼₂𝛽₂₁𝑀_𝑛𝑀^′_𝑛Δ𝑡²_𝑛𝐮_𝑛+𝑂Δ𝑡³_𝑛,𝑘₃=𝑀_𝑛+𝛼₃𝑀^′_𝑛Δ𝑡_𝑛+𝛽₃₁+𝛽₃₂𝑀²_𝑛Δ𝑡_𝑛+𝛼²₃2𝑀_𝑛^′′Δ𝑡²_𝑛+𝛽₃₁𝛼₃+𝛽₃₂𝛼₂+𝛼₃𝑀^′_𝑛𝑀_𝑛Δ𝑡²_𝑛+𝛽₃₂𝛽₂₁𝑀³_𝑛Δ𝑡²_𝑛𝐮_𝑛+𝑂Δ𝑡³_𝑛. (2.20)
Since by consistency, ³_𝑖=1𝛾_𝑖𝛼=1,_𝑞=^𝑞−1_𝑗=1𝛽_𝑞𝑗𝛾,𝑞=2,3,₂𝛼₂+𝛾₃𝛼₃=12,𝛾₂𝛼²₂+𝛾₃𝛼²₃=13,𝛾₃𝛼₂𝛽₃₂=16. (2.21)
Followed by (2.20) and the consistency condition, (2.19) can be simplified as 𝐮_𝑛+1=𝐮_𝑛+𝑀_𝑛Δ𝑡_𝑛+12𝑀^′_𝑛Δ𝑡²_𝑛+16𝑀_𝑛^′′Δ𝑡³_𝑛+12𝑀²_𝑛Δ𝑡²_𝑛+12𝑀^′_𝑛𝑀_𝑛Δ𝑡³_𝑛+16𝑀³_𝑛Δ𝑡³_𝑛𝐮_𝑛+𝑂Δ𝑡⁴_𝑛. (2.22)
Now, we assume ³_𝑖=1𝜏_𝑖𝜏=1,₂𝜇+𝜏₃=12,𝜏₂𝜇²+𝜏₃=13. (2.23)
Then, we can rewrite (2.22) as 𝐮_𝑛+1=𝐮_𝑛+𝜏₁+𝜏₂+𝜏₃𝑀_𝑛Δ𝑡_𝑛+𝜏₂𝜇+𝜏₃𝑀^′_𝑛Δ𝑡²_𝑛+12𝜏₂𝜇²+𝜏₃𝑀_𝑛^′′Δ𝑡³_𝑛+12𝜏₁+𝜏₂+𝜏₃²𝑀²_𝑛Δ𝑡²_𝑛+𝜏₁+𝜏₂+𝜏₃𝜏₂𝜇+𝜏₃𝑀^′_𝑛𝑀_𝑛Δ𝑡³_𝑛+16𝜏₁+𝜏₂+𝜏₃³𝑀³_𝑛Δ𝑡³_𝑛𝐮_𝑛+𝑂Δ𝑡⁴_𝑛. (2.24)
This way, we arrive the following desirable form through rearrangement: 𝐮_𝑛+1=𝐮_𝑛+𝜏₁𝑀_𝑛+𝜏₂𝑀_𝑛+𝜇Δ𝑡_𝑛𝑀^′_𝑛+𝜇Δ𝑡_𝑛²2𝑀_𝑛^′′+𝜏₃𝑀_𝑛+Δ𝑡_𝑛𝑀^′_𝑛+Δ𝑡²_𝑛2𝑀_𝑛^′′+12𝜏₁+𝜏₂+𝜏₃𝑀_𝑛+𝜏₂+𝜏₃𝜇𝑀^′_𝑛Δ𝑡_𝑛+𝑂Δ𝑡²_𝑛²Δ𝑡²_𝑛+16𝜏₁+𝜏₂+𝜏₃³𝑀³_𝑛Δ𝑡³_𝑛𝐮_𝑛+𝑂Δ𝑡⁴_𝑛. (2.25)
Let 𝐽_𝑛=𝜏₁𝑀_𝑛+𝜏₂𝑀_𝑛+𝜇+𝜏₃𝑀_𝑛+1, (2.26) where 𝜏₂1=−,𝜏6𝜇𝜇−1₃=3𝜇−26.𝜇−1 (2.27)
Then the previous equation implies 𝐮_𝑛+1=1+𝐽_𝑛Δ𝑡_𝑛+12𝐽_𝑛Δ𝑡_𝑛²+1𝐽3!^𝑛Δ𝑡_𝑛³𝐮_𝑛+𝑂Δ𝑡⁴_𝑛. (2.28)
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 𝐽_𝑛 is replaced by 𝑀. 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𝐮_𝑛+𝛽₂₁Δ𝑡𝑘₁,𝑘₃=𝑀_𝑛+𝛼3𝐮_𝑛+Δ𝑡𝛽₃₁𝑘₁+𝛽₃₂𝑘₂,𝑘₄=𝑀_𝑛+𝛼4𝐮_𝑛+Δ𝑡𝛽₄₁𝑘₁+𝛽₄₂𝑘₂+𝛽₄₃𝑘₃,𝐮_𝑛+1=𝐮_𝑛+Δ𝑡𝑟₁𝑘₁+𝑟₂𝑘₂+𝑟₃𝑘₃+𝑟₄𝑘₄. (2.29)
Since 𝑀_{𝑛+𝛼𝑖}=𝑀_𝑛+𝛼_𝑖Δ𝑡_𝑛𝑀^′_𝑛+((𝛼_𝑖Δ𝑡_𝑛)²/2!)𝑀_𝑛^′′+((𝛼_𝑖Δ𝑡_𝑛)³/3!)𝑀_𝑛^′′′+𝑂(Δ𝑡⁴), 𝑘₂=𝑀_𝑛+𝛼₂𝑀^′_𝑛+𝛽₂₁𝑀²_𝑛Δ𝑡_𝑛+𝛼²₂2𝑀_𝑛^′′+𝛼₂𝛽₂₁𝑀_𝑛𝑀^′_𝑛Δ𝑡²_𝑛+𝛼³₂𝑀3!^{𝑛′′′}+𝛼²₂𝛽₂₁2𝑀_𝑛^′′𝑀_𝑛Δ𝑡³_𝑛𝐮_𝑛+𝑂Δ𝑡⁴_𝑛,𝑘₃=𝑀_𝑛+𝛼₃𝑀^′_𝑛Δ𝑡_𝑛+𝛽₃₁+𝛽₃₂𝑀²_𝑛Δ𝑡_𝑛+𝛼²₃2𝑀_𝑛^′′Δ𝑡²_𝑛+𝛽₃₁𝛼₃+𝛽₃₂𝛼₂+𝛼₃𝑀^′_𝑛𝑀_𝑛Δ𝑡²_𝑛+𝛽₃₂𝛽₂₁𝑀³_𝑛Δ𝑡²_𝑛+𝛼³₃𝑀3!^{𝑛′′′}Δ𝑡³_𝑛+12𝛽₃₁𝛼²₃+𝛽₃₂𝛼²₂+𝛼²₃𝑀_𝑛^′′𝑀_𝑛Δ𝑡³_𝑛+𝛽₃₂𝛼₂𝛼₃𝑀^′_𝑛²Δ𝑡³_𝑛+𝛽₃₂𝛽₂₁𝛼₂+𝛼₃𝑀^′_𝑛𝑀²_𝑛Δ𝑡³_𝑛𝐮_𝑛+𝑂Δ𝑡⁴_𝑛,𝑘₄=𝑀_𝑛+𝛼₄𝑀^′_𝑛Δ𝑡_𝑛+𝛽₄₁+𝛽₄₂+𝛽₄₃𝑀²_𝑛Δ𝑡_𝑛+12𝛼²₄𝑀_𝑛^′′Δ𝑡²_𝑛+𝛽₄₁𝛼₄+𝛽₄₂𝛼₂+𝛼₄+𝛽₄₃𝛼₃+𝛼₄𝑀^′_𝑛𝑀_𝑛Δ𝑡²_𝑛+𝛽₄₂𝛽₂₁+𝛽₄₃𝛽₃₁+𝛽₃₂𝑀³_𝑛Δ𝑡²_𝑛+1𝛼3!³⁴𝑀_𝑛^′′′Δ𝑡³_𝑛+12𝛽₄₁𝛼²₄+𝛽₄₂𝛼²₂+𝛼²₄+𝛽₄₃𝛼²₃+𝛼²₄𝑀_𝑛^′′𝑀_𝑛Δ𝑡³_𝑛+𝛽₄₂𝛽₂₁𝛼₂+𝛼₄+𝛽₄₃𝛽₃₁𝛼₃+𝛼₄+𝛽₄₃𝛽₃₂𝛼₂+𝛼₃+𝛼₄𝑀^′_𝑛𝑀²_𝑛Δ𝑡³_𝑛+𝛽₄₂𝛼₂𝛼₄+𝛽₄₃𝛼₃𝛼₄𝑀^′_𝑛²Δ𝑡³_𝑛+𝛽₄₃𝛽₃₂𝛽₂₁𝑀⁴_𝑛Δ𝑡³_𝑛𝐮_𝑛+𝑂Δ𝑡⁴_𝑛. (2.30)
By consistency of the method, ⁴_𝑖=1𝛾_𝑖𝛼=1,_𝑞=^𝑞−1_𝑗=1𝛽_𝑞𝑗𝛾,𝑞=2,…,4,₂𝛼₂+𝛾₃𝛼₃+𝛾₄𝛼₄=12,𝛾₂𝛼²₂+𝛾₃𝛼²₃+𝛾₄𝛼²₄=13,𝛾₂𝛼³₂+𝛾₃𝛼³₃+𝛾₄𝛼³₄=14,𝛾₃𝛼₂𝛽₃₂+𝛾₄𝛼₂𝛽₄₂+𝛽₄₃𝛼₃=16,𝛾₃𝛽₃₂𝛼₂𝛼₃+𝛾₄𝛽₄₂𝛼₂𝛼₄+𝛾₄𝛽₄₃𝛼₃𝛼₄=18,𝛾₃𝛼²₂𝛽₃₂+𝛾₄𝛼²₂𝛽₄₂+𝛾₄𝛼²₃𝛽₄₃=1,𝛾12₄𝛼₂𝛽₄₃𝛽₃₂=1.24 (2.31)
Thus, 𝐮_𝑛+1=𝐮_𝑛+𝑀_𝑛Δ𝑡_𝑛+12𝑀^′_𝑛Δ𝑡²_𝑛+16𝑀_𝑛^′′Δ𝑡³_𝑛+1𝑀24^{𝑛′′′}Δ𝑡⁴_𝑛+12𝑀²_𝑛Δ𝑡²_𝑛+12𝑀^′_𝑛𝑀_𝑛Δ𝑡³_𝑛+16𝑀_𝑛^′′𝑀_𝑛Δ𝑡⁴_𝑛+18𝑀^′_𝑛²Δ𝑡⁴_𝑛+16𝑀³_𝑛Δ𝑡³_𝑛+14𝑀^′_𝑛𝑀²_𝑛Δ𝑡⁴_𝑛+1𝑀24^4𝑛Δ𝑡⁴_𝑛𝐮_𝑛+𝑂Δ𝑡⁵_𝑛. (2.32)
We further assume ⁴_𝑖=1𝜏_𝑖𝜏=1,₂𝜇₁+𝜏₃𝜇₂+𝜏₄=12,𝜏₂𝜇²₁+𝜏₃𝜇²₂+𝜏₄=13,𝜏₂𝜇³₁+𝜏₃𝜇³₂+𝜏₄=14. (2.33)
The previous equation hence can be written as 𝐮_𝑛+1=𝐮_𝑛+𝜏₁+𝜏₂+𝜏₃+𝜏₄𝑀_𝑛Δ𝑡_𝑛+𝜏₂𝜇₁+𝜏₃𝜇₂+𝜏₄𝑀^′_𝑛Δ𝑡²_𝑛+12𝜏₂𝜇²₁+𝜏₃𝜇²₂+𝜏₄𝑀_𝑛^′′Δ𝑡³_𝑛+16𝜏₂𝜇³₁+𝜏₃𝜇³₂+𝜏₄𝑀_𝑛^′′′Δ𝑡⁴_𝑛+12𝜏₁+𝜏₂+𝜏₃+𝜏₄²𝑀²_𝑛Δ𝑡²_𝑛+𝜏₁+𝜏₂+𝜏₃+𝜏₄∗𝜏₂𝜇₁+𝜏₃𝜇₂+𝜏₄𝑀^′_𝑛𝑀_𝑛Δ𝑡³_𝑛+12𝜏₁+𝜏₂+𝜏₃+𝜏₄𝜏₂𝜇²₁+𝜏₃𝜇²₂+𝜏₄𝑀_𝑛^′′𝑀_𝑛Δ𝑡⁴_𝑛+12𝜏₂𝜇₁+𝜏₃𝜇₂+𝜏₄²𝑀^′_𝑛²Δ𝑡⁴_𝑛+16𝜏₁+𝜏₂+𝜏₃+𝜏₄³𝑀³_𝑛Δ𝑡³_𝑛+12𝜏₁+𝜏₂+𝜏₃+𝜏₄𝜏₂𝜇₁+𝜏₃𝜇₂+𝜏₄𝑀^′_𝑛𝑀²_𝑛Δ𝑡⁴_𝑛+1𝑀24^4𝑛Δ𝑡⁴_𝑛𝐮_𝑛+𝑂Δ𝑡⁵_𝑛. (2.34)
That is, 𝐮_𝑛+1=1+𝐽_𝑛Δ𝑡_𝑛+12𝐽_𝑛Δ𝑡_𝑛²+1𝐽3!^𝑛Δ𝑡_𝑛³+1𝐽4!^𝑛Δ𝑡_𝑛⁴𝐮_𝑛+𝑂Δ𝑡⁵_𝑛, (2.35) where 𝐽_𝑛=𝜏₁𝑀_𝑛+𝜏₂𝑀_𝑛+𝜇1+𝜏₃𝑀_𝑛+𝜇2+𝜏₄𝑀_𝑛+1,𝜏₂=2𝜇₂−112−𝜇₂+𝜇₁𝜇₁−1𝜇₁,𝜏₃=−−1+2𝜇₁12𝜇₂𝜇₁𝜇₂+𝜇₂−𝜇²₂−𝜇₁,𝜏₄=6𝜇₁𝜇₂−4𝜇₂+3−4𝜇₁12𝜇₂−1𝜇₁.−1 (2.36)
We now turn to examine the more general situation in which the reaction term is not zero, that is, 𝐛≠0. In this case, we have𝐮_𝑛+1=^𝑚_𝑖=0Δ𝑡_𝑛𝐽_𝑛^𝑖𝐮𝑖!^𝑛+𝐅𝐛_𝑛,𝐛_{𝑛+𝝁𝟏},…,𝐛_𝑛+1+𝑂Δ𝑡_𝑛^𝑖+1. (2.37)
Similar to our discussion in the previous section, the accumulation error for (2.37) is 𝐮_𝑛+1−𝐯_𝑛+1=^𝑚_𝑖=0Δ𝑡_𝑛𝐽_𝑛^𝑖×𝐮𝑖!^𝑛−𝐯_𝑛+𝑂Δ𝑡_𝑛^𝑖+1. (2.38)
We remark again that the amplification factor ∑^𝑚_𝑖=0((Δ𝑡_𝑛𝐽_𝑛)^𝑖/𝑖!) in (2.38) at the 𝑛th local step has the same form with (2.5) except 𝑀 being replaced by 𝐽_𝑛. 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 𝑎_max and 𝑎_min appeared in the optimal step size formulas then, respectively, represent the greatest and smallest real parts of the eigenvalues of the matrix 𝐽_𝑛. 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 condition 𝐮_𝟎 and initial step size Δ𝑡₀. (2)Set 𝐽_𝑛 in accordance with the stage number 𝑚, and apply 𝑚-stage Runge-Kutta schemes. (3)Use the Gerschgorin theorem to estimate the interval, [−𝑎^𝑛_max,−𝑎^𝑛_min] for negative real parts of the eigenvalues of 𝐽_𝑛.(4)Place Δ𝑡_𝑛+1, (i)𝑚=2:Δ𝑡_𝑛+1=2/(𝑎^𝑛_max+𝑎^𝑛_min), (ii)𝑚=3:Δ𝑡_𝑛+1=2.512745327/𝑎^𝑛_max, (iii)𝑚=4:Δ𝑡_𝑛+1=13172𝑎⁶_max+192𝑎⁵_max𝑎_min+276𝑎⁴_max𝑎²_min+448𝑎³_max𝑎³_min+276𝑎⁴_min𝑎²_max+192𝑎_max𝑎⁵_min+172𝑎⁶_min+3629𝑎⁶_max+30𝑎⁵_max𝑎_min+51𝑎⁴_max𝑎²_min+68𝑎³_max𝑎³_min+51𝑎⁴_min𝑎²_max+30𝑎_max𝑎⁵_min+29𝑎⁶_min^1/2𝑎³_max+3629𝑎⁶_max+30𝑎⁵_max𝑎_min+51𝑎⁴_max𝑎²_min+68𝑎³_max𝑎³_min+51𝑎⁴_min𝑎²_max+30𝑎_max𝑎⁵_min+29𝑎⁶_min^1/2𝑎_min𝑎²_max+3629𝑎⁶_max+30𝑎⁵_max+51𝑎⁴_max𝑎²_min+68𝑎³_max𝑎³_min+51𝑎⁴_min𝑎²_max+30𝑎_max𝑎⁵_min+29𝑎⁶_min^1/2𝑎_max𝑎²_min+3629𝑎⁶_max+30𝑎⁵_max𝑎_min+51𝑎⁴_max𝑎²_min+68𝑎³_max𝑎³_min+51𝑎⁴_min𝑎²_max+30𝑎_max𝑎⁵_min+29𝑎⁶_min^1/2𝑎³_min^1/3𝑎³_max+𝑎_min𝑎²_max+𝑎_max𝑎²_min+𝑎³_min−435𝑎⁴_max+10𝑎³_max𝑎_min+6𝑎²_max𝑎²_min+10𝑎_max𝑎³_min+5𝑎⁴_min𝑎³_max+𝑎_min𝑎²_max+𝑎_max𝑎²_min+𝑎³_min×172𝑎⁶_max+192𝑎⁵_max𝑎_min+276𝑎⁴_max𝑎²_min+448𝑎³_max𝑎³_min+276𝑎⁴_min𝑎²_max+192𝑎_max𝑎⁵_min+172𝑎⁶_min+3629𝑎⁶_max+30𝑎⁵_max𝑎_min+51𝑎⁴_max𝑎²_min+68𝑎³_max𝑎³_min+51𝑎⁴_min𝑎²_max+30𝑎_max𝑎⁵_min+29𝑎⁶_min^1/2𝑎𝑟³+3629𝑎⁶_max+30𝑎⁵_max𝑎_min+51𝑎⁴_max𝑎²_min+68𝑎³_max𝑎³_min+51𝑎⁴_min𝑎²_max+30𝑎_max𝑎⁵_min+29𝑎⁶_min^1/2𝑎_min𝑎²_max+3629𝑎⁶_max+30𝑎⁵_max𝑎_min+51𝑎⁴_max𝑎²_min+68𝑎³_max𝑎³_min+51𝑎⁴_min𝑎²_max+30𝑎_max𝑎⁵_min+29𝑎⁶_min^1/2𝑎_max𝑎²_min+3629𝑎⁶_max+30𝑎⁵_max𝑎_min+51𝑎⁴_max𝑎²_min+68𝑎³_max𝑎³_min+51𝑎⁴_min𝑎²_max+30𝑎_max𝑎⁵_min+29𝑎⁶_min^1/2𝑎³_min^1/3+43𝑎²_max+𝑎_max𝑎_min+𝑎²_min𝑎³_max+𝑎_min𝑎²_max+𝑎_max𝑎²_min+𝑎³_min, (2.39) and then let 𝑡_𝑛+1=𝑡_𝑛+Δ𝑡_𝑛+1. When an 𝑎^𝑛_max is equal to or close to zero, simply choose a reasonable step size for Δ𝑡_𝑛+1.
Note that the time step for 𝑚=4 can be solved by mathematical software such as MAPLE from the following equation: ⁴_𝑖=11𝑎𝑖!^𝑖max−𝑎^𝑖_minΔ𝑡^𝑖−1_𝑛+1=0. (2.40)

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, 𝑎^𝑛_min=0. In this way, only 𝑎^𝑛_max are needed in forming optimal time step formulas and time steps Δ𝑡_𝑛+1 in step 4 are switched as(i)𝑚=2:Δ𝑡_𝑛+1=2/𝑎^𝑛_max, (ii)𝑚=3:Δ𝑡_𝑛+1=2.512745327/𝑎^𝑛_max, (iii)𝑚=4:Δ𝑡_𝑛+1=2.785293563/𝑎^𝑛_max,

where 𝑎^𝑛_max is defined as before in step 2. We remark that when Δ𝑡_𝑛+1 is specified as above values, 𝑝_𝑚(𝑎,Δ𝑡) 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, 𝑝_𝑚(𝑎,Δ𝑡) 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 the 𝑚-stage Runge-kutta method as some equivalent 𝑚-stage formulas while the order of accuracy remains invariant. For example, 𝑚=3, we have𝐮_𝑛^(𝟎)=𝐮_𝑛+Δ𝑡_𝑛3𝐽_𝑛𝐮_𝑛,𝐮_𝑛^(𝟏)=𝐮_𝑛+Δ𝑡_𝑛2𝐽_𝑛𝐮_𝑛^(𝟎),𝐮_𝑛+1=𝐮_𝑛+Δ𝑡𝐽_𝑛𝐮_𝑛^(𝟏). (2.41)

For 𝑚=4, 𝐮_𝑛^(𝟎)=𝐮_𝑛+Δ𝑡_𝑛4𝐽_𝑛𝐮_𝑛,𝐮_𝑛^(𝟏)=𝐮_𝑛+Δ𝑡_𝑛3𝐽_𝑛𝐮_𝑛^(𝟎),𝐮_𝑛^(𝟐)=𝐮_𝑛+Δ𝑡_𝑛2𝐽_𝑛𝐮_𝑛^(𝟏),𝐮_𝑛+1=𝐮_𝑛+Δ𝑡_𝑛𝐽_𝑛𝐮_𝑛^(𝟐). (2.42)

Of course, 𝐽_𝑛 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 standard 𝑚-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𝑢_𝑡=𝑢_𝑥𝑥−(𝑡𝑥𝑢)_𝑥+𝑓,𝑥∈(0,2𝜋),𝑡>0, (3.1) where 𝑢=𝑢(𝑥,𝑡) and 𝑓=(𝑡𝑥−1)𝑒^−𝑡cos(𝑥−𝑡)+𝑡𝑒^−𝑡sin(𝑥−𝑡). The three terms on the right-hand side represent diffusion, chemotaxis, and reaction, respectively. Assume the periodic boundary condition: 𝑢(𝑜,𝑡)=𝑢(2𝜋,𝑡),𝑡≥0 and the initial condition: 𝑢(𝑥,0)=sin𝑥. The analytic solution of (3.1) is𝑢(𝑥,𝑡)=𝑒^−𝑡sin(𝑥−𝑡). (3.2)

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 𝑚-stage Runge-Kutta method with the adaptive optimal time stepping scheme as the time integration strategy. By denoting ℎ=2𝜋/𝑁, 𝑥_𝑖=𝑖ℎ, 𝑖=0,1,…,𝑁, and the numerical solution as 𝑢_𝑖, the resulting ODE system is 𝑑𝐮(𝑡)(𝐭),𝑑𝑡=𝑀𝐮+𝐟 (3.3) where 𝐮=𝑢₀,𝑢₁,…,𝑢_𝑁^𝑇,𝐟(𝐭)=𝑓𝑥₀,𝑡,𝑓𝑥₁,𝑡,…,𝑓𝑥_𝑁,𝑡^𝑇,𝑀𝑡_𝑛=−2ℎ²−𝑥₀𝑡_𝑛−1ℎ1ℎ²1ℎ²−𝑡_𝑛−101ℎ²+𝑥₀𝑡_𝑛−1ℎ−2ℎ²−𝑥₁𝑡_𝑛−1ℎ1ℎ²100ℎ²+𝑥₁𝑡_𝑛−1ℎ−2ℎ²−𝑥₂𝑡_𝑛−1ℎ1ℎ²0100⋱⋱⋱⋱⋱⋱⋱⋱⋱ℎ²01ℎ²+𝑥_𝑁−1𝑡_𝑛−1ℎ−2ℎ²−𝑥_𝑁𝑡_𝑛−1ℎ. (3.4) Generally, the actual eigenvalues of the spatial matrix 𝑀(𝑡) 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 𝐾_𝑖∶=𝜇∈𝜇−𝑎_𝑖𝑖≤^𝑛_{𝑘=1,𝑘≠𝑖}𝑎_𝑖𝑘 (3.5) contains all eigenvalues of the 𝑛×𝑛 matrix 𝐴=[𝑎_𝑖𝑘].

The theorem states that if 𝐴 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 −𝑎_max≤Re(𝜆)<0, (3.6) where 𝑎_max=4Δ𝑥²+𝑡_𝑛4𝜋.Δ𝑥+1 (3.7)

Thus, for 𝑚-stages, where 𝑚=2,…,4, the optimal time steps at the current step are𝑚=2:Δ𝑡_opt=24/Δ𝑥²+𝑡_𝑛)),((4𝜋/Δ𝑥+1𝑚=3:Δ𝑡_opt=2.5127453274/Δ𝑥²+𝑡_𝑛)),((4𝜋/Δ𝑥+1𝑚=4:Δ𝑡_opt=2.7852935634/Δ𝑥²+𝑡_𝑛)).((4𝜋/Δ𝑥+1 (3.8)

Figure 3 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 𝑡=5.2, but the optimal time-stepping 𝑚-stage Runge-Kutta methods still converge to the analytic solution.

Figure 3 

Relative error between numerical simulations and the analytic solution.

Tables 1, 2, and 3 demonstrate the number of iteration for our scheme (Optimal 𝑚-stage RK) with that of two other methods: standard 𝑚-stage Runge-Kutta method with fixed step size (Std 𝑚-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 𝑚-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.

Table 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.

Table 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.

Table 4 

Comparison of CPU time among the optimal two-stage RK, RKF 23, and standard two-stage RK with Δ𝑡=0.001.

Table 5 

Comparison of CPU time among the optimal three-stage RK, RKF 34, and standard three-stage RK with Δ𝑡=0.0005 and Δ𝑡=0.001, respectively.

Table 6 

Comparison of CPU time among the optimal four-stage RK, RKF 45, and standard four-stage RK with Δ𝑡=0.001.

Tables 4–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 Sections 2.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Δ𝑡=2.785293563𝑎_max. (3.9)

To get the “worst" case, set 𝑎_max=𝑎_min. Then Δ𝑡=1.596071638𝑎_max. (3.10)

Now, let us assign Δ𝑡 to be the average of (3.9) and (3.10). We obtainΔ𝑡=2.1906826005𝑎_max. (3.11)

Table 7 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.

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𝜕𝑢𝐷𝜕𝑡=∇⋅_𝑢∇𝑢−𝛼𝑢(1+𝑣)²∇𝑣,𝜕𝑣𝑢𝜕𝑡=Δ𝑣+𝜔²𝜇+𝑢², (3.12) where 𝑢,𝑣, and 𝜔 represent cell density, concentration of chemotactic substance, and succinate, respectively. Simulations were done with zero flux boundary condition for both 𝑢 and 𝑣. Initial conditions are set to be zero for 𝑣 and the random perturbation at 𝑢₀=1 for the initial condition of 𝑢. The noise amplitude is of order 10⁻¹. They were created by the 𝐶⁺⁺ random number generator. The initial, boundary conditions for 𝑢 and 𝑣, 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 Figure 4. Simulations did capture the characteristics as described in (2.28) qualitatively. The aggregates first are divided into small regions and then cojoin with others.

(a)

(b)

(c)

(d)

Figure 4 

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.

4. Conclusions and Discussions

We investigate an optimal adaptive time stepping strategy for fully explicit 𝑚-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 𝑚-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.

Acknowledgment

This work was supported by National Science Council under Grant no. NSC 97-2115-M-126-003 in Taiwan.