# A Strongly A-Stable Time Integration Method for Solving the Nonlinear Reaction-Diffusion Equation

Authors: Wenyuan Liao
Journal: Abstract and Applied Analysis
 (2015)

---

## Abstract

The semidiscrete ordinary differential equation (ODE) system resulting from compact higher-order finite difference spatial discretization of a nonlinear parabolic partial differential equation, for instance, the reaction-diffusion equation, is highly stiff. Therefore numerical time integration methods with stiff stability such as implicit Runge-Kutta methods and implicit multistep methods are required to solve the large-scale stiff ODE system. However those methods are computationally expensive, especially for nonlinear cases. Rosenbrock method is efficient since it is iteration-free; however it suffers from order reduction when it is used for nonlinear parabolic partial differential equation. In this work we construct a new fourth-order Rosenbrock method to solve the nonlinear parabolic partial differential equation supplemented with Dirichlet or Neumann boundary condition. We successfully resolved the phenomena of order reduction, so the new method is fourth-order in time when it is used for nonlinear parabolic partial differential equations. Moreover, it has been shown that the Rosenbrock method is strongly A-stable hence suitable for the stiff ODE system obtained from compact finite difference discretization of the nonlinear parabolic partial differential equation. Several numerical experiments have been conducted to demonstrate the efficiency, stability, and accuracy of the new method.


---

## Body

## 1. Introduction

Let us consider the following parabolic partial differential equation:(1)
u
t
=
D
u
x
x
+
f
u
,
x
,
t
,
x
,
t
∈
a
,
b
×
0
,
T
,with the initial condition:(2)
u
x
,
0
=
u
0
x
,
x
∈
a
,
b
,where D is a positive constant describing the diffusion property and f
(
u
,
x
,
t
) is a function representing the reaction term, which is nonlinear on u. The unknown function u represents, depending on the applications, variables such as mass concentration in chemical reaction process, temperature in heat conduction, neutron flux in nuclear reactors, and population density in population dynamics. On the boundary, either Dirichlet condition(3)
u
a
,
t
=
g
1
t
,
u
b
,
t
=
g
2
t
,
0
≤
t
≤
Tor Neumann condition(4)
u
x
a
,
t
=
g
3
t
,
u
x
b
,
t
=
g
4
t
,
0
≤
t
≤
Tis specified, where g
1, g
2, g
3, and g
4 are sufficiently smooth functions. Here in this paper we restrict our attention on Dirichlet and Neumann boundary conditions, while the developed techniques can be easily extended to Robin boundary condition.Efficient and accurate numerical methods for solving (1) had attracted great attentions from scientists and engineers, as for many application problems in science and engineering, it is preferable to use high-order compact numerical algorithms to compute accurate solutions. In the past several decades a great deal of work has been done in the development of efficient, accurate, and robust numerical algorithm for solving such problem. For more details, the reader is referred to [1–4].Since both temporal and spatial derivatives are involved in the equation, we discuss the numerical treatments in time and space separately. Here we first apply the high-order compact finite difference approximation to the spatial derivative, so a semidiscrete ODE system is obtained, which is then solved by a fourth-order Rosenbrock method that will be discussed later.Recently there have been attempts to develop high-order compact scheme for the spatial derivative. In [5], a three-point combined compact difference scheme was proposed to approximate the first and second derivatives for problem with periodic boundary condition. The resulting scheme has up to sixth-order accuracy at all grid points including the boundary nodes for periodic boundaries; however it is only fourth-order accurate for nonperiodic boundary condition.Because the semidiscrete ODE system obtained from spatial discretization, such as method of lines, of the nonlinear parabolic partial differential equation is highly stiff, the choices of time integration methods are limited to implicit methods only. Explicit algorithm is efficient in a single time step but suffers from strict step size restriction, which makes it less efficient. Implicit method, on the other side, is less efficient in a single step but the unconditional stability allows the use of larger time step; hence the overall computational efficiency can be significantly improved. One issue, however, is that the iteration is usually slow when large step sizeΔ
t is used. Also, due to the stiffness of the ODE system, only Newton-type iterative methods are applicable to solve the nonlinear algebraic system. Furthermore, strong A-stability or L-stability of the time integration method is necessary for error damping. A great deal of work has been done in the development of efficient time stepping methods for the stiff ODE system. In [6] explicit exponential Rosenbrock methods of order five have been constructed to solve the large-scale stiff ODE system. Through the derivation of stiff order condition, new pairs of embedded methods of higher-order can be obtained. Similarly, fifth-order explicit exponential Runge-Kutta methods were constructed to efficiently integrate the semilinear stiff problems in [7]. The authors have also shown that there does not exist an explicit exponential Runge-Kutta method of order 5 with less than or equal to 6 stages; therefore the resultant methods are 8-stage methods. In [8] a fourth-order time stepping method, which is a modification of the exponential time-differencing fourth-order Runge-Kutta method, has been developed for stiff ODEs. These methods are efficient and accurate. However, A-stability of the time stepping method is not sufficient for highly stiff problem. To overcome these difficulties, it is desirable to construct new algorithms with strong A-stability or L-stability that are free of solving nonlinear equations. It turns out that the Rosenbrock method, which was firstly reported by Rosenbrock [9] and then improved by Haines [10], responded to these issues with considerably satisfying and promising results.The objective herein is to develop a strongly A-stable Rosenbrock method to solve the semidiscrete stiff ODEs resulting from compact high-order finite difference approximation of a semilinear parabolic partial differential equation. The rest of the paper is organized as the following. In Section2, we discretize the spatial derivatives of the semilinear parabolic partial differential equation using a fourth-order compact finite difference scheme, which is then combined with a newly proposed compact fourth-order boundary condition treatment to form the semidiscrete ODE system. In Section 3 we focus on the development of a fourth-order strongly A-stable Rosenbrock method and the stability analysis. Several numerical examples are used to demonstrate the accuracy and efficiency of the new algorithm in Section 4, which is followed by conclusions and possible future work.


## 2. Compact Fourth-Order Spatial Discretization

For the sake of simplicity, we assume that the 1D spatial domainΩ
=
[
a
,
b
] is divided into M subintervals with equal length h
=
(
b
-
a
)
/
M. Let x
i
=
a
+
i
·
h, i
=
0,1
,
…
,
M, be the grid points. A variety of compact high-order discretizations can be utilized to approximate the second derivative u
x
x in (1).Here we introduce a compact finite difference scheme to approximateu
x
x, such that the resulting semidiscrete ODE system is an accurate and compact approximation to the original semilinear parabolic partial differential equation. This operator-approximation based method has been widely used to solve various multidimensional problems. We first define the central finite difference operator δ
x
2 as(5)
δ
x
2
u
i
=
u
i
+
1
-
2
u
i
+
u
i
-
1
;then δ
x
2
/
h
2 gives second-order accurate approximation to u
x
x. Using Taylor series to expand all terms on the right-hand side of (5), under the assumption that u
(
x
) is sufficiently smooth, we have(6)
δ
x
2
h
2
u
i
=
u
x
x
x
i
+
h
2
12
u
x
x
x
x
x
i
+
O
h
4
,
1
≤
i
≤
M
-
1
.To improve the above finite difference approximation to fourth-order accurate, one just needs to eliminate the second-order error term. Applying δ
x
2
/
12 to both sides of (6), we have(7)
δ
x
2
h
2
δ
x
2
12
u
i
=
u
x
x
x
x
x
i
+
O
h
4
,
1
≤
i
≤
M
-
1
.Combining (6) with (7), neglecting O
(
h
4
), we obtain the following fourth-order accurate approximation to u
x
x at node x
i:(8)
u
x
x
x
i
≈
1
h
2
δ
x
2
1
-
δ
x
2
12
u
i
,
1
≤
i
≤
M
-
1
.The drawback is that a five-point stencil is required; therefore the compactness is destroyed so the method becomes less efficient. Further investigation shows that the difference between 1
-
δ
x
2
/
12 and 1
+
δ
x
2
/
12
-
1 is O
(
h
4
), so a natural way is to approximate u
x
x
(
x
i
) as δ
x
2
1
+
δ
x
2
/
12
-
1
u
i, which is fourth-order accurate and compact.Applying the fourth-order Padé approximation tou
x
x in (1), we obtain the following ODE system(9)
u
i
′
t
=
D
h
2
δ
x
2
1
+
δ
x
2
/
12
u
i
t
+
f
u
i
t
,
x
i
,
t
,
hhhhhhhhhhhhhhhhiihhhhh
1
≤
i
≤
M
-
1
,which is a fourth-order accurate approximation (in space) to the original semilinear parabolic partial differential equation defined in (1).However, the above algorithm is difficult to implement, so we multiply1
+
δ
x
2
/
12 to both sides to obtain the following implicit ODE system:(10)
1
+
δ
x
2
12
u
i
′
t
=
D
h
2
δ
x
2
u
i
t
+
1
+
δ
x
2
12
f
u
i
t
,
x
i
,
t
,
hhhhhhhhhhhhhhhhhhhhh
1
≤
i
≤
M
-
1
,
0
<
t
≤
T
,which can be written in vector form as(11)
A
U
′
t
=
F
U
,
X
,
t
,
0
<
t
≤
T
,where U
(
t
)
=
(
u
1
(
t
)
,
u
2
(
t
)
,
…
,
u
M
-
1
(
t
)
) is the discrete solution of (1) at time t, with u
i
(
t
)
=
u
(
x
i
,
t
), A is an (
M
-
1
)
×
(
M
-
1
) tridiagonal matrix, and F is a vector-valued function defined through (10). To complete the ODE system, we need the boundary conditions at x
=
x
0 and x
=
x
M, which can be derived from the original boundary conditions defined in (3) or (4).First, if the Dirichlet boundary condition (3) is specified, one can add the following two ODEs to (9):(12)
u
0
′
t
=
g
1
′
t
,
u
M
′
t
=
g
2
′
t
.Consequently the matrix is modified as(13)
A
=
1
0
0
0
⋯
0
0
0
1
12
5
6
1
12
0
⋯
0
0
0
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
0
0
0
0
⋯
1
12
5
6
1
12
0
0
0
0
⋯
⋯
⋯
1
,while the vector-valued function F
(
U
,
X
,
t
), after modifications, is defined as(14)
h
F
0
=
g
1
′
t
,
h
F
i
=
D
h
2
δ
x
2
u
i
t
+
1
+
δ
x
2
12
f
u
i
t
,
x
i
,
t
,
hhhhhhhhhhhhhhhhhhhhhh
1
≤
i
≤
M
-
1
,
F
M
=
g
2
′
t
.Alternatively, we can incorporate the boundary condition by replacingu
0
(
t
) and u
M
(
t
) in (10) with g
1
(
t
) and g
2
(
t
), respectively, so the ODE system equation (11) has only M
-
2 equations.As one can imagine, the situation is more complicated when the Neumann boundary condition (4) is specified. To complete the ODE system and maintain the higher-order overall accuracy, a compact fourth-order approximation of the Neumann boundary condition is needed. Let us use the boundary condition at x
=
a as the example to demonstrate the idea of the new algorithm.Unlike the Dirichlet boundary condition, which specifies the solutionu on the boundary point explicitly, the Neumann boundary condition defines u
x at the boundary points; thus u
0
(
t
) and u
M
(
t
) need to be calculated along with solution at the interior grid points. Consequently, the range for i in (10) should be changed to 0
≤
i
≤
M, so A is an (
M
+
1
)
×
(
M
+
1
) matrix. To approximate the derivative at x
0, we introduce a ghost point x
-
1
=
a
-
h and assume that (1) holds and the solution u is sufficiently smooth on the extended domain [
a
-
h
,
b
]. Let u
-
1
(
t
) denote the solution at x
-
1
=
a
-
h and then apply the second-order central finite difference approximation to u
x
(
a
,
t
),(15)
u
1
t
-
u
-
1
t
2
h
=
u
x
a
,
t
+
h
2
6
u
x
x
x
a
,
t
+
O
h
4
.Taking partial derivative with respect tox on both sides of (1), we have(16)
u
x
x
x
=
1
D
u
x
t
-
f
x
-
f
u
·
u
x
.Letting x
→
a in (16) and then applying the Neumann boundary condition (4), we obtain(17)
u
x
x
x
a
=
1
D
g
3
′
t
-
f
x
u
0
t
,
a
,
t
-
f
u
u
0
t
,
a
,
t
·
g
3
t
.Combining (15) with (17), we obtain the following fourth-order compact approximation for u
-
1
(
t
):(18)
u
-
1
t
=
u
1
t
-
2
h
g
3
t
-
h
3
3
D
g
3
′
t
f
u
u
0
t
,
a
,
t
·
g
3
t
g
3
′
t
-
f
x
u
0
t
,
a
,
t
hhhhhhhhhhhhhhhhhh
-
g
3
′
t
f
u
u
0
t
,
a
,
t
·
g
3
t
,which involves u
0
(
t
) and u
1
(
t
) only, so the compact structure is preserved.Similarly, the fourth-order compact approximation foru
M
+
1
(
t
) can be derived as(19)
u
M
+
1
t
=
u
M
-
1
t
+
2
h
g
4
t
+
h
3
3
D
g
4
′
t
-
f
x
u
M
t
,
b
,
t
hhhhhhh
-
f
u
u
M
t
,
b
,
t
·
g
4
t

g
4
′
t
.Now the matrixA involves t and U, and the first and last rows are modified as(20)
A
0,0
⋯
M
=
10
12
+
h
3
36
D
f
x
u
u
0
,
a
,
t
hhhhhhhhhihh
+
f
u
u
u
0
,
a
,
t
·
g
3
t
f
x
u
u
0
,
a
,
t
,
1
6
,
…
,
0
h
3
36
D
,
A
M
,
0
⋯
M
=
0
,
…
,
1
6
,
10
12
-
h
3
36
D
f
x
u
1
f
x
u
(
u
M
,
b
,
t
)
hhhhhhhhhhhhhiihhhh
+
f
u
u
u
M
,
b
,
t
·
g
4
t
f
x
u
1
h
3
36
D
.Consequently the first and last components ofF are(21)
F
0
=
D
h
2
2
u
1
-
2
u
0
-
2
h
g
3
+
1
12
10
f
u
0
,
a
,
t
+
f
u
1
,
a
+
h
,
t
-
h
3
g
3
′
-
f
x
u
0
,
a
,
t
-
f
u
u
0
,
a
,
t
·
g
3
+
h
6
g
3
′
+
h
3
36
D
g
3
′
′
-
f
x
t
u
0
,
a
,
t
-
f
u
u
0
,
a
,
t
·
g
3
′
hhhhhhhh
-
f
u
t
u
0
,
a
,
t
·
g
3
g
3
′′
-
f
x
t
u
0
,
a
,
t
+
1
12
f
u
1
-
2
h
g
3
-
h
3
3
D
hhhhhhhih
×
g
3
′
-
f
x
u
0
,
a
,
t
-
f
u
u
0
,
a
,
t
·
g
3
,
hhhhhhhih
a
-
h
,
t
h
3
3
D
,
F
M
=
2
D
h
2
u
M
-
1
-
u
M
+
h
g
4
+
1
12
10
f
u
M
,
b
,
t
+
f
u
M
-
1
,
b
-
h
,
t
-
h
3
(
g
4
′
-
f
x
(
u
M
,
b
,
t
)
-
f
u
(
u
M
,
b
,
t
)
·
g
4
)
-
h
6
g
4
′
-
h
3
36
D
g
4
′
′
-
f
x
t
u
M
,
b
,
t
-
f
u
u
M
,
b
,
t
·
g
4
′
hhhhhhhh
-
f
u
t
u
M
,
b
,
t
·
g
4
g
4
′′
-
f
x
t
u
M
,
b
,
t
+
1
12
f
u
M
-
1
+
2
h
g
4
+
h
3
3
D
hhhhhhhih
×
g
4
′
-
f
x
u
M
,
b
,
t
-
f
u
u
M
,
b
,
t
·
g
4
,
hhhhhhhih
b
+
h
,
t
h
3
3
D
.Finally, the ODE system is written in the form of A
(
t
,
U
)
U
′
=
F
(
U
,
t
). Apparently, the matrix A preserves tridiagonal structure, but it depends on t and U; hence the development of Rosenbrock method becomes difficult. Fortunately, we notice that t and U are involved in two entries: A
0,0 and A
M
,
M only. Further investigation shows that the extra terms in (20) are(22)
h
3
36
D
f
x
u
u
0
,
a
,
t
+
f
u
u
u
0
,
a
,
t
·
g
1
t
,
(23)
-
h
3
36
D
f
x
u
u
M
,
b
,
t
+
f
u
u
u
M
,
b
,
t
·
g
4
t
,respectively. We can eliminate these two extra terms by incorporating them into vector F; therefore F
0 and F
M are modified as(24)
F
¯
0
=
F
0
-
h
3
36
D
f
x
u
u
0
,
a
,
t
+
f
u
u
u
0
,
a
,
t
·
g
3
t
·
u
0
′
t
,
(25)
F
¯
M
=
F
M
+
h
3
36
D
f
x
u
u
M
,
b
,
t
+
f
u
u
u
M
,
b
,
t
·
g
4
t
·
u
M
′
t
.Using (1), we have(26)
u
0
′
t
=
2
D
h
2
u
1
t
-
u
0
t
-
h
g
3
t
+
f
u
0
t
,
a
,
t
+
O
h
,
(27)
u
M
′
t
=
2
D
h
2
u
M
-
1
t
-
u
M
t
+
h
g
4
t
+
f
u
M
t
,
b
,
t
+
O
h
.Inserting (26) into (24) and then ignoring the fourth-order error term O
(
h
4
), we obtain(28)
F
¯
0
=
F
0
-
h
3
36
D
f
x
u
u
0
,
a
,
t
+
f
u
u
u
0
,
a
,
t
·
g
3
t
×
2
D
h
2
u
1
-
u
0
-
h
g
3
t
+
f
u
0
,
a
,
t
.Similarly, inserting (27) into (25) and then ignoring the fourth-order error term, we obtain(29)
F
¯
M
=
F
M
-
h
3
36
D
f
x
u
u
M
,
b
,
t
+
f
u
u
u
M
,
b
,
t
·
g
4
t
×
2
D
h
2
u
M
-
1
t
-
u
M
t
+
h
g
4
t
hhhhh
+
f
u
M
t
,
b
,
t
2
D
h
2
.We then obtain the closed ODE system A
U
′
(
t
)
=
F
(
U
,
t
), where the vector-valued function is given as F
=
(
F
¯
0
,
F
1
,
…
,
F
M
-
1
,
F
¯
M
), and A is a nonsingular (
M
+
1
)
×
(
M
+
1
) constant matrix given as(30)
A
=
5
6
1
12
0
0
⋯
0
0
0
1
12
5
6
1
12
0
⋯
0
0
0
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
0
0
0
0
⋯
1
12
5
6
1
12
0
0
0
0
⋯
⋯
1
12
5
6
.If the Robin boundary condition is specified, a similar numerical technique can be used to derive the semidiscrete ODE system.Here we mention, without theoretical proof, that the resulting semidiscrete ODE system (9) is a fourth-order accurate approximation to the original semilinear parabolic partial differential equation given in (1), supplemented with either Dirichlet or Neumann boundary conditions. Interested readers can find a similar theorem and proof in [11].


## 3. Fourth-Order Strongly A-Stable Rosenbrock Method

Various numerical methods can be used to solve the ODE system in (11). However, due to the stiffness of the problem, only stiffly stable methods are applicable; thus the choices are limited to the subclass of implicit methods such as implicit linear multistep methods and implicit Runge-Kutta methods. It is known that A-stability is necessary for stiff problem, and in general strong A-stable or even L-stable methods are preferred. The A-stability was firstly introduced and defined by Dahlquist [12] as the following.Definition 1.
A numerical method is called A-stable if there is no restriction on the step size, when it is applied to solve the test equationy
′
=
λ
y, where Re
(
λ
)
<
0.For a single-step method such as Runge-Kutta method that can be written asy
n
+
1
=
R
(
y
n
), the A-stability is equivalent to the condition R
(
z
)
≤
1 for any z
∈
C
-
1, where R
(
z
) is called the stability function of the method. Although successfully used in various applications, an A-stable linear multistep method has the highest order of 2. In fact, the second-order A-stable linear multistep method with optimal error constant is the Trapezium rule [12]. Further, Gourlay [13] pointed out that an A-stable method is necessary but not sufficient for very stiff system as it has the incorrect damping rate. For example, the widely used Trapezium rule has stability function R
(
z
)
=
1
+
z
/
1
-
z satisfying R
(
z
)
<
1 for any z
∈
C
-
1, but its damping rate converges to -
1 when z
→
-
∞. To overcome this difficulty, strong A-stability was introduced.Definition 2.
A numerical method is called strongly A-stable if it is A-stable andR
(
-
∞
)
<
1.It has been shown that a numerical method with strong A-stability is effective in damping numerical oscillations for highly stiff system. For more details about the description and comparison of A-stability and strong A-stability, the readers are referred to [14].Implicit Runge-Kutta method is usually unconditionally stable but suffers from the issue of high computational complexity, especially for nonlinear ODE system. For example, during each time step, an algebraic system withs
×
M unknown variables needs to be solved, if an s-stage implicit Runge-Kutta method is used to solve an ODE system with M equations. Therefore, fully implicit Runge-Kutta methods are too computationally expensive to be useful for large-scale problems. In the past several decades, efforts have been made to reduce the computational cost, which results in various modified implicit Runge-Kutta methods, such as diagonally implicit Runge-Kutta method, singly diagonally implicit Runge-Kutta method, explicit-implicit Runge-Kutta method, to name a few. For more details of these methods, the reader is referred to [15–19].To completely avoid solving a nonlinear algebraic system, Rosenbrock method, which is a special class of Runge-Kutta method, had been proposed. Since the spatial discretization is fourth-order, our aim herein is to develop a strongly A-stable fourth-order Rosenbrock method for solving the ODE system (11), so the new algorithm is fourth-order accurate in both temporal and spatial dimensions.


### 3.1. Rosenbrock Method for Scalar Equation

We first derive the Rosenbrock method based on an autonomous scalar equationy
′
=
f
(
y
), for which the initial condition is given as y
(
t
0
)
=
y
0. Nonautonomous equations can be converted to autonomous form by adding an extra equation to the system. Some previous research [20] suggested that it is unlikely to find a 3-stage fourth-order Rosenbrock method with strong A-stability or L-stability, so herein we focus on the development of a 4-stage Rosenbrock method. Suppose the numerical solution at time t
n is known as y
n, the 4-stage Rosenbrock method calculates the numerical solution at t
n
+
1 as(31)
y
n
+
1
=
y
n
+
b
1
k
1
+
b
2
k
2
+
b
3
k
3
+
b
4
k
4
,
(32)
k
i
=
Δ
t
f
y
n
+
∑
j
=
1
i
-
1
‍
α
i
j
k
j
+
Δ
t
J
y
n
∑
j
=
1
i
‍
γ
i
j
k
j
,
hhhhhhhhhhhhhhhhihhhh
i
=
1,2
,
3,4
,where b
i
,
α
i
j, and γ
i
j are coefficients to be determined and J
(
y
n
)
=
f
y
(
y
n
).To extend the above algorithm to the nonautonomous problemsy
′
=
f
(
y
,
t
), we first convert it to autonomous form by adding a new equation t
′
=
1 and then apply the algorithm (31) to the augmented system. Note that the components corresponding to the last variable can be computed explicitly; thus we can derive the modified algorithm as the following:(33)
y
n
+
1
=
y
n
+
b
1
k
1
+
b
2
k
2
+
b
3
k
3
+
b
4
k
4
,
hh
k
i
=
Δ
t
f
y
n
+
∑
j
=
1
i
-
1
‍
α
i
j
k
j
,
t
n
+
α
i
Δ
t
hhhhh
+
γ
i
Δ
t
2
f
t
y
n
,
t
n
+
Δ
t
f
y
y
n
,
t
n
∑
j
=
1
i
‍
γ
i
j
k
j
,
hhhhhhhhhhhhhhhhhhhhhhhhhh
i
=
1,2
,
3,4
,where(34)
α
i
=
∑
j
=
1
i
-
1
‍
α
i
j
,
γ
i
=
∑
j
=
1
i
‍
γ
i
j
.A Rosenbrock method of orderp is obtained through choosing coefficients in (31)-(32) so that the local error satisfies y
(
t
n
+
Δ
t
)
-
y
n
+
1
=
O
(
Δ
t
p
+
1
). This can be done either by solving the so-calledButcher Series [14] or by straightforward differentiation. Here we derive the order conditions for the fourth-order Rosenbrock method in a different way using Taylor series. First, both y
(
t
n
+
Δ
t
) and k
i are expanded as Taylor series so that the difference y
(
t
n
+
Δ
t
)
-
y
n
+
1 can be expressed as a Taylor series and its coefficients of the terms up to O
(
Δ
t
4
) are set to 0, which results in a set of equations involving these coefficients.Giveny
(
t
n
)
=
y
n, the Taylor series of y
(
t
n
+
Δ
t
) at t
n is expanded as(35)
y
t
n
+
Δ
t
=
y
n
+
Δ
t
f
n
+
Δ
t
2
2
J
n
f
n
+
Δ
t
3
6
J
n
2
f
n
+
J
n
′
f
n
2
+
Δ
t
4
24
J
n
′
′
f
n
3
+
2
J
n
′
J
n
f
n
2
+
2
J
n
J
n
′
f
n
2
+
J
n
3
f
n
+
O
Δ
t
5
,where f
n
=
f
(
y
n
), J
n
=
∂
f
/
∂
y
(
y
n
), J
n
′
=
∂
2
f
/
∂
y
2
y
n, and J
n
′′
=
∂
3
f
/
∂
y
3
y
n. For the sake of simplicity, let β
i
j
=
α
i
j
+
γ
i
j and α
i
i
=
0.Lettingi
=
1 in (32) we have(36)
k
1
=
Δ
t
1
-
Δ
t
γ
11
J
n
-
1
f
n
=
Δ
t
f
n
+
Δ
t
2
γ
11
J
n
f
n
+
Δ
t
3
γ
11
2
J
n
2
f
n
+
Δ
t
4
γ
11
3
J
n
3
f
n
+
O
Δ
t
5
.Similarly, lettingi
=
2 in (32) we have(37)
k
2
=
Δ
t
1
-
Δ
t
γ
22
J
n
-
1
f
y
n
+
α
21
k
1
+
γ
21
J
n
k
1
.Combining it with the Taylor series ofk
1 in (36), we obtain(38)
k
2
=
Δ
t
f
n
+
Δ
t
2
α
21
+
γ
21
+
γ
22
J
n
f
n
+
Δ
t
3
γ
11
+
γ
22
α
21
+
γ
21
+
γ
22
2
J
n
2
f
n
+
Δ
t
3
2
α
21
2
J
n
′
f
n
2
+
Δ
t
4
γ
11
γ
22
+
γ
11
2
+
γ
22
2
α
21
+
γ
21
+
γ
22
3
J
n
3
f
n
+
Δ
t
4
γ
11
α
21
2
J
n
′
J
n
f
n
2
+
Δ
t
4
2
γ
22
α
21
2
J
n
J
n
′
f
n
2
+
Δ
t
4
6
α
21
3
J
n
′
′
f
n
3
+
O
Δ
t
5
.Following the same way, we have(39)
k
3
=
Δ
t
f
n
+
Δ
t
2
(
β
31
+
β
32
+
γ
33
)
J
n
f
n
+
Δ
t
3
∑
j
=
1
3
‍
α
3
j
+
γ
3
j
·
∑
l
=
1
j
‍
α
j
l
+
γ
j
l
J
n
2
f
n
+
Δ
t
3
2
α
3
2
J
n
′
f
n
2
+
Δ
t
4
∑
i
=
1
3
‍
α
3
i
+
γ
3
i
·
∑
j
=
1
i
‍
∑
l
=
1
j
‍

α
i
j
+
γ
i
j
hhhhhhhhhhhhhhhhhhhhhhh
·
∑
l
=
1
j
‍
α
j
l
+
γ
j
l
∑
l
=
1
j
‍

α
i
j
+
γ
i
j
α
3
i
+
γ
3
i
·
∑
j
=
1
i
‍
∑
i
=
1
3
‍
α
3
i
+
γ
3
i
·
∑
j
=
1
i
‍
∑
l
=
1
j
‍

α
i
j
+
γ
i
j
J
n
3
f
n
+
Δ
t
4
2
[
α
2
2
β
32
+
γ
33
α
3
2
]
J
n
J
n
′
f
n
2
+
Δ
t
4
α
3
β
21
+
γ
22
α
32
+
γ
11
α
31
J
n
′
J
n
f
n
2
+
Δ
t
4
6
α
31
+
α
32
3
J
n
′
′
f
n
3
+
O
Δ
t
5
,
k
4
=
Δ
t
f
n
+
Δ
t
2
β
41
+
β
42
+
β
42
+
γ
44
J
n
f
n
+
Δ
t
3
∑
j
=
1
4
‍
α
4
j
+
γ
4
j
·
∑
l
=
1
j
‍
α
j
l
+
γ
j
l
J
n
2
f
n
+
Δ
t
3
2
α
4
2
J
n
′
f
n
2
+
Δ
t
4
∑
i
=
1
4
‍
α
4
i
+
γ
4
i
·
∑
j
=
1
i
‍
∑
l
=
1
j
‍
α
i
j
+
γ
i
j
hhhhhhhhhhhhhhhhhhhhhhh
·
∑
l
=
1
j
‍
α
j
l
+
γ
j
l
∑
l
=
1
j
‍
α
i
j
+
γ
i
j
α
4
i
+
γ
4
i
·
∑
j
=
1
i
‍
∑
l
=
1
j
‍

α
i
j
+
γ
i
j
∑
i
=
1
4
‍
α
4
i
+
γ
4
i
·
∑
j
=
1
i
‍
∑
l
=
1
j
‍

α
i
j
+
γ
i
j
J
n
3
f
n
+
Δ
t
4
α
4
γ
33
+
β
32
+
β
31
α
43
hhhhhhhhh
+
γ
22
+
β
21
α
42
+
γ
11
α
41
J
n
′
J
n
f
n
2
+
Δ
t
4
2
γ
44
α
4
2
+
β
43
α
3
2
+
β
42
α
2
2
J
n
J
n
′
f
n
2
+
Δ
t
4
6
α
41
+
α
42
+
α
43
3
J
n
′
′
f
n
3
+
O
Δ
t
5
.Inserting the four Taylor series into (35) and matching the coefficients of Δ
t
p for p
=
0,1
,
2,3
,
4 on both sides of (31), we obtain the following order conditions:(40)
1
=
∑
i
=
1
4
‍
b
i
,
(41)
1
2
=
∑
i
=
1
4
‍
b
i
·
γ
i
i
+
∑
j
=
1
i
-
1
‍
β
i
j
,
(42)
1
6
=
∑
i
=
1
4
‍
b
i
·
∑
j
=
1
i
‍
β
i
j
·
∑
l
=
1
j
‍
β
j
l
,
(43)
1
3
=
∑
i
=
1
4
‍
b
i
·
∑
j
=
1
i
-
1
‍
α
i
j
2
,
(44)
1
4
=
∑
i
=
1
4
‍
b
i
·
∑
j
=
1
i
-
1
‍
α
i
j
3
,
(45)
1
24
=
∑
i
=
1
4
‍
b
i
·
α
i
∑
j
=
1
i
-
1
‍
α
i
j
·
γ
j
j
+
∑
l
=
1
j
-
1
‍
β
j
l
,
(46)
1
12
=
∑
i
=
1
4
‍
b
i
·
∑
j
=
1
i
‍
β
i
j
·
∑
l
=
1
j
-
1
‍
α
j
l
2
,
(47)
1
24
=
∑
i
=
1
4
‍
b
i
·
∑
j
=
1
i
‍
β
i
j
·
∑
l
=
1
j
‍
β
j
l
·
∑
k
=
1
l
‍
β
l
k
.Note that the set of order conditions obtained by using Butcher series [14, page 108] is a special case here when γ
i
i
=
γ.Also it is worthy to point out that iff
(
y
) is a scalar function of a single variable, J
n
J
n
′
=
J
n
′
J
n, thus conditions (45) and (46) can be combined to one single condition. However here since the method is developed for ODE system, both conditions should be satisfied individually.Apparently there are 12 degrees of freedom to determine the method defined in (32) since there are 20 parameters while only 8 constraints are given by (40)–(47). To simplify the procedure of determining these parameters and reduce the number of matrix inversion, we assume that β
i
i
=
γ, for i
=
1,2
,
3,4, so that only one matrix inversion is required to solve all k
i during each time step. After this simplification, there are 17 parameters left, and the order conditions are simplified as well.It is known that the Rosenbrock method suffers from order reduction when applied to nonlinear parabolic partial differential equations; see [20, 21]. In order to avoid such reduction of accuracy, the following two extra order conditions, simplified after taking into account the previous order conditions, should be satisfied [22]:(48)
b
4
β
43
β
32
α
2
2
=
-
2
γ
4
+
4
γ
3
-
5
3
γ
2
+
γ
6
,
(49)
b
4
β
43
β
32
α
2
2
=
-
8
γ
4
3
+
5
γ
3
-
2
γ
2
+
7
γ
36
,which imply(50)
γ
3
-
3
2
γ
2
+
γ
2
-
1
24
=
0
.Solving (50) results in three distinct real roots, while stability analysis shows that only one (
γ
=
1.06857902130162885
) can ensure A-stability for the Rosenbrock method. More details regarding stability will be provided later.To perform step size control and error estimation, we need a third-order embedded formula [23] defined as y
n
+
1
=
y
n
+
∑
j
=
1
s
‍
b
~
j
k
j which uses the same values of k
j’s as the Rosenbrock method in (31) but coefficients of k
j’s are different. To ensure the existence of such embedded formula, we need Det
(
A
^
)
=
0, where A
^ is a 4
×
4 square matrix defined as(51)
A
^
=
1
1
1
1
β
1
β
2
β
3
β
4
β
1
2
∑
i
=
1
2
‍
β
2
i
·
β
i
∑
i
=
1
3
‍
β
3
i
·
β
i
∑
i
=
1
4
‍
β
4
i
·
β
i
0
α
2
2
α
3
2
α
4
2and β
i
=
∑
j
=
1
i
β
i
j.Taking into account the order conditions, we obtain the simplified condition asβ
¯
3
2
α
2
2
β
43
=
0, which unfortunately contradicts the order condition given in (47). Thus one can either sacrifices efficiency by relaxing the condition γ
i
i
=
γ for the existence of a third-order embedded formula or simply skip the condition and then ignore step size control and error estimation. In this work, we are interested in the efficiency and order reduction of the Rosenbrock method, so the condition for step size control is bypassed.Upon the determination ofγ, there are 17 free parameters remaining while the order conditions are simplified as the following:(52)
b
1
+
b
2
+
b
3
+
b
4
=
1
,
(53)
b
2
β
¯
2
+
b
3
β
¯
3
+
b
4
β
¯
4
=
1
2
-
γ
,
(54)
b
3
β
32
β
¯
2
+
b
4
(
β
42
β
¯
2
+
β
43
β
¯
3
)
=
1
6
-
γ
+
γ
2
,
(55)
b
2
α
21
2
+
b
3
α
3
2
+
b
4
α
4
2
=
1
3
,
(56)
b
2
α
21
3
+
b
3
α
3
3
+
b
4
α
4
3
=
1
4
,
(57)
b
3
α
3
α
32
β
¯
2
+
b
4
α
4
(
α
42
β
¯
2
+
α
43
β
¯
3
)
=
1
8
-
γ
3
,
(58)
b
3
β
32
α
2
2
+
b
4
(
β
42
α
2
2
+
β
43
α
3
2
)
=
1
12
-
γ
3
,
(59)
b
4
β
43
β
32
β
¯
2
=
1
24
-
γ
2
+
3
2
γ
2
-
γ
3
,which are equivalent to those given in [14, page 108], where β
¯
i
=
∑
j
=
1
i
-
1
β
i
j.Sinceγ is a root of (50), b
4
β
43
β
32
β
¯
2
=
0, but b
4
β
43
β
32
α
2
2
≠
0, thus β
21
=
0. Then the order conditions in (52)–(58) can be further simplified.We now choose the coefficientsα
i
j with the purpose of reducing the number of functional evaluations. Particularly, letting α
2
=
α
3
=
α
4, (55)-(56) imply that α
21
=
α
2
=
3
/
4, which further implies b
2
+
b
3
+
b
4
=
16
/
27, so b
1
=
11
/
27. Letting α
31
=
α
21 and α
32
=
0, we can reduce the number of functional evaluation by 1.We chooseβ
43 as a free parameter, so(60)
b
4
β
¯
3
=
1
-
6
γ
+
6
γ
2
6
β
43
,
α
43
=
(
1
-
8
γ
/
3
)
β
43
1
-
6
γ
+
6
γ
2
.By choosingα
42 as a free parameter, we have(61)
α
41
=
3
4
-
α
42
-
(
1
-
8
γ
/
3
)
β
43
1
-
6
γ
+
6
γ
2
.Lettingβ
32 be a free parameter, we obtain(62)
b
4
=
16
(
-
2
γ
4
+
4
γ
3
-
5
γ
2
/
3
+
γ
/
6
)
9
β
43
β
32
,
β
31
=
1
-
6
γ
+
6
γ
2
6
β
43
b
4
-
β
32
.Finally, choose b
3 as free parameter, so we have b
2
=
16
/
27
-
b
3
-
b
4 and(63)
β
42
=
4
-
16
γ
-
27
b
3
β
32
27
b
4
-
β
43
,
β
41
=
1
/
2
-
γ
-
b
3
β
31
+
β
32
b
4
-
β
42
-
β
43
.In summary, the coefficients of the Rosenbrock method are listed in Table1.Table 1
Coefficients of the 4th-order Rosenbrock method.
b
1
=
0.4074074074074
b
2
=
-
0.2568608534470
b
3
=
0.2
b
4
=
0.6494534460396
α
21
=
0.75
α
31
=
0.75
α
32
=
0
α
41
=
2.9193596398302
α
42
=
0.4
α
43
=
-
2.5693596398302
γ
21
=
-
0.75
γ
31
=
-
1.3152686912402
γ
32
=
0.75
γ
41
=
-
2.8738466294648
γ
42
=
-
3.3778743470341
γ
43
=
4.5693596398302
γ
=
1.068579021301629Note that one can also use these free parameters to eliminate some fifth-order truncation error terms, so the truncation error constant can be further optimized.


### 3.2. Extend the Rosenbrock Method to ODE System

We now extend the Rosenbrock algorithm to the ODE systemA
U
′
(
t
)
=
F
(
U
,
t
). Suppose that U
n
=
(
u
0
n
,
u
1
n
,
…
,
u
M
n
) is the numerical solution of (10) at time t
n; then U
n
+
1 is calculated as(64)
U
n
+
1
=
U
n
+
b
1
K
1
+
b
2
K
2
+
b
3
K
3
+
b
4
K
4
,
(65)
A
-
γ
Δ
t
J
n
K
i
=
F
U
n
+
∑
j
=
1
i
-
1
‍
α
i
j
K
j
,
t
n
+
σ
i
Δ
t
+
μ
i
Δ
t
2
∂
F
∂
t
U
n
,
t
n
+
Δ
t
J
n
∑
j
=
1
i
-
1
‍
γ
i
j
K
j
,
hhhhhhhhhhhhhhhhhhhhhhh
i
=
1,2
,
3,4
,where J
n
=
∂
F
/
∂
U
(
U
n
,
t
n
) is the Jacobi matrix σ
i and μ
i are defined through (34).Obviously, each stage of this method consists of solving a linear system ofM
+
1 equations. Since γ
11
=
⋯
=
γ
44
=
γ, only oneLU-decomposition is required during each time step.


### 3.3. Stability of the Rosenbrock Method

We now study the stability of the Rosenbrock method defined by (31)-(32) with coefficients given in Table 1. Applying the Rosenbrock method to the linear test equation y
′
=
λ
y, λ
<
0, and taking into account the order condition equation (49), we obtain the stability function of the Rosenbrock method as(66)
R
z
=
5
γ
2
4
-
5
γ
8
+
1
16
z
4
1
+
1
-
4
γ
z
+
1
2
-
4
γ
+
6
γ
2
z
2
hhh
-
5
γ
2
4
-
5
γ
8
+
1
16
z
4
×
(
1
-
γ
z
)
4
-
1
,where z
=
λ
Δ
t.In order for the Rosenbrock method to be strongly A-stable, we need(67)
R
(
z
)
≤
1
,
∀
z
∈
C
-
,
(68)
lim
z
→
-
∞
R
(
z
)
<
1
.Apparently,(69)
lim
z
→
-
∞
R
(
z
)
=
5
γ
2
/
4
-
5
γ
/
8
+
1
/
16
γ
4
=
0.6304149382
<
1
,so (68) is satisfied.To prove (67), we can substitute z with x
+
i
y and then show that R
(
x
+
i
y
)
≤
1 for any x
<
0. However such proof is technically tedious although possible. Here we use the Boundary-Locus method to draw the region of absolute stability of the method, which is shown in Figure 1. Obviously, the region of absolute stability contains the whole half-plane C
-, so the Rosenbrock method is A-stable. Combined with (69), we can prove that the Rosenbrock method is strongly A-stable.Figure 1
Absolute stable region of the Rosenbrock method.


## 3.1. Rosenbrock Method for Scalar Equation

We first derive the Rosenbrock method based on an autonomous scalar equationy
′
=
f
(
y
), for which the initial condition is given as y
(
t
0
)
=
y
0. Nonautonomous equations can be converted to autonomous form by adding an extra equation to the system. Some previous research [20] suggested that it is unlikely to find a 3-stage fourth-order Rosenbrock method with strong A-stability or L-stability, so herein we focus on the development of a 4-stage Rosenbrock method. Suppose the numerical solution at time t
n is known as y
n, the 4-stage Rosenbrock method calculates the numerical solution at t
n
+
1 as(31)
y
n
+
1
=
y
n
+
b
1
k
1
+
b
2
k
2
+
b
3
k
3
+
b
4
k
4
,
(32)
k
i
=
Δ
t
f
y
n
+
∑
j
=
1
i
-
1
‍
α
i
j
k
j
+
Δ
t
J
y
n
∑
j
=
1
i
‍
γ
i
j
k
j
,
hhhhhhhhhhhhhhhhihhhh
i
=
1,2
,
3,4
,where b
i
,
α
i
j, and γ
i
j are coefficients to be determined and J
(
y
n
)
=
f
y
(
y
n
).To extend the above algorithm to the nonautonomous problemsy
′
=
f
(
y
,
t
), we first convert it to autonomous form by adding a new equation t
′
=
1 and then apply the algorithm (31) to the augmented system. Note that the components corresponding to the last variable can be computed explicitly; thus we can derive the modified algorithm as the following:(33)
y
n
+
1
=
y
n
+
b
1
k
1
+
b
2
k
2
+
b
3
k
3
+
b
4
k
4
,
hh
k
i
=
Δ
t
f
y
n
+
∑
j
=
1
i
-
1
‍
α
i
j
k
j
,
t
n
+
α
i
Δ
t
hhhhh
+
γ
i
Δ
t
2
f
t
y
n
,
t
n
+
Δ
t
f
y
y
n
,
t
n
∑
j
=
1
i
‍
γ
i
j
k
j
,
hhhhhhhhhhhhhhhhhhhhhhhhhh
i
=
1,2
,
3,4
,where(34)
α
i
=
∑
j
=
1
i
-
1
‍
α
i
j
,
γ
i
=
∑
j
=
1
i
‍
γ
i
j
.A Rosenbrock method of orderp is obtained through choosing coefficients in (31)-(32) so that the local error satisfies y
(
t
n
+
Δ
t
)
-
y
n
+
1
=
O
(
Δ
t
p
+
1
). This can be done either by solving the so-calledButcher Series [14] or by straightforward differentiation. Here we derive the order conditions for the fourth-order Rosenbrock method in a different way using Taylor series. First, both y
(
t
n
+
Δ
t
) and k
i are expanded as Taylor series so that the difference y
(
t
n
+
Δ
t
)
-
y
n
+
1 can be expressed as a Taylor series and its coefficients of the terms up to O
(
Δ
t
4
) are set to 0, which results in a set of equations involving these coefficients.Giveny
(
t
n
)
=
y
n, the Taylor series of y
(
t
n
+
Δ
t
) at t
n is expanded as(35)
y
t
n
+
Δ
t
=
y
n
+
Δ
t
f
n
+
Δ
t
2
2
J
n
f
n
+
Δ
t
3
6
J
n
2
f
n
+
J
n
′
f
n
2
+
Δ
t
4
24
J
n
′
′
f
n
3
+
2
J
n
′
J
n
f
n
2
+
2
J
n
J
n
′
f
n
2
+
J
n
3
f
n
+
O
Δ
t
5
,where f
n
=
f
(
y
n
), J
n
=
∂
f
/
∂
y
(
y
n
), J
n
′
=
∂
2
f
/
∂
y
2
y
n, and J
n
′′
=
∂
3
f
/
∂
y
3
y
n. For the sake of simplicity, let β
i
j
=
α
i
j
+
γ
i
j and α
i
i
=
0.Lettingi
=
1 in (32) we have(36)
k
1
=
Δ
t
1
-
Δ
t
γ
11
J
n
-
1
f
n
=
Δ
t
f
n
+
Δ
t
2
γ
11
J
n
f
n
+
Δ
t
3
γ
11
2
J
n
2
f
n
+
Δ
t
4
γ
11
3
J
n
3
f
n
+
O
Δ
t
5
.Similarly, lettingi
=
2 in (32) we have(37)
k
2
=
Δ
t
1
-
Δ
t
γ
22
J
n
-
1
f
y
n
+
α
21
k
1
+
γ
21
J
n
k
1
.Combining it with the Taylor series ofk
1 in (36), we obtain(38)
k
2
=
Δ
t
f
n
+
Δ
t
2
α
21
+
γ
21
+
γ
22
J
n
f
n
+
Δ
t
3
γ
11
+
γ
22
α
21
+
γ
21
+
γ
22
2
J
n
2
f
n
+
Δ
t
3
2
α
21
2
J
n
′
f
n
2
+
Δ
t
4
γ
11
γ
22
+
γ
11
2
+
γ
22
2
α
21
+
γ
21
+
γ
22
3
J
n
3
f
n
+
Δ
t
4
γ
11
α
21
2
J
n
′
J
n
f
n
2
+
Δ
t
4
2
γ
22
α
21
2
J
n
J
n
′
f
n
2
+
Δ
t
4
6
α
21
3
J
n
′
′
f
n
3
+
O
Δ
t
5
.Following the same way, we have(39)
k
3
=
Δ
t
f
n
+
Δ
t
2
(
β
31
+
β
32
+
γ
33
)
J
n
f
n
+
Δ
t
3
∑
j
=
1
3
‍
α
3
j
+
γ
3
j
·
∑
l
=
1
j
‍
α
j
l
+
γ
j
l
J
n
2
f
n
+
Δ
t
3
2
α
3
2
J
n
′
f
n
2
+
Δ
t
4
∑
i
=
1
3
‍
α
3
i
+
γ
3
i
·
∑
j
=
1
i
‍
∑
l
=
1
j
‍

α
i
j
+
γ
i
j
hhhhhhhhhhhhhhhhhhhhhhh
·
∑
l
=
1
j
‍
α
j
l
+
γ
j
l
∑
l
=
1
j
‍

α
i
j
+
γ
i
j
α
3
i
+
γ
3
i
·
∑
j
=
1
i
‍
∑
i
=
1
3
‍
α
3
i
+
γ
3
i
·
∑
j
=
1
i
‍
∑
l
=
1
j
‍

α
i
j
+
γ
i
j
J
n
3
f
n
+
Δ
t
4
2
[
α
2
2
β
32
+
γ
33
α
3
2
]
J
n
J
n
′
f
n
2
+
Δ
t
4
α
3
β
21
+
γ
22
α
32
+
γ
11
α
31
J
n
′
J
n
f
n
2
+
Δ
t
4
6
α
31
+
α
32
3
J
n
′
′
f
n
3
+
O
Δ
t
5
,
k
4
=
Δ
t
f
n
+
Δ
t
2
β
41
+
β
42
+
β
42
+
γ
44
J
n
f
n
+
Δ
t
3
∑
j
=
1
4
‍
α
4
j
+
γ
4
j
·
∑
l
=
1
j
‍
α
j
l
+
γ
j
l
J
n
2
f
n
+
Δ
t
3
2
α
4
2
J
n
′
f
n
2
+
Δ
t
4
∑
i
=
1
4
‍
α
4
i
+
γ
4
i
·
∑
j
=
1
i
‍
∑
l
=
1
j
‍
α
i
j
+
γ
i
j
hhhhhhhhhhhhhhhhhhhhhhh
·
∑
l
=
1
j
‍
α
j
l
+
γ
j
l
∑
l
=
1
j
‍
α
i
j
+
γ
i
j
α
4
i
+
γ
4
i
·
∑
j
=
1
i
‍
∑
l
=
1
j
‍

α
i
j
+
γ
i
j
∑
i
=
1
4
‍
α
4
i
+
γ
4
i
·
∑
j
=
1
i
‍
∑
l
=
1
j
‍

α
i
j
+
γ
i
j
J
n
3
f
n
+
Δ
t
4
α
4
γ
33
+
β
32
+
β
31
α
43
hhhhhhhhh
+
γ
22
+
β
21
α
42
+
γ
11
α
41
J
n
′
J
n
f
n
2
+
Δ
t
4
2
γ
44
α
4
2
+
β
43
α
3
2
+
β
42
α
2
2
J
n
J
n
′
f
n
2
+
Δ
t
4
6
α
41
+
α
42
+
α
43
3
J
n
′
′
f
n
3
+
O
Δ
t
5
.Inserting the four Taylor series into (35) and matching the coefficients of Δ
t
p for p
=
0,1
,
2,3
,
4 on both sides of (31), we obtain the following order conditions:(40)
1
=
∑
i
=
1
4
‍
b
i
,
(41)
1
2
=
∑
i
=
1
4
‍
b
i
·
γ
i
i
+
∑
j
=
1
i
-
1
‍
β
i
j
,
(42)
1
6
=
∑
i
=
1
4
‍
b
i
·
∑
j
=
1
i
‍
β
i
j
·
∑
l
=
1
j
‍
β
j
l
,
(43)
1
3
=
∑
i
=
1
4
‍
b
i
·
∑
j
=
1
i
-
1
‍
α
i
j
2
,
(44)
1
4
=
∑
i
=
1
4
‍
b
i
·
∑
j
=
1
i
-
1
‍
α
i
j
3
,
(45)
1
24
=
∑
i
=
1
4
‍
b
i
·
α
i
∑
j
=
1
i
-
1
‍
α
i
j
·
γ
j
j
+
∑
l
=
1
j
-
1
‍
β
j
l
,
(46)
1
12
=
∑
i
=
1
4
‍
b
i
·
∑
j
=
1
i
‍
β
i
j
·
∑
l
=
1
j
-
1
‍
α
j
l
2
,
(47)
1
24
=
∑
i
=
1
4
‍
b
i
·
∑
j
=
1
i
‍
β
i
j
·
∑
l
=
1
j
‍
β
j
l
·
∑
k
=
1
l
‍
β
l
k
.Note that the set of order conditions obtained by using Butcher series [14, page 108] is a special case here when γ
i
i
=
γ.Also it is worthy to point out that iff
(
y
) is a scalar function of a single variable, J
n
J
n
′
=
J
n
′
J
n, thus conditions (45) and (46) can be combined to one single condition. However here since the method is developed for ODE system, both conditions should be satisfied individually.Apparently there are 12 degrees of freedom to determine the method defined in (32) since there are 20 parameters while only 8 constraints are given by (40)–(47). To simplify the procedure of determining these parameters and reduce the number of matrix inversion, we assume that β
i
i
=
γ, for i
=
1,2
,
3,4, so that only one matrix inversion is required to solve all k
i during each time step. After this simplification, there are 17 parameters left, and the order conditions are simplified as well.It is known that the Rosenbrock method suffers from order reduction when applied to nonlinear parabolic partial differential equations; see [20, 21]. In order to avoid such reduction of accuracy, the following two extra order conditions, simplified after taking into account the previous order conditions, should be satisfied [22]:(48)
b
4
β
43
β
32
α
2
2
=
-
2
γ
4
+
4
γ
3
-
5
3
γ
2
+
γ
6
,
(49)
b
4
β
43
β
32
α
2
2
=
-
8
γ
4
3
+
5
γ
3
-
2
γ
2
+
7
γ
36
,which imply(50)
γ
3
-
3
2
γ
2
+
γ
2
-
1
24
=
0
.Solving (50) results in three distinct real roots, while stability analysis shows that only one (
γ
=
1.06857902130162885
) can ensure A-stability for the Rosenbrock method. More details regarding stability will be provided later.To perform step size control and error estimation, we need a third-order embedded formula [23] defined as y
n
+
1
=
y
n
+
∑
j
=
1
s
‍
b
~
j
k
j which uses the same values of k
j’s as the Rosenbrock method in (31) but coefficients of k
j’s are different. To ensure the existence of such embedded formula, we need Det
(
A
^
)
=
0, where A
^ is a 4
×
4 square matrix defined as(51)
A
^
=
1
1
1
1
β
1
β
2
β
3
β
4
β
1
2
∑
i
=
1
2
‍
β
2
i
·
β
i
∑
i
=
1
3
‍
β
3
i
·
β
i
∑
i
=
1
4
‍
β
4
i
·
β
i
0
α
2
2
α
3
2
α
4
2and β
i
=
∑
j
=
1
i
β
i
j.Taking into account the order conditions, we obtain the simplified condition asβ
¯
3
2
α
2
2
β
43
=
0, which unfortunately contradicts the order condition given in (47). Thus one can either sacrifices efficiency by relaxing the condition γ
i
i
=
γ for the existence of a third-order embedded formula or simply skip the condition and then ignore step size control and error estimation. In this work, we are interested in the efficiency and order reduction of the Rosenbrock method, so the condition for step size control is bypassed.Upon the determination ofγ, there are 17 free parameters remaining while the order conditions are simplified as the following:(52)
b
1
+
b
2
+
b
3
+
b
4
=
1
,
(53)
b
2
β
¯
2
+
b
3
β
¯
3
+
b
4
β
¯
4
=
1
2
-
γ
,
(54)
b
3
β
32
β
¯
2
+
b
4
(
β
42
β
¯
2
+
β
43
β
¯
3
)
=
1
6
-
γ
+
γ
2
,
(55)
b
2
α
21
2
+
b
3
α
3
2
+
b
4
α
4
2
=
1
3
,
(56)
b
2
α
21
3
+
b
3
α
3
3
+
b
4
α
4
3
=
1
4
,
(57)
b
3
α
3
α
32
β
¯
2
+
b
4
α
4
(
α
42
β
¯
2
+
α
43
β
¯
3
)
=
1
8
-
γ
3
,
(58)
b
3
β
32
α
2
2
+
b
4
(
β
42
α
2
2
+
β
43
α
3
2
)
=
1
12
-
γ
3
,
(59)
b
4
β
43
β
32
β
¯
2
=
1
24
-
γ
2
+
3
2
γ
2
-
γ
3
,which are equivalent to those given in [14, page 108], where β
¯
i
=
∑
j
=
1
i
-
1
β
i
j.Sinceγ is a root of (50), b
4
β
43
β
32
β
¯
2
=
0, but b
4
β
43
β
32
α
2
2
≠
0, thus β
21
=
0. Then the order conditions in (52)–(58) can be further simplified.We now choose the coefficientsα
i
j with the purpose of reducing the number of functional evaluations. Particularly, letting α
2
=
α
3
=
α
4, (55)-(56) imply that α
21
=
α
2
=
3
/
4, which further implies b
2
+
b
3
+
b
4
=
16
/
27, so b
1
=
11
/
27. Letting α
31
=
α
21 and α
32
=
0, we can reduce the number of functional evaluation by 1.We chooseβ
43 as a free parameter, so(60)
b
4
β
¯
3
=
1
-
6
γ
+
6
γ
2
6
β
43
,
α
43
=
(
1
-
8
γ
/
3
)
β
43
1
-
6
γ
+
6
γ
2
.By choosingα
42 as a free parameter, we have(61)
α
41
=
3
4
-
α
42
-
(
1
-
8
γ
/
3
)
β
43
1
-
6
γ
+
6
γ
2
.Lettingβ
32 be a free parameter, we obtain(62)
b
4
=
16
(
-
2
γ
4
+
4
γ
3
-
5
γ
2
/
3
+
γ
/
6
)
9
β
43
β
32
,
β
31
=
1
-
6
γ
+
6
γ
2
6
β
43
b
4
-
β
32
.Finally, choose b
3 as free parameter, so we have b
2
=
16
/
27
-
b
3
-
b
4 and(63)
β
42
=
4
-
16
γ
-
27
b
3
β
32
27
b
4
-
β
43
,
β
41
=
1
/
2
-
γ
-
b
3
β
31
+
β
32
b
4
-
β
42
-
β
43
.In summary, the coefficients of the Rosenbrock method are listed in Table1.Table 1
Coefficients of the 4th-order Rosenbrock method.
b
1
=
0.4074074074074
b
2
=
-
0.2568608534470
b
3
=
0.2
b
4
=
0.6494534460396
α
21
=
0.75
α
31
=
0.75
α
32
=
0
α
41
=
2.9193596398302
α
42
=
0.4
α
43
=
-
2.5693596398302
γ
21
=
-
0.75
γ
31
=
-
1.3152686912402
γ
32
=
0.75
γ
41
=
-
2.8738466294648
γ
42
=
-
3.3778743470341
γ
43
=
4.5693596398302
γ
=
1.068579021301629Note that one can also use these free parameters to eliminate some fifth-order truncation error terms, so the truncation error constant can be further optimized.


## 3.2. Extend the Rosenbrock Method to ODE System

We now extend the Rosenbrock algorithm to the ODE systemA
U
′
(
t
)
=
F
(
U
,
t
). Suppose that U
n
=
(
u
0
n
,
u
1
n
,
…
,
u
M
n
) is the numerical solution of (10) at time t
n; then U
n
+
1 is calculated as(64)
U
n
+
1
=
U
n
+
b
1
K
1
+
b
2
K
2
+
b
3
K
3
+
b
4
K
4
,
(65)
A
-
γ
Δ
t
J
n
K
i
=
F
U
n
+
∑
j
=
1
i
-
1
‍
α
i
j
K
j
,
t
n
+
σ
i
Δ
t
+
μ
i
Δ
t
2
∂
F
∂
t
U
n
,
t
n
+
Δ
t
J
n
∑
j
=
1
i
-
1
‍
γ
i
j
K
j
,
hhhhhhhhhhhhhhhhhhhhhhh
i
=
1,2
,
3,4
,where J
n
=
∂
F
/
∂
U
(
U
n
,
t
n
) is the Jacobi matrix σ
i and μ
i are defined through (34).Obviously, each stage of this method consists of solving a linear system ofM
+
1 equations. Since γ
11
=
⋯
=
γ
44
=
γ, only oneLU-decomposition is required during each time step.


## 3.3. Stability of the Rosenbrock Method

We now study the stability of the Rosenbrock method defined by (31)-(32) with coefficients given in Table 1. Applying the Rosenbrock method to the linear test equation y
′
=
λ
y, λ
<
0, and taking into account the order condition equation (49), we obtain the stability function of the Rosenbrock method as(66)
R
z
=
5
γ
2
4
-
5
γ
8
+
1
16
z
4
1
+
1
-
4
γ
z
+
1
2
-
4
γ
+
6
γ
2
z
2
hhh
-
5
γ
2
4
-
5
γ
8
+
1
16
z
4
×
(
1
-
γ
z
)
4
-
1
,where z
=
λ
Δ
t.In order for the Rosenbrock method to be strongly A-stable, we need(67)
R
(
z
)
≤
1
,
∀
z
∈
C
-
,
(68)
lim
z
→
-
∞
R
(
z
)
<
1
.Apparently,(69)
lim
z
→
-
∞
R
(
z
)
=
5
γ
2
/
4
-
5
γ
/
8
+
1
/
16
γ
4
=
0.6304149382
<
1
,so (68) is satisfied.To prove (67), we can substitute z with x
+
i
y and then show that R
(
x
+
i
y
)
≤
1 for any x
<
0. However such proof is technically tedious although possible. Here we use the Boundary-Locus method to draw the region of absolute stability of the method, which is shown in Figure 1. Obviously, the region of absolute stability contains the whole half-plane C
-, so the Rosenbrock method is A-stable. Combined with (69), we can prove that the Rosenbrock method is strongly A-stable.Figure 1
Absolute stable region of the Rosenbrock method.


## 4. Numerical Examples and Discussions

We solve three numerical examples to demonstrate the accuracy and efficiency of the new algorithm. The first and second examples are reaction-diffusion equations supplemented with Dirichlet and Neumann boundary conditions, respectively. The two examples are solved by the new method to demonstrate the fourth-order accuracy in time and space. The third example is used to demonstrate that the new algorithm is free of order reduction and more efficient than several other existing fourth-order Rosenbrock methods. For comparison, the rate of convergence is compared with several other fourth-order methods which suffer from order reduction. In what follows, we useHOC-ROSB4 to represent the new algorithm developed in this paper andHOC-GRK4A to represent the fourth-order method that combines high-order compact finite difference scheme and fourth-orderGRK4A [23].HOC-SHAMP represents the algorithm that combines high-order compact finite difference scheme with the fourth-order Rosenbrock methodShampine [24].HOC-LSTAB represents the algorithm that combines the high-order compact finite difference scheme in space and the L-stable fourth-order Rosenbrock method [14] in time. Finally, we useHOC-VELD to represent the algorithm that combines the high-order compact finite difference scheme in space with fourth-order Rosenbrock method proposed in [25].


### 4.1. Example  1

Consider the following reaction-diffusion equation:(70)
u
t
=
u
x
x
+
cos
⁡
⁡
u
-
cos
⁡
⁡
e
-
t
cos
⁡
⁡
x
,
hhhhhhhhihhh
(
x
,
t
)
∈
(
0,2
)
×
(
0,1
]
,
u
(
x
,
0
)
=
cos
⁡
(
x
)
,
x
∈
[
0,2
]
,
u
(
0
,
t
)
=
e
-
t
,
u
(
2
,
t
)
=
cos
⁡
(
2
)
e
-
t
,
t
∈
[
0,1
]
,for which the analytical solution is u
(
x
,
t
)
=
e
-
t
cos
⁡
(
x
).We first show that the compact finite difference scheme for spatial discretization and the compact boundary condition treatment are fourth-order accurate. In order to do so, we fixedΔ
t
=
0.0001 so the truncation error from time discretization is negligible. The results included in Table 2 clearly demonstrate the fourth-order convergence in space, as the maximal error reduced by a factor of 16 (roughly) when the grid size h is halved.Table 2
Numerical results ofexample 1 by HOC-ROSB4 with Δ
t
=
0.0001.
h
1/10
1/20
1/40
1/80
1/160
E
(
h
)
7.38
e
-
08
4.62
e
-
09
2.89
e
-
10
1.80
e
-
11
1.06
e
-
12
E
(
h
)
/
E
(
h
/
2
)
—
15.9562
16.0049
16.0607
17.0153
Convergence rate
—
3.9960
4.0004
4.0055
4.0888We then show the fourth-order accuracy in time by fixingh
=
0.001. The data in Table 3 confirmed that the Rosenbrock method is fourth-order accurate in time, so the new algorithm is free of order reduction.Table 3
Numerical results ofexample 1 by HOC-ROSB4 with h
=
0.001.
Δ
t
1/10
1/20
1/40
1/80
1/160
E
(
Δ
t
)
9.03
e
-
06
6.16
e
-
07
3.96
e
-
08
2.45
e
-
09
1.49
e
-
10
E
(
Δ
t
)/E(Δ
t
/
2)
—
14.6679
15.5402
16.1607
16.4631
Convergence rate
—
3.8746
3.9579
4.0144
4.0412Since the method is fourth-order accurate in both temporal and spatial dimensions, we can write the leading term of the truncation error asE
(
h
,
Δ
t
)
=
C
1
h
4
+
C
2
Δ
t
4, where C
1 and C
2 are two constants. To obtain optimal performance, we can adjust the ratio Δ
t
/
h to balance the two error terms. To do so, we adjust Δ
t and h so that the two error terms are balanced, which is utilized by letting C
2
Δ
t
4
/
C
1
h
4
≈
1, from which we can estimate the optimal ratio as (
C
1
/
C
2
)
1
/
4. To estimate C
1, we solve the example using small Δ
t and then C
1
≈
E
(
h
)
/
h
4. Similarly, we can estimate C
2 using the same method.Simple calculation based on numerical results from Tables2 and 3 suggests that the optimal ratio for this example is h
/
Δ
t
≈
3.2. We then solve the reaction-diffusion equation using the optimal ratio to show that the new algorithm is fourth-order accurate in both temporal and spatial dimensions. The numerical results in Table 4 indicate that the new algorithm is fourth-order accurate in both temporal and spatial dimensions, as the maximal error is reduced by a factor of 16 (roughly) when Δ
t and h are halved simultaneously. We also notice that the algorithm is very efficient as the maximal error drops to the level of 1
0
-
12 when h and Δ
t are still reasonably large.Table 4
Numerical results ofexample 1 by HOC-ROSB4 with h
/
Δ
t
=
3.2.
h
1/10
1/20
1/40
1/80
1/160
Δ
t
1/32
1/64
1/128
1/256
1/512
E
(
Δ
t
,
h
)
5.94
e
-
08
4.09
e
-
09
2.73
e
-
10
1.78
e
-
11
1.15
e
-
12
E
(
Δ
t
,
h
)
/
E
(
Δ
t
/
2
,
h
/
2
)
—
14.4970
14.9793
15.3189
15.4959
Convergence rate
—
3.8577
3.9049
3.9372
3.9538


### 4.2. Example  2

We solve the following semilinear parabolic partial differential equation with the Neumann boundary conditions:(71)
u
t
=
2
u
x
x
+
u
+
u
2
-
e
-
2
t
cos
⁡
2
(
x
)
,
(
x
,
t
)
∈
(
0,2
)
×
(
0,1
]
,
u
(
x
,
0
)
=
cos
⁡
(
x
)
,
x
∈
[
0,2
]
,
u
x
(
0
,
t
)
=
0
,
u
x
(
2
,
t
)
=
-
sin
(
2
)
e
-
t
,
t
∈
[
0,1
]
,for which the analytical solution is u
(
x
,
t
)
=
e
-
t
cos
⁡
(
x
).Notice that the boundary conditions are approximated by the compact fourth-order boundary scheme given in (18)-(19), so the compactness of the resulting linear system is preserved; consequently, in each time step, only a tridiagonal system is solved for each k
i, so the solution procedure is very efficient. However we point out that some tedious work is needed to form the semidiscrete ODE system, for instance, the first and last rows of the matrix A in (20) and the first and last components of F in (28)-(29).To obtain the optimal performance, we choose the optimal ratio asΔ
t as h
/
Δ
t
=
2.5 to balance the two error terms. The data in Table 5 shows that the method is fourth-order accurate in both time and space, as the maximum error is reduced by a factor of 16 (roughly), when h and Δ
t are halved.Table 5
Numerical results ofexample 2 by HOC-ROSB4 with h
/
Δ
t
=
2.5.
h
1/10
1/20
1/40
1/80
1/160
Δ
t
1/25
1/50
1/100
1/200
1/400
E
(
h
,
Δ
t
)
6.27
e
-
06
4.40
e
-
07
2.95
e
-
08
1.96
e
-
09
1.31
e
-
10
E
(
h
,
Δ
t
)
/
E
(
h
/
2
,
Δ
t
/
2
)
—
14.2495
14.9178
15.0510
15.0135
Convergence rate
—
3.8328
3.8990
3.9118
3.9082It is worthy to mention that one can also use any other one-sided formula to approximate the Neumann boundary condition, which apparently reduces the effort to derive the semidiscrete ODE system, but will destroy the tridiagonal structure of the matrixA and consequently affect the efficiency of the algorithm. Another side-effect of using one-sided approximation on the boundary is the stability issue that may arise because the modifications to matrix A may result in positive eigenvalues of A.


### 4.3. Example  3

In this example we compare the new algorithm with several other fourth-order Rosenbrock methods in terms of the rate of convergence in time and efficiency. The nonlinear parabolic partial differential equation to be solved is defined as(72)
u
t
=
u
x
x
+
u
3
-
e
-
3
t
cos
⁡
3
(
x
)
,
(
x
,
t
)
∈
(
0,1
)
×
(
0,1
]
u
(
x
,
0
)
=
cos
⁡
(
x
)
,
x
∈
[
0,1
]
,
u
(
0
,
t
)
=
e
-
t
,
u
(
1
,
t
)
=
cos
⁡
(
1
)
e
-
t
,
t
∈
[
0,1
]
,for which the exact solution is u
(
x
,
t
)
=
e
-
t
cos
⁡
(
x
).The data in Table6 shows thatHOC-ROSB4 is fourth-order accurate and thus is free of order reduction, while the other four fourth-order Rosenbrock methods, which obviously suffer from order reduction, show rates of convergence ranging from 3.0 to 3.25. Note that h
=
0.001 is fixed and the same spatial discretization method is used for all five methods, so the demonstrated rate of convergence is in time. One can see that theHOC-LSTAB method has the most optimal error constant, while theHOC-ROSB4 method has the largest error constant. However, due to the freedom in determining those coefficients in Rosenbrock method, better error constant can be accomplished by eliminating several fifth-order truncation error terms.Table 6
Comparison of the rate of convergence among the five fourth-order Rosenbrock methods forexample 3 with h
=
0.001. E(Δ
t) represents the maximum error of the numerical solution obtained by using Δ
t.
Δ
t
1/10
1/20
1/40
1/80
HOC-ROSB4
E(Δ
t)
9.59
e
-
06
6.94
e
-
0
4.58
e
-
08
2.88
e
-
09
Rate of Conv.
—
3.7875
7
3.9221
3.9896
HOC-GRK4A
E(Δ
t)
5.88
e
-
06
6.16
e
-
07
6.83
e
-
08
7.56
e
-
09
Rate of Conv.
—
3.2552
3.1728
3.1757
HOC-LSTAB
E(Δ
t)
3.96
e
-
06
4.56
e
-
07
5.23
e
-
08
6.19
e
-
09
Rate of Conv.
—
3.1205
3.1223
3.0795
HOC-VELD
E(Δ
t)
8.41
e
-
06
8.77
e
-
07
9.77
e
-
08
1.14
e
-
08
Rate of Conv.
—
3.2613
3.1653
3.0996
HOC-SHAMP
E(Δ
t)
6.35
e
-
06
6.34
e
-
07
6.92
e
-
08
8.01
e
-
09
Rate of Conv.
—
3.3230
3.1965
3.1099Finally, we show thatHOC-ROSB4 is the most efficient method. For consistency, we adjust Δ
t and h for each method to reach the same error level and record the average CPU time of 5 simulation runs. Since all of the five methods use the same spatial discretization, the same h is used for all methods; hence we adjust Δ
t only. The comparison result included in Table 7 clearly indicates that the new method is the most efficient one. The higher efficiency apparently is obtained from the fact that the new method is free of order reduction; hence large time step can be used to reach the same error level.Table 7
Comparison of the efficiency among the five fourth-order Rosenbrock methods forexample 3.
Method
(
Δ
t
,
h
)
Max. error
CPU time (seconds)
HOC-ROSB4
(1/180, 1/40)
7.72
e
-
011
0.107
HOC-GRK4A
(1/256, 1/40)
7.28
e
-
011
0.134
HOC-LSTAB
(1/360, 1/40)
7.83
e
-
011
0.174
HOC-VELD
(1/400, 1/40)
7.77
e
-
011
0.188
HOC-SHAMP
(1/360, 1/40)
7.60
e
-
011
0.170


## 4.1. Example  1

Consider the following reaction-diffusion equation:(70)
u
t
=
u
x
x
+
cos
⁡
⁡
u
-
cos
⁡
⁡
e
-
t
cos
⁡
⁡
x
,
hhhhhhhhihhh
(
x
,
t
)
∈
(
0,2
)
×
(
0,1
]
,
u
(
x
,
0
)
=
cos
⁡
(
x
)
,
x
∈
[
0,2
]
,
u
(
0
,
t
)
=
e
-
t
,
u
(
2
,
t
)
=
cos
⁡
(
2
)
e
-
t
,
t
∈
[
0,1
]
,for which the analytical solution is u
(
x
,
t
)
=
e
-
t
cos
⁡
(
x
).We first show that the compact finite difference scheme for spatial discretization and the compact boundary condition treatment are fourth-order accurate. In order to do so, we fixedΔ
t
=
0.0001 so the truncation error from time discretization is negligible. The results included in Table 2 clearly demonstrate the fourth-order convergence in space, as the maximal error reduced by a factor of 16 (roughly) when the grid size h is halved.Table 2
Numerical results ofexample 1 by HOC-ROSB4 with Δ
t
=
0.0001.
h
1/10
1/20
1/40
1/80
1/160
E
(
h
)
7.38
e
-
08
4.62
e
-
09
2.89
e
-
10
1.80
e
-
11
1.06
e
-
12
E
(
h
)
/
E
(
h
/
2
)
—
15.9562
16.0049
16.0607
17.0153
Convergence rate
—
3.9960
4.0004
4.0055
4.0888We then show the fourth-order accuracy in time by fixingh
=
0.001. The data in Table 3 confirmed that the Rosenbrock method is fourth-order accurate in time, so the new algorithm is free of order reduction.Table 3
Numerical results ofexample 1 by HOC-ROSB4 with h
=
0.001.
Δ
t
1/10
1/20
1/40
1/80
1/160
E
(
Δ
t
)
9.03
e
-
06
6.16
e
-
07
3.96
e
-
08
2.45
e
-
09
1.49
e
-
10
E
(
Δ
t
)/E(Δ
t
/
2)
—
14.6679
15.5402
16.1607
16.4631
Convergence rate
—
3.8746
3.9579
4.0144
4.0412Since the method is fourth-order accurate in both temporal and spatial dimensions, we can write the leading term of the truncation error asE
(
h
,
Δ
t
)
=
C
1
h
4
+
C
2
Δ
t
4, where C
1 and C
2 are two constants. To obtain optimal performance, we can adjust the ratio Δ
t
/
h to balance the two error terms. To do so, we adjust Δ
t and h so that the two error terms are balanced, which is utilized by letting C
2
Δ
t
4
/
C
1
h
4
≈
1, from which we can estimate the optimal ratio as (
C
1
/
C
2
)
1
/
4. To estimate C
1, we solve the example using small Δ
t and then C
1
≈
E
(
h
)
/
h
4. Similarly, we can estimate C
2 using the same method.Simple calculation based on numerical results from Tables2 and 3 suggests that the optimal ratio for this example is h
/
Δ
t
≈
3.2. We then solve the reaction-diffusion equation using the optimal ratio to show that the new algorithm is fourth-order accurate in both temporal and spatial dimensions. The numerical results in Table 4 indicate that the new algorithm is fourth-order accurate in both temporal and spatial dimensions, as the maximal error is reduced by a factor of 16 (roughly) when Δ
t and h are halved simultaneously. We also notice that the algorithm is very efficient as the maximal error drops to the level of 1
0
-
12 when h and Δ
t are still reasonably large.Table 4
Numerical results ofexample 1 by HOC-ROSB4 with h
/
Δ
t
=
3.2.
h
1/10
1/20
1/40
1/80
1/160
Δ
t
1/32
1/64
1/128
1/256
1/512
E
(
Δ
t
,
h
)
5.94
e
-
08
4.09
e
-
09
2.73
e
-
10
1.78
e
-
11
1.15
e
-
12
E
(
Δ
t
,
h
)
/
E
(
Δ
t
/
2
,
h
/
2
)
—
14.4970
14.9793
15.3189
15.4959
Convergence rate
—
3.8577
3.9049
3.9372
3.9538


## 4.2. Example  2

We solve the following semilinear parabolic partial differential equation with the Neumann boundary conditions:(71)
u
t
=
2
u
x
x
+
u
+
u
2
-
e
-
2
t
cos
⁡
2
(
x
)
,
(
x
,
t
)
∈
(
0,2
)
×
(
0,1
]
,
u
(
x
,
0
)
=
cos
⁡
(
x
)
,
x
∈
[
0,2
]
,
u
x
(
0
,
t
)
=
0
,
u
x
(
2
,
t
)
=
-
sin
(
2
)
e
-
t
,
t
∈
[
0,1
]
,for which the analytical solution is u
(
x
,
t
)
=
e
-
t
cos
⁡
(
x
).Notice that the boundary conditions are approximated by the compact fourth-order boundary scheme given in (18)-(19), so the compactness of the resulting linear system is preserved; consequently, in each time step, only a tridiagonal system is solved for each k
i, so the solution procedure is very efficient. However we point out that some tedious work is needed to form the semidiscrete ODE system, for instance, the first and last rows of the matrix A in (20) and the first and last components of F in (28)-(29).To obtain the optimal performance, we choose the optimal ratio asΔ
t as h
/
Δ
t
=
2.5 to balance the two error terms. The data in Table 5 shows that the method is fourth-order accurate in both time and space, as the maximum error is reduced by a factor of 16 (roughly), when h and Δ
t are halved.Table 5
Numerical results ofexample 2 by HOC-ROSB4 with h
/
Δ
t
=
2.5.
h
1/10
1/20
1/40
1/80
1/160
Δ
t
1/25
1/50
1/100
1/200
1/400
E
(
h
,
Δ
t
)
6.27
e
-
06
4.40
e
-
07
2.95
e
-
08
1.96
e
-
09
1.31
e
-
10
E
(
h
,
Δ
t
)
/
E
(
h
/
2
,
Δ
t
/
2
)
—
14.2495
14.9178
15.0510
15.0135
Convergence rate
—
3.8328
3.8990
3.9118
3.9082It is worthy to mention that one can also use any other one-sided formula to approximate the Neumann boundary condition, which apparently reduces the effort to derive the semidiscrete ODE system, but will destroy the tridiagonal structure of the matrixA and consequently affect the efficiency of the algorithm. Another side-effect of using one-sided approximation on the boundary is the stability issue that may arise because the modifications to matrix A may result in positive eigenvalues of A.


## 4.3. Example  3

In this example we compare the new algorithm with several other fourth-order Rosenbrock methods in terms of the rate of convergence in time and efficiency. The nonlinear parabolic partial differential equation to be solved is defined as(72)
u
t
=
u
x
x
+
u
3
-
e
-
3
t
cos
⁡
3
(
x
)
,
(
x
,
t
)
∈
(
0,1
)
×
(
0,1
]
u
(
x
,
0
)
=
cos
⁡
(
x
)
,
x
∈
[
0,1
]
,
u
(
0
,
t
)
=
e
-
t
,
u
(
1
,
t
)
=
cos
⁡
(
1
)
e
-
t
,
t
∈
[
0,1
]
,for which the exact solution is u
(
x
,
t
)
=
e
-
t
cos
⁡
(
x
).The data in Table6 shows thatHOC-ROSB4 is fourth-order accurate and thus is free of order reduction, while the other four fourth-order Rosenbrock methods, which obviously suffer from order reduction, show rates of convergence ranging from 3.0 to 3.25. Note that h
=
0.001 is fixed and the same spatial discretization method is used for all five methods, so the demonstrated rate of convergence is in time. One can see that theHOC-LSTAB method has the most optimal error constant, while theHOC-ROSB4 method has the largest error constant. However, due to the freedom in determining those coefficients in Rosenbrock method, better error constant can be accomplished by eliminating several fifth-order truncation error terms.Table 6
Comparison of the rate of convergence among the five fourth-order Rosenbrock methods forexample 3 with h
=
0.001. E(Δ
t) represents the maximum error of the numerical solution obtained by using Δ
t.
Δ
t
1/10
1/20
1/40
1/80
HOC-ROSB4
E(Δ
t)
9.59
e
-
06
6.94
e
-
0
4.58
e
-
08
2.88
e
-
09
Rate of Conv.
—
3.7875
7
3.9221
3.9896
HOC-GRK4A
E(Δ
t)
5.88
e
-
06
6.16
e
-
07
6.83
e
-
08
7.56
e
-
09
Rate of Conv.
—
3.2552
3.1728
3.1757
HOC-LSTAB
E(Δ
t)
3.96
e
-
06
4.56
e
-
07
5.23
e
-
08
6.19
e
-
09
Rate of Conv.
—
3.1205
3.1223
3.0795
HOC-VELD
E(Δ
t)
8.41
e
-
06
8.77
e
-
07
9.77
e
-
08
1.14
e
-
08
Rate of Conv.
—
3.2613
3.1653
3.0996
HOC-SHAMP
E(Δ
t)
6.35
e
-
06
6.34
e
-
07
6.92
e
-
08
8.01
e
-
09
Rate of Conv.
—
3.3230
3.1965
3.1099Finally, we show thatHOC-ROSB4 is the most efficient method. For consistency, we adjust Δ
t and h for each method to reach the same error level and record the average CPU time of 5 simulation runs. Since all of the five methods use the same spatial discretization, the same h is used for all methods; hence we adjust Δ
t only. The comparison result included in Table 7 clearly indicates that the new method is the most efficient one. The higher efficiency apparently is obtained from the fact that the new method is free of order reduction; hence large time step can be used to reach the same error level.Table 7
Comparison of the efficiency among the five fourth-order Rosenbrock methods forexample 3.
Method
(
Δ
t
,
h
)
Max. error
CPU time (seconds)
HOC-ROSB4
(1/180, 1/40)
7.72
e
-
011
0.107
HOC-GRK4A
(1/256, 1/40)
7.28
e
-
011
0.134
HOC-LSTAB
(1/360, 1/40)
7.83
e
-
011
0.174
HOC-VELD
(1/400, 1/40)
7.77
e
-
011
0.188
HOC-SHAMP
(1/360, 1/40)
7.60
e
-
011
0.170


## 5. Conclusion

An efficient fourth-order numerical algorithm that combines the Padé approximation in space and fourth-order accurate Rosenbrock method in time is proposed in this paper. Our investigation shows that many widely used fourth-order Rosenbrock methods (A-stable or L-stable) suffer from order reduction when they are used to solve nonlinear parabolic partial differential equation. To avoid order reduction, extra order conditions are required, which are implemented in this paper to develop the new algorithm. Also, it has been shown [26] that the extra condition to resolve order reduction contradicts the L-stability of the Rosenbrock method, so there is no L-stable Rosenbrock method that is also free of order reduction. The new method can be used to solve nonlinear parabolic partial differential equation with all types of boundary conditions; however, for Neumann or Robin boundary condition, extra efforts are needed to form the semidiscrete ODE system. Two numerical examples are solved to demonstrate that the new method is fourth-order accurate in both time and space, while the third example shows that the new method is free of order reduction and is very efficient. In the future, we plan to extend the new method to multidimensional problems supplemented with various types of boundary conditions.



---

*Source: 539652-2015-03-25.xml*