Some numerical algorithms for solving the highly ... - Semantic Scholar
Journal of Computational Physics 276 (2014) 235–251
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
Some numerical algorithms for solving the highly oscillatory second-order initial value problems ✩ Wenjie Liu, Boying Wu, Jiebao Sun ∗ Department of Mathematics, Harbin Institute of Technology, Harbin 150001, PR China
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 9 January 2014 Received in revised form 16 May 2014 Accepted 15 July 2014 Available online 28 July 2014 Keywords: Highly oscillatory problems Second-order initial value problems Block spectral collocation method Block boundary value method Implicit Runge–Kutta method Diagonally implicit Runge–Kutta method Total variation diminishing Runge–Kutta method
In this paper, some numerical algorithms (spectral collocation method, block spectral collocation method, boundary value method, block boundary value method, implicit Runge–Kutta method, diagonally implicit Runge–Kutta method and total variation diminishing Runge–Kutta method) are used to solve the highly oscillatory second-order initial value problems. We first derive these methods for the first-order initial value problems, and then extend these methods to the highly oscillatory nonlinear systems by matrix analysis methods. These new methods preserve the accuracy of the original methods and the main advantages of these new methods are low storage requirements and high efficiency. Extensive numerical results are presented to demonstrate the convergence properties of these methods. © 2014 Elsevier Inc. All rights reserved.
1. Introduction In this paper, we study the spectral collocation method (or the block spectral collocation method) [1–5], the boundary value method (or the block boundary value method) [6–8], the implicit Runge–Kutta method [9–11], the diagonally implicit Runge–Kutta method [9,11,12] and the total variation diminishing (TVD) Runge–Kutta method [13,14] for the following highly oscillatory second-order initial value problem
y (t ) = f t , y (t ), y (t ) , t 0 < t ≤ T , y (t 0 ) = y 0 , y (t 0 ) = y 0 .
(1)
The problem (1) is the fundamental model in physics for describing mechanical vibrations behavior [15]. This model also appears in circuit simulations, flexible body dynamics and various quantum dynamics calculations [16–18]. Therefore, the construction of efficient numerical schemes for solving (1) is an important task. Recently, several numerical schemes have been developed for (1), and some useful approaches to construct Runge–Kutta– Nystro¨ m (RKN) type methods have also been proposed. For details we refer to the monograph [16]. The block spectral collocation method (or the spectral collocation method), the block boundary value method (or the boundary value method), the implicit Runge–Kutta method, the diagonally implicit Runge–Kutta method and the TVD
✩ This work is partially supported by the National Science Foundation of China (11271100, 11301113) and the Fundamental Research Funds for the Central Universities (Grant No. HIT. IBRSEM. A. 201412), Harbin Science and Technology Innovative Talents Project of Special Fund (2013RFXYJ044). Corresponding author. E-mail address: [email protected] (J. Sun).
*
http://dx.doi.org/10.1016/j.jcp.2014.07.033 0021-9991/© 2014 Elsevier Inc. All rights reserved.
236
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
Runge–Kutta method are high-accuracy schemes for the first-order nonlinear initial value problems. So, we extend these methods for solving the second-order initial value problems by matrix analysis methods [19]. The main advantages of these new methods are low storage requirements and high efficiency. The outline of this paper is as follows: in Section 2, we describe the spectral collocation method (SCM) and the block spectral collocation method (BSCM). The boundary value method (BVM) and the block boundary value method (BBVM) are briefly discussed in Section 3. In Section 4, we present the implicit Runge–Kutta (IRK) method and the diagonally implicit Runge–Kutta (DIRK) method. The TVD Runge–Kutta (TVDRK) method is described in Section 5. In Section 6, we discuss the stability of our methods. In Section 7, we provide some extensive numerical results to assess the convergence and accuracy of our methods. Finally, we summarize the main features of our methods and briefly comment on the extension in Section 8. 2. The spectral collocation method In this section, we briefly introduce the spectral collocation method for solving the nonlinear system of second-order ODEs. We introduce the Chebyshev–Gauss–Lobatto points in Λ = [−1, 1],
x˜ j = cos
jπ
j = 0, 1 , · · · , N .
,
N
By differentiating the polynomial and evaluating the polynomial at the same collocation points with ˜f k = ˜f (˜xk ), we have N
F˜ N (˜x) =
˜f k L˜ k (˜x),
(2)
j =0
where L˜ k (˜x) are the Lagrange basis polynomials given by (see [5])
L˜ k (˜x) =
(−1)k+1 (1 − x˜ 2 ) T N (˜x) , ck N 2 (˜x − x˜ k )
k = 0, 1 , · · · , N ,
(3)
where T N (˜x) = cos( N cos−1 x˜ ) is the Chebyshev polynomial and
ck =
k = 0 or N , otherwise.
2, 1,
Let F˜ = [ ˜f (˜x0 ), · · · , ˜f (˜x N )] T , F˜ (m) = [ ˜f (m) (˜x0 ), · · · , ˜f (m) (˜x N )] T and approximate the derivative of F˜ at x˜ j by differentiating and evaluating (2), we get
F˜ N (˜x) = (m)
N
˜f k L˜ (m) (˜x), k
m = 1, 2, · · · .
(4)
j =0
Then (4) is equivalent to the following matrix equation
˜ (m) F˜ , F˜ (m) = D
m = 1, 2, · · · ,
˜ (m) is the (m + 1) × (m + 1) matrix whose entries are given by where D
˜ D = L˜ k (˜x j ), jk (m)
(m)
j , k = 0, 1 , · · · , N .
˜ (1) = D˜ = (d˜ kj ) is given by (see [1–5]) The first-order Chebyshev differentiation matrix D d˜ kj =
⎧ ck ⎪ ⎪ ⎨ − 2c j
(−1) j+k
π
π ,
sin((k+ j ) ) sin((k− j ) 2N ) 2N 1 kπ − cos ( )( 1 + cot2 ( kNπ )), ⎪ 2 N ⎪ 2 ⎩˜ d00 = −d˜ N N = 2N 6+1 .
k = j , k = j , k = 0, N ,
Higher derivative matrices can be obtained as matrix powers, i.e.,
˜ (m) = D˜ (1) D
m
.
a Let x j = a − b− (˜x j − 1) be the Chebyshev–Gauss–Lobatto points in [a, b], such that 2
x˜ j = 1 −
2 b−a
F (m) = D (m) F ,
(x j − a), m = 1, 2, · · · ,
(5)
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
where
T
F = f (x0 ), · · · , f (x N )
b−a
f (x j ) = f a −
2
,
˜ ˜ (˜x j − 1) = f (˜x j ) = f 1 −
T F (m) = f (m) (x0 ), · · · , f (m) (x N ) , and
D jk = L˜ k (m)
(m)
Then we have
1−
b−a
m
2
(m)
D jk = −
2
(x j − a) ,
L˜ k (˜x j ), (m)
b−a
2 b−a
237
(x j − a) ,
j , k = 0, 1 , · · · , N .
j , k = 0, 1 , · · · , N .
2 m ˜ (m) Therefore, D (m) = (− b− ) D is the Chebyshev differentiation matrix for the Chebyshev–Gauss–Lobatto points in [a, b]. a We now turn to the spectral collocation methods for solving the nonlinear system of the second-order ODEs. We first give the definition of the Kronecker product of two matrices.
Definition 2.1. (See [19].) Let A = (ai j )m×n and B be arbitrary matrices. Then the matrix
⎛
⎞ . . . a1n B . . . a2n B ⎟ ⎟ A⊗B =⎜ . . .. ⎟ , .. .. ⎝ .. . . ⎠ am1 B am2 B . . . amn B a11 B ⎜ a21 B ⎜
a12 B a22 B
is called the Kronecker product of A and B. Definition 2.2. (See [19].) Let A = (ai j )m×n . Then the vec( A ) is defined to be a column vector of size m × n made of the row of A stacked atop one another from left to right
vec( A ) = (a11 , a12 , . . . , a1n , a21 , a22 , . . . , am1 , . . . , amn ) T . Lemma 2.1. (See [19].) Let A = (ai j )m×n , B = (b i j )n× p , C = (c i j ) p ×q . Then
vec( A BC ) = A ⊗ C T vec( B ). Lemma 2.2. (See [19].) Let A = (ai j )m×1 , B = (b i j )n×1 . Then
vec A B T = A ⊗ B . Now we consider the following first-order nonlinear initial value problem
y = f t , y (t ) , y (t 0 ) = y 0 .
t0 < t ≤ T ,
(6) jπ
Let N t be a nonnegative integer and denote z j = cos( N ), 0 ≤ j ≤ N t the Chebyshev–Gauss–Lobatto points in the interval
[−1, 1], then t j = t 0 −
T −t 0 (z j 2
t
− 1) are the Chebyshev–Gauss–Lobatto points in the interval [t 0 , T ]. Discretizing (6) by the
spectral collocation method, we obtain
˜ Y˜ = F (Y ), A
(7)
˜ is the first-order Chebyshev differentiation matrix in the interval [t 0 , T ] given by where A
˜ =− A
2 T − t0
˜ (1) (2 : N t + 1, 1 : N t + 1), D
and
Y˜ = [ y 0 , Y ] T ,
T
F (Y ) = f t 1 , y (t 1 ) , f t 2 , y (t 2 ) , . . . , f t Nt , y (t Nt )
.
238
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
˜ = [a0 , A ] and Y˜ = [ y 0 , Y ] T , we can rewrite (7) by Using the partitions A AY + a0 y 0 = F (Y ).
(8)
If we consider the system of nonlinear ordinary differential equations
y = f t , y (t ) , y (t 0 ) = y 0 ,
where
t0 < t ≤ T ,
(9)
y (t ) = y 1 (t ), y 2 (t ), · · · , ym (t ) ,
f (t , y ) = f 1 (t , y ), f 2 (t , y ), · · · , f m (t , y ) . Then the spectral collocation method for (9) can be written as
AY + a0 y 0 = F (Y ), where
(10)
T
y = y (t 1 ), y (t 2 ), . . . , y (t Nt )
,
y (t j ) = y 1 (t j ), y 2 (t j ), · · · , ym (t j ) ,
T
F (Y ) = f t 1 , y (t 1 ) , f t 2 , y (t 2 ) , . . . , f t Nt , y (t Nt )
, f t j , y (t j ) = f 1 t j , y (t j ) , f 2 t j , y (t j ) , · · · , f m t j , y (t j ) .
Let us now consider solving the second-order initial value problems
y (t ) = f t , y (t ), y (t ) , t 0 < t ≤ T , y (t 0 ) = y 0 , y (t 0 ) = y 0 ,
where
y (t ) = y 1 (t ), y 2 (t ), · · · , ym (t ) ,
(11)
f t , y, y = f 1 t , y, y , f 2 t , y, y , · · · , fm t , y, y . By introducing the auxiliary variable q = y , we can rewrite (11) in the following form
y = q, q (t ) = f t , y (t ), q(t ) ,
(12)
with initial conditions y (t 0 ) = y 0 , q(t 0 ) = y 0 . By Eq. (10), we have
AY + a0 y 0 = Q ,
(13)
A Q + a0 q0 = F (Y , Q ).
(14)
Insert (13) into (14), we obtain
A ( AY + a0 y 0 ) + a0 y 0 = F (Y , AY + a0 y 0 ).
(15)
Using Lemma 2.1 and Lemma 2.2, we have
T A 2 ⊗ I m vec(Y ) = vec F (Y , AY + a0 y 0 ) − a0 ⊗ y 0 − ( Aa0 ) ⊗ y 0T ,
(16)
where I m is the m × m unitary matrix. Thus, we can calculate vec(Y ) from Eq. (16) with the Newton–Raphson iteration [20]. If we divide [t 0 , T ] to subintervals t 0 < t 1 < · · · < t Nt = T and apply the spectral collocation method (16) on each subinterval [t i −1 , t i ], then the resulted method is called the block spectral collocation method (BSCM). Let
˜ = D˜ (1) (2 : nt + 1, 2 : nt + 1), A ˜ (1) (2 : nt + 1, 1), a˜ 0 = D where nt + 1 is the number of the Chebyshev–Gauss–Lobatto points in the interval [tn , tn+1 ]. Then the spectral collocation method for solving (11) in [tn , tn+1 ] can be written as
2 2 ˜ − 2 a˜ 0 yn − 2 a˜ 0 yn = F Y , − 2 AY ˜ − 2 a˜ 0 yn , − A˜ − AY h
h
h
h
h
h
(17)
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
239
where h = tn+1 − tn . Eq. (17) is equivalent to
2 2 2 2 ˜ ˜ ˜ ˜ ˜ ˜ 4 A Y + 4 Aa0 yn − 2ha0 yn = h F Y , − AY − a0 yn . h
(18)
h
Using Lemma 2.1 and Lemma 2.2, we have
2 2 2 2 ˜ ˜ + 2ha˜ 0 ⊗ ynT −4 ( A˜ a˜ 0 ) ⊗ ynT . 4 A ⊗ I m vec(Y ) = h vec F Y , − AY − a˜ 0 yn h
(19)
h
Thus we can calculate vec(Y ) from Eq. (19) with Eq. (13) and the Newton–Raphson iteration [20]. Then both yn+1 and yn +1 can be obtained from the following form
yn+1 = Y (nt , :), 2 ˜ (nt , :)Y − 2 a˜ 0 (nt ) yn . yn +1 = − A h h
(20) (21)
3. The boundary value method In this section, we briefly introduce the boundary value methods (BVMs) for solving the nonlinear system of second-order ODEs which are based on the linear multistep methods (LMMs). We now consider the following first-order initial value problem
y = f t , y (t ) , y (t 0 ) = y 0 .
t0 < t ≤ T ,
(22)
A k-step formula of LMMs for approximating (22) can be written as k
αr yr + j = h
r =0
k
βr f r + j ,
j = 0, 1 , · · · , N t − k ,
(23)
r =0
where y r = y (tr ), f r = f (tr , y r ), tr = t 0 + rh and h = ( T − t 0 )/ N t . The BVMs (23) requires that final conditions are of the form k
αr yr + j = h
r =0
k
βr f r + j ,
γ initial conditions and k − γ
j = 0, 1 , · · · , N t − k ,
(24)
r =0
and k
αr( j) y Nt −k+r = h
r =0
k
( j)
βr f Nt −k+r ,
j = N t − k + γ + 1, · · · , N t ,
(25)
r =0
( j)
( j)
where the coefficients αr and βr are chosen such that truncation errors for the initial conditions and final conditions are of the same order as those for the basic formula (23). Eqs. (23)–(25) can be written as
AY + a0 y 0 = hB F (Y ) + hb0 f (t 0 , y 0 ), where [a0 | A ] = A e , [b0 | B ] = B e , A e , B e ∈ R f (t Nt , y (t Nt ))] T and
⎛
α0(1)
α1(1)
...
αk(1)
(26) N ×( N +1)
, Y = [ y (t 1 ), y (t 2 ), . . . , y (t Nt )] , F (Y ) = [ f (t 1 , y (t 1 )), f (t 2 , y (t 2 )), . . . , T
⎞
⎜ . ⎟ .. .. ⎜ . ⎟ . . ⎜ . ⎟ ⎜ ( γ −1 ) ⎟ ( γ − 1 ) ( γ − 1 ) ⎜α ⎟ α . . . α 1 k ⎜ 0 ⎟ ⎜ α0 ⎟ α . . . α 1 k ⎜ ⎟ ⎜ ⎟ α α . . . α 0 1 k ⎜ ⎟ Ae = ⎜ ⎟ .. .. ⎜ ⎟ . . ⎜ ⎟ ⎜ ⎟ α α . . . α 0 1 k ⎜ ⎟ ( N t −γ +1 ) ( N t −γ +1 ) ( N t −γ +1 ) ⎟ ⎜ α0 α1 . . . αk ⎜ ⎟ ⎜ ⎟ .. .. .. ⎜ ⎟ ⎝ ⎠ . . . Nt Nt Nt α0 α1 ... αk N Replacing
α by β yields the structure for B e .
.
t ×( N t +1)
240
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
We choose a particular fourth-order BVM approximation for (22) with k = 3 and please see [6,7]. The A, B, a0 and b0 in (26) can be written as follows
⎞
⎛
9/12 9/24 −1/24 ⎜ −9/12 9/12 1/12 ⎜ ⎜ −1/12 −9/12 9/12 ⎜
A=⎜
..
⎜ ⎜ ⎝
γ = 2. For the detailed procedure,
⎟ ⎟ ⎟ 1/12 ⎟ ⎟, .. .. .. ⎟ . . . ⎟ 1/12 ⎠ −1/12 −9/12 9/12 1/24 −9/24 −9/24 17/24 ⎞
.
⎛
3 /4 ⎟ ⎜ 1 /2 1 /2 ⎟ ⎜ ⎟ ⎜ . . .. .. B =⎜ ⎟, ⎟ ⎜ ⎠ ⎝ 1 /2 1/2 3 /4 1 /4 T T 17 1 1 a 0 = − , − , 0, · · · , 0 , b0 = , 0, · · · , 0 . 24
12
4
Let us now consider the following second-order initial value problem
y = f t , y (t ), y (t ) , t 0 < t ≤ T , y (t 0 ) = y 0 , y (t 0 ) = y 0 ,
where
y (t ) = y 1 (t ), y 2 (t ), · · · , ym (t ) ,
(27)
f t , y, y = f 1 t , y, y , f 2 t , y, y , · · · , fm t , y, y . We introduce a new variable q = y and rewrite (27) in the following form
y = q, q (t ) = f t , y (t ), q(t ) ,
(28)
with initial conditions y (t 0 ) = y 0 , q(t 0 ) = y 0 . The BVM for (28) can be written as
AY + a0 y 0 = hB Q + hb0 q0 ,
(29)
A Q + a0 q0 = hB F (Y , Q ) + hb0 f 0 .
(30)
Insert (29) into (30), we obtain
A B −1 AY + A B −1 a0 y 0 − h A B −1 b0 q0 + ha0 q0 = h2 B F Y , (hB )−1 ( AY + a0 y 0 − hb0 q0 ) + h2 b0 f 0 . Using Lemma 2.1 and Lemma 2.2, we obtain
(31)
A B −1 A ⊗ I m vec(Y ) − h2 B ⊗ I m vec F Y , (hB )−1 ( AY + a0 y 0 − hb0 q0 )
= − A B −1 a0 ⊗ y 0T + h A B −1 b0 − a0 ⊗ y 0T + h2 b0 ⊗ f 0T ,
(32)
where I m is the m × m unitary matrix and f 0 = f (t 0 , y (t 0 )). Thus, we can calculate vec(Y ) from Eq. (32) by the Newton–Raphson iteration [20]. If we divide [t 0 , T ] to subintervals t 0 < t 1 < · · · < t Nt = T and apply the boundary value method (32) on each subinterval [t i −1 , t i ], then the resulted method is called the block boundary value method (BBVM). Thus, the boundary value method for (27) in [tn , tn+1 ] can be written as
A B −1 AY + A B −1 a0 yn − h A B −1 b0 yn + ha0 yn = h2 B F Y , (hB )−1 AY + a0 yn − hb0 yn
+ h2 b0 f n ,
(33)
where h = (tn+1 , tn )/nt and nt + 1 is the number of grid points in the interval [tn , tn+1 ]. Using Lemma 2.1 and Lemma 2.2, we have
A B −1 A ⊗ I m vec(Y ) − h2 B ⊗ I m vec F Y , (hB )−1 ( AY + a0 yn − hb0 qn )
= − A B −1 a0 ⊗ ynT + h A B −1 b0 − a0 ⊗ ynT + h2 b0 ⊗ f nT .
(34)
Thus we can calculate vec(Y ) from Eq. (34) by Eq. (29) and the Newton–Raphson iteration [20]. Then yn+1 and yn +1 can be obtained from the following form
yn+1 = Y (nt , :),
−1
yn+1 = (hB )
AY + a0 yn − hb0 yn (nt , :).
(35) (36)
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
241
4. The implicit Runge–Kutta method and the diagonally implicit Runge–Kutta method In this section, we briefly introduce the implicit Runge–Kutta method and the diagonally implicit Runge–Kutta method for solving the second-order initial value problems. We now consider the following first-order initial value problem
y = f t , y (t ) , y (t 0 ) = y 0 .
t0 < t ≤ T ,
(37)
An s-stage formula of the implicit Runge–Kutta method for approximating (37) can be written as
y n +1 = y n + h
s
b i f (tn + c i h, ki ),
(38)
i =1
ki = yn + h
s
ai j f (tn + c j h, k j ),
i = 1, 2, · · · , s ,
(39)
j =1
where h = ( T − t 0 )/( N t − 1) is the integration step, ki (i = 1, 2, · · · , s) are the internal stages and tn = t 0 + hn. If ai j = 0 (i > j ), the method is the diagonally implicit Runge–Kutta method. A convenient way to represent the Runge–Kutta method is the Butcher table [9]
· · · a1s . .. . .. · · · ass · · · bs
c 1 a11
.. .
.. .
c s a11 b1
Then (38) and (39) can be written by the following form
K = Y n + h A F ( K ), yn+1 = yn + hB T F ( K ),
(40)
where B = [b1 , b2 , · · · , b s ] , K = [k1 , k2 , · · · , k s ] , Y n = [ yn , yn , · · · , yn ] and F ( K ) = [ f (tn + c 1 h, k1 ), f (tn + c 2 h, k2 ), · · · , f (tn + c s h, k s )] T . Let us now consider the following second-order initial value problem T
T
T
y = f t , y (t ), y (t ) , t 0 < t ≤ T , y (t 0 ) = y 0 , y (t 0 ) = y 0 ,
where
y (t ) = y 1 (t ), y 2 (t ), · · · , ym (t ) ,
(41)
f t , y, y = f 1 t , y, y , f 2 t , y, y , · · · , fm t , y, y . We introduce the auxiliary variable q = y and rewrite (41) in the following form
y = q, q (t ) = f t , y (t ), q(t ) ,
(42)
with initial conditions y (t 0 ) = y 0 , q(t 0 ) = y 0 . The s-stage formula of the Runge–Kutta method (40) for approximating (42) can be written as
yn+1 = yn + hB T P ,
(43)
T
qn+1 = qn + hB F ( K , P ),
(44)
K = Yn + h A P ,
(45)
P = Q n + h A F ( K , P ),
(46)
where Q n = [qn , qn , · · · , qn ] . By Eqs. (43)–(46), we have T
yn+1 = yn + hB T P , T
(47)
qn+1 = qn + hB F (Y n + h A P , P ),
(48)
P = Q n + h A F (Y n + h A P , P ).
(49)
242
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
Table 1 Implicit Runge–Kutta method. 1 2 1 2 1 2
− +
√
15 10
5 36 5 36 5 36 5 18
√
15 10
− +
2 9 2 9 2 9 4 9
√
15 √24 15 30
− +
√
5 36 5 36 5 36 5 18
15 15
√
15 15
− −
√
15
√30
15 24
Table 2 Diagonally implicit Runge–Kutta method. 1+λ 2 1 2 1−λ 2
1+λ 2
0
0
− λ2
1+λ 2
0
1+λ
−1 − 2λ
1
1+λ 2 1 6λ2
1
1 − 3λ2
6λ2
Using Lemma 2.1 and Lemma 2.2, we obtain
vec( P ) = vec( Q n ) + h A ⊗ I m vec F (Y n + h A P , P ) .
(50)
Thus, we can calculate vec( P ) from Eq. (50) with the Newton–Raphson iteration [20]. In Table 1, we choose a 3-stage implicit Runge–Kutta method. For details, please see [9–11]. In Table 2, we choose a 3-stage diagonally implicit Runge–Kutta method [9,11,12], where λ = √2 cos(10◦ ). 3
5. The TVD Runge–Kutta method In this section, we briefly introduce the TVD Runge–Kutta method for solving the second-order initial value problems. We now consider the following first-order initial value problem
y = f t , y (t ) , y (t 0 ) = y 0 .
t0 < t ≤ T ,
(51)
An s + 1-stage formula of the TVD Runge–Kutta method for approximating (51) can be written as
k0 = yn , ki =
i −1
(52)
ai j k j + hb i j f (tn + d j h, k j ) ,
i = 1, 2, · · · , s + 1,
(53)
j =0
y n +1 = k s +1 ,
(54)
where h = ( T − t 0 )/( N t − 1) is the integration step, ki (i = 1, 2, · · · , s + 1) are the internal stages and tn = t 0 + hn. Eqs. (52)–(54) are equivalent to
ki = ai0 yn + b i0 hf (tn + d0 h, yn ) +
i −1
ai j k j + b i j hf (tn + d j h, k j ) ,
i = 1, 2, · · · , s ,
(55)
j =1
yn+1 = as+1,0 yn + b s+1,0 hf (tn + d0 h, yn ) +
s
as+1, j k j + b s+1, j hf (tn + d j h, k j ) .
(56)
j =1
Then (55) and (56) are rewritten by the following matrix form
˜ yn + h B˜ f (tn + d0 h, yn ) + A K + hB F ( K ), K=A
(57)
ˆ K + h Bˆ F ( K ), yn+1 = as+1,0 yn + b s+1,0 hf (tn + d0 h, yn ) + A
(58)
˜ = [a10 , a20 , · · · , as0 ] T , B˜ = [b10 , b20 , · · · , b s0 ] T , where K = [k1 , k2 , · · · , k s ] T , F ( K ) = [ f (tn + d1 h, k1 ), · · · , f (tn + d s h, k s )] T , A ˆA = [as+1,1 , as+1,2 , · · · , as+1,s ], Bˆ = [b s+1,1 , b s+1,2 , · · · , b s+1,s ] and
⎛
0 ⎜ a21 ⎜
A=⎜ . .
⎝ .
as1
0 0
..
.
as2
⎞ ··· 0 ··· 0 ⎟ ⎟ , . . .. ⎟ . . ⎠ ··· 0 s×s
⎛
0 ⎜ b21 ⎜
B =⎜ . .
⎝ .
b s1
0 0
..
.
b s2
⎞ ··· 0 ··· 0 ⎟ ⎟ . . . .. ⎟ . . ⎠ ··· 0 s×s
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
243
Let us now consider the following second-order initial value problem
y = f t , y (t ), y (t ) , t 0 < t ≤ T , y (t 0 ) = y 0 , y (t 0 ) = y 0 .
(59)
We introduce the auxiliary variable q = y and rewrite (59) in the following form
y = q, q (t ) = f t , y (t ), q(t ) ,
(60)
with initial conditions y (t 0 ) = y 0 , q(t 0 ) = y 0 . The TVD Runge–Kutta method for (60) can be written as
˜ yn + h Bq ˜ n + A K + hB P , K=A
(61)
˜ n + h B˜ f (tn + d0 h, yn ) + A K + hB F ( K , P ), P = Aq
(62)
ˆ P + h Bˆ P , yn+1 = as+1,0 yn + b s+1,0 hqn + A
(63)
ˆ P + h Bˆ F ( K , P ). qn+1 = as+1,0 qn + b s+1,0 hf (tn + d0 h, yn ) + A
(64)
Eqs. (61)–(64) are equivalent to
˜ n + h B˜ f (tn + d0 h, yn , qn ) + A P + hB F ( I s − A )−1 ( A˜ yn + h Bq ˜ n + hB P ), P , P = Aq
ˆ ( I s − A )−1 A˜ yn + b s+1,0 + Aˆ ( I s − A )−1 B˜ hqn + h Bˆ + Aˆ ( I s − A )−1 B P , yn+1 = as+1,0 + A ˆ P + h Bˆ F ( I s − A )−1 ( A˜ yn + h Bq ˜ n + hB P ), P . qn+1 = as+1,0 qn + b s+1,0 hf (tn + d0 h, yn , qn ) + A
Denote
A=
˜ A as+1,0
A ˆ A
,
B=
B˜ b s+1,0
B Bˆ
(65) (66) (67)
.
We choose a third-order TVDRK approximation for (59) with s = 2. For details, please see [13]. Then A, B and d can be written as follows
⎞
⎛
1 ⎠, A = ⎝ 3 /4 1 /4 1 /3 0 2/3
⎞
⎛
1 B = ⎝ 0 1 /4 0 0
⎠,
d = [0, 1, 1/2] T .
2/ 3
Then Eqs. (65)–(67) become
p 1 = qn + hf (tn , yn , qn ), 1 3 1 p 2 = qn + p 1 + h f (tn + h, yn + hqn , p 1 ), 4 4 4 1 1 2 y n + 1 = y n + h qn + h p 1 + h p 2 , 6 6 3 2 1 1 1 2 1 qn+1 = qn + p 2 + h f tn + h, yn + h qn + h p 1 , p 2 . 3 3 3 2 4 4
(68)
6. Stability In this section, we study the stability regions of the methods discussed in this paper, that is the five methods including BSCM, BBVM, IRK method, DIRK method and TVDRK method. The stability properties of our methods are analyzed by employing the second-order homogeneous linear test equation
y (t ) = −λ2 y (t ),
λ > 0.
(69)
Using the BSCM for Eq. (69), we have
˜ 2 Y + 4 A˜ a˜ 0 yn − 2ha˜ 0 yn = −λ2 h2 Y . 4A By Eqs. (20) and (21), we obtain
(nt , :) 2ha˜ 0 yn − 4 A˜ a˜ 0 yn , 2 ˜ (nt , :) 4 A˜ 2 + λ2 h2 I nt −1 2ha˜ 0 yn − 4 A˜ a˜ 0 yn − 2 a˜ 0 (nt ) yn , yn +1 = − A
˜ 2 + λ2 h2 I nt y n +1 = 4 A h
− 1
h
244
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
so that
˜ (nt , :) 4 A˜ 2 + λ2 h2 I nt hyn +1 = −2 A Then we have
y n +1 hyn +1
= S ( z)
yn hyn
− 1
˜ a˜ 0 yn − 2a˜ 0 (nt ) yn . 2ha˜ 0 yn − 4 A
(70)
,
where the stability matrix S ( z) is determined by
S ( z) =
˜ 2 + zI nt )−1 (nt , :)˜a0 2(4 A −4(4 A˜ 2 + zI nt )−1 (nt , :) A˜ a˜ 0 2 − 1˜˜ ˜ ˜ ˜ ˜ 2 + zI nt )−1 a˜ 0 8 A (nt , :)(4 A + zI nt ) A a0 − 2a˜ 0 (nt ) −4 A (nt , :)(4 A
,
z = λ2 h2 .
Using the BBVM for Eq. (69), we have
A B −1 AY + A B −1 a0 yn − h A B −1 b0 yn + ha0 yn = −λ2 h2 B Y − λ2 h2 b0 yn . Thus we obtain
Y = A B −1 A + zB
−1
− A B −1 a0 yn + h A B −1 b0 yn − ha0 yn − zb0 yn ,
where z = λ2 h2 . By Eq. (29), we have
Q = (hB )−1 AY + a0 yn − hb0 yn . Then
hyn +1 = h Q (nt ).
yn+1 = Y (nt ), Therefore
y n +1 hyn +1
= S ( z)
yn hyn
(71)
,
where the stability matrix S ( z) is determined by
S 11 = − R 2 (nt , :)a0 − zR 1 (nt , :)b0 , S 12 = R 2 (nt , :)b0 − R 1 (nt , :)a0 , S 21 = − R 4 (nt , :)a0 − zR 3 (nt , :)b0 + B −1 (nt , :)a0 , S 22 = R 4 (nt , :)b0 − R 3 (nt , :)a0 − B −1 (nt , :)b0 ,
R 1 = A B −1 A + zB R 3 = B −1 A R 1 ,
− 1
R 2 = R 1 A B −1 ,
,
R 4 = B −1 A R 2 .
Using the IRK method or the DIRK method for Eq. (69), we have
yn+1 = yn + hB T P , qn+1 = qn − λ2 hB T ( ˆI s yn + h A P ), P = ˆI s qn − λ2 h A ( ˆI s yn + h A P ), where ˆI s = [1, 1, · · · , 1] T . Then we have
s
yn+1 = 1 − zB T R −1 A ˆI s yn + hB T R −1 ˆI s qn ,
hqn+1 = − zB T ˆI s − z A R −1 A ˆI s yn + 1 − zB T A R −1 ˆI s hqn , where R = ( I s + λ2 h2 A 2 ) and z = λ2 h2 . Therefore, we obtain
y n +1 hyn +1
= S ( z)
yn hyn
(72)
,
where the stability matrix S ( z) is determined by
S ( z) =
1 − zB T R −1 A ˆI s − zB T ( ˆI s − z A R −1 A ˆI s )
B T R −1 ˆI s 1 − zB T A R −1 ˆI s
.
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
245
Using the TVDRK method for Eq. (69), we have
˜ n − hλ2 B˜ yn + A P − hλ2 B ( I s − A )−1 ( A˜ yn + h Bq ˜ n + hB P ), P = Aq
ˆ ( I s − A )−1 A˜ yn + b s+1,0 + Aˆ ( I s − A )−1 B˜ hqn + h Bˆ + Aˆ ( I s − A )−1 B P , yn+1 = as+1,0 + A ˆ P − hλ2 Bˆ ( I s − A )−1 ( A˜ yn + h Bq ˜ n + hB P ). qn+1 = as+1,0 qn − b s+1,0 hλ2 yn + A Thus, we obtain
˜ n − z B˜ yn − zB ( I s − A )−1 ( A˜ yn + h Bq ˜ n) , h P = R −1 Ahq
ˆ ( I s − A )−1 A˜ yn + b s+1,0 + Aˆ ( I s − A )−1 B˜ hqn + Bˆ + Aˆ ( I s − A )−1 B h P , yn+1 = as+1,0 + A ˆ P − z Bˆ ( I s − A )−1 ( A˜ yn + h Bq ˜ n + hB P ), hqn+1 = as+1,0 hqn − b s+1,0 zyn + Ah where R = I s − A + h2 λ2 B ( I s − A )−1 B and z = λ2 h2 . Therefore, we have
y n +1 hyn +1
= S ( z)
yn hyn
(73)
,
where the stability matrix S ( z) is determined by
ˆ ( I s − A )−1 A˜ − z Bˆ + Aˆ ( I s − A )−1 B R −1 B˜ + B ( I s − A )−1 A˜ , S 11 = as+1,0 + A
ˆ ( I s − A )−1 B˜ + Bˆ + Aˆ ( I s − A )−1 B R −1 A˜ − zB ( I s − A )−1 B˜ , S 12 = b s+1,0 + A
˜ − z Aˆ − z Bˆ ( I s − A )−1 B R −1 B˜ + B ( I s − A )−1 A˜ , S 21 = − zb s+1,0 − z Bˆ ( I s − A )−1 A
ˆ − z Bˆ ( I s − A )−1 B R −1 A˜ − zB ( I s − A )−1 B˜ . S 22 = as+1,0 − z Bˆ ( I s − A )−1 B˜ + A The stability is characterized by the spectral radius
ρ ( S (z)). We now give the definition of the stability.
Definition 6.1. (See [16,21,22,24].) (i) (ii) (iii) (iv) (v)
R s = { z ∈ R| z > 0 and ρ ( S ) < 1} is called the stability region. R p = { z ∈ R| z > 0, ρ ( S ) = 1 and tr( S )2 < 4 det( S )} is called the periodicity region. R q = { z ∈ R| z > 0 and det( S ) < 1, | tr( S )| < det( S ) + 1} is called the strong stability region. If R s = { z ∈ (0, +∞)}, the method is called A-stable. If R p = { z ∈ (0, +∞)}, the method is called P -stable.
Fig. 1 shows the spectral radius ρ ( S ) for different z ∈ [0, 90] and the comparison with the method of [23] and the method of [24, Eq. (4.4)]. Fig. 2 shows the spectral radius ρ ( S ) for different z ∈ [0, 400] and different nt (left panel for BSCM and right panel for BBVM). In Fig. 3, we show the tr( S )2 − 4 det( S ) for different z ∈ [0, 400] and different nt (left panel for BSCM and right panel for BBVM). Fig. 4 shows det( S ) and | tr( S )| − det( S ) − 1 for different z ∈ [0, 400] and different nt for BSCM. In Fig. 5, we show ρ ( S ), tr( S )2 − 4 det( S ), det( S ) and | tr( S )| − det( S ) − 1 for different z ∈ [0, 400] for IRK and DIRK methods. By these figures, we can determine the stability of the methods. 7. Numerical results In this section, we present some numerical examples to demonstrate the convergence and accuracy of the methods. The implementation is written in Matlab Version 7.1. 7.1. Test problem 1 We consider the problem (1) with f (t , y , y ) = −100 y + 1+ y 21+ y 2 + Φ(t ), t 0 = 0, T = 4,
Φ(t ) = 100 sin(t ) − sin(t ) −
1
(cos(10t ) + sin(10t ) + sin(t ))2 + (cos(t ) + 10 cos(10t ) − 10 sin(10t ))2 + 1
,
and initial conditions
y (0) = 1,
y (0) = ν + 1.
The exact solution is given by
y (t ) = cos(10t ) + sin(10t ) + sin(t ). This test problem is given in [23]. In Table 3, we show the errors and computational orders obtained for Test problem 1. In Table 4, we show the errors and computational orders obtained for Test problem 1 by using the spectral collocation method (SCM).
246
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
Fig. 1. Spectral radius
Fig. 2. Spectral radius
ρ ( S ) for different z ∈ [0, 90].
ρ ( S ) for different z ∈ [0, 400] and different nt (left panel for BSCM and right panel for BBVM).
Fig. 3. tr( S )2 − 4 det( S ) for different z ∈ [0, 400] and different nt (left panel for BSCM and right panel for BBVM).
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
247
Fig. 4. det( S ) and | tr( S )| − det( S ) − 1 for different z ∈ [0, 400] and different nt for BSCM.
Fig. 5.
ρ ( S ), tr( S )2 − 4 det( S ), det( S ) and | tr( S )| − det( S ) − 1 for different z ∈ [0, 400] for IRK and DIRK methods.
Table 3 History of convergence for Test problem 1. Nt
40 60 80 100 120 140
BVM
DIRK
IRK
TVDRK
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
4.59e−01 9.24e−02 2.97e−02 1.17e−02 5.71e−03 3.14e−03
−
1.20e+00 7.09e−01 3.44e−01 1.74e−01 9.34e−02 5.27e−02
−
6.38e−04 5.44e−05 9.43e−06 2.38e−06 8.24e−07 3.12e−07
−
1.36e+00 5.53e−01 2.68e−01 1.43e−01 8.35e−02 5.29e−02
2.21 2.51 2.82 2.95 2.96
3.95 3.94 4.18 3.94 3.88
1.29 2.51 3.07 3.40 3.71
6.07 6.09 6.17 5.82 6.30
Table 4 History of convergence for Test problem 1 by using the spectral collocation method. Nt
20
25
30
35
40
L ∞ error Order
1.10e+00
3.32e−02 15.68
7.08e−04 21.10
2.03e−06 37.99
4.33e−09 46.04
−
−
248
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
Table 5 History of convergence for Test problem 2. BSCM (nt = 6)
Nt
20 40 60 80 100 120
BBVM (nt = 6)
L ∞ error
order
L ∞ error
order
1.09e−04 7.55e−07 4.28e−08 5.68e−09 1.20e−09 3.35e−10
−
2.79e−04 8.34e−06 1.35e−06 3.92e−07 1.54e−07 7.26e−08
5.06 4.50 4.29 4.19 4.13
7.18 7.08 7.02 6.96 7.01
−
Table 6 History of convergence for Test problem 2. Nt
DIRK
200 300 400 500 600 700
IRK
TVDRK
Method of [23]
Method of [24]
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
8.22e−03 9.99e−04 2.29e−04 7.34e−05 2.90e−05 1.32e−05
−
3.85e−09 3.35e−10 5.94e−11 1.56e−11 5.20e−12 2.02e−12
−
1.54e−01 4.78e−02 2.04e−02 1.05e−02 6.07e−03 3.82e−03
−
1.04e−01 4.27e−02 2.30e−02 1.43e−02 9.74e−03 7.06e−03
−
3.20e−03 7.76e−04 3.46e−04 1.95e−04 1.26e−04 8.81e−05
−
5.20 5.12 5.10 5.09 5.10
6.02 6.01 5.98 6.04 6.13
2.89 2.96 2.98 2.99 3.00
2.20 2.15 2.12 2.10 2.09
3.49 2.81 2.57 2.41 2.31
7.2. Test problem 2 We consider the problem (1) with f (t , y , y ) = [−
0 0
y (0) =
y (0) =
,
0 1
y1
( y 21 + y 22 )3/2
,−
y2
( y 21 + y 22 )3/2
] T , t 0 = 0, T = 10 and initial conditions
.
The exact solution is given by
y (t ) =
cos(t ) sin(t )
.
This test problem is given in [23]. In Table 5 and Table 6, we show the errors and computational orders obtained for Test problem 2. 7.3. Test problem 3 We consider the problem (1) with t 0 = 0, T = 10,
f t , y, y = −
13 −12 −12 13
and initial conditions
1 0
y (0) =
y (0) =
,
y+
−4
y (t ) =
,
8
The exact solution is given by
9 cos(2t ) − 12 sin(2t ) −12 cos(2t ) + 9 sin(2t )
sin(t ) − sin(5t ) + cos(2t ) cos(t ) + sin(5t ) + sin(2t )
. .
This test problem is given in [25]. In Table 7 and Table 8, we show the errors and computational orders obtained for Test problem 3. 7.4. Test problem 4 We consider the problem (1) with t 0 = 0, T = 10,
f t , y, y
=−
101 2 − 99 2
− 99 2 101 2
y+ε
93 2 93 2
cos(2t ) − sin(2t ) −
99 sin(2t ) 2 99 cos(2t ) 2
,
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
249
Table 7 History of convergence for Test problem 3. BSCM (nt = 6)
Nt
100 200 300 400 500 600
BBVM (nt = 6)
L ∞ error
order
L ∞ error
order
1.21e−05 1.81e−07 1.57e−08 2.80e−09 7.33e−10 2.54e−10
−
4.05e−04 2.83e−05 5.61e−06 1.77e−06 7.28e−07 3.50e−07
3.84 3.99 4.01 3.98 4.02
6.06 6.04 5.99 6.01 5.81
−
Table 8 History of convergence for Test problem 3. Nt
DIRK
400 600 800 1000 1200 1400
IRK
TVDRK
Method of [23]
Method of [24]
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
5.68e−02 1.21e−02 3.84e−03 1.57e−03 7.64e−04 4.13e−04
−
2.38e−07 2.08e−08 3.71e−09 9.72e−10 3.25e−10 1.29e−10
−
6.26e−02 1.90e−02 8.03e−03 4.12e−03 2.39e−03 1.50e−03
−
3.24e−01 1.43e−01 8.05e−02 5.15e−02 3.57e−02 2.62e−02
−
3.15e−02 6.87e−03 2.53e−03 1.29e−03 7.97e−04 5.49e−04
3.76 3.47 3.01 2.65 2.41
3.82 3.98 4.00 3.97 3.98
6.00 6.00 6.00 6.01 6.01
2.95 2.98 2.99 3.00 3.00
2.01 2.00 2.00 2.00 2.00
−
Table 9 History of convergence for Test problem 4. BSCM (nt = 6)
Nt
100 200 300 400 500 600
BBVM (nt = 6)
L ∞ error
order
L ∞ error
order
2.03e−03 3.16e−05 2.90e−06 5.19e−07 1.36e−07 4.53e−08
−
1.58e−02 1.22e−03 2.48e−04 7.96e−05 3.30e−05 1.58e−05
3.69 3.93 3.95 3.94 4.04
6.01 5.89 5.98 6.01 6.02
−
Table 10 History of convergence for Test problem 4. Nt
DIRK
400 600 800 1000 1200 1400
IRK
TVDRK
Method of [23]
Method of [24]
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
1.08e+00 4.22e−01 1.58e−01 6.86e−02 3.41e−02 1.87e−02
−
4.30e−05 3.78e−06 6.75e−07 1.78e−07 5.95e−08 2.36e−08
−
9.40e−01 3.76e−01 1.72e−01 9.08e−02 5.33e−02 3.38e−02
−
3.76e+00 1.68e+00 9.46e−01 5.95e−01 4.09e−01 3.01e−01
−
1.15e+00 2.58e−01 8.40e−02 3.43e−02 1.62e−02 8.58e−03
3.68 3.90 4.01 4.10 4.14
2.32 3.42 3.73 3.83 3.90
and initial conditions
y (0) =
−1 + ε 1
,
y (0) =
6.00 5.99 5.98 6.00 6.00
−10 10 + 2ε
The exact solution is given by
y (t ) =
− cos(10t ) − sin(10t ) + ε cos(2t ) cos(10t ) + sin(10t ) + ε sin(2t )
2.26 2.73 2.85 2.92 2.96
1.99 2.00 2.08 2.05 2.00
−
. ,
where ε = 0.01. This test problem is given in [25]. In Table 9 and Table 10, we show the errors and computational orders obtained for Test problem 4. 7.5. Test problem 5 We consider the following one-dimensional sine-Gordon equation
d2 u dt 2
=
d2 u dx2
− sin(u ),
−1 < x < 1, 0 < t ≤ 30,
(74)
250
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
Fig. 6. Curve plot of the absolute errors for Test problem 5 at t = 30.
with initial condition −1
√
u (x, 0) = 4 tan
1 − ω2
1
√
cosh(x 1 − ω2 )
ω
,
and Dirichlet boundary conditions −1
√
u (−1, t ) = 4 tan
u (1, t ) = 4 tan−1
√
1 − ω2
cos ωt
√
cosh(−1 1 − ω2 )
ω
cos ωt
ω
√
1 − ω2
cos ωt
√
cosh(x 1 − ω2 )
ω
,
. √ cosh( 1 − ω2 )
1 − ω2
The exact solution is given by
u (x, t ) = 4 tan−1
ut (x, 0) = 0,
.
We discretize this problem in space by the standard finite differences with N x + 1 grid points.
⎧ 2 ⎨ d uj
=
u j +1 − 2u j + u j −1
2 x2 ⎩ udt(x, 0) = u (x , 0 ), j j
− sin(u j ),
0
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
Some numerical algorithms for solving the highly oscillatory second-order initial value problems ✩ Wenjie Liu, Boying Wu, Jiebao Sun ∗ Department of Mathematics, Harbin Institute of Technology, Harbin 150001, PR China
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 9 January 2014 Received in revised form 16 May 2014 Accepted 15 July 2014 Available online 28 July 2014 Keywords: Highly oscillatory problems Second-order initial value problems Block spectral collocation method Block boundary value method Implicit Runge–Kutta method Diagonally implicit Runge–Kutta method Total variation diminishing Runge–Kutta method
In this paper, some numerical algorithms (spectral collocation method, block spectral collocation method, boundary value method, block boundary value method, implicit Runge–Kutta method, diagonally implicit Runge–Kutta method and total variation diminishing Runge–Kutta method) are used to solve the highly oscillatory second-order initial value problems. We first derive these methods for the first-order initial value problems, and then extend these methods to the highly oscillatory nonlinear systems by matrix analysis methods. These new methods preserve the accuracy of the original methods and the main advantages of these new methods are low storage requirements and high efficiency. Extensive numerical results are presented to demonstrate the convergence properties of these methods. © 2014 Elsevier Inc. All rights reserved.
1. Introduction In this paper, we study the spectral collocation method (or the block spectral collocation method) [1–5], the boundary value method (or the block boundary value method) [6–8], the implicit Runge–Kutta method [9–11], the diagonally implicit Runge–Kutta method [9,11,12] and the total variation diminishing (TVD) Runge–Kutta method [13,14] for the following highly oscillatory second-order initial value problem
y (t ) = f t , y (t ), y (t ) , t 0 < t ≤ T , y (t 0 ) = y 0 , y (t 0 ) = y 0 .
(1)
The problem (1) is the fundamental model in physics for describing mechanical vibrations behavior [15]. This model also appears in circuit simulations, flexible body dynamics and various quantum dynamics calculations [16–18]. Therefore, the construction of efficient numerical schemes for solving (1) is an important task. Recently, several numerical schemes have been developed for (1), and some useful approaches to construct Runge–Kutta– Nystro¨ m (RKN) type methods have also been proposed. For details we refer to the monograph [16]. The block spectral collocation method (or the spectral collocation method), the block boundary value method (or the boundary value method), the implicit Runge–Kutta method, the diagonally implicit Runge–Kutta method and the TVD
✩ This work is partially supported by the National Science Foundation of China (11271100, 11301113) and the Fundamental Research Funds for the Central Universities (Grant No. HIT. IBRSEM. A. 201412), Harbin Science and Technology Innovative Talents Project of Special Fund (2013RFXYJ044). Corresponding author. E-mail address: [email protected] (J. Sun).
*
http://dx.doi.org/10.1016/j.jcp.2014.07.033 0021-9991/© 2014 Elsevier Inc. All rights reserved.
236
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
Runge–Kutta method are high-accuracy schemes for the first-order nonlinear initial value problems. So, we extend these methods for solving the second-order initial value problems by matrix analysis methods [19]. The main advantages of these new methods are low storage requirements and high efficiency. The outline of this paper is as follows: in Section 2, we describe the spectral collocation method (SCM) and the block spectral collocation method (BSCM). The boundary value method (BVM) and the block boundary value method (BBVM) are briefly discussed in Section 3. In Section 4, we present the implicit Runge–Kutta (IRK) method and the diagonally implicit Runge–Kutta (DIRK) method. The TVD Runge–Kutta (TVDRK) method is described in Section 5. In Section 6, we discuss the stability of our methods. In Section 7, we provide some extensive numerical results to assess the convergence and accuracy of our methods. Finally, we summarize the main features of our methods and briefly comment on the extension in Section 8. 2. The spectral collocation method In this section, we briefly introduce the spectral collocation method for solving the nonlinear system of second-order ODEs. We introduce the Chebyshev–Gauss–Lobatto points in Λ = [−1, 1],
x˜ j = cos
jπ
j = 0, 1 , · · · , N .
,
N
By differentiating the polynomial and evaluating the polynomial at the same collocation points with ˜f k = ˜f (˜xk ), we have N
F˜ N (˜x) =
˜f k L˜ k (˜x),
(2)
j =0
where L˜ k (˜x) are the Lagrange basis polynomials given by (see [5])
L˜ k (˜x) =
(−1)k+1 (1 − x˜ 2 ) T N (˜x) , ck N 2 (˜x − x˜ k )
k = 0, 1 , · · · , N ,
(3)
where T N (˜x) = cos( N cos−1 x˜ ) is the Chebyshev polynomial and
ck =
k = 0 or N , otherwise.
2, 1,
Let F˜ = [ ˜f (˜x0 ), · · · , ˜f (˜x N )] T , F˜ (m) = [ ˜f (m) (˜x0 ), · · · , ˜f (m) (˜x N )] T and approximate the derivative of F˜ at x˜ j by differentiating and evaluating (2), we get
F˜ N (˜x) = (m)
N
˜f k L˜ (m) (˜x), k
m = 1, 2, · · · .
(4)
j =0
Then (4) is equivalent to the following matrix equation
˜ (m) F˜ , F˜ (m) = D
m = 1, 2, · · · ,
˜ (m) is the (m + 1) × (m + 1) matrix whose entries are given by where D
˜ D = L˜ k (˜x j ), jk (m)
(m)
j , k = 0, 1 , · · · , N .
˜ (1) = D˜ = (d˜ kj ) is given by (see [1–5]) The first-order Chebyshev differentiation matrix D d˜ kj =
⎧ ck ⎪ ⎪ ⎨ − 2c j
(−1) j+k
π
π ,
sin((k+ j ) ) sin((k− j ) 2N ) 2N 1 kπ − cos ( )( 1 + cot2 ( kNπ )), ⎪ 2 N ⎪ 2 ⎩˜ d00 = −d˜ N N = 2N 6+1 .
k = j , k = j , k = 0, N ,
Higher derivative matrices can be obtained as matrix powers, i.e.,
˜ (m) = D˜ (1) D
m
.
a Let x j = a − b− (˜x j − 1) be the Chebyshev–Gauss–Lobatto points in [a, b], such that 2
x˜ j = 1 −
2 b−a
F (m) = D (m) F ,
(x j − a), m = 1, 2, · · · ,
(5)
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
where
T
F = f (x0 ), · · · , f (x N )
b−a
f (x j ) = f a −
2
,
˜ ˜ (˜x j − 1) = f (˜x j ) = f 1 −
T F (m) = f (m) (x0 ), · · · , f (m) (x N ) , and
D jk = L˜ k (m)
(m)
Then we have
1−
b−a
m
2
(m)
D jk = −
2
(x j − a) ,
L˜ k (˜x j ), (m)
b−a
2 b−a
237
(x j − a) ,
j , k = 0, 1 , · · · , N .
j , k = 0, 1 , · · · , N .
2 m ˜ (m) Therefore, D (m) = (− b− ) D is the Chebyshev differentiation matrix for the Chebyshev–Gauss–Lobatto points in [a, b]. a We now turn to the spectral collocation methods for solving the nonlinear system of the second-order ODEs. We first give the definition of the Kronecker product of two matrices.
Definition 2.1. (See [19].) Let A = (ai j )m×n and B be arbitrary matrices. Then the matrix
⎛
⎞ . . . a1n B . . . a2n B ⎟ ⎟ A⊗B =⎜ . . .. ⎟ , .. .. ⎝ .. . . ⎠ am1 B am2 B . . . amn B a11 B ⎜ a21 B ⎜
a12 B a22 B
is called the Kronecker product of A and B. Definition 2.2. (See [19].) Let A = (ai j )m×n . Then the vec( A ) is defined to be a column vector of size m × n made of the row of A stacked atop one another from left to right
vec( A ) = (a11 , a12 , . . . , a1n , a21 , a22 , . . . , am1 , . . . , amn ) T . Lemma 2.1. (See [19].) Let A = (ai j )m×n , B = (b i j )n× p , C = (c i j ) p ×q . Then
vec( A BC ) = A ⊗ C T vec( B ). Lemma 2.2. (See [19].) Let A = (ai j )m×1 , B = (b i j )n×1 . Then
vec A B T = A ⊗ B . Now we consider the following first-order nonlinear initial value problem
y = f t , y (t ) , y (t 0 ) = y 0 .
t0 < t ≤ T ,
(6) jπ
Let N t be a nonnegative integer and denote z j = cos( N ), 0 ≤ j ≤ N t the Chebyshev–Gauss–Lobatto points in the interval
[−1, 1], then t j = t 0 −
T −t 0 (z j 2
t
− 1) are the Chebyshev–Gauss–Lobatto points in the interval [t 0 , T ]. Discretizing (6) by the
spectral collocation method, we obtain
˜ Y˜ = F (Y ), A
(7)
˜ is the first-order Chebyshev differentiation matrix in the interval [t 0 , T ] given by where A
˜ =− A
2 T − t0
˜ (1) (2 : N t + 1, 1 : N t + 1), D
and
Y˜ = [ y 0 , Y ] T ,
T
F (Y ) = f t 1 , y (t 1 ) , f t 2 , y (t 2 ) , . . . , f t Nt , y (t Nt )
.
238
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
˜ = [a0 , A ] and Y˜ = [ y 0 , Y ] T , we can rewrite (7) by Using the partitions A AY + a0 y 0 = F (Y ).
(8)
If we consider the system of nonlinear ordinary differential equations
y = f t , y (t ) , y (t 0 ) = y 0 ,
where
t0 < t ≤ T ,
(9)
y (t ) = y 1 (t ), y 2 (t ), · · · , ym (t ) ,
f (t , y ) = f 1 (t , y ), f 2 (t , y ), · · · , f m (t , y ) . Then the spectral collocation method for (9) can be written as
AY + a0 y 0 = F (Y ), where
(10)
T
y = y (t 1 ), y (t 2 ), . . . , y (t Nt )
,
y (t j ) = y 1 (t j ), y 2 (t j ), · · · , ym (t j ) ,
T
F (Y ) = f t 1 , y (t 1 ) , f t 2 , y (t 2 ) , . . . , f t Nt , y (t Nt )
, f t j , y (t j ) = f 1 t j , y (t j ) , f 2 t j , y (t j ) , · · · , f m t j , y (t j ) .
Let us now consider solving the second-order initial value problems
y (t ) = f t , y (t ), y (t ) , t 0 < t ≤ T , y (t 0 ) = y 0 , y (t 0 ) = y 0 ,
where
y (t ) = y 1 (t ), y 2 (t ), · · · , ym (t ) ,
(11)
f t , y, y = f 1 t , y, y , f 2 t , y, y , · · · , fm t , y, y . By introducing the auxiliary variable q = y , we can rewrite (11) in the following form
y = q, q (t ) = f t , y (t ), q(t ) ,
(12)
with initial conditions y (t 0 ) = y 0 , q(t 0 ) = y 0 . By Eq. (10), we have
AY + a0 y 0 = Q ,
(13)
A Q + a0 q0 = F (Y , Q ).
(14)
Insert (13) into (14), we obtain
A ( AY + a0 y 0 ) + a0 y 0 = F (Y , AY + a0 y 0 ).
(15)
Using Lemma 2.1 and Lemma 2.2, we have
T A 2 ⊗ I m vec(Y ) = vec F (Y , AY + a0 y 0 ) − a0 ⊗ y 0 − ( Aa0 ) ⊗ y 0T ,
(16)
where I m is the m × m unitary matrix. Thus, we can calculate vec(Y ) from Eq. (16) with the Newton–Raphson iteration [20]. If we divide [t 0 , T ] to subintervals t 0 < t 1 < · · · < t Nt = T and apply the spectral collocation method (16) on each subinterval [t i −1 , t i ], then the resulted method is called the block spectral collocation method (BSCM). Let
˜ = D˜ (1) (2 : nt + 1, 2 : nt + 1), A ˜ (1) (2 : nt + 1, 1), a˜ 0 = D where nt + 1 is the number of the Chebyshev–Gauss–Lobatto points in the interval [tn , tn+1 ]. Then the spectral collocation method for solving (11) in [tn , tn+1 ] can be written as
2 2 ˜ − 2 a˜ 0 yn − 2 a˜ 0 yn = F Y , − 2 AY ˜ − 2 a˜ 0 yn , − A˜ − AY h
h
h
h
h
h
(17)
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
239
where h = tn+1 − tn . Eq. (17) is equivalent to
2 2 2 2 ˜ ˜ ˜ ˜ ˜ ˜ 4 A Y + 4 Aa0 yn − 2ha0 yn = h F Y , − AY − a0 yn . h
(18)
h
Using Lemma 2.1 and Lemma 2.2, we have
2 2 2 2 ˜ ˜ + 2ha˜ 0 ⊗ ynT −4 ( A˜ a˜ 0 ) ⊗ ynT . 4 A ⊗ I m vec(Y ) = h vec F Y , − AY − a˜ 0 yn h
(19)
h
Thus we can calculate vec(Y ) from Eq. (19) with Eq. (13) and the Newton–Raphson iteration [20]. Then both yn+1 and yn +1 can be obtained from the following form
yn+1 = Y (nt , :), 2 ˜ (nt , :)Y − 2 a˜ 0 (nt ) yn . yn +1 = − A h h
(20) (21)
3. The boundary value method In this section, we briefly introduce the boundary value methods (BVMs) for solving the nonlinear system of second-order ODEs which are based on the linear multistep methods (LMMs). We now consider the following first-order initial value problem
y = f t , y (t ) , y (t 0 ) = y 0 .
t0 < t ≤ T ,
(22)
A k-step formula of LMMs for approximating (22) can be written as k
αr yr + j = h
r =0
k
βr f r + j ,
j = 0, 1 , · · · , N t − k ,
(23)
r =0
where y r = y (tr ), f r = f (tr , y r ), tr = t 0 + rh and h = ( T − t 0 )/ N t . The BVMs (23) requires that final conditions are of the form k
αr yr + j = h
r =0
k
βr f r + j ,
γ initial conditions and k − γ
j = 0, 1 , · · · , N t − k ,
(24)
r =0
and k
αr( j) y Nt −k+r = h
r =0
k
( j)
βr f Nt −k+r ,
j = N t − k + γ + 1, · · · , N t ,
(25)
r =0
( j)
( j)
where the coefficients αr and βr are chosen such that truncation errors for the initial conditions and final conditions are of the same order as those for the basic formula (23). Eqs. (23)–(25) can be written as
AY + a0 y 0 = hB F (Y ) + hb0 f (t 0 , y 0 ), where [a0 | A ] = A e , [b0 | B ] = B e , A e , B e ∈ R f (t Nt , y (t Nt ))] T and
⎛
α0(1)
α1(1)
...
αk(1)
(26) N ×( N +1)
, Y = [ y (t 1 ), y (t 2 ), . . . , y (t Nt )] , F (Y ) = [ f (t 1 , y (t 1 )), f (t 2 , y (t 2 )), . . . , T
⎞
⎜ . ⎟ .. .. ⎜ . ⎟ . . ⎜ . ⎟ ⎜ ( γ −1 ) ⎟ ( γ − 1 ) ( γ − 1 ) ⎜α ⎟ α . . . α 1 k ⎜ 0 ⎟ ⎜ α0 ⎟ α . . . α 1 k ⎜ ⎟ ⎜ ⎟ α α . . . α 0 1 k ⎜ ⎟ Ae = ⎜ ⎟ .. .. ⎜ ⎟ . . ⎜ ⎟ ⎜ ⎟ α α . . . α 0 1 k ⎜ ⎟ ( N t −γ +1 ) ( N t −γ +1 ) ( N t −γ +1 ) ⎟ ⎜ α0 α1 . . . αk ⎜ ⎟ ⎜ ⎟ .. .. .. ⎜ ⎟ ⎝ ⎠ . . . Nt Nt Nt α0 α1 ... αk N Replacing
α by β yields the structure for B e .
.
t ×( N t +1)
240
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
We choose a particular fourth-order BVM approximation for (22) with k = 3 and please see [6,7]. The A, B, a0 and b0 in (26) can be written as follows
⎞
⎛
9/12 9/24 −1/24 ⎜ −9/12 9/12 1/12 ⎜ ⎜ −1/12 −9/12 9/12 ⎜
A=⎜
..
⎜ ⎜ ⎝
γ = 2. For the detailed procedure,
⎟ ⎟ ⎟ 1/12 ⎟ ⎟, .. .. .. ⎟ . . . ⎟ 1/12 ⎠ −1/12 −9/12 9/12 1/24 −9/24 −9/24 17/24 ⎞
.
⎛
3 /4 ⎟ ⎜ 1 /2 1 /2 ⎟ ⎜ ⎟ ⎜ . . .. .. B =⎜ ⎟, ⎟ ⎜ ⎠ ⎝ 1 /2 1/2 3 /4 1 /4 T T 17 1 1 a 0 = − , − , 0, · · · , 0 , b0 = , 0, · · · , 0 . 24
12
4
Let us now consider the following second-order initial value problem
y = f t , y (t ), y (t ) , t 0 < t ≤ T , y (t 0 ) = y 0 , y (t 0 ) = y 0 ,
where
y (t ) = y 1 (t ), y 2 (t ), · · · , ym (t ) ,
(27)
f t , y, y = f 1 t , y, y , f 2 t , y, y , · · · , fm t , y, y . We introduce a new variable q = y and rewrite (27) in the following form
y = q, q (t ) = f t , y (t ), q(t ) ,
(28)
with initial conditions y (t 0 ) = y 0 , q(t 0 ) = y 0 . The BVM for (28) can be written as
AY + a0 y 0 = hB Q + hb0 q0 ,
(29)
A Q + a0 q0 = hB F (Y , Q ) + hb0 f 0 .
(30)
Insert (29) into (30), we obtain
A B −1 AY + A B −1 a0 y 0 − h A B −1 b0 q0 + ha0 q0 = h2 B F Y , (hB )−1 ( AY + a0 y 0 − hb0 q0 ) + h2 b0 f 0 . Using Lemma 2.1 and Lemma 2.2, we obtain
(31)
A B −1 A ⊗ I m vec(Y ) − h2 B ⊗ I m vec F Y , (hB )−1 ( AY + a0 y 0 − hb0 q0 )
= − A B −1 a0 ⊗ y 0T + h A B −1 b0 − a0 ⊗ y 0T + h2 b0 ⊗ f 0T ,
(32)
where I m is the m × m unitary matrix and f 0 = f (t 0 , y (t 0 )). Thus, we can calculate vec(Y ) from Eq. (32) by the Newton–Raphson iteration [20]. If we divide [t 0 , T ] to subintervals t 0 < t 1 < · · · < t Nt = T and apply the boundary value method (32) on each subinterval [t i −1 , t i ], then the resulted method is called the block boundary value method (BBVM). Thus, the boundary value method for (27) in [tn , tn+1 ] can be written as
A B −1 AY + A B −1 a0 yn − h A B −1 b0 yn + ha0 yn = h2 B F Y , (hB )−1 AY + a0 yn − hb0 yn
+ h2 b0 f n ,
(33)
where h = (tn+1 , tn )/nt and nt + 1 is the number of grid points in the interval [tn , tn+1 ]. Using Lemma 2.1 and Lemma 2.2, we have
A B −1 A ⊗ I m vec(Y ) − h2 B ⊗ I m vec F Y , (hB )−1 ( AY + a0 yn − hb0 qn )
= − A B −1 a0 ⊗ ynT + h A B −1 b0 − a0 ⊗ ynT + h2 b0 ⊗ f nT .
(34)
Thus we can calculate vec(Y ) from Eq. (34) by Eq. (29) and the Newton–Raphson iteration [20]. Then yn+1 and yn +1 can be obtained from the following form
yn+1 = Y (nt , :),
−1
yn+1 = (hB )
AY + a0 yn − hb0 yn (nt , :).
(35) (36)
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
241
4. The implicit Runge–Kutta method and the diagonally implicit Runge–Kutta method In this section, we briefly introduce the implicit Runge–Kutta method and the diagonally implicit Runge–Kutta method for solving the second-order initial value problems. We now consider the following first-order initial value problem
y = f t , y (t ) , y (t 0 ) = y 0 .
t0 < t ≤ T ,
(37)
An s-stage formula of the implicit Runge–Kutta method for approximating (37) can be written as
y n +1 = y n + h
s
b i f (tn + c i h, ki ),
(38)
i =1
ki = yn + h
s
ai j f (tn + c j h, k j ),
i = 1, 2, · · · , s ,
(39)
j =1
where h = ( T − t 0 )/( N t − 1) is the integration step, ki (i = 1, 2, · · · , s) are the internal stages and tn = t 0 + hn. If ai j = 0 (i > j ), the method is the diagonally implicit Runge–Kutta method. A convenient way to represent the Runge–Kutta method is the Butcher table [9]
· · · a1s . .. . .. · · · ass · · · bs
c 1 a11
.. .
.. .
c s a11 b1
Then (38) and (39) can be written by the following form
K = Y n + h A F ( K ), yn+1 = yn + hB T F ( K ),
(40)
where B = [b1 , b2 , · · · , b s ] , K = [k1 , k2 , · · · , k s ] , Y n = [ yn , yn , · · · , yn ] and F ( K ) = [ f (tn + c 1 h, k1 ), f (tn + c 2 h, k2 ), · · · , f (tn + c s h, k s )] T . Let us now consider the following second-order initial value problem T
T
T
y = f t , y (t ), y (t ) , t 0 < t ≤ T , y (t 0 ) = y 0 , y (t 0 ) = y 0 ,
where
y (t ) = y 1 (t ), y 2 (t ), · · · , ym (t ) ,
(41)
f t , y, y = f 1 t , y, y , f 2 t , y, y , · · · , fm t , y, y . We introduce the auxiliary variable q = y and rewrite (41) in the following form
y = q, q (t ) = f t , y (t ), q(t ) ,
(42)
with initial conditions y (t 0 ) = y 0 , q(t 0 ) = y 0 . The s-stage formula of the Runge–Kutta method (40) for approximating (42) can be written as
yn+1 = yn + hB T P ,
(43)
T
qn+1 = qn + hB F ( K , P ),
(44)
K = Yn + h A P ,
(45)
P = Q n + h A F ( K , P ),
(46)
where Q n = [qn , qn , · · · , qn ] . By Eqs. (43)–(46), we have T
yn+1 = yn + hB T P , T
(47)
qn+1 = qn + hB F (Y n + h A P , P ),
(48)
P = Q n + h A F (Y n + h A P , P ).
(49)
242
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
Table 1 Implicit Runge–Kutta method. 1 2 1 2 1 2
− +
√
15 10
5 36 5 36 5 36 5 18
√
15 10
− +
2 9 2 9 2 9 4 9
√
15 √24 15 30
− +
√
5 36 5 36 5 36 5 18
15 15
√
15 15
− −
√
15
√30
15 24
Table 2 Diagonally implicit Runge–Kutta method. 1+λ 2 1 2 1−λ 2
1+λ 2
0
0
− λ2
1+λ 2
0
1+λ
−1 − 2λ
1
1+λ 2 1 6λ2
1
1 − 3λ2
6λ2
Using Lemma 2.1 and Lemma 2.2, we obtain
vec( P ) = vec( Q n ) + h A ⊗ I m vec F (Y n + h A P , P ) .
(50)
Thus, we can calculate vec( P ) from Eq. (50) with the Newton–Raphson iteration [20]. In Table 1, we choose a 3-stage implicit Runge–Kutta method. For details, please see [9–11]. In Table 2, we choose a 3-stage diagonally implicit Runge–Kutta method [9,11,12], where λ = √2 cos(10◦ ). 3
5. The TVD Runge–Kutta method In this section, we briefly introduce the TVD Runge–Kutta method for solving the second-order initial value problems. We now consider the following first-order initial value problem
y = f t , y (t ) , y (t 0 ) = y 0 .
t0 < t ≤ T ,
(51)
An s + 1-stage formula of the TVD Runge–Kutta method for approximating (51) can be written as
k0 = yn , ki =
i −1
(52)
ai j k j + hb i j f (tn + d j h, k j ) ,
i = 1, 2, · · · , s + 1,
(53)
j =0
y n +1 = k s +1 ,
(54)
where h = ( T − t 0 )/( N t − 1) is the integration step, ki (i = 1, 2, · · · , s + 1) are the internal stages and tn = t 0 + hn. Eqs. (52)–(54) are equivalent to
ki = ai0 yn + b i0 hf (tn + d0 h, yn ) +
i −1
ai j k j + b i j hf (tn + d j h, k j ) ,
i = 1, 2, · · · , s ,
(55)
j =1
yn+1 = as+1,0 yn + b s+1,0 hf (tn + d0 h, yn ) +
s
as+1, j k j + b s+1, j hf (tn + d j h, k j ) .
(56)
j =1
Then (55) and (56) are rewritten by the following matrix form
˜ yn + h B˜ f (tn + d0 h, yn ) + A K + hB F ( K ), K=A
(57)
ˆ K + h Bˆ F ( K ), yn+1 = as+1,0 yn + b s+1,0 hf (tn + d0 h, yn ) + A
(58)
˜ = [a10 , a20 , · · · , as0 ] T , B˜ = [b10 , b20 , · · · , b s0 ] T , where K = [k1 , k2 , · · · , k s ] T , F ( K ) = [ f (tn + d1 h, k1 ), · · · , f (tn + d s h, k s )] T , A ˆA = [as+1,1 , as+1,2 , · · · , as+1,s ], Bˆ = [b s+1,1 , b s+1,2 , · · · , b s+1,s ] and
⎛
0 ⎜ a21 ⎜
A=⎜ . .
⎝ .
as1
0 0
..
.
as2
⎞ ··· 0 ··· 0 ⎟ ⎟ , . . .. ⎟ . . ⎠ ··· 0 s×s
⎛
0 ⎜ b21 ⎜
B =⎜ . .
⎝ .
b s1
0 0
..
.
b s2
⎞ ··· 0 ··· 0 ⎟ ⎟ . . . .. ⎟ . . ⎠ ··· 0 s×s
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
243
Let us now consider the following second-order initial value problem
y = f t , y (t ), y (t ) , t 0 < t ≤ T , y (t 0 ) = y 0 , y (t 0 ) = y 0 .
(59)
We introduce the auxiliary variable q = y and rewrite (59) in the following form
y = q, q (t ) = f t , y (t ), q(t ) ,
(60)
with initial conditions y (t 0 ) = y 0 , q(t 0 ) = y 0 . The TVD Runge–Kutta method for (60) can be written as
˜ yn + h Bq ˜ n + A K + hB P , K=A
(61)
˜ n + h B˜ f (tn + d0 h, yn ) + A K + hB F ( K , P ), P = Aq
(62)
ˆ P + h Bˆ P , yn+1 = as+1,0 yn + b s+1,0 hqn + A
(63)
ˆ P + h Bˆ F ( K , P ). qn+1 = as+1,0 qn + b s+1,0 hf (tn + d0 h, yn ) + A
(64)
Eqs. (61)–(64) are equivalent to
˜ n + h B˜ f (tn + d0 h, yn , qn ) + A P + hB F ( I s − A )−1 ( A˜ yn + h Bq ˜ n + hB P ), P , P = Aq
ˆ ( I s − A )−1 A˜ yn + b s+1,0 + Aˆ ( I s − A )−1 B˜ hqn + h Bˆ + Aˆ ( I s − A )−1 B P , yn+1 = as+1,0 + A ˆ P + h Bˆ F ( I s − A )−1 ( A˜ yn + h Bq ˜ n + hB P ), P . qn+1 = as+1,0 qn + b s+1,0 hf (tn + d0 h, yn , qn ) + A
Denote
A=
˜ A as+1,0
A ˆ A
,
B=
B˜ b s+1,0
B Bˆ
(65) (66) (67)
.
We choose a third-order TVDRK approximation for (59) with s = 2. For details, please see [13]. Then A, B and d can be written as follows
⎞
⎛
1 ⎠, A = ⎝ 3 /4 1 /4 1 /3 0 2/3
⎞
⎛
1 B = ⎝ 0 1 /4 0 0
⎠,
d = [0, 1, 1/2] T .
2/ 3
Then Eqs. (65)–(67) become
p 1 = qn + hf (tn , yn , qn ), 1 3 1 p 2 = qn + p 1 + h f (tn + h, yn + hqn , p 1 ), 4 4 4 1 1 2 y n + 1 = y n + h qn + h p 1 + h p 2 , 6 6 3 2 1 1 1 2 1 qn+1 = qn + p 2 + h f tn + h, yn + h qn + h p 1 , p 2 . 3 3 3 2 4 4
(68)
6. Stability In this section, we study the stability regions of the methods discussed in this paper, that is the five methods including BSCM, BBVM, IRK method, DIRK method and TVDRK method. The stability properties of our methods are analyzed by employing the second-order homogeneous linear test equation
y (t ) = −λ2 y (t ),
λ > 0.
(69)
Using the BSCM for Eq. (69), we have
˜ 2 Y + 4 A˜ a˜ 0 yn − 2ha˜ 0 yn = −λ2 h2 Y . 4A By Eqs. (20) and (21), we obtain
(nt , :) 2ha˜ 0 yn − 4 A˜ a˜ 0 yn , 2 ˜ (nt , :) 4 A˜ 2 + λ2 h2 I nt −1 2ha˜ 0 yn − 4 A˜ a˜ 0 yn − 2 a˜ 0 (nt ) yn , yn +1 = − A
˜ 2 + λ2 h2 I nt y n +1 = 4 A h
− 1
h
244
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
so that
˜ (nt , :) 4 A˜ 2 + λ2 h2 I nt hyn +1 = −2 A Then we have
y n +1 hyn +1
= S ( z)
yn hyn
− 1
˜ a˜ 0 yn − 2a˜ 0 (nt ) yn . 2ha˜ 0 yn − 4 A
(70)
,
where the stability matrix S ( z) is determined by
S ( z) =
˜ 2 + zI nt )−1 (nt , :)˜a0 2(4 A −4(4 A˜ 2 + zI nt )−1 (nt , :) A˜ a˜ 0 2 − 1˜˜ ˜ ˜ ˜ ˜ 2 + zI nt )−1 a˜ 0 8 A (nt , :)(4 A + zI nt ) A a0 − 2a˜ 0 (nt ) −4 A (nt , :)(4 A
,
z = λ2 h2 .
Using the BBVM for Eq. (69), we have
A B −1 AY + A B −1 a0 yn − h A B −1 b0 yn + ha0 yn = −λ2 h2 B Y − λ2 h2 b0 yn . Thus we obtain
Y = A B −1 A + zB
−1
− A B −1 a0 yn + h A B −1 b0 yn − ha0 yn − zb0 yn ,
where z = λ2 h2 . By Eq. (29), we have
Q = (hB )−1 AY + a0 yn − hb0 yn . Then
hyn +1 = h Q (nt ).
yn+1 = Y (nt ), Therefore
y n +1 hyn +1
= S ( z)
yn hyn
(71)
,
where the stability matrix S ( z) is determined by
S 11 = − R 2 (nt , :)a0 − zR 1 (nt , :)b0 , S 12 = R 2 (nt , :)b0 − R 1 (nt , :)a0 , S 21 = − R 4 (nt , :)a0 − zR 3 (nt , :)b0 + B −1 (nt , :)a0 , S 22 = R 4 (nt , :)b0 − R 3 (nt , :)a0 − B −1 (nt , :)b0 ,
R 1 = A B −1 A + zB R 3 = B −1 A R 1 ,
− 1
R 2 = R 1 A B −1 ,
,
R 4 = B −1 A R 2 .
Using the IRK method or the DIRK method for Eq. (69), we have
yn+1 = yn + hB T P , qn+1 = qn − λ2 hB T ( ˆI s yn + h A P ), P = ˆI s qn − λ2 h A ( ˆI s yn + h A P ), where ˆI s = [1, 1, · · · , 1] T . Then we have
s
yn+1 = 1 − zB T R −1 A ˆI s yn + hB T R −1 ˆI s qn ,
hqn+1 = − zB T ˆI s − z A R −1 A ˆI s yn + 1 − zB T A R −1 ˆI s hqn , where R = ( I s + λ2 h2 A 2 ) and z = λ2 h2 . Therefore, we obtain
y n +1 hyn +1
= S ( z)
yn hyn
(72)
,
where the stability matrix S ( z) is determined by
S ( z) =
1 − zB T R −1 A ˆI s − zB T ( ˆI s − z A R −1 A ˆI s )
B T R −1 ˆI s 1 − zB T A R −1 ˆI s
.
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
245
Using the TVDRK method for Eq. (69), we have
˜ n − hλ2 B˜ yn + A P − hλ2 B ( I s − A )−1 ( A˜ yn + h Bq ˜ n + hB P ), P = Aq
ˆ ( I s − A )−1 A˜ yn + b s+1,0 + Aˆ ( I s − A )−1 B˜ hqn + h Bˆ + Aˆ ( I s − A )−1 B P , yn+1 = as+1,0 + A ˆ P − hλ2 Bˆ ( I s − A )−1 ( A˜ yn + h Bq ˜ n + hB P ). qn+1 = as+1,0 qn − b s+1,0 hλ2 yn + A Thus, we obtain
˜ n − z B˜ yn − zB ( I s − A )−1 ( A˜ yn + h Bq ˜ n) , h P = R −1 Ahq
ˆ ( I s − A )−1 A˜ yn + b s+1,0 + Aˆ ( I s − A )−1 B˜ hqn + Bˆ + Aˆ ( I s − A )−1 B h P , yn+1 = as+1,0 + A ˆ P − z Bˆ ( I s − A )−1 ( A˜ yn + h Bq ˜ n + hB P ), hqn+1 = as+1,0 hqn − b s+1,0 zyn + Ah where R = I s − A + h2 λ2 B ( I s − A )−1 B and z = λ2 h2 . Therefore, we have
y n +1 hyn +1
= S ( z)
yn hyn
(73)
,
where the stability matrix S ( z) is determined by
ˆ ( I s − A )−1 A˜ − z Bˆ + Aˆ ( I s − A )−1 B R −1 B˜ + B ( I s − A )−1 A˜ , S 11 = as+1,0 + A
ˆ ( I s − A )−1 B˜ + Bˆ + Aˆ ( I s − A )−1 B R −1 A˜ − zB ( I s − A )−1 B˜ , S 12 = b s+1,0 + A
˜ − z Aˆ − z Bˆ ( I s − A )−1 B R −1 B˜ + B ( I s − A )−1 A˜ , S 21 = − zb s+1,0 − z Bˆ ( I s − A )−1 A
ˆ − z Bˆ ( I s − A )−1 B R −1 A˜ − zB ( I s − A )−1 B˜ . S 22 = as+1,0 − z Bˆ ( I s − A )−1 B˜ + A The stability is characterized by the spectral radius
ρ ( S (z)). We now give the definition of the stability.
Definition 6.1. (See [16,21,22,24].) (i) (ii) (iii) (iv) (v)
R s = { z ∈ R| z > 0 and ρ ( S ) < 1} is called the stability region. R p = { z ∈ R| z > 0, ρ ( S ) = 1 and tr( S )2 < 4 det( S )} is called the periodicity region. R q = { z ∈ R| z > 0 and det( S ) < 1, | tr( S )| < det( S ) + 1} is called the strong stability region. If R s = { z ∈ (0, +∞)}, the method is called A-stable. If R p = { z ∈ (0, +∞)}, the method is called P -stable.
Fig. 1 shows the spectral radius ρ ( S ) for different z ∈ [0, 90] and the comparison with the method of [23] and the method of [24, Eq. (4.4)]. Fig. 2 shows the spectral radius ρ ( S ) for different z ∈ [0, 400] and different nt (left panel for BSCM and right panel for BBVM). In Fig. 3, we show the tr( S )2 − 4 det( S ) for different z ∈ [0, 400] and different nt (left panel for BSCM and right panel for BBVM). Fig. 4 shows det( S ) and | tr( S )| − det( S ) − 1 for different z ∈ [0, 400] and different nt for BSCM. In Fig. 5, we show ρ ( S ), tr( S )2 − 4 det( S ), det( S ) and | tr( S )| − det( S ) − 1 for different z ∈ [0, 400] for IRK and DIRK methods. By these figures, we can determine the stability of the methods. 7. Numerical results In this section, we present some numerical examples to demonstrate the convergence and accuracy of the methods. The implementation is written in Matlab Version 7.1. 7.1. Test problem 1 We consider the problem (1) with f (t , y , y ) = −100 y + 1+ y 21+ y 2 + Φ(t ), t 0 = 0, T = 4,
Φ(t ) = 100 sin(t ) − sin(t ) −
1
(cos(10t ) + sin(10t ) + sin(t ))2 + (cos(t ) + 10 cos(10t ) − 10 sin(10t ))2 + 1
,
and initial conditions
y (0) = 1,
y (0) = ν + 1.
The exact solution is given by
y (t ) = cos(10t ) + sin(10t ) + sin(t ). This test problem is given in [23]. In Table 3, we show the errors and computational orders obtained for Test problem 1. In Table 4, we show the errors and computational orders obtained for Test problem 1 by using the spectral collocation method (SCM).
246
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
Fig. 1. Spectral radius
Fig. 2. Spectral radius
ρ ( S ) for different z ∈ [0, 90].
ρ ( S ) for different z ∈ [0, 400] and different nt (left panel for BSCM and right panel for BBVM).
Fig. 3. tr( S )2 − 4 det( S ) for different z ∈ [0, 400] and different nt (left panel for BSCM and right panel for BBVM).
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
247
Fig. 4. det( S ) and | tr( S )| − det( S ) − 1 for different z ∈ [0, 400] and different nt for BSCM.
Fig. 5.
ρ ( S ), tr( S )2 − 4 det( S ), det( S ) and | tr( S )| − det( S ) − 1 for different z ∈ [0, 400] for IRK and DIRK methods.
Table 3 History of convergence for Test problem 1. Nt
40 60 80 100 120 140
BVM
DIRK
IRK
TVDRK
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
4.59e−01 9.24e−02 2.97e−02 1.17e−02 5.71e−03 3.14e−03
−
1.20e+00 7.09e−01 3.44e−01 1.74e−01 9.34e−02 5.27e−02
−
6.38e−04 5.44e−05 9.43e−06 2.38e−06 8.24e−07 3.12e−07
−
1.36e+00 5.53e−01 2.68e−01 1.43e−01 8.35e−02 5.29e−02
2.21 2.51 2.82 2.95 2.96
3.95 3.94 4.18 3.94 3.88
1.29 2.51 3.07 3.40 3.71
6.07 6.09 6.17 5.82 6.30
Table 4 History of convergence for Test problem 1 by using the spectral collocation method. Nt
20
25
30
35
40
L ∞ error Order
1.10e+00
3.32e−02 15.68
7.08e−04 21.10
2.03e−06 37.99
4.33e−09 46.04
−
−
248
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
Table 5 History of convergence for Test problem 2. BSCM (nt = 6)
Nt
20 40 60 80 100 120
BBVM (nt = 6)
L ∞ error
order
L ∞ error
order
1.09e−04 7.55e−07 4.28e−08 5.68e−09 1.20e−09 3.35e−10
−
2.79e−04 8.34e−06 1.35e−06 3.92e−07 1.54e−07 7.26e−08
5.06 4.50 4.29 4.19 4.13
7.18 7.08 7.02 6.96 7.01
−
Table 6 History of convergence for Test problem 2. Nt
DIRK
200 300 400 500 600 700
IRK
TVDRK
Method of [23]
Method of [24]
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
8.22e−03 9.99e−04 2.29e−04 7.34e−05 2.90e−05 1.32e−05
−
3.85e−09 3.35e−10 5.94e−11 1.56e−11 5.20e−12 2.02e−12
−
1.54e−01 4.78e−02 2.04e−02 1.05e−02 6.07e−03 3.82e−03
−
1.04e−01 4.27e−02 2.30e−02 1.43e−02 9.74e−03 7.06e−03
−
3.20e−03 7.76e−04 3.46e−04 1.95e−04 1.26e−04 8.81e−05
−
5.20 5.12 5.10 5.09 5.10
6.02 6.01 5.98 6.04 6.13
2.89 2.96 2.98 2.99 3.00
2.20 2.15 2.12 2.10 2.09
3.49 2.81 2.57 2.41 2.31
7.2. Test problem 2 We consider the problem (1) with f (t , y , y ) = [−
0 0
y (0) =
y (0) =
,
0 1
y1
( y 21 + y 22 )3/2
,−
y2
( y 21 + y 22 )3/2
] T , t 0 = 0, T = 10 and initial conditions
.
The exact solution is given by
y (t ) =
cos(t ) sin(t )
.
This test problem is given in [23]. In Table 5 and Table 6, we show the errors and computational orders obtained for Test problem 2. 7.3. Test problem 3 We consider the problem (1) with t 0 = 0, T = 10,
f t , y, y = −
13 −12 −12 13
and initial conditions
1 0
y (0) =
y (0) =
,
y+
−4
y (t ) =
,
8
The exact solution is given by
9 cos(2t ) − 12 sin(2t ) −12 cos(2t ) + 9 sin(2t )
sin(t ) − sin(5t ) + cos(2t ) cos(t ) + sin(5t ) + sin(2t )
. .
This test problem is given in [25]. In Table 7 and Table 8, we show the errors and computational orders obtained for Test problem 3. 7.4. Test problem 4 We consider the problem (1) with t 0 = 0, T = 10,
f t , y, y
=−
101 2 − 99 2
− 99 2 101 2
y+ε
93 2 93 2
cos(2t ) − sin(2t ) −
99 sin(2t ) 2 99 cos(2t ) 2
,
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
249
Table 7 History of convergence for Test problem 3. BSCM (nt = 6)
Nt
100 200 300 400 500 600
BBVM (nt = 6)
L ∞ error
order
L ∞ error
order
1.21e−05 1.81e−07 1.57e−08 2.80e−09 7.33e−10 2.54e−10
−
4.05e−04 2.83e−05 5.61e−06 1.77e−06 7.28e−07 3.50e−07
3.84 3.99 4.01 3.98 4.02
6.06 6.04 5.99 6.01 5.81
−
Table 8 History of convergence for Test problem 3. Nt
DIRK
400 600 800 1000 1200 1400
IRK
TVDRK
Method of [23]
Method of [24]
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
5.68e−02 1.21e−02 3.84e−03 1.57e−03 7.64e−04 4.13e−04
−
2.38e−07 2.08e−08 3.71e−09 9.72e−10 3.25e−10 1.29e−10
−
6.26e−02 1.90e−02 8.03e−03 4.12e−03 2.39e−03 1.50e−03
−
3.24e−01 1.43e−01 8.05e−02 5.15e−02 3.57e−02 2.62e−02
−
3.15e−02 6.87e−03 2.53e−03 1.29e−03 7.97e−04 5.49e−04
3.76 3.47 3.01 2.65 2.41
3.82 3.98 4.00 3.97 3.98
6.00 6.00 6.00 6.01 6.01
2.95 2.98 2.99 3.00 3.00
2.01 2.00 2.00 2.00 2.00
−
Table 9 History of convergence for Test problem 4. BSCM (nt = 6)
Nt
100 200 300 400 500 600
BBVM (nt = 6)
L ∞ error
order
L ∞ error
order
2.03e−03 3.16e−05 2.90e−06 5.19e−07 1.36e−07 4.53e−08
−
1.58e−02 1.22e−03 2.48e−04 7.96e−05 3.30e−05 1.58e−05
3.69 3.93 3.95 3.94 4.04
6.01 5.89 5.98 6.01 6.02
−
Table 10 History of convergence for Test problem 4. Nt
DIRK
400 600 800 1000 1200 1400
IRK
TVDRK
Method of [23]
Method of [24]
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
L ∞ error
order
1.08e+00 4.22e−01 1.58e−01 6.86e−02 3.41e−02 1.87e−02
−
4.30e−05 3.78e−06 6.75e−07 1.78e−07 5.95e−08 2.36e−08
−
9.40e−01 3.76e−01 1.72e−01 9.08e−02 5.33e−02 3.38e−02
−
3.76e+00 1.68e+00 9.46e−01 5.95e−01 4.09e−01 3.01e−01
−
1.15e+00 2.58e−01 8.40e−02 3.43e−02 1.62e−02 8.58e−03
3.68 3.90 4.01 4.10 4.14
2.32 3.42 3.73 3.83 3.90
and initial conditions
y (0) =
−1 + ε 1
,
y (0) =
6.00 5.99 5.98 6.00 6.00
−10 10 + 2ε
The exact solution is given by
y (t ) =
− cos(10t ) − sin(10t ) + ε cos(2t ) cos(10t ) + sin(10t ) + ε sin(2t )
2.26 2.73 2.85 2.92 2.96
1.99 2.00 2.08 2.05 2.00
−
. ,
where ε = 0.01. This test problem is given in [25]. In Table 9 and Table 10, we show the errors and computational orders obtained for Test problem 4. 7.5. Test problem 5 We consider the following one-dimensional sine-Gordon equation
d2 u dt 2
=
d2 u dx2
− sin(u ),
−1 < x < 1, 0 < t ≤ 30,
(74)
250
W. Liu et al. / Journal of Computational Physics 276 (2014) 235–251
Fig. 6. Curve plot of the absolute errors for Test problem 5 at t = 30.
with initial condition −1
√
u (x, 0) = 4 tan
1 − ω2
1
√
cosh(x 1 − ω2 )
ω
,
and Dirichlet boundary conditions −1
√
u (−1, t ) = 4 tan
u (1, t ) = 4 tan−1
√
1 − ω2
cos ωt
√
cosh(−1 1 − ω2 )
ω
cos ωt
ω
√
1 − ω2
cos ωt
√
cosh(x 1 − ω2 )
ω
,
. √ cosh( 1 − ω2 )
1 − ω2
The exact solution is given by
u (x, t ) = 4 tan−1
ut (x, 0) = 0,
.
We discretize this problem in space by the standard finite differences with N x + 1 grid points.
⎧ 2 ⎨ d uj
=
u j +1 − 2u j + u j −1
2 x2 ⎩ udt(x, 0) = u (x , 0 ), j j
− sin(u j ),
0
