Pseudospectral vs Finite Di erence Methods for Initial Value Problems ...
Pseudospectral vs Finite Di erence Methods for Initial Value Problems with Discontinuous Coecients Erding Luo
Heinz-Otto Kreissy
June 9, 1995
Abstract
An initial value problem with piecewise constant coecients is considered. The accuracies for both nite di erence methods and the pseudospectral method are analyzed, and a modi cation of the initial value problem is suggested. The modi ed problem is shown to have the same temporal period as the original problem does, and a second order accuracy is obtained for the pseudospectral method at integral multiples of the temporal period.
1 Introduction The pseudospectral method is very powerful for periodic initial value problems with constant coecients. Well known results include its convergence at a spectral rate, [2, 3, 4, 8, 9]. This is largely because the eigenfunctions for the pseudospectral method are exactly those of fourier modes. Hence, the error for the pseudospectral method is not bigger than the truncated error of the initial value. For equations with jump coecients, however, the pseudospectral method is far from being clear, even for one dimensional case, [1, 4, 5, 6, 7]. For such a problem, a large number of fourier eigenmodes to well represent all the eigenfunctions of the problem is needed in order to maintain a fast convergence rate, [3, 4]. The purpose of this paper is to compare the pseudospectral method with nite di erence methods. For this, we will focus on an one dimensional periodic boundary wave equation, with a speci c initial condition. The coecient of the problem is assumed to be a piecewise constant function. We study the accuracies both for nite di erence methods and the pseudospectral method using eigenfunction analysis. Our approach is to determine how many fourier eigenmodes are needed to well represent the initial condition with a given accuracy, and to analyze the accuracies of these eigenmodes. For nite di erence methods, local truncation error analysis is used to obtain a global error estimate. For the pseudospectral method, the error estimate is generated by numerical computations. Since the coecient is discontinuous, with an arbitrary discretization in the spatial domain, the accuracies of numerical methods may be a ected by the values around the discontinuous positions. For y
IMA, University of Minnesota, Mpls, MN 55455. Research supported by NSF through IMA. Department of Mathematics, University of California at Los Angeles, LA, CA90024.
1
the case in which the discontinuous points align with grid points, the harmonic values of the adjacent points are suggested for the coecient at these discontinuous points. p Under this construction, nite di erence schemes give us an accuracy with order O(h h) in l2 norm and O(h) in l1 norm. In general, we do not know the exact locations of the discontinuous points. Both nite di erence methods and the pseudospectral method give us an approximation with a shift. First order accuracy is always what we can expect, [1]. To improve this accuracy, we suggest a modi ed problem by truncating the fourier expansion of the inverse of the coecient. This modi ed problem retains the same temporal period as the original one.pIt turns out that the pseudospectral method still has O(h) accuracy in l1 norm, but O(h h) in l2 norm. Moreover, at integral multiples of the temporal period, the pseudospectral method always has O(h2) accuracy in l1 norm. The rest of the paper is organized as follows: our model problem is presented in section 2. In section 3, Numerical schemes are introduced. The accuracies of both numerical methods are demonstrated in section 4 and some numerical results are presented in section 5.
2 Model Problem Consider the following equation @u(x; t) = @t
(x; t) ; x 2 ( 1; 1) a(x) @u@x
(2.1)
with periodic boundary condition, where the coecient a(x) is also periodic with two discontinuous points at x = 1=2; 1=2, such that ( a(x) = 01:5; ; xx 22 ([ 10;:51]; =0[:5)0;:5; 0:5]: (2.2) The initial condition takes the following form,
f (x) = exp( Cx ); x 2 [ 1; 1]: 2
2.1 Analytic Solution
(2.3)
Denote the kth (k 2 Z ) eigenfunction of (2.1) by k , which satis es the periodic boundary condition. Then, k (x) a(x) @ @x = k (x): Thus, Zx 1 k (x) = exp( a() d); where
1
Z
= 2ik; = (
2
1 d ) 1 = 1 : 3 1 a( )
1
Let Then, k can be rewritten as
k = 2ik; k 2 Z:
Zx 1 (2.4) a() d); x 2 [ 1; 1]: The sequence f k (x)g forms an orthogonal base with the weighing function a x in the domain [ 1; 1]. The solution of (2.1) can thus be expanded in terms of the eigenfunctions as Zx 1 X u(x; t) ck exp(k(t (2.5) a() d)); k where ck (k 2 Z ) is a constant that is determined by the initial condition (2.3). Remark. Using the characteristic functions, it can be shown that the analytic solution of equation (2.1) will not be a ected by the values of the coecient a(x) at discontinuous points.
k (x) = exp( k
1
1 ( )
1
2.2 Expansion of Initial Condition
Denote the fourier eigenmodes by 'k (x) = exp( kx); k 2 Z: (2.6) The initial condition (2.3) can also be expanded in terms of a certain number of the fourier eigenmodes (2.6) within endurance accuracy. Theorem 2.1 Let function f be as in (2.3) with the parameter C = 60. Then, a few number of fourier eigenmodes are sucient to well represent the initial condition with 1% accuracy. More precisely, let M = 26, k<M= X2 1 ; 8x 2 [ 1; 1]; ck'k( 23x + 12 )j 100 jf (x) 31 k M=2 where
ck =
Z
1 f (x) (x)dx; k 1 a(x)
1
and k (x) is as in (2.4). Proof. De ne an auxiliary 1-periodic function g(y) = 3f ( 32 y 34 ); y 2 [0; 1]: Then, its fourier series is as follows,
g(y) =
X
k<M=2 k M=2
Gk exp( ky) + M ; 8y; 3
(2.7)
where
Gk =
Z
1 0
g(y) exp(ky)dy = 3
and M (y ) =
Z
1
0
kM=2
X
k< M=2
f ( 23 y 34 ) exp(ky)dy;
Gk exp( ky)
denotes the residual of the truncated fourier series for g (y ). Since 12 < jxj 1, it follows f (x) for some negligibly small number . From (2.7), it follows Z1 1 Zx 1 ck = f (x) exp(k d)dx 1 a(x) 1 a( ) Z 1=2 = 2 f (x) exp(k(2x + 32 ))dx + O() 1=2 Z 5=6 3 3 = 3 f ( 2 y 4 ) exp(ky)dy + O() 1=6 Z11 3 3 f ( 2 y 4 ) exp(ky)dy + O() = 0 = Gk + O( ): Therefore,
g(y) = which implies,
X
k<M=2 k M=2
ck exp( ky) + O(M ) + M ; 8y;
k<M= X2 1 f (x) = 3 ck exp( k(2x + 23 )) + O(M ) + M ; 8x 2 [ 3=4; 3=4]: k M=2
Since g (y ) is 1-periodic, for all y 2 [0; 1]=(1=6; 5=6);
g(y) = O(): Equivalently, for x 2 [ 1; 1]=( 3=4; 3=4); X2 1 k<M= 3 c k exp( k (2x + )) + O(M ) + M = O( ): 3 k M=2 2 Thus, for all x 2 [ 1; 1], k<M= X2 1 f (x) 3 ck'k( 23x + 12 ) = O(M ) + M : k M=2
4
Because is negligible, setting s = 2=3k it follows from integration by parts that Z 3=4 jG j j 4C exp( =2) exp(sx Cx2 )(1 2Cx2)dxj k
k s = Z = 8 C j s exp( k=2) exp(sx Cx )x(3 2Cx )dxj = Z = 8 C exp(sx Cx )(3 12Cx + 4C x )dxj j s exp( k=2) = Z = 16 C j s exp( k=2) exp(sx Cx )(15 90Cx + 60C x 8C x )dxj = Z = 64 C j s exp( k=2) exp(sx Cx )(945 9450Cx + 12600C x = 5040C x + 720C x 32C x )dxj: 2
3 4
3 4
2
3
2
2
3 4
3 4
2
4
2
2
2
4
3 4
3
3 4
6
2
2
2
4
3
6
3 4
5
3 4
10
2
2
2
4
3 4
3
6
4
8
5
We can thus establish
10
C: jGkj 90720 s 5
(2.8)
10
Consequently,
jM j
kM=2
X
k< M=2
jGkj
kM=2
X 90720C s k< M=
C 3 90720 2 M : 5
10
In order that
5
10
2
10
10
9
1 ; jM j 100
it suces to set M = 26. Therefore, k<M= X2 1 ; 8x 2 [ 1; 1]: jf (x) 13 ck'k ( 23x + 21 )j 100 k M=2 In summary, for the given initial condition (2.3) with the parameter C = 60, the number M of fourier eigenmodes needed to well represent the initial condition (2.3) is 26. The accuracy of the approximation of (2.1) is 1%. Since 'k ( x + ) is orthogonal in [ ; ], a sucient number of fourier eigenmodes in the entire domain [ 1; 1] is 32 in order to 2 3
1 2
3 4
3 4
maintain the same accuracy.
2.3 Modi ed Problem
We modify equation (2.1) as follows @v(x; t) =
@t
1 @v (x; t) b(x) @x ; x 2 ( 1; 1); 5
(2.9)
where the coecient b(x) is the truncated fourier series of a(1x) . Note rst the fourier expansion of a(1x) is 1 1 X ^ a(x) != 1 b! exp(i!(x + 1)); where ^b! is the fourier coecient of a(1x) with an order of O( !1 ) given by ( 1 ! ; ! 6= 0; 2 ^b! = R 1! sin 1 dx; ! = 0: 1 a(x) If we let
RN (x) = then b(x) can be written as
!N=2
X ^ b! exp(i!(x + 1));
!< N=2
!
Heinz-Otto Kreissy
June 9, 1995
Abstract
An initial value problem with piecewise constant coecients is considered. The accuracies for both nite di erence methods and the pseudospectral method are analyzed, and a modi cation of the initial value problem is suggested. The modi ed problem is shown to have the same temporal period as the original problem does, and a second order accuracy is obtained for the pseudospectral method at integral multiples of the temporal period.
1 Introduction The pseudospectral method is very powerful for periodic initial value problems with constant coecients. Well known results include its convergence at a spectral rate, [2, 3, 4, 8, 9]. This is largely because the eigenfunctions for the pseudospectral method are exactly those of fourier modes. Hence, the error for the pseudospectral method is not bigger than the truncated error of the initial value. For equations with jump coecients, however, the pseudospectral method is far from being clear, even for one dimensional case, [1, 4, 5, 6, 7]. For such a problem, a large number of fourier eigenmodes to well represent all the eigenfunctions of the problem is needed in order to maintain a fast convergence rate, [3, 4]. The purpose of this paper is to compare the pseudospectral method with nite di erence methods. For this, we will focus on an one dimensional periodic boundary wave equation, with a speci c initial condition. The coecient of the problem is assumed to be a piecewise constant function. We study the accuracies both for nite di erence methods and the pseudospectral method using eigenfunction analysis. Our approach is to determine how many fourier eigenmodes are needed to well represent the initial condition with a given accuracy, and to analyze the accuracies of these eigenmodes. For nite di erence methods, local truncation error analysis is used to obtain a global error estimate. For the pseudospectral method, the error estimate is generated by numerical computations. Since the coecient is discontinuous, with an arbitrary discretization in the spatial domain, the accuracies of numerical methods may be a ected by the values around the discontinuous positions. For y
IMA, University of Minnesota, Mpls, MN 55455. Research supported by NSF through IMA. Department of Mathematics, University of California at Los Angeles, LA, CA90024.
1
the case in which the discontinuous points align with grid points, the harmonic values of the adjacent points are suggested for the coecient at these discontinuous points. p Under this construction, nite di erence schemes give us an accuracy with order O(h h) in l2 norm and O(h) in l1 norm. In general, we do not know the exact locations of the discontinuous points. Both nite di erence methods and the pseudospectral method give us an approximation with a shift. First order accuracy is always what we can expect, [1]. To improve this accuracy, we suggest a modi ed problem by truncating the fourier expansion of the inverse of the coecient. This modi ed problem retains the same temporal period as the original one.pIt turns out that the pseudospectral method still has O(h) accuracy in l1 norm, but O(h h) in l2 norm. Moreover, at integral multiples of the temporal period, the pseudospectral method always has O(h2) accuracy in l1 norm. The rest of the paper is organized as follows: our model problem is presented in section 2. In section 3, Numerical schemes are introduced. The accuracies of both numerical methods are demonstrated in section 4 and some numerical results are presented in section 5.
2 Model Problem Consider the following equation @u(x; t) = @t
(x; t) ; x 2 ( 1; 1) a(x) @u@x
(2.1)
with periodic boundary condition, where the coecient a(x) is also periodic with two discontinuous points at x = 1=2; 1=2, such that ( a(x) = 01:5; ; xx 22 ([ 10;:51]; =0[:5)0;:5; 0:5]: (2.2) The initial condition takes the following form,
f (x) = exp( Cx ); x 2 [ 1; 1]: 2
2.1 Analytic Solution
(2.3)
Denote the kth (k 2 Z ) eigenfunction of (2.1) by k , which satis es the periodic boundary condition. Then, k (x) a(x) @ @x = k (x): Thus, Zx 1 k (x) = exp( a() d); where
1
Z
= 2ik; = (
2
1 d ) 1 = 1 : 3 1 a( )
1
Let Then, k can be rewritten as
k = 2ik; k 2 Z:
Zx 1 (2.4) a() d); x 2 [ 1; 1]: The sequence f k (x)g forms an orthogonal base with the weighing function a x in the domain [ 1; 1]. The solution of (2.1) can thus be expanded in terms of the eigenfunctions as Zx 1 X u(x; t) ck exp(k(t (2.5) a() d)); k where ck (k 2 Z ) is a constant that is determined by the initial condition (2.3). Remark. Using the characteristic functions, it can be shown that the analytic solution of equation (2.1) will not be a ected by the values of the coecient a(x) at discontinuous points.
k (x) = exp( k
1
1 ( )
1
2.2 Expansion of Initial Condition
Denote the fourier eigenmodes by 'k (x) = exp( kx); k 2 Z: (2.6) The initial condition (2.3) can also be expanded in terms of a certain number of the fourier eigenmodes (2.6) within endurance accuracy. Theorem 2.1 Let function f be as in (2.3) with the parameter C = 60. Then, a few number of fourier eigenmodes are sucient to well represent the initial condition with 1% accuracy. More precisely, let M = 26, k<M= X2 1 ; 8x 2 [ 1; 1]; ck'k( 23x + 12 )j 100 jf (x) 31 k M=2 where
ck =
Z
1 f (x) (x)dx; k 1 a(x)
1
and k (x) is as in (2.4). Proof. De ne an auxiliary 1-periodic function g(y) = 3f ( 32 y 34 ); y 2 [0; 1]: Then, its fourier series is as follows,
g(y) =
X
k<M=2 k M=2
Gk exp( ky) + M ; 8y; 3
(2.7)
where
Gk =
Z
1 0
g(y) exp(ky)dy = 3
and M (y ) =
Z
1
0
kM=2
X
k< M=2
f ( 23 y 34 ) exp(ky)dy;
Gk exp( ky)
denotes the residual of the truncated fourier series for g (y ). Since 12 < jxj 1, it follows f (x) for some negligibly small number . From (2.7), it follows Z1 1 Zx 1 ck = f (x) exp(k d)dx 1 a(x) 1 a( ) Z 1=2 = 2 f (x) exp(k(2x + 32 ))dx + O() 1=2 Z 5=6 3 3 = 3 f ( 2 y 4 ) exp(ky)dy + O() 1=6 Z11 3 3 f ( 2 y 4 ) exp(ky)dy + O() = 0 = Gk + O( ): Therefore,
g(y) = which implies,
X
k<M=2 k M=2
ck exp( ky) + O(M ) + M ; 8y;
k<M= X2 1 f (x) = 3 ck exp( k(2x + 23 )) + O(M ) + M ; 8x 2 [ 3=4; 3=4]: k M=2
Since g (y ) is 1-periodic, for all y 2 [0; 1]=(1=6; 5=6);
g(y) = O(): Equivalently, for x 2 [ 1; 1]=( 3=4; 3=4); X2 1 k<M= 3 c k exp( k (2x + )) + O(M ) + M = O( ): 3 k M=2 2 Thus, for all x 2 [ 1; 1], k<M= X2 1 f (x) 3 ck'k( 23x + 12 ) = O(M ) + M : k M=2
4
Because is negligible, setting s = 2=3k it follows from integration by parts that Z 3=4 jG j j 4C exp( =2) exp(sx Cx2 )(1 2Cx2)dxj k
k s = Z = 8 C j s exp( k=2) exp(sx Cx )x(3 2Cx )dxj = Z = 8 C exp(sx Cx )(3 12Cx + 4C x )dxj j s exp( k=2) = Z = 16 C j s exp( k=2) exp(sx Cx )(15 90Cx + 60C x 8C x )dxj = Z = 64 C j s exp( k=2) exp(sx Cx )(945 9450Cx + 12600C x = 5040C x + 720C x 32C x )dxj: 2
3 4
3 4
2
3
2
2
3 4
3 4
2
4
2
2
2
4
3 4
3
3 4
6
2
2
2
4
3
6
3 4
5
3 4
10
2
2
2
4
3 4
3
6
4
8
5
We can thus establish
10
C: jGkj 90720 s 5
(2.8)
10
Consequently,
jM j
kM=2
X
k< M=2
jGkj
kM=2
X 90720C s k< M=
C 3 90720 2 M : 5
10
In order that
5
10
2
10
10
9
1 ; jM j 100
it suces to set M = 26. Therefore, k<M= X2 1 ; 8x 2 [ 1; 1]: jf (x) 13 ck'k ( 23x + 21 )j 100 k M=2 In summary, for the given initial condition (2.3) with the parameter C = 60, the number M of fourier eigenmodes needed to well represent the initial condition (2.3) is 26. The accuracy of the approximation of (2.1) is 1%. Since 'k ( x + ) is orthogonal in [ ; ], a sucient number of fourier eigenmodes in the entire domain [ 1; 1] is 32 in order to 2 3
1 2
3 4
3 4
maintain the same accuracy.
2.3 Modi ed Problem
We modify equation (2.1) as follows @v(x; t) =
@t
1 @v (x; t) b(x) @x ; x 2 ( 1; 1); 5
(2.9)
where the coecient b(x) is the truncated fourier series of a(1x) . Note rst the fourier expansion of a(1x) is 1 1 X ^ a(x) != 1 b! exp(i!(x + 1)); where ^b! is the fourier coecient of a(1x) with an order of O( !1 ) given by ( 1 ! ; ! 6= 0; 2 ^b! = R 1! sin 1 dx; ! = 0: 1 a(x) If we let
RN (x) = then b(x) can be written as
!N=2
X ^ b! exp(i!(x + 1));
!< N=2
!
