- Email: [email protected]
Unconditionally stable difference schemes for a one-space-dimensional linear hyperbolic equation
Applied Mathematics and Computation 187 (2007) 1272–1276 www.elsevier.com/locate/amc
Unconditionally stable difference schemes for a one-space-dimensional linear hyperbolic equation Feng Gao
a,*
, Chunmei Chi
b
a b
Science School, Qingdao Technological University, Qingdao 266033, PR China Computer School, Qingdao Technological University, Qingdao 266033, PR China
Abstract A few explicit difference schemes are discussed for the numerical solution of the linear hyperbolic equation utt + 2a ut + b2u = uxx + f(x, t), a > 0, b > 0, in the region X = {(x, t)ja < x < b, t > 0} subject to appropriate initial and Dirichlet boundary conditions, where a and b are real numbers. The proposed scheme is showed to be unconditionally stable, and numerical result is presented. 2006 Elsevier Inc. All rights reserved. Keywords: Second-order linear hyperbolic equation; Damped wave equation; Telegraph equation; Explicit scheme; Unconditionally stable; Pade approximation
1. Introduction In recent years, much attention has been given in the literature to the development, analysis, and implementation of stable methods for the numerical solution of second-order hyperbolic equations, see, for example, [1– 3]. These methods are conditionally stable. Recently, Mohanty et al. [4,5] have developed new three-level implicit unconditionally stable alternating direction implicit schemes for the two and three-space-dimensional linear hyperbolic equations. These schemes are second-order accurate both in space and time. In [6], Mohanty carried over a new technique to the following linear one-space-dimensional hyperbolic equation with given initial and Dirichlet boundary conditions o2 u ou o2 u þ 2a þ b2 u ¼ 2 þ f ðx; tÞ; 2 ot ot ox
a > 0; b P 0;
over a region X = [a < x < b] · [t > 0],with the initial conditions uðx; 0Þ ¼ /ðxÞ;
*
ouðx; 0Þ ¼ wðxÞ; ot
Corresponding author. E-mail address: [email protected] (F. Gao).
0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.09.057
F. Gao, C. Chi / Applied Mathematics and Computation 187 (2007) 1272–1276
1273
and an unconditionally stable implicit three level scheme was obtained. The above equation represents a damped wave equation and a telegraph equation respectively, see [1]. In this paper, however, we will show that if the above problem can be reduced to the following one, o2 u ou o2 u 2 þ b þ 2a u ¼ þ f ðxÞ; ot2 ot ox2
ð1Þ
a > 0; b P 0;
over a region X = [a < x < b] · [t > 0], with the initial conditions uðx; 0Þ ¼ /ðxÞ;
ouðx; 0Þ ¼ wðxÞ ot
and boundary conditions u(a, t) = 0,u(b, t) = 0 where a and b are constants, We can simply solve it by a semidiscretion method. We assume that /(x), w(x), and their derivatives are continuous functions of x. For a difference solution of the above initial boundary value problem, we divide the interval a 6 x 6 b into (N + 1) subintervals each of width h > 0, so that (N + 1)h = b a. Let xi = a + ih, i = 1, 2, . . . , N. 2. Two difference equations We assume that u(x, t) is the exact solution to (1). Let ui and vi be approximations to u(xi, t) and ut(xi, t) 2 i þui1 respectively, and denote f(xi) by fi. Let ooxu2 be replaced by uiþ1 2u . Then Eq. (1) can be semi-discretized at h2 the point (xi,t) by dvi uiþ1 2ui þ ui1 þ 2avi þ b2 ui ¼ þ fi : dt h2
ð2Þ
Let U = (u1,u2, . . . ,u N)t, V = (v1,v2, . . . ,vN)t. Eq. (2) becomes ( dU ¼V; dt dV dt
ð3Þ
¼ AU 2aV þ b;
where 8 2 þ h2 b2 > > > > > 1 > > > > 1< A¼ 2 h > > > > > > > > > :
9 > > > > > > > > > =
1 2 þ h2 b2
1
1
2 þ h2 b2 .. .
1 .. . 1
..
.
2 þ h2 b2 1
> > > > > > 1 > > > ; 2 2 2þh b
and b = (f1,f2, . . . ,f N)t. Denote (u1,u2, . . . ,un,v1,v2, . . . ,vN)t with X, we can obtain dX
¼ BX þ d; X ð0Þ ¼ X 0 ; dt
ð4Þ
where B¼
0 A
I n ; 2aI n
d¼
0 ; b
and t
X 0 ¼ ð/ðx1 Þ; /ðx2 Þ; . . . ; /ðxN Þ; wðx1 Þ; wðx2 Þ; . . . ; wðxN ÞÞ : The exact solution of (4) is X ¼ etB ðB1 d þ X 0 Þ þ B1 d:
ð5Þ
1274
F. Gao, C. Chi / Applied Mathematics and Computation 187 (2007) 1272–1276
Denote B1d + X0 by Y0 and X B1d by Y, solution (5) becomes Y ¼ etB Y 0 :
ð6Þ
Discretize solution (6), let the step size be s, ti = is, i = 0,1, 2, . . . , then Y ðtkþ1 Þ ¼ eðkþ1ÞsB Y 0 ¼ esB Y ðtk Þ: Thus we get a difference scheme for solution (6), Y kþ1 ¼ esB Y k : Then the problem is how to approximate the matrix series esB to get the numerical solution of (6). To do so, we study the function ez. It has the property ea+b = eaeb. Therefore, it satisfies the following equation ez ez ¼ 1:
ð7Þ
A good approximation to ez should preserve the property (7), for example, a [n, n] diagonal Pade approximation to ez always satisfies the property (7). The [1, 1] and [2, 2] Pade approximations to ez are ez
2þz ; 2z
ð8Þ
ez
12 þ 6z þ z2 ; 12 6z þ z2
ð9Þ
and
respectively. One can easily see that the [1, 1] and [2, 2] Pade approximations to ez are ez
2z ; 2þz
ð10Þ
ez
12 6z þ z2 ; 12 þ 6z þ z2
ð11Þ
and
respectively. Approximations (10) and (11) satisfies 2z ez ¼ Oðz3 Þ: 2þz 12 6z þ z2 ez ¼ Oðz5 Þ: 12 þ 6z þ z2
ð12Þ ð13Þ
and it is easy to verify 2 2 z < 1; 12 6z þ z < 1; 2 þ z 12 þ 6z þ z2
ð14Þ
if the real part of z is positive. We use (2I sB)(2I + sB)1 and (12I 6sB + s2 B2)(12 + 6sB + s2B2)1 to approximate the matrix series sB e successively, Then we get the following difference schemes for the numerical solution of (6), 1
Y kþ1 ¼ ð2I sBÞð2I þ sBÞ Y k : 2
2
ð15Þ 2
2 1
Y kþ1 ¼ ð12I 6sB þ s B Þð12 þ 6sB þ s B Þ Y k ;
ð16Þ 3
where Yk is the numerical solution of Y(tk). The accuracy order of difference scheme (15) is O(s ) in time direction because of (12), and O(h2) in space direction, so that is of O(s3 + h2) accurate. In the same way, the accuracy order of difference scheme (16) is O(s5 + h2) because of (13).
F. Gao, C. Chi / Applied Mathematics and Computation 187 (2007) 1272–1276
1275
3. Stability analysis 0 I n Let B ¼ , and k be an eigenvalue of B. We first prove that the real part of k is positive. In fact, A 2aI n since k is an eigenvalue of B, we have n1 n1 0 I n ¼k ; ð17Þ A 2aI n n2 n2 where n1, n2 are some vectors in Rn. From (17), we have n2 ¼ kn1 ;
An1 þ 2an2 ¼ kn2 ;
therefore we have An1 = (2a k)kn1. Thus (2a k)k is an eigenvalue of A. But one can easily see that all eigenvalues of A are positive real numbers, so we have (2a k)k is a positive real number. Let k = a + bi, where a, b are real numbers,then we have ð2a a biÞða þ biÞ > 0:
ð18Þ
From (18), we obtain a > 0. Combine with the condition a > 0, we can finally get the conclusion that the real part of k is positive. We can now show that difference schemes (15) and (16) are unconditionally stable. Let ki, i = 1, 2, . . . , n be eigenvalues of B, one can see that eigenvalues of (2I sB)(2I + sB)1 are (2 ski)(2 + sk i)1, i = 1, 2, . . . , n. We have j(2 sk i)(2 + ski)1j < 1, i = 1,2, . . . ,n, if we note (14) and the result that the real part of ki is positive. Thus, the difference scheme (15) is unconditionally stable. In the same way, we can prove difference scheme (16) is also unconditionally stable. 4. Numerical results and conclusions We apply scheme (16) to the test function uðt; xÞ ¼ et sin x: pffiffiffi e.g., let a ¼ 2; b ¼ 2 in (1). We divide the interval [0,p] into 30 subintervals, e.g., let h = p/30, xk = kh, k = 1,2, . . . ,30, and set the time step size s = 0.1. The errors of numerical results are listed in the following table: x1 t = 0.5 t = 1.0 t = 1.5 t = 2.0
x8 4
0.0483 · 10 0.0904 · 104 0.0990 · 104 0.0884 · 104
x15 4
0.3462 · 10 0.6479 · 104 0.7095 · 104 0.6337 · 104
x22 4
0.4771 · 10 0.8928 · 104 0.9776 · 104 0.8731 · 104
x29 4
0.3777 · 10 0.7069 · 104 0.7740 · 104 0.6913 · 104
0.1430 · 104 0.0904 · 104 0.0990 · 104 0.0884 · 104
The results are satisfactory. To the author’s knowledge, it seems difficult to construct unconditionally stable difference schemes linear hyperbolic equation if by only use of finite difference method. But in some special cases, for example, under the conditions of (1), the technique of semi-discretion and approximating matrix series showed in this paper seems effective. References [1] E.H. Twizell, An explicit difference method for the wave equation with extended stability range, BIT 19 (1979) 378–383. [2] R.K. Mohanty, M.K. Jain, K. George, On the use of high order difference methods for the system of one space second order non-linear hyperbolic equations with variable coefficients, I. Comp. Appl. Math. 72 (1996) 421–431. [3] M. Ciment, S.H. Leventhal, A note on the operator compact implicit method for the wave equation, Math. Comput. 32 (1978) 143–147.
1276
F. Gao, C. Chi / Applied Mathematics and Computation 187 (2007) 1272–1276
[4] R.K. Mohanty, M.K. Jam, An unconditionally stable alternating direction implicit scheme for the two space dimensional linear hyperbolic equation, Namer. Meth. Partial D. Eq. 17 (2001) 684–688. [5] R.K. Mohanty, M.K. Jain, U. Arora, An unconditionally stable AD1 method for the linear hyperbolic equation in three space dimensions, Int. J. Comput. Math. 79 (2002) 133–142. [6] R.K. Mohanty, An unconditionally stable difference scheme for the one-space-dimensional linear hyperbolic equation, Appl. Math. Lett. 17 (2004) 101–105.
Unconditionally stable difference schemes for a one-space-dimensional linear hyperbolic equation Feng Gao
a,*
, Chunmei Chi
b
a b
Science School, Qingdao Technological University, Qingdao 266033, PR China Computer School, Qingdao Technological University, Qingdao 266033, PR China
Abstract A few explicit difference schemes are discussed for the numerical solution of the linear hyperbolic equation utt + 2a ut + b2u = uxx + f(x, t), a > 0, b > 0, in the region X = {(x, t)ja < x < b, t > 0} subject to appropriate initial and Dirichlet boundary conditions, where a and b are real numbers. The proposed scheme is showed to be unconditionally stable, and numerical result is presented. 2006 Elsevier Inc. All rights reserved. Keywords: Second-order linear hyperbolic equation; Damped wave equation; Telegraph equation; Explicit scheme; Unconditionally stable; Pade approximation
1. Introduction In recent years, much attention has been given in the literature to the development, analysis, and implementation of stable methods for the numerical solution of second-order hyperbolic equations, see, for example, [1– 3]. These methods are conditionally stable. Recently, Mohanty et al. [4,5] have developed new three-level implicit unconditionally stable alternating direction implicit schemes for the two and three-space-dimensional linear hyperbolic equations. These schemes are second-order accurate both in space and time. In [6], Mohanty carried over a new technique to the following linear one-space-dimensional hyperbolic equation with given initial and Dirichlet boundary conditions o2 u ou o2 u þ 2a þ b2 u ¼ 2 þ f ðx; tÞ; 2 ot ot ox
a > 0; b P 0;
over a region X = [a < x < b] · [t > 0],with the initial conditions uðx; 0Þ ¼ /ðxÞ;
*
ouðx; 0Þ ¼ wðxÞ; ot
Corresponding author. E-mail address: [email protected] (F. Gao).
0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.09.057
F. Gao, C. Chi / Applied Mathematics and Computation 187 (2007) 1272–1276
1273
and an unconditionally stable implicit three level scheme was obtained. The above equation represents a damped wave equation and a telegraph equation respectively, see [1]. In this paper, however, we will show that if the above problem can be reduced to the following one, o2 u ou o2 u 2 þ b þ 2a u ¼ þ f ðxÞ; ot2 ot ox2
ð1Þ
a > 0; b P 0;
over a region X = [a < x < b] · [t > 0], with the initial conditions uðx; 0Þ ¼ /ðxÞ;
ouðx; 0Þ ¼ wðxÞ ot
and boundary conditions u(a, t) = 0,u(b, t) = 0 where a and b are constants, We can simply solve it by a semidiscretion method. We assume that /(x), w(x), and their derivatives are continuous functions of x. For a difference solution of the above initial boundary value problem, we divide the interval a 6 x 6 b into (N + 1) subintervals each of width h > 0, so that (N + 1)h = b a. Let xi = a + ih, i = 1, 2, . . . , N. 2. Two difference equations We assume that u(x, t) is the exact solution to (1). Let ui and vi be approximations to u(xi, t) and ut(xi, t) 2 i þui1 respectively, and denote f(xi) by fi. Let ooxu2 be replaced by uiþ1 2u . Then Eq. (1) can be semi-discretized at h2 the point (xi,t) by dvi uiþ1 2ui þ ui1 þ 2avi þ b2 ui ¼ þ fi : dt h2
ð2Þ
Let U = (u1,u2, . . . ,u N)t, V = (v1,v2, . . . ,vN)t. Eq. (2) becomes ( dU ¼V; dt dV dt
ð3Þ
¼ AU 2aV þ b;
where 8 2 þ h2 b2 > > > > > 1 > > > > 1< A¼ 2 h > > > > > > > > > :
9 > > > > > > > > > =
1 2 þ h2 b2
1
1
2 þ h2 b2 .. .
1 .. . 1
..
.
2 þ h2 b2 1
> > > > > > 1 > > > ; 2 2 2þh b
and b = (f1,f2, . . . ,f N)t. Denote (u1,u2, . . . ,un,v1,v2, . . . ,vN)t with X, we can obtain dX
¼ BX þ d; X ð0Þ ¼ X 0 ; dt
ð4Þ
where B¼
0 A
I n ; 2aI n
d¼
0 ; b
and t
X 0 ¼ ð/ðx1 Þ; /ðx2 Þ; . . . ; /ðxN Þ; wðx1 Þ; wðx2 Þ; . . . ; wðxN ÞÞ : The exact solution of (4) is X ¼ etB ðB1 d þ X 0 Þ þ B1 d:
ð5Þ
1274
F. Gao, C. Chi / Applied Mathematics and Computation 187 (2007) 1272–1276
Denote B1d + X0 by Y0 and X B1d by Y, solution (5) becomes Y ¼ etB Y 0 :
ð6Þ
Discretize solution (6), let the step size be s, ti = is, i = 0,1, 2, . . . , then Y ðtkþ1 Þ ¼ eðkþ1ÞsB Y 0 ¼ esB Y ðtk Þ: Thus we get a difference scheme for solution (6), Y kþ1 ¼ esB Y k : Then the problem is how to approximate the matrix series esB to get the numerical solution of (6). To do so, we study the function ez. It has the property ea+b = eaeb. Therefore, it satisfies the following equation ez ez ¼ 1:
ð7Þ
A good approximation to ez should preserve the property (7), for example, a [n, n] diagonal Pade approximation to ez always satisfies the property (7). The [1, 1] and [2, 2] Pade approximations to ez are ez
2þz ; 2z
ð8Þ
ez
12 þ 6z þ z2 ; 12 6z þ z2
ð9Þ
and
respectively. One can easily see that the [1, 1] and [2, 2] Pade approximations to ez are ez
2z ; 2þz
ð10Þ
ez
12 6z þ z2 ; 12 þ 6z þ z2
ð11Þ
and
respectively. Approximations (10) and (11) satisfies 2z ez ¼ Oðz3 Þ: 2þz 12 6z þ z2 ez ¼ Oðz5 Þ: 12 þ 6z þ z2
ð12Þ ð13Þ
and it is easy to verify 2 2 z < 1; 12 6z þ z < 1; 2 þ z 12 þ 6z þ z2
ð14Þ
if the real part of z is positive. We use (2I sB)(2I + sB)1 and (12I 6sB + s2 B2)(12 + 6sB + s2B2)1 to approximate the matrix series sB e successively, Then we get the following difference schemes for the numerical solution of (6), 1
Y kþ1 ¼ ð2I sBÞð2I þ sBÞ Y k : 2
2
ð15Þ 2
2 1
Y kþ1 ¼ ð12I 6sB þ s B Þð12 þ 6sB þ s B Þ Y k ;
ð16Þ 3
where Yk is the numerical solution of Y(tk). The accuracy order of difference scheme (15) is O(s ) in time direction because of (12), and O(h2) in space direction, so that is of O(s3 + h2) accurate. In the same way, the accuracy order of difference scheme (16) is O(s5 + h2) because of (13).
F. Gao, C. Chi / Applied Mathematics and Computation 187 (2007) 1272–1276
1275
3. Stability analysis 0 I n Let B ¼ , and k be an eigenvalue of B. We first prove that the real part of k is positive. In fact, A 2aI n since k is an eigenvalue of B, we have n1 n1 0 I n ¼k ; ð17Þ A 2aI n n2 n2 where n1, n2 are some vectors in Rn. From (17), we have n2 ¼ kn1 ;
An1 þ 2an2 ¼ kn2 ;
therefore we have An1 = (2a k)kn1. Thus (2a k)k is an eigenvalue of A. But one can easily see that all eigenvalues of A are positive real numbers, so we have (2a k)k is a positive real number. Let k = a + bi, where a, b are real numbers,then we have ð2a a biÞða þ biÞ > 0:
ð18Þ
From (18), we obtain a > 0. Combine with the condition a > 0, we can finally get the conclusion that the real part of k is positive. We can now show that difference schemes (15) and (16) are unconditionally stable. Let ki, i = 1, 2, . . . , n be eigenvalues of B, one can see that eigenvalues of (2I sB)(2I + sB)1 are (2 ski)(2 + sk i)1, i = 1, 2, . . . , n. We have j(2 sk i)(2 + ski)1j < 1, i = 1,2, . . . ,n, if we note (14) and the result that the real part of ki is positive. Thus, the difference scheme (15) is unconditionally stable. In the same way, we can prove difference scheme (16) is also unconditionally stable. 4. Numerical results and conclusions We apply scheme (16) to the test function uðt; xÞ ¼ et sin x: pffiffiffi e.g., let a ¼ 2; b ¼ 2 in (1). We divide the interval [0,p] into 30 subintervals, e.g., let h = p/30, xk = kh, k = 1,2, . . . ,30, and set the time step size s = 0.1. The errors of numerical results are listed in the following table: x1 t = 0.5 t = 1.0 t = 1.5 t = 2.0
x8 4
0.0483 · 10 0.0904 · 104 0.0990 · 104 0.0884 · 104
x15 4
0.3462 · 10 0.6479 · 104 0.7095 · 104 0.6337 · 104
x22 4
0.4771 · 10 0.8928 · 104 0.9776 · 104 0.8731 · 104
x29 4
0.3777 · 10 0.7069 · 104 0.7740 · 104 0.6913 · 104
0.1430 · 104 0.0904 · 104 0.0990 · 104 0.0884 · 104
The results are satisfactory. To the author’s knowledge, it seems difficult to construct unconditionally stable difference schemes linear hyperbolic equation if by only use of finite difference method. But in some special cases, for example, under the conditions of (1), the technique of semi-discretion and approximating matrix series showed in this paper seems effective. References [1] E.H. Twizell, An explicit difference method for the wave equation with extended stability range, BIT 19 (1979) 378–383. [2] R.K. Mohanty, M.K. Jain, K. George, On the use of high order difference methods for the system of one space second order non-linear hyperbolic equations with variable coefficients, I. Comp. Appl. Math. 72 (1996) 421–431. [3] M. Ciment, S.H. Leventhal, A note on the operator compact implicit method for the wave equation, Math. Comput. 32 (1978) 143–147.
1276
F. Gao, C. Chi / Applied Mathematics and Computation 187 (2007) 1272–1276
[4] R.K. Mohanty, M.K. Jam, An unconditionally stable alternating direction implicit scheme for the two space dimensional linear hyperbolic equation, Namer. Meth. Partial D. Eq. 17 (2001) 684–688. [5] R.K. Mohanty, M.K. Jain, U. Arora, An unconditionally stable AD1 method for the linear hyperbolic equation in three space dimensions, Int. J. Comput. Math. 79 (2002) 133–142. [6] R.K. Mohanty, An unconditionally stable difference scheme for the one-space-dimensional linear hyperbolic equation, Appl. Math. Lett. 17 (2004) 101–105.