THE CONVERGENCE OF SPECTRAL AND ... - Semantic Scholar
SIAM J. SCI. COMPUT. Vol. 23, No. 5, pp. 1731–1751
c 2002 Society for Industrial and Applied Mathematics
THE CONVERGENCE OF SPECTRAL AND FINITE DIFFERENCE METHODS FOR INITIAL-BOUNDARY VALUE PROBLEMS∗ NATASHA FLYER† ‡ AND PAUL N. SWARZTRAUBER† Abstract. The general theory of compatibility conditions for the differentiability of solutions to initial-boundary value problems is well known. This paper introduces the application of that theory to numerical solutions of partial differential equations and its ramifications on the performance of high-order methods. Explicit application of boundary conditions (BCs) that are independent of the initial condition (IC) results in the compatibility conditions not being satisfied. Since this is the case in most science and engineering applications, it is shown that not only does the error in a spectral method, as measured in the maximum norm, converge algebraically, but the accuracy of finite differences is also reduced. For the heat equation with a parabolic IC and Dirichlet BCs, we prove that the Fourier method converges quadratically in the neighborhood of t = 0 and the boundaries and quartically for large t when the first-order compatibility conditions are violated. For the same problem, the Chebyshev method initially yields quartic convergence and exponential convergence for t > 0. In contrast, the wave equation subject to the same conditions results in inferior convergence rates with all spectral methods yielding quadratic convergence for all t. These results naturally direct attention to finite difference methods that are also algebraically convergent. In the case of the wave equation, we prove that a second-order finite difference method is reduced to 4/3-order convergence and numerically show that a fourth-order finite difference scheme is apparently reduced to 3/2-order. Finally, for the wave equation subject to general ICs and zero BCs, we give a conjecture on the error for a second-order finite difference scheme, showing that an O(N −2 log N ) convergence is possible. Key words. compatibility conditions, convergence theory, spectral and finite difference methods AMS subject classifications. 65M06, 65M12, 65M70 PII. S1064827500374169
1. Introduction. The regularity of solutions for initial-boundary value problems (IBVPs) is determined by the compatibility of the initial condition (IC) and the boundary conditions (BCs). IBVPs will have singular solutions at the corners of the space-time domain if the BCs do not satisfy the partial differential equation (PDE) or any of its higher-order derivatives, even if the IC is C ∞ . The Cauchy–Kowalesky theorem [19] provides a solution in the neighborhood of t = 0 to the same problem without BCs. If compatibility exists, then the BCs can always be computed by the solution obtained from the Cauchy–Kowalesky theorem and the IC. In such cases, the IBVP could in fact be posed as an initial-value problem (IVP). The set of compatibility conditions are derived by equating the time derivatives of the BCs to those of the solution to the IVP given by the Cauchy–Kowalesky theorem. Literature on the theory of compatibility conditions and the regularity of solutions for IBVPs is extensive, starting in the 1950s with the work of Ladyzenskaja [10], [11]. In the 1960s, Ladyzenskaja, Solonnikov, and Ural’ceva [12] and Friedman [3] discussed the regularity of the solution for linear and quasi-linear parabolic systems. During this same time, Kreiss [4] and Hersh [7] extended the theory for constant coefficient hyperbolic systems by providing the necessary and sufficient conditions for ∗ Received by the editors June 20, 2000; accepted for publication (in revised form) August 6, 2001; published electronically January 30, 2002. http://www.siam.org/journals/sisc/23-5/37416.html † National Center for Atmospheric Research, Advanced Study Program, P.O. Box 3000, Boulder, CO 80307-3000 ([email protected]). NCAR is sponsored by the National Science Foundation. ‡ Current address: Department of Applied Mathematics, University of Colorado, Boulder, CO 80309 (fl[email protected]).
1731
1732
NATASHA FLYER AND PAUL N. SWARZTRAUBER
the problem to be well-posed. It was further developed in the 1970s by Rauch and Massey [16] and Sakamoto [17] for hyperbolic systems with time-dependent coefficients. In the 1980s, Temam [20] extended the theory by obtaining the compatibility conditions for semilinear evolution equations while Smale [18] gave a rigorous proof of the compatibility conditions for the heat and wave equation on a bounded domain. These accomplishments have been confined to the realm of theoretical mathematics. Yet, not satisfying the compatibility conditions of a system can have a significant computational impact, particularly for high-order accurate methods, such as spectral methods, whose performance intrinsically depends on the smoothness of the solution. For nonsmooth initial data, there is a wealth of information on the convergence theory of numerical schemes in a variety of norms for parabolic and hyperbolic IVPs [13], [5], [21], [22]. These analyses can be applied to the IBVPs if the IBVP can be posed as an IVP with discontinuous ICs through periodic extensions of the ICs to the entire real axis [6], [9], [14]. However, such periodic extensions may induce singularities that do not exist in a smooth nonperiodic solution of an IBVP. More recently, the impact of noncompatible BCs has been explored by Boyd and Flyer [2] in a computational framework. This paper appeared in a special issue on spectral methods. In the prologue, Karniadakis states “This is the first such effort in bringing the theory of compatibility conditions from the mathematical literature to the numerical community.” This work illustrated the ubiquity of compatibility conditions and analyzed the connection between incompatibility and the rate of convergence of Chebyshev spectral series. However, the convergence of the error was not discussed or compared to alternatives such as the finite difference method. The focus of the paper, as has been in the mathematical literature, is on smoothing the initial condition so as to satisfy the compatibility conditions [2], [5], [18], [15]. The difficulty with this approach is that it leaves the fundamental question unanswered: namely, “What is the rate of convergence to the solution of the original unperturbed problem?” If the IC is smoothed, spectral convergence might be achieved but to a solution that differs algebraically from the original problem. The focus of the current work is on the convergence of the approximate solution to the exact solution. As discussed in the first paragraph, the temporal derivatives of the solution as defined by the IC, differential operator, and Cauchy–Kowalesky theorem will not equal those determined by the independent BCs. The resulting singularities in the corners of the temporal-spatial domain, which are independent of the smoothness of the IC, disrupt the convergence of spectral and finite difference methods in a manner that differs significantly between parabolic and hyperbolic systems. A detailed analysis of the induced singularities for two model problems and their effect on the convergence of both numerical methods is determined. Since an incompatibility exists any time the BC is independent of the IC, algebraic convergence in the maximum norm is the expectation for IBVPs. In most science and engineering applications, this is the usual scenario. In section 2, the infinite set of compatibility conditions for the generalized heat and wave equations with time-dependent BCs is presented. The purpose is to lay the groundwork for describing the inherent singular nature of IBVPs and how it affects the convergence rate of the error. In section 3, we examine the prototype parabolic case, namely the heat equation subject to Dirichlet BCs and a smooth parabolic IC. We use the Fourier method to determine an exact expression for the approximate solution, represented by a truncated Fourier series. Convergence rates are then developed for both finite difference and Chebyshev spectral methods, the latter of which is normally the basis set of choice for a problem on a finite domain. For the parabolic example
INITIAL-BOUNDARY VALUE PROBLEMS
1733
in section 3, ut is discontinuous at t = 0. Nevertheless, the Chebyshev method is shown to yield spectral convergence for t > 0, but only quartic convergence in the neighborhood of the discontinuity, which would be the relevant consideration for the computation of transient heat flow. Furthermore, the usual Fourier method is only quadratically convergent near the discontinuity and quartically convergent elsewhere. Thus, the convergence rate of the error is shown to be nonuniform and algebraic as measured in the maximum norm. The hyperbolic case is discussed in section 4. Discontinuities that are induced by the incompatibilities now propagate throughout the region. It is shown that all spectral methods will yield algebraic convergence for t > 0. The Fourier method is used again to derive an exact expression for the error. Similar to the parabolic case, the results are compared with the centered finite difference and Chebyshev methods which yield algebraic convergence for smooth ICs. The performance of finite difference methods is also significantly diminished when the compatibility conditions are not satisfied. For this reason, the appendix contains the rates of convergence for a second-order finite difference method for a variety of smooth ICs and shows, for example, that second-order convergence is reduced to 4/3-order. 2. Compatibility conditions for the generalized heat and wave equations. Here compatibility conditions are reviewed to illustrate their fundamental impact on the smoothness of the solutions to IBVPs. We simply state the necessary conditions for the solution to be C ∞ on the domain for the heat and wave equation with time-dependent BCs. The complete proofs of necessity and sufficiency are given for parabolic systems in [11], [12], [3], [18] and for hyperbolic systems in [10], [16], [17], [4]. Theorem 2.1. The domain is [0, T ] × Ω, where Ω is a d-dimensional spatial domain with a boundary ∂Ω which is a C ∞ manifold of dimension (d − 1). Let L be an elliptic operator of the form (2.1)
L=
d d
d
Aij (x)
i=1 j=1
d
∂ ∂ ∂ + Bj (x) + Cj (x). ∂xi ∂xj j=1 ∂xj j=1
The generalized linear diffusion problem is ut = Lu,
(2.2)
u(x, 0) = u0 (x) ∈ Ω.
The generalized wave equation is utt = Lu,
(2.3)
u(x, 0) = u0 (x),
ut (x, 0) = v0 (x) ∈ Ω.
Both are subject to the BC (2.4)
u = f (t) ∈ ∂Ω ∀t.
Then, the necessary and sufficient compatibility conditions for u(x,t) to be C 2k for the heat equation are Lk u0 =
(2.5)
∂kf , ∂tk
k = 0, 1, 2, . . . ∈ [0, 0] × ∂Ω
and for the wave equation are (2.6)
Lk u0 =
∂ 2k f ∂t2k
and
Lk v0 =
∂ 2k+1 f , ∂t2k+1
k = 0, 1, 2, . . . ∈ [0, 0] × ∂Ω.
1734
NATASHA FLYER AND PAUL N. SWARZTRAUBER
In essence (2.5) and (2.6) state that at time t = 0, the temporal derivatives of the solution determined by the application of the differential operator to the IC must equal the temporal derivatives of the BCs. When these conditions are satisfied, we do not need to impose BCs because they can be derived from the Taylor series expansion of the IC. When (2.5) and (2.6) do not hold, we induce a singularity in the 2kth-order derivative of the solution at t = 0 on ∂Ω. Thus, imposing BCs that are independent of the IC places an infinite set of constraints on the solution which differ from the constraints enforced by the PDE, resulting in discontinuities in the corners of the temporal-spatial domain. The manner in which these incompatibilities affect the accuracy of the numerical method depends on whether the system is parabolic or hyperbolic. In the next section we will discuss the parabolic case, followed in section 4 by the hyperbolic case. 3. A parabolic example. Our goal is to study the convergence of the error for approximate spectral and finite difference solutions to IBVPs when the compatibility conditions of the system are not satisfied. With this in mind, we first study the effect on parabolic systems, taking as our example the one-dimensional heat equation ut = uxx
(3.1) subject to BCs (3.2)
u(0, t) = u(π, t) = 0
and IC (3.3)
u0 = u(x, 0) =
π x(π − x). 8
The first-order compatibility condition, (k = 1) in (2.5), is not satisfied because 2 π ∂ u0 =− Lu0 = (3.4) ∀x 2 ∂x 4 and ∂u |x=0,π = 0 ∂t
(3.5)
∀t.
Thus from (3.1), uxx |x=0,π = ut |x=0,π = 0 yet uxx = −2 in the interior. To see the impact of the jump discontinuity at the boundary in the second derivative at t = 0 on the decay rate of the error we will solve the problem exactly at time t = 0 via a truncated Fourier series, then compare this together with finite difference and Chebyshev methods against the exact solution for all time. Time integration is assumed exact throughout, which yields the approximate spectral representation (3.6)
N −1
uN (xi , t) =
2
cN (n)e−n t sin nxi
n=1,odd
of the exact solution, given by (3.7)
u(xi , t) =
∞ n=1,odd
2
e−n t sin nxi . n3
1735
INITIAL-BOUNDARY VALUE PROBLEMS
Gottlieb and Orszag [8] note that for the exact coefficients cN (n), uN (x, t) converges spectrally to u(x, t). This is also true if u(x, 0) is known for all x because cN (n) can be computed to any accuracy. In this work, we assume, as with most scientific data, that u(x, 0) is given in tabular form ui = u(xi , t), where xi = iπ/N . In science and engineering, we cannot assume we know the function exactly and therefore must approach the problem from a numerical standpoint. The coefficients cN (n) are then computed in the traditional manner, using trigonometric interpolation, as cN (n) =
(3.8)
N −1 2 ui sin nxi . N i=1
With this implementation of the spectral method, it will be shown that the approximate spectral representation (3.6) converges algebraically to the exact solution in the maximum norm, and, furthermore, convergence is nonuniform. This characteristic algebraic convergence is method-independent due to the inherent singularity induced by violating the compatibility conditions from explicit application of the BCs. In the next section, we derive an exact expression for cN (n) which will provide the initial error in the coefficients for the Fourier method. 3.1. The approximate spectral solution. In what follows, we derive an analytical representation of the coefficients cN (n) in (3.6). From (3.7) ∞ 1 iπ u(xi , 0) = (3.9) , sin n n3 N n=1,odd
(3.10)
N −1
u(xi , 0) =
n=1,odd ∞
1 iπ sin(n ) n3 N
1 iπ sin (2mN + n)i (2mN + n)3 N m=1 1 iπ + , sin (2mN − n)i 3 (2mN − n) N
+
or (3.11)
u(xi , 0) =
N −1 n=1,odd
∞ 1 1 1 + − 3 3 n (2mN + n) (2mN − n)3 m=1
iπ sin n N
and (3.12)
u(xi , 0) =
N −1 n=1,odd
1 n iπ , − 4 Ψ(n/N ) sin n n3 N N
where (3.13)
Ψ(s) = 2
∞ 3m2 + s2 . [m2 − s2 ]3 m=1
Ψ(s) can be called the alias function since it represents the error in the coefficients due to truncating the series after N terms. The initial error in the coefficients is due
1736
NATASHA FLYER AND PAUL N. SWARZTRAUBER
to the fact that a sine series cannot represent a nonperiodic polynomial with a finite number of terms. For use later, we observe that Ψ(s) is positive and monotonically increasing on the interval [0, 1] with extremes Ψ(0) ≈ 0.4058712 and Ψ(1) = 1. Equating (3.6) and (3.12), we obtain the desired result 1 nΨ(n) − . 3 n N4 We study convergence only as it relates to spatial discretization, and therefore temporal integrations are represented analytically. The spectral approximation is then given by N −1 2 1 nΨ(n/N ) uN (xi , t) = sin nxi . (3.15) e−n t 3 − n N4 (3.14)
cN (n) =
n=1,odd
In the next two sections it will be proved that although the spectral series for both the exact and approximate solutions decay exponentially, the error decays algebraically and nonuniformly over the domain. Section 3.2 explores convergence in the neighborhood of the discontinuity induced by the incompatibility of imposing BCs. Section 3.3 looks at convergence in the interior of the domain, bounded away from the singularities at (0, 0) and (0, π). 3.2. Quadratic convergence near (0, 0) and (0, π). The discontinuity in ut at (0, 0) induces an error in ∂uN /∂t(π/N, 0) that is bounded away from zero for all N. In what follows, we show that this produces an O(N −2 ) error in the neighborhood of (0, 0) and likewise near (π, 0). With a derivation similar to that preceding (3.12), the exact solution is given by
2 2 2 N −1 ∞ e−n t e−(mN +n) t e−(mN −n) t sin nxi . (3.16) u(xi , t) = − − n3 [mN + n]3 [mN − n]3 m=2,even n=1,odd
Subtracting uN (xi , t), given by (3.15), defines the error (3.17) eN (xi , t) =
N −1 n=1,odd
∞
m=2,even
2
e−(mN +n) t − e−n [mN + n]3
2
t
2
e−(mN −n) t − e−n − [mN − n]3
2
t
sin nxi .
Then, on the discrete trajectory x1 = π/N and tN = 1/N 2 π 1 1 = 2 IN , (3.18) eN , N N2 N where (3.19) IN
1 = N
N −1 n=1,odd
∞
m=2,even
2
2
2
e−(m+n/N ) − e−(n/N ) e−(m−n/N ) − e−(n/N ) − [m + n/N ]3 [m − n/N ]3
2
sin n
IN defines the midpoint quadrature of the smooth positive function, which is the integrand in (3.20) below. Therefore,
1 2 2 2 2 ∞ e−(m+s) − e−s e−(m−s) − e−s (3.20) I = lim IN = ds. sin πs − N →∞ [m + s]3 [m − s]3 0 m=2,even
π . N
1737
INITIAL-BOUNDARY VALUE PROBLEMS
−4
10
−5
N=32
|u(x=π/N,t) − uN (x=π/N,t)|
10
−6
10
N=64 −7
10
N=128
−8
10
−9
10
−10
10
−11
10
−7
10
−6
10
−5
−4
10
−3
10
t
−2
10
−1
10
0
10
10
Fig. 3.1. |eN | for the Fourier method on the interval t(0, 1] for x = π/N and fixed N .
Numerical evaluation of the midpoint rule (3.20) with N = 512 provides the approximate value I ≈ .0876323127. This completes the proof that convergence is at best quadratic. For each N there exists a point (π/N, 1/N 2 ) at which the error is I/N 2 . The error as a function of time for x = π/N is plotted in Figure 3.1, where the maximum occurs at t = 1/N 2 . Likewise, we can plot the error as a function of space for t = 1/N 2 and see that the maximum error occurs at the endpoints as shown in Figure 3.2(a). Numerical experiments, in the following sections, will confirm that indeed convergence is quadratic. 3.3. Quartic convergence for a fixed t > > 0. For any fixed t > > 0 in the interior of the domain, bounded away from the singularities at (0, 0) and (0, π), eN (xi , t)2 converges at best quartically. As was shown in (3.15), the coefficients of the approximate trigonometric series have an initial error of O(1/N 4 ). This, in turn, induces quartic convergence for any fixed t. The difference between (3.7) and (3.15) provides the error in the spectral approximation (3.21)
eN (xi , t) = u(xi , t) − uN (xi , t) = SN (xi , t) + EN (xi , t),
where (3.22)
SN (xi , t) =
1 N4
N −1
2
nΨ(n/N )e−n t sin n
n=1,odd
and (3.23)
EN (xi , t) =
∞ n=N +1,odd
2
e−n t iπ sin n . n3 N
iπ N
1738
NATASHA FLYER AND PAUL N. SWARZTRAUBER
We will rigorously establish quartic convergence for fixed t > 0 by demonstrating spectral convergence for EN (xi , t) and, at best, quartic convergence for SN (xi , t). Starting with (3.23), (3.24)
∞
EN (xi , t) =
2
n=1,odd
(3.25)
= e−N
2
e−(N +n) t sin(N + n)xi (N + n)3 ∞
t
n=1,odd
2
e−(2N +n )t sin(N + n)xi . (N + n)3
Therefore, (3.26)
∞
−N 2 t
|EN (xi , t)| < e
n=1,odd
2 1 < e−N t (N + n)3
∞
N
ds s3
and finally 2
e−N t |EN (xi , t)| < . 2N 2
(3.27)
Now consider the asymptotic behavior of SN (xi , t). From (3.22), (3.28)
1/2 N −1 2 1 e−t e−t SN (xi , t)2 = 4 [nΨ(n/N )e−n t ]2 > Ψ(0) 4 ≈ 0.4058712 4 , N N N n=1,odd
where the last inequality is from the earlier observation that Ψ(s) is monotone increasing on the interval [0,1] and consequently Ψ(1/N ) > Ψ(0) ≈ 0.4058712. This demonstrates that convergence in the l2 norm is at best quartic. In contrast to Figure 3.2(a), Figure 3.2(b) shows that for any fixed t, the maximum error, as measured by |eN |, occurs in the middle of the domain at x = π/2 as opposed to the endpoint, x = π/N , for small t. The reason the error curves change as t becomes large 2 is that the exponential term e−n t dominates, and the error can be approximated by the lowest wavenumber, n = 1. Thus, for large t, the error is essentially proportional to sin x which has a maximum at π/2. As will be seen in the next section, numerical experiments demonstrate that convergence in the l∞ norm is also quartic for t bounded away from (0, 0) and (0, π). 3.4. Comparison with finite difference methods. Given the nonuniform convergence of the Fourier method, as illustrated by quadratic convergence for small t and quartic convergence for large t, it is natural to compare its performance with fourth and second-order finite difference methods that require less computation. The finite difference approximations are derived by substituting (3.29)
N −1
ui =
an (t) sin nxi
n=1,odd
into (3.30)
u it =
ui−1 − 2ui + ui+1 , δx2
1739
INITIAL-BOUNDARY VALUE PROBLEMS
−6
10
−4
10
N=32
−5
10
−7
10 −6
|u(x,t=1)−u N (x,t=1)|
|u(x,t=1/N 2 ) − u N (x,t=1/N 2 )|
10
N=32
−7
10
N=64
−8
10
−8
10
N=64
−9
10
−9
10
N=128
N=128
−10
10
−10
10 −11
10
N=256 −11
−12
10
0
0.5
1
1.5
x
2
2.5
3
3.5
10
0
0.5
1
1.5
(a)
x
2
2.5
3
3.5
(b)
Fig. 3.2. (a) |eN (x, N −2 )| on the interval x(0, π) for t = 1/N 2 . (b) |eN (x, 1)| on the interval x(0, π) for t = 1.
where δx is π/N , yielding (3.31)
2
an (t) = cN (n)eλn t ,
where
λ2n =
nπ ) −4N 2 sin2 ( 2N . 2 π
Likewise, the solution to the fourth-order finite difference approximation (3.32)
u it =
−ui−2 + 16ui−1 − 30ui + 16ui+1 − ui+2 12δx2
is given by (3.31), except now (3.33)
λ2n =
nπ ) −4N 2 sin2 ( 2N 2 nπ . 3 + sin 3π 2 2N
The values of λn are the only difference between solutions to the finite difference and Fourier method for which λn = n. However, if we compare the performance of the FD4 to the Fourier method for large t, it can be seen that there is little difference, as shown in Figure 3.3(a). For small t, a FD2 method performs just as well as the Fourier method and FD4 with only a slightly larger constant of proportionality, as seen in Figure 3.3(b). This marginal error reduction may not be sufficient to justify the expense of the Fourier method. However, it is reasonable to speculate that a FD4 scheme on a Chebyshev grid could give comparable results to the Chebyshev method that yields an O(N −4 ) convergence error near the endpoints, as will be seen in the next section. 3.5. Comparison with the Chebyshev method. The method of choice for the interpolation of nonperiodic functions on a finite interval is the Chebyshev method. Unlike the Fourier method, the Chebyshev method does yield spectral convergence
1740
NATASHA FLYER AND PAUL N. SWARZTRAUBER
−3
−3
10
10
−4
10
−5
10
−4
−6
10
|e N ( π/2,1)|
|e N (pi/N,1/N 2 )|
10
1/N 4
−7
10
−5
10 −8
10
fd4 Fourier
−9
10
fd2 fd4 Fourier
−10
10
−6
0
10
1
10
2
N
(a)
10
3
10
10
0
10
1
10
2
N
10
3
10
(b)
Fig. 3.3. (a) The error of a fourth-order centered finite difference method versus the error in the spectral approximation for x = π/2, t = 1. (b) The error of a fourth-order and second-order centered finite difference method versus the error in the spectral approximation for x = π/N and t = 1/N 2 .
for fixed t. The reason is simply that the cosine series representation of a smooth IC does not alias. Nevertheless, convergence in the neighborhood of (0, 0) and (π, 0) is algebraic. To see why this is the case, we implement the Chebyshev method by expanding the solution in terms of Chebyshev cardinal functions (see [1]). The error for the Chebyshev method, as a function of space, is qualitatively similar to Figures 3.2(a) and 3.2(b), with the maximum error occurring at the endpoints for small t and in the middle of the domain for large t. The error as a function of time is similar to Figure 3.1 in that it reaches a maximum and then decays, but it differs in the fact that the decay becomes exponential after the maximum as seen in Figure 3.4. However, there exists a trajectory along which the max error converges algebraically as shown in Figure 3.5(a). The difference is that the max error is converging quartically as opposed to quadratically, as seen for the Fourier case. The increase in accuracy occurs because the Chebyshev grid is quadratically clustered near the endpoints where the singularities are occurring, giving an extra factor of O(1/N 2 ). In contrast to the Fourier method for large t, the Chebyshev method does indeed give the expected exponential rate of convergence for the error, as indicated in Figure 3.5(b). This is a direct result of the Chebyshev method’s ability to represent the IC exactly with only three polynomials. Thus unlike the Fourier method, at time t = 0, there is no initial error in the coefficients due to truncation. For the parabolic case, we can see that the convergence rate of the error in the maximum norm is nonuniform and algebraic. In general, as long as we are interested in solutions that are bounded away from singularities induced by violating the compatibility conditions, the spectral method will yield spectral convergence. However, if we are interested in heat transients such as those that occur on a microprocessor,
1741
INITIAL-BOUNDARY VALUE PROBLEMS −3
10
N=8 −4
10
N=16 −5
10
−6
|e N (π /N,t)|
10
N=32
−7
10
−8
10
−9
10
−10
10
−6
10
−5
10
−4
−3
10
−2
10
−1
10
t
0
10
10
Fig. 3.4. |eN | for the Chebyshev method on the interval t(0, 1] for x = π/N and fixed N .
CONVERGENCE OF THE MAX ERROR AT THE ENDPOINT − CHEBYSHEV
−3
CONVERGENCE FOR THE MAX ERROR FOR T=0.1
−4
10
10
−6
10 −4
10
−8
10
−10
10
−5
||e N|| ∞
10
Chebyshev 1/N 4
−12
|e N|
10
−6
10
e −2N
−14
10
−16
10 −7
10
−18
10
−8
10
−20
0
10
1
10
N
(a)
2
10
10
4
6
8
10
12
N
14
16
18
20
22
(b)
Fig. 3.5. (a) The max error in the Chebyshev approximation at the endpoint x = π/N as a function of resolution N . (b) The error in the Chebyshev approximation for x = π/2 and t = 0.1.
a computationally more efficient method with comparable algebraic convergence may be more attractive. 4. The hyperbolic case. In contrast to the parabolic case, a hyperbolic operator will propagate the solution through the time-space domain. As a result, singularities induced by incompatible BCs and the IC are not smoothed as t increases.
1742
NATASHA FLYER AND PAUL N. SWARZTRAUBER
To observe the effect on the convergence of the error for spectral methods, we will consider the one-dimensional wave equation utt = uxx
(4.1) subject to BCs (4.2)
u(0, t) = u(π, t) = 0
and ICs (4.3)
u(x, 0) =
π x(π − x), 8
ut (x, 0) = 0.
The first-order compatibility condition, (2.6) for u0 and k = 1, is not satisfied. The difference in the analysis of section 3.1 is that the time dependence of the exact and approximate spectral solution is given by an oscillatory function. Thus, (4.4)
u(xi , t) =
∞ n=1,odd
cos(nt) sin nxi n3
with the analogue to (3.15) being (4.5)
uN (xi , t) =
N −1 n=1,odd
1 n n sin nxi . cos(nt) 3 − 4 Ψ n N N
4.1. Global quadratic convergence. The decay rate of the error for hyperbolic cases is more dramatic than for parabolic problems. The error induced by violating the first-order compatibility condition for a hyperbolic problem not only produces an O(N −2 ) error for the Fourier method in the neighborhood of (0, 0) and (π, 0) but propagates throughout the domain, resulting in global quadratic convergence in the maximum norm. The propagation of the error through the x−t plane is along the characteristic x = t and x = π − t, as shown in Figures 4.1(a) and 4.1(b) for x ∈ [0, π] and t ∈ [0, π]. The proof of quadratic convergence of the error for small t follows that given in section 3.2. For large t, the proof is the same except the error is defined along the discrete trajectory xN/2 = π/2 and tN = π/2 − 1/N . Numerical experiments, given in Figure 4.2, show that convergence is quadratic along both of these trajectories, demonstrating the global nature of the error propagation. 4.2. The Chebyshev method. Next, we compare the performance of a Chebyshev method with the Fourier method that yields quadratic convergence. As noted in section 3.5, we use Chebyshev cardinal functions. Figure 4.3 shows the propagation of the error through the x−t plane for the Chebyshev method. As expected, we see that the error propagates along the characteristics t = 1 ± x. The convergence of the error along the characteristics for any value of t differs significantly from the parabolic case, where spectral convergence was obtained for t > 0. The convergence of the error, as shown in Figure 4.4, is quadratic, no better than the Fourier series expansion and requiring O(N 2 ) operations. At first, one might suspect that the convergence rate would be O(N 4 ), the same as was observed for the heat equation in the neighborhood
1743
INITIAL-BOUNDARY VALUE PROBLEMS
3 −4
x 10 2
2.5
2
1
t
|e N=32 (x,t)|
1.5
1.5
0.5
0 3.5
1
3 2.5 3.5
2
0.5
3 2.5
1.5 2
1
1.5 1
0.5
t
0
0
0.5 0
x
0.5
1
(a)
1.5
2
x
2.5
3
(b)
Fig. 4.1. (a) Propagation of the error through the x−t plane. (b) Contour plot of (a).
−1
10
−2
10
−3
|e N|
10
−4
10
1/N2
|eN(x=π/2,t=π/2−1/N)| −5
10
|eN(x=π/N,t=1/N)| −6
10
0
10
1
10
2
N
10
3
10
Fig. 4.2. The convergence of the error in the Fourier method on the discrete trajectories x = π/2 and t = π/2 − 1/N , and x = π/N and t = 1/N .
of the discontinuities (x = 0, t = 0) and (x = π, t = 0). The difference, however, is that in the heat equation there is only a jump discontinuity in the second derivative of the solution at t = 0. For any t > 0, the discontinuity becomes smoothed out due to the diffusive operator. Thus, in the neighborhood of the singularities for 0 < t < , the original jump discontinuities become steep gradients that are analytically continuous.
1744
NATASHA FLYER AND PAUL N. SWARZTRAUBER 2
1.8
1.6
1.4
t
1.2
1
0.8
0.6
0.4
0.2
−0.8
−0.6
−0.4
−0.2
0
x
0.2
0.4
0.6
0.8
Fig. 4.3. The propagation of the error in the Chebyshev approximation for N = 128.
−1
10
−2
|eN(x=t)|
10
1/N2
−3
10
Chebyshev
−4
10
0
10
1
10
N
2
10
Fig. 4.4. The convergence of the error in the Chebyshev approximation along the characteristic t = x.
Spectral methods have difficulty resolving these gradients, but the Chebyshev method gives quartic convergence due to its quadratically clustered grid near the endpoints, as opposed to the uniform grid of the Fourier method. However, in the hyperbolic case, the discontinuities at t = 0 in the second derivative of the analytical solution are present for all time. Thus, the type of grid used is irrelevant, and any spectral method will converge quadratically.
1745
INITIAL-BOUNDARY VALUE PROBLEMS
4.3. Comparison with finite difference methods. The quadratic convergence of spectral methods naturally refocuses attention on second-order finite difference methods. The second-order centered finite difference method given by (3.30) yields the solution uN (xi , t) =
(4.6)
N
cos(λn t)
n=1,odd
where (4.7)
λn =
1 n − 4 Ψ(n/N ) sin nxi , n3 N
nπ 2N sin π 2N
is a second-order approximation to the eigenvalues n. The error is then given by (4.8)
eN (xi , t) = u(xi , t) − uN (xi , t) = TN (xi , t) + SN (xi , t) + EN (xi , t),
where N
TN (xi , t) =
(4.9)
n=1,odd
SN (xi , t) =
(4.10) (4.11)
cos nt − cos λn t sin nxi , n3
N
1 N4
nΨ(n/N ) cos λn t sin nxi ,
n=1,odd ∞
EN (xi , t) =
n=N +1,odd
cos nt sin nxi . n3
Using these formulas, we will analytically prove as well as computationally show that not only does a second-order finite difference method not give quadratic convergence but rather converges only slightly better than linearly. The dominant error term is TN (xi , t) which we demonstrate by first showing that |SN (xi , t)| < O(N −2 )
(4.12)
and
|EN (xi , t)| < O(N −2 ).
For t = x, we have |SN (t)| =
(4.13)
N 1 | nΨ(n/N ) cos λn t sin nt| 2N 4 n=1,odd
≤
4/3 , π 0 N where C is a constant. These bounds hold for any second-order finite difference operator in which there is a discontinuity in the second derivative of the solution. The fact that the convergence of the second-order finite difference scheme is not quadratic may not be surprising since the approximation is based on Taylor series expansions that assume continuity of the function at least up to the fourth derivative. However, the fact that it is 4/3 is surprising. For a variety of ICs that would correspond to violating different order compatibility conditions we propose the following conjecture for the upper bounds of the error term, TN (xi , t), due to a second-order finite difference approximation. Let (4.17)
TN (t) =
N n=1,odd
cos λn,N t − cos nt sin nt = n3
N
bn,N (t).
n=1,odd
When we enforce BCs, we violate some order compatibility condition and the coefficients of the error, |bn,N (t)|, will be algebraically decreasing due to the resulting
INITIAL-BOUNDARY VALUE PROBLEMS
1747
discontinuity in a derivative of the Taylor expansion of the IC about the boundary. Thus, we can always impose the bound |bn,N (t)| ≤
(4.18)
c , nβ
β ≥ 1,
where c is a constant. From (4.18) and generalizing the analysis in Appendix A we make the following conjecture whose derivation is relegated to Appendix B. Conjecture. The error term TN (xi , t) is bounded from above by N
n=1
|bn,N (t)|
≤ ≤ ≤
O(N (1−β)2/3 ), O(N −2 log N ), O(N −2 ),
β < 4, β = 4, β > 4.
For β < 4, the error from the IC dominates, while for β = 4 both the error from the IC and second-order difference approximation contribute equally. In the last case the error from the finite difference operator dominates. β > 4 implies that at least the fourth derivative of the IC is continuous. Thus, the best a second-order method can perform is O(N −2 ). β < 4 implies that the function is not C 4 , and therefore we would expect less than quadratic convergence since the error in a second-order method depends on the continuity of the fourth derivative of the function. However, the fact that β = 4 yields a convergence rate that is proportional to log(N ) is not only new but quite surprising. The above conjecture has been supported by the numerical evidence. For completeness, the convergence rate of the FD4 method that was used in section 3.4 is considered. However, we still see that convergence is less than quadratic and in fact appears to be only O(N −3/2 ), as demonstrated in Figure 4.5(b). The reduction in the convergence rates of the finite difference methods demonstrates the significant impact of incompatibilities which occur anytime the BCs are independent of the IC. In general, violation of the compatibility conditions will likely lead to algebraic convergence of spectral methods as well as reduce the convergence rate of finite difference schemes. 5. Conclusions. In this paper, we have applied the theory of compatibility conditions to the numerical solutions of PDEs and have demonstrated how the subtle singularities inherent in IBVPs impact the convergence properties of spectral and finite difference methods. The temporal derivatives of the solution as defined by the initial condition, differential operator, and Cauchy–Kowalesky theorem will not equal those determined by the independent BCs. The resulting singularities in the corners of the temporal-spatial domain, which are independent of the smoothness of the IC, disrupt the spectral convergence of the error normally associated with spectral methods, in a manner that differs significantly between parabolic and hyperbolic systems. For parabolic systems, the Chebyshev spectral method may be considered “selfhealing,” leading to a spectral rate of convergence for the error. However, this is not true for the Fourier method where aliasing is likely to induce an algebraic error in the initial Fourier representation. In the neighborhood of the singularities, all methods yield convergence rates that are algebraic. Thus, as measured in the maximum norm, the convergence of the approximate solution to the exact solution for spectral methods is both algebraic and nonuniform. Therefore, the usefulness of spectral methods may be diminished if we are interested in transient solutions for small t.
1748
NATASHA FLYER AND PAUL N. SWARZTRAUBER
For hyperbolic systems, the situation is more interesting but worse because the singularities at t = 0 on the boundary propagate throughout the temporal-spatial domain along the characteristic lines x ± ct. Thus, the algebraic rate of convergence for the error that existed in the neighborhood of the singularities for the parabolic case exists for all time in hyperbolic problems. As a result, the spectral method does not yield spectral convergence for hyperbolic IBVPs. This naturally draws attention to finite difference methods which require only O(N ) operations. However, it was shown that violation of the compatibility conditions can also reduce the order of convergence for finite difference schemes. The above study details the impact of the theory of compatibility conditions on numerical calculations for IBVPs. Application is broad because many scientific problems are calculated on finite domains where independent BCs are imposed. As a result, singularities in the corners of the temporal-spatial domain are likely to disrupt the accuracy of numerical methods. It is in this regard that the above study has analyzed the convergence rate of the error for spectral and finite difference methods and explored their application to solving IBVPs. Appendix A: Proof of 4/3 convergence. The term TN (xi , t), which accounts for the error in the eigenvalue approximation, dominates the error of the second-order finite difference scheme. Our goal is therefore to bound the error from above and below to prove the 4/3 convergence rate that was computationally observed in section 4.3. We will first prove the upper bound for the error term TN (x, t) along the characteristic t = x for x (0, π) |TN (t)| ≤ O(N −4/3 ).
(5.1) Proof. (5.2)
TN (t) =
N n=1,odd
cos λn,N t − cos nt sin nt = n3
N
bn,N (t),
n=1,odd
where (5.3)
λn,N = n
nπ ) sin( 2N nπ 2N
.
Since the lim supn→N cos(c(n)) sin(d(n)) = 1, we can immediately achieve an upper bound on the coefficients bn,N (5.4)
|bn,N | ≤
2 . n3
This bound is due to the discontinuity in the second derivative of the IC from violating the compatibility conditions. However, because the series is decreasing slowly as 1/n3 , this bound is only tight for the tail end of the series, when n is large. To describe the behavior for the initial terms (i.e., small n), we need to derive a secondary bound. This bound will come from the error in the second-order finite difference operator. Using trigonometric identities, we can rewrite (5.5)
cos λn,N t − cos nt = 2 sin
n − λn,N n + λn,N t sin t. 2 2
1749
INITIAL-BOUNDARY VALUE PROBLEMS
To find an upper bound on the product of the sine functions, we use the known result derived from the Taylor expansion of sin y/y for y [0, π] which gives (5.6)
| sin
n − λn,N π 2 n3 n − λn,N t| < t< t. 2 2 24 N 2
Therefore, |bn,N | ≤
(5.7)
|(cos λn,N t − cos nt) sin nt| π2 t ≤ . 3 n 24 N 2
The trick now is to divide the summation in (5.2) into two parts, γ
N
(5.8)
N
bn,N (t) =
n=1,odd
N
bn,N (t) +
bn,N (t),
n=N γ +1,odd
n=1,odd
where 0 < γ < 1. Equation (5.7) gives the tightest bound for the initial terms, while (5.4) gives a tighter bound for the tail end of the series. The bounds are minimized, giving equal contributions, for γ = 2/3. Therefore, we have N
(5.9) |TN (t)| = |
bn,N (t)| ≤
n=1,odd
2/3 N
n=1,odd
(5.10)
≤ N 2/3
(5.11)
≤
N
|bn,N (t)| +
|bn,N (t)|
n=N 2/3 +1,odd
π2 t + 24 N 2
N
2 3 n n=N 2/3 +1,odd 1 π2 t N −4/3 . + 4/3 = 1 + 24 N
π2 t 24 N 4/3
For the second part of the proof, we show that TN (x, t) for any average value of t along the characteristic t = x on x (0, π) is bounded from below by 1 π C (5.12) TN (t)dt > 4/3 . π 0 N Proof. Beginning with 1 π
(5.13)
0
π
TN (t)dt =
1 π
N n=1,odd
0
π
cos nt − cos λn t sin nt, n3
then by using trigonometric identities and integrating, we have (5.14)
(5.15)
1 π
0
π
TN (t)dt =
1 2
1 ≥ 2
N n=1,odd 2/3 N
n=1,odd
f ( n(1+λ2 n /n) π) + f ( n(1−λ2 n /n) π) n3 f ( n(1−λ2 n /n) π) , n3
where (5.16)
f (y) =
1 − cos(y) . y
1750
NATASHA FLYER AND PAUL N. SWARZTRAUBER
Next, we will bound the argument of f which is of the form 1 − siny y as shown below. By inspection, we can show that on y (0, π), y 2 /12 < 1 − siny y < y 2 /6, which yields the following bounds: nπ sin 2N π n(1 − λn /n) π = n 1 − nπ (5.17) 2 2 2N nπ π 3 n3 nπ 2N = 2 6 48 N 2 nπ π 3 n3 nπ 2N = (5.19) . > 2 12 96 N 2 In the range 0 < y < 1, f (y) can be bounded by its argument 1/4y < f (y) < 1/2y. Thus, for n < N 2/3 n(1 − λn /n) 1 n(1 − λn /n) (5.20) π 4.
(5.24)
|bn,N (t)| ≤
Similarly, for the tail end of the series as N becomes large, we have, using (5.23), (5.25)
∞ n=N γ +1,odd
∼
c N γ(1−β) . β−1
From these approximations, the conjecture in section 4.3 follows for the upper bound on the error term, TN (xi , t), due to the second-order finite difference scheme.
INITIAL-BOUNDARY VALUE PROBLEMS
1751
Acknowledgments. The authors gratefully thank Dr. Keith Lindsay at the National Center for Atmospheric Research for his discussions, help, and contributions, particularly with the convergence proofs in the appendices and the factorizations (5.8) and (5.9). REFERENCES [1] J.P. Boyd, Chebyshev and Fourier Spectral Methods, Springer-Verlag, New York, 1989. [2] J.P. Boyd and N. Flyer, Compatibility conditions for time-dependent partial differential equations and the rate of convergence of Chebyshev and Fourier spectral methods, Comput. Methods Appl. Mech. Engrg., 175 (1999), pp. 281–309. [3] A. Friedman, Partial Differential Equations of Parabolic Type, Prentice-Hall, Englewood Cliffs, NJ, 1964. [4] H.O. Kreiss, Initial boundary value problems for hyperbolic systems, Comm. Pure Appl. Math, 23 (1970), pp. 277–298. [5] H.O. Kreiss, V. Thomee, and O. Widlund, Smoothing of the initial data and rates of convergence for parabolic difference equations, Comm. Pure Appl. Math., 23 (1970), pp. 241–259. [6] H.O. Kreiss and J. Lorenz, Initial-Boundary Value Problems and the Navier-Stokes Equations, Academic Press, Boston, 1989. [7] R. Hersh, Mixed problems in several variables, J. Math. Mech., 12 (1963), pp. 317–334. [8] D. Gottlieb and S.A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, SIAM, Philadelphia, 1977. [9] B. Gustafsson, H.O. Kriess, and J. Oliger, Time-Dependent Problems and Difference Methods, Wiley-Interscience, New York, 1995. [10] O. Ladyzenskaja, On the convergence of Fourier series defining a solution of a mixed problem for hyperbolic equations, Dokl. Akad. Nauk SSSR, 85 (1952), pp. 481–484. [11] O. Ladyzenskaja, On the solvability of the fundamental boundary problems for equations of parabolic and hyperbolic type, Dokl. Akad. Nauk SSSR, 87 (1954), pp. 395–398. [12] O. Ladyzenskaja, V. Solonnikov, and N. Ural’ceva, Linear and Quasi-Linear Equations of Parabolic Type, Transl. Math. Monogr. 23, AMS, Providence, RI, 1968. [13] A. Majda, J. McDonough, and S. Osher, The Fourier method for nonsmooth initial data, Math. Comp., 32 (1978), pp. 1041–1081. [14] A. Majda and S. Osher, Propagation of error into regions of smoothness for accurate difference approximations to hyperbolic equations, Comm. Pure Appl. Math., 30 (1977), pp. 671–705. [15] S. Osher, Smoothing for spectral methods, in Spectral Methods for Partial Differential Equations, SIAM, Philadelphia, 1984, pp. 209–216. [16] J.B. Rauch and F.J. Massey, Differentiability of solutions to hyperbolic initial-boundary value problems, Trans. Amer. Math. Soc., 189 (1974), pp. 303–318. [17] R. Sakamoto, Hyperbolic Boundary Value Problems, Cambridge University Press, Cambridge, UK, 1982. [18] S. Smale, Smooth solutions of the heat and diffusion equations, Comment. Math. Helv., 55 (1980), pp. 1–12. [19] V.I. Smirnov, A Course of Higher Mathematics. Integral Equations and Partial Differential Equations, Vol. 4, Pergamon Press, London, 1964. [20] R. Temam, Behaviour at time t = 0 of the solutions of semi-linear evolution equations, J. Differential Equations, 43 (1982), pp. 73–92. [21] V. Thomee, Numerical Solution of Partial Differential Equations II, Academic Press, 1970, pp. 585–622. [22] V. Thomee and L. Wahlbin, Convergence rates of parabolic difference schemes for non-smooth data, Math. Comp., 28 (1974), pp. 1–13.
c 2002 Society for Industrial and Applied Mathematics
THE CONVERGENCE OF SPECTRAL AND FINITE DIFFERENCE METHODS FOR INITIAL-BOUNDARY VALUE PROBLEMS∗ NATASHA FLYER† ‡ AND PAUL N. SWARZTRAUBER† Abstract. The general theory of compatibility conditions for the differentiability of solutions to initial-boundary value problems is well known. This paper introduces the application of that theory to numerical solutions of partial differential equations and its ramifications on the performance of high-order methods. Explicit application of boundary conditions (BCs) that are independent of the initial condition (IC) results in the compatibility conditions not being satisfied. Since this is the case in most science and engineering applications, it is shown that not only does the error in a spectral method, as measured in the maximum norm, converge algebraically, but the accuracy of finite differences is also reduced. For the heat equation with a parabolic IC and Dirichlet BCs, we prove that the Fourier method converges quadratically in the neighborhood of t = 0 and the boundaries and quartically for large t when the first-order compatibility conditions are violated. For the same problem, the Chebyshev method initially yields quartic convergence and exponential convergence for t > 0. In contrast, the wave equation subject to the same conditions results in inferior convergence rates with all spectral methods yielding quadratic convergence for all t. These results naturally direct attention to finite difference methods that are also algebraically convergent. In the case of the wave equation, we prove that a second-order finite difference method is reduced to 4/3-order convergence and numerically show that a fourth-order finite difference scheme is apparently reduced to 3/2-order. Finally, for the wave equation subject to general ICs and zero BCs, we give a conjecture on the error for a second-order finite difference scheme, showing that an O(N −2 log N ) convergence is possible. Key words. compatibility conditions, convergence theory, spectral and finite difference methods AMS subject classifications. 65M06, 65M12, 65M70 PII. S1064827500374169
1. Introduction. The regularity of solutions for initial-boundary value problems (IBVPs) is determined by the compatibility of the initial condition (IC) and the boundary conditions (BCs). IBVPs will have singular solutions at the corners of the space-time domain if the BCs do not satisfy the partial differential equation (PDE) or any of its higher-order derivatives, even if the IC is C ∞ . The Cauchy–Kowalesky theorem [19] provides a solution in the neighborhood of t = 0 to the same problem without BCs. If compatibility exists, then the BCs can always be computed by the solution obtained from the Cauchy–Kowalesky theorem and the IC. In such cases, the IBVP could in fact be posed as an initial-value problem (IVP). The set of compatibility conditions are derived by equating the time derivatives of the BCs to those of the solution to the IVP given by the Cauchy–Kowalesky theorem. Literature on the theory of compatibility conditions and the regularity of solutions for IBVPs is extensive, starting in the 1950s with the work of Ladyzenskaja [10], [11]. In the 1960s, Ladyzenskaja, Solonnikov, and Ural’ceva [12] and Friedman [3] discussed the regularity of the solution for linear and quasi-linear parabolic systems. During this same time, Kreiss [4] and Hersh [7] extended the theory for constant coefficient hyperbolic systems by providing the necessary and sufficient conditions for ∗ Received by the editors June 20, 2000; accepted for publication (in revised form) August 6, 2001; published electronically January 30, 2002. http://www.siam.org/journals/sisc/23-5/37416.html † National Center for Atmospheric Research, Advanced Study Program, P.O. Box 3000, Boulder, CO 80307-3000 ([email protected]). NCAR is sponsored by the National Science Foundation. ‡ Current address: Department of Applied Mathematics, University of Colorado, Boulder, CO 80309 (fl[email protected]).
1731
1732
NATASHA FLYER AND PAUL N. SWARZTRAUBER
the problem to be well-posed. It was further developed in the 1970s by Rauch and Massey [16] and Sakamoto [17] for hyperbolic systems with time-dependent coefficients. In the 1980s, Temam [20] extended the theory by obtaining the compatibility conditions for semilinear evolution equations while Smale [18] gave a rigorous proof of the compatibility conditions for the heat and wave equation on a bounded domain. These accomplishments have been confined to the realm of theoretical mathematics. Yet, not satisfying the compatibility conditions of a system can have a significant computational impact, particularly for high-order accurate methods, such as spectral methods, whose performance intrinsically depends on the smoothness of the solution. For nonsmooth initial data, there is a wealth of information on the convergence theory of numerical schemes in a variety of norms for parabolic and hyperbolic IVPs [13], [5], [21], [22]. These analyses can be applied to the IBVPs if the IBVP can be posed as an IVP with discontinuous ICs through periodic extensions of the ICs to the entire real axis [6], [9], [14]. However, such periodic extensions may induce singularities that do not exist in a smooth nonperiodic solution of an IBVP. More recently, the impact of noncompatible BCs has been explored by Boyd and Flyer [2] in a computational framework. This paper appeared in a special issue on spectral methods. In the prologue, Karniadakis states “This is the first such effort in bringing the theory of compatibility conditions from the mathematical literature to the numerical community.” This work illustrated the ubiquity of compatibility conditions and analyzed the connection between incompatibility and the rate of convergence of Chebyshev spectral series. However, the convergence of the error was not discussed or compared to alternatives such as the finite difference method. The focus of the paper, as has been in the mathematical literature, is on smoothing the initial condition so as to satisfy the compatibility conditions [2], [5], [18], [15]. The difficulty with this approach is that it leaves the fundamental question unanswered: namely, “What is the rate of convergence to the solution of the original unperturbed problem?” If the IC is smoothed, spectral convergence might be achieved but to a solution that differs algebraically from the original problem. The focus of the current work is on the convergence of the approximate solution to the exact solution. As discussed in the first paragraph, the temporal derivatives of the solution as defined by the IC, differential operator, and Cauchy–Kowalesky theorem will not equal those determined by the independent BCs. The resulting singularities in the corners of the temporal-spatial domain, which are independent of the smoothness of the IC, disrupt the convergence of spectral and finite difference methods in a manner that differs significantly between parabolic and hyperbolic systems. A detailed analysis of the induced singularities for two model problems and their effect on the convergence of both numerical methods is determined. Since an incompatibility exists any time the BC is independent of the IC, algebraic convergence in the maximum norm is the expectation for IBVPs. In most science and engineering applications, this is the usual scenario. In section 2, the infinite set of compatibility conditions for the generalized heat and wave equations with time-dependent BCs is presented. The purpose is to lay the groundwork for describing the inherent singular nature of IBVPs and how it affects the convergence rate of the error. In section 3, we examine the prototype parabolic case, namely the heat equation subject to Dirichlet BCs and a smooth parabolic IC. We use the Fourier method to determine an exact expression for the approximate solution, represented by a truncated Fourier series. Convergence rates are then developed for both finite difference and Chebyshev spectral methods, the latter of which is normally the basis set of choice for a problem on a finite domain. For the parabolic example
INITIAL-BOUNDARY VALUE PROBLEMS
1733
in section 3, ut is discontinuous at t = 0. Nevertheless, the Chebyshev method is shown to yield spectral convergence for t > 0, but only quartic convergence in the neighborhood of the discontinuity, which would be the relevant consideration for the computation of transient heat flow. Furthermore, the usual Fourier method is only quadratically convergent near the discontinuity and quartically convergent elsewhere. Thus, the convergence rate of the error is shown to be nonuniform and algebraic as measured in the maximum norm. The hyperbolic case is discussed in section 4. Discontinuities that are induced by the incompatibilities now propagate throughout the region. It is shown that all spectral methods will yield algebraic convergence for t > 0. The Fourier method is used again to derive an exact expression for the error. Similar to the parabolic case, the results are compared with the centered finite difference and Chebyshev methods which yield algebraic convergence for smooth ICs. The performance of finite difference methods is also significantly diminished when the compatibility conditions are not satisfied. For this reason, the appendix contains the rates of convergence for a second-order finite difference method for a variety of smooth ICs and shows, for example, that second-order convergence is reduced to 4/3-order. 2. Compatibility conditions for the generalized heat and wave equations. Here compatibility conditions are reviewed to illustrate their fundamental impact on the smoothness of the solutions to IBVPs. We simply state the necessary conditions for the solution to be C ∞ on the domain for the heat and wave equation with time-dependent BCs. The complete proofs of necessity and sufficiency are given for parabolic systems in [11], [12], [3], [18] and for hyperbolic systems in [10], [16], [17], [4]. Theorem 2.1. The domain is [0, T ] × Ω, where Ω is a d-dimensional spatial domain with a boundary ∂Ω which is a C ∞ manifold of dimension (d − 1). Let L be an elliptic operator of the form (2.1)
L=
d d
d
Aij (x)
i=1 j=1
d
∂ ∂ ∂ + Bj (x) + Cj (x). ∂xi ∂xj j=1 ∂xj j=1
The generalized linear diffusion problem is ut = Lu,
(2.2)
u(x, 0) = u0 (x) ∈ Ω.
The generalized wave equation is utt = Lu,
(2.3)
u(x, 0) = u0 (x),
ut (x, 0) = v0 (x) ∈ Ω.
Both are subject to the BC (2.4)
u = f (t) ∈ ∂Ω ∀t.
Then, the necessary and sufficient compatibility conditions for u(x,t) to be C 2k for the heat equation are Lk u0 =
(2.5)
∂kf , ∂tk
k = 0, 1, 2, . . . ∈ [0, 0] × ∂Ω
and for the wave equation are (2.6)
Lk u0 =
∂ 2k f ∂t2k
and
Lk v0 =
∂ 2k+1 f , ∂t2k+1
k = 0, 1, 2, . . . ∈ [0, 0] × ∂Ω.
1734
NATASHA FLYER AND PAUL N. SWARZTRAUBER
In essence (2.5) and (2.6) state that at time t = 0, the temporal derivatives of the solution determined by the application of the differential operator to the IC must equal the temporal derivatives of the BCs. When these conditions are satisfied, we do not need to impose BCs because they can be derived from the Taylor series expansion of the IC. When (2.5) and (2.6) do not hold, we induce a singularity in the 2kth-order derivative of the solution at t = 0 on ∂Ω. Thus, imposing BCs that are independent of the IC places an infinite set of constraints on the solution which differ from the constraints enforced by the PDE, resulting in discontinuities in the corners of the temporal-spatial domain. The manner in which these incompatibilities affect the accuracy of the numerical method depends on whether the system is parabolic or hyperbolic. In the next section we will discuss the parabolic case, followed in section 4 by the hyperbolic case. 3. A parabolic example. Our goal is to study the convergence of the error for approximate spectral and finite difference solutions to IBVPs when the compatibility conditions of the system are not satisfied. With this in mind, we first study the effect on parabolic systems, taking as our example the one-dimensional heat equation ut = uxx
(3.1) subject to BCs (3.2)
u(0, t) = u(π, t) = 0
and IC (3.3)
u0 = u(x, 0) =
π x(π − x). 8
The first-order compatibility condition, (k = 1) in (2.5), is not satisfied because 2 π ∂ u0 =− Lu0 = (3.4) ∀x 2 ∂x 4 and ∂u |x=0,π = 0 ∂t
(3.5)
∀t.
Thus from (3.1), uxx |x=0,π = ut |x=0,π = 0 yet uxx = −2 in the interior. To see the impact of the jump discontinuity at the boundary in the second derivative at t = 0 on the decay rate of the error we will solve the problem exactly at time t = 0 via a truncated Fourier series, then compare this together with finite difference and Chebyshev methods against the exact solution for all time. Time integration is assumed exact throughout, which yields the approximate spectral representation (3.6)
N −1
uN (xi , t) =
2
cN (n)e−n t sin nxi
n=1,odd
of the exact solution, given by (3.7)
u(xi , t) =
∞ n=1,odd
2
e−n t sin nxi . n3
1735
INITIAL-BOUNDARY VALUE PROBLEMS
Gottlieb and Orszag [8] note that for the exact coefficients cN (n), uN (x, t) converges spectrally to u(x, t). This is also true if u(x, 0) is known for all x because cN (n) can be computed to any accuracy. In this work, we assume, as with most scientific data, that u(x, 0) is given in tabular form ui = u(xi , t), where xi = iπ/N . In science and engineering, we cannot assume we know the function exactly and therefore must approach the problem from a numerical standpoint. The coefficients cN (n) are then computed in the traditional manner, using trigonometric interpolation, as cN (n) =
(3.8)
N −1 2 ui sin nxi . N i=1
With this implementation of the spectral method, it will be shown that the approximate spectral representation (3.6) converges algebraically to the exact solution in the maximum norm, and, furthermore, convergence is nonuniform. This characteristic algebraic convergence is method-independent due to the inherent singularity induced by violating the compatibility conditions from explicit application of the BCs. In the next section, we derive an exact expression for cN (n) which will provide the initial error in the coefficients for the Fourier method. 3.1. The approximate spectral solution. In what follows, we derive an analytical representation of the coefficients cN (n) in (3.6). From (3.7) ∞ 1 iπ u(xi , 0) = (3.9) , sin n n3 N n=1,odd
(3.10)
N −1
u(xi , 0) =
n=1,odd ∞
1 iπ sin(n ) n3 N
1 iπ sin (2mN + n)i (2mN + n)3 N m=1 1 iπ + , sin (2mN − n)i 3 (2mN − n) N
+
or (3.11)
u(xi , 0) =
N −1 n=1,odd
∞ 1 1 1 + − 3 3 n (2mN + n) (2mN − n)3 m=1
iπ sin n N
and (3.12)
u(xi , 0) =
N −1 n=1,odd
1 n iπ , − 4 Ψ(n/N ) sin n n3 N N
where (3.13)
Ψ(s) = 2
∞ 3m2 + s2 . [m2 − s2 ]3 m=1
Ψ(s) can be called the alias function since it represents the error in the coefficients due to truncating the series after N terms. The initial error in the coefficients is due
1736
NATASHA FLYER AND PAUL N. SWARZTRAUBER
to the fact that a sine series cannot represent a nonperiodic polynomial with a finite number of terms. For use later, we observe that Ψ(s) is positive and monotonically increasing on the interval [0, 1] with extremes Ψ(0) ≈ 0.4058712 and Ψ(1) = 1. Equating (3.6) and (3.12), we obtain the desired result 1 nΨ(n) − . 3 n N4 We study convergence only as it relates to spatial discretization, and therefore temporal integrations are represented analytically. The spectral approximation is then given by N −1 2 1 nΨ(n/N ) uN (xi , t) = sin nxi . (3.15) e−n t 3 − n N4 (3.14)
cN (n) =
n=1,odd
In the next two sections it will be proved that although the spectral series for both the exact and approximate solutions decay exponentially, the error decays algebraically and nonuniformly over the domain. Section 3.2 explores convergence in the neighborhood of the discontinuity induced by the incompatibility of imposing BCs. Section 3.3 looks at convergence in the interior of the domain, bounded away from the singularities at (0, 0) and (0, π). 3.2. Quadratic convergence near (0, 0) and (0, π). The discontinuity in ut at (0, 0) induces an error in ∂uN /∂t(π/N, 0) that is bounded away from zero for all N. In what follows, we show that this produces an O(N −2 ) error in the neighborhood of (0, 0) and likewise near (π, 0). With a derivation similar to that preceding (3.12), the exact solution is given by
2 2 2 N −1 ∞ e−n t e−(mN +n) t e−(mN −n) t sin nxi . (3.16) u(xi , t) = − − n3 [mN + n]3 [mN − n]3 m=2,even n=1,odd
Subtracting uN (xi , t), given by (3.15), defines the error (3.17) eN (xi , t) =
N −1 n=1,odd
∞
m=2,even
2
e−(mN +n) t − e−n [mN + n]3
2
t
2
e−(mN −n) t − e−n − [mN − n]3
2
t
sin nxi .
Then, on the discrete trajectory x1 = π/N and tN = 1/N 2 π 1 1 = 2 IN , (3.18) eN , N N2 N where (3.19) IN
1 = N
N −1 n=1,odd
∞
m=2,even
2
2
2
e−(m+n/N ) − e−(n/N ) e−(m−n/N ) − e−(n/N ) − [m + n/N ]3 [m − n/N ]3
2
sin n
IN defines the midpoint quadrature of the smooth positive function, which is the integrand in (3.20) below. Therefore,
1 2 2 2 2 ∞ e−(m+s) − e−s e−(m−s) − e−s (3.20) I = lim IN = ds. sin πs − N →∞ [m + s]3 [m − s]3 0 m=2,even
π . N
1737
INITIAL-BOUNDARY VALUE PROBLEMS
−4
10
−5
N=32
|u(x=π/N,t) − uN (x=π/N,t)|
10
−6
10
N=64 −7
10
N=128
−8
10
−9
10
−10
10
−11
10
−7
10
−6
10
−5
−4
10
−3
10
t
−2
10
−1
10
0
10
10
Fig. 3.1. |eN | for the Fourier method on the interval t(0, 1] for x = π/N and fixed N .
Numerical evaluation of the midpoint rule (3.20) with N = 512 provides the approximate value I ≈ .0876323127. This completes the proof that convergence is at best quadratic. For each N there exists a point (π/N, 1/N 2 ) at which the error is I/N 2 . The error as a function of time for x = π/N is plotted in Figure 3.1, where the maximum occurs at t = 1/N 2 . Likewise, we can plot the error as a function of space for t = 1/N 2 and see that the maximum error occurs at the endpoints as shown in Figure 3.2(a). Numerical experiments, in the following sections, will confirm that indeed convergence is quadratic. 3.3. Quartic convergence for a fixed t > > 0. For any fixed t > > 0 in the interior of the domain, bounded away from the singularities at (0, 0) and (0, π), eN (xi , t)2 converges at best quartically. As was shown in (3.15), the coefficients of the approximate trigonometric series have an initial error of O(1/N 4 ). This, in turn, induces quartic convergence for any fixed t. The difference between (3.7) and (3.15) provides the error in the spectral approximation (3.21)
eN (xi , t) = u(xi , t) − uN (xi , t) = SN (xi , t) + EN (xi , t),
where (3.22)
SN (xi , t) =
1 N4
N −1
2
nΨ(n/N )e−n t sin n
n=1,odd
and (3.23)
EN (xi , t) =
∞ n=N +1,odd
2
e−n t iπ sin n . n3 N
iπ N
1738
NATASHA FLYER AND PAUL N. SWARZTRAUBER
We will rigorously establish quartic convergence for fixed t > 0 by demonstrating spectral convergence for EN (xi , t) and, at best, quartic convergence for SN (xi , t). Starting with (3.23), (3.24)
∞
EN (xi , t) =
2
n=1,odd
(3.25)
= e−N
2
e−(N +n) t sin(N + n)xi (N + n)3 ∞
t
n=1,odd
2
e−(2N +n )t sin(N + n)xi . (N + n)3
Therefore, (3.26)
∞
−N 2 t
|EN (xi , t)| < e
n=1,odd
2 1 < e−N t (N + n)3
∞
N
ds s3
and finally 2
e−N t |EN (xi , t)| < . 2N 2
(3.27)
Now consider the asymptotic behavior of SN (xi , t). From (3.22), (3.28)
1/2 N −1 2 1 e−t e−t SN (xi , t)2 = 4 [nΨ(n/N )e−n t ]2 > Ψ(0) 4 ≈ 0.4058712 4 , N N N n=1,odd
where the last inequality is from the earlier observation that Ψ(s) is monotone increasing on the interval [0,1] and consequently Ψ(1/N ) > Ψ(0) ≈ 0.4058712. This demonstrates that convergence in the l2 norm is at best quartic. In contrast to Figure 3.2(a), Figure 3.2(b) shows that for any fixed t, the maximum error, as measured by |eN |, occurs in the middle of the domain at x = π/2 as opposed to the endpoint, x = π/N , for small t. The reason the error curves change as t becomes large 2 is that the exponential term e−n t dominates, and the error can be approximated by the lowest wavenumber, n = 1. Thus, for large t, the error is essentially proportional to sin x which has a maximum at π/2. As will be seen in the next section, numerical experiments demonstrate that convergence in the l∞ norm is also quartic for t bounded away from (0, 0) and (0, π). 3.4. Comparison with finite difference methods. Given the nonuniform convergence of the Fourier method, as illustrated by quadratic convergence for small t and quartic convergence for large t, it is natural to compare its performance with fourth and second-order finite difference methods that require less computation. The finite difference approximations are derived by substituting (3.29)
N −1
ui =
an (t) sin nxi
n=1,odd
into (3.30)
u it =
ui−1 − 2ui + ui+1 , δx2
1739
INITIAL-BOUNDARY VALUE PROBLEMS
−6
10
−4
10
N=32
−5
10
−7
10 −6
|u(x,t=1)−u N (x,t=1)|
|u(x,t=1/N 2 ) − u N (x,t=1/N 2 )|
10
N=32
−7
10
N=64
−8
10
−8
10
N=64
−9
10
−9
10
N=128
N=128
−10
10
−10
10 −11
10
N=256 −11
−12
10
0
0.5
1
1.5
x
2
2.5
3
3.5
10
0
0.5
1
1.5
(a)
x
2
2.5
3
3.5
(b)
Fig. 3.2. (a) |eN (x, N −2 )| on the interval x(0, π) for t = 1/N 2 . (b) |eN (x, 1)| on the interval x(0, π) for t = 1.
where δx is π/N , yielding (3.31)
2
an (t) = cN (n)eλn t ,
where
λ2n =
nπ ) −4N 2 sin2 ( 2N . 2 π
Likewise, the solution to the fourth-order finite difference approximation (3.32)
u it =
−ui−2 + 16ui−1 − 30ui + 16ui+1 − ui+2 12δx2
is given by (3.31), except now (3.33)
λ2n =
nπ ) −4N 2 sin2 ( 2N 2 nπ . 3 + sin 3π 2 2N
The values of λn are the only difference between solutions to the finite difference and Fourier method for which λn = n. However, if we compare the performance of the FD4 to the Fourier method for large t, it can be seen that there is little difference, as shown in Figure 3.3(a). For small t, a FD2 method performs just as well as the Fourier method and FD4 with only a slightly larger constant of proportionality, as seen in Figure 3.3(b). This marginal error reduction may not be sufficient to justify the expense of the Fourier method. However, it is reasonable to speculate that a FD4 scheme on a Chebyshev grid could give comparable results to the Chebyshev method that yields an O(N −4 ) convergence error near the endpoints, as will be seen in the next section. 3.5. Comparison with the Chebyshev method. The method of choice for the interpolation of nonperiodic functions on a finite interval is the Chebyshev method. Unlike the Fourier method, the Chebyshev method does yield spectral convergence
1740
NATASHA FLYER AND PAUL N. SWARZTRAUBER
−3
−3
10
10
−4
10
−5
10
−4
−6
10
|e N ( π/2,1)|
|e N (pi/N,1/N 2 )|
10
1/N 4
−7
10
−5
10 −8
10
fd4 Fourier
−9
10
fd2 fd4 Fourier
−10
10
−6
0
10
1
10
2
N
(a)
10
3
10
10
0
10
1
10
2
N
10
3
10
(b)
Fig. 3.3. (a) The error of a fourth-order centered finite difference method versus the error in the spectral approximation for x = π/2, t = 1. (b) The error of a fourth-order and second-order centered finite difference method versus the error in the spectral approximation for x = π/N and t = 1/N 2 .
for fixed t. The reason is simply that the cosine series representation of a smooth IC does not alias. Nevertheless, convergence in the neighborhood of (0, 0) and (π, 0) is algebraic. To see why this is the case, we implement the Chebyshev method by expanding the solution in terms of Chebyshev cardinal functions (see [1]). The error for the Chebyshev method, as a function of space, is qualitatively similar to Figures 3.2(a) and 3.2(b), with the maximum error occurring at the endpoints for small t and in the middle of the domain for large t. The error as a function of time is similar to Figure 3.1 in that it reaches a maximum and then decays, but it differs in the fact that the decay becomes exponential after the maximum as seen in Figure 3.4. However, there exists a trajectory along which the max error converges algebraically as shown in Figure 3.5(a). The difference is that the max error is converging quartically as opposed to quadratically, as seen for the Fourier case. The increase in accuracy occurs because the Chebyshev grid is quadratically clustered near the endpoints where the singularities are occurring, giving an extra factor of O(1/N 2 ). In contrast to the Fourier method for large t, the Chebyshev method does indeed give the expected exponential rate of convergence for the error, as indicated in Figure 3.5(b). This is a direct result of the Chebyshev method’s ability to represent the IC exactly with only three polynomials. Thus unlike the Fourier method, at time t = 0, there is no initial error in the coefficients due to truncation. For the parabolic case, we can see that the convergence rate of the error in the maximum norm is nonuniform and algebraic. In general, as long as we are interested in solutions that are bounded away from singularities induced by violating the compatibility conditions, the spectral method will yield spectral convergence. However, if we are interested in heat transients such as those that occur on a microprocessor,
1741
INITIAL-BOUNDARY VALUE PROBLEMS −3
10
N=8 −4
10
N=16 −5
10
−6
|e N (π /N,t)|
10
N=32
−7
10
−8
10
−9
10
−10
10
−6
10
−5
10
−4
−3
10
−2
10
−1
10
t
0
10
10
Fig. 3.4. |eN | for the Chebyshev method on the interval t(0, 1] for x = π/N and fixed N .
CONVERGENCE OF THE MAX ERROR AT THE ENDPOINT − CHEBYSHEV
−3
CONVERGENCE FOR THE MAX ERROR FOR T=0.1
−4
10
10
−6
10 −4
10
−8
10
−10
10
−5
||e N|| ∞
10
Chebyshev 1/N 4
−12
|e N|
10
−6
10
e −2N
−14
10
−16
10 −7
10
−18
10
−8
10
−20
0
10
1
10
N
(a)
2
10
10
4
6
8
10
12
N
14
16
18
20
22
(b)
Fig. 3.5. (a) The max error in the Chebyshev approximation at the endpoint x = π/N as a function of resolution N . (b) The error in the Chebyshev approximation for x = π/2 and t = 0.1.
a computationally more efficient method with comparable algebraic convergence may be more attractive. 4. The hyperbolic case. In contrast to the parabolic case, a hyperbolic operator will propagate the solution through the time-space domain. As a result, singularities induced by incompatible BCs and the IC are not smoothed as t increases.
1742
NATASHA FLYER AND PAUL N. SWARZTRAUBER
To observe the effect on the convergence of the error for spectral methods, we will consider the one-dimensional wave equation utt = uxx
(4.1) subject to BCs (4.2)
u(0, t) = u(π, t) = 0
and ICs (4.3)
u(x, 0) =
π x(π − x), 8
ut (x, 0) = 0.
The first-order compatibility condition, (2.6) for u0 and k = 1, is not satisfied. The difference in the analysis of section 3.1 is that the time dependence of the exact and approximate spectral solution is given by an oscillatory function. Thus, (4.4)
u(xi , t) =
∞ n=1,odd
cos(nt) sin nxi n3
with the analogue to (3.15) being (4.5)
uN (xi , t) =
N −1 n=1,odd
1 n n sin nxi . cos(nt) 3 − 4 Ψ n N N
4.1. Global quadratic convergence. The decay rate of the error for hyperbolic cases is more dramatic than for parabolic problems. The error induced by violating the first-order compatibility condition for a hyperbolic problem not only produces an O(N −2 ) error for the Fourier method in the neighborhood of (0, 0) and (π, 0) but propagates throughout the domain, resulting in global quadratic convergence in the maximum norm. The propagation of the error through the x−t plane is along the characteristic x = t and x = π − t, as shown in Figures 4.1(a) and 4.1(b) for x ∈ [0, π] and t ∈ [0, π]. The proof of quadratic convergence of the error for small t follows that given in section 3.2. For large t, the proof is the same except the error is defined along the discrete trajectory xN/2 = π/2 and tN = π/2 − 1/N . Numerical experiments, given in Figure 4.2, show that convergence is quadratic along both of these trajectories, demonstrating the global nature of the error propagation. 4.2. The Chebyshev method. Next, we compare the performance of a Chebyshev method with the Fourier method that yields quadratic convergence. As noted in section 3.5, we use Chebyshev cardinal functions. Figure 4.3 shows the propagation of the error through the x−t plane for the Chebyshev method. As expected, we see that the error propagates along the characteristics t = 1 ± x. The convergence of the error along the characteristics for any value of t differs significantly from the parabolic case, where spectral convergence was obtained for t > 0. The convergence of the error, as shown in Figure 4.4, is quadratic, no better than the Fourier series expansion and requiring O(N 2 ) operations. At first, one might suspect that the convergence rate would be O(N 4 ), the same as was observed for the heat equation in the neighborhood
1743
INITIAL-BOUNDARY VALUE PROBLEMS
3 −4
x 10 2
2.5
2
1
t
|e N=32 (x,t)|
1.5
1.5
0.5
0 3.5
1
3 2.5 3.5
2
0.5
3 2.5
1.5 2
1
1.5 1
0.5
t
0
0
0.5 0
x
0.5
1
(a)
1.5
2
x
2.5
3
(b)
Fig. 4.1. (a) Propagation of the error through the x−t plane. (b) Contour plot of (a).
−1
10
−2
10
−3
|e N|
10
−4
10
1/N2
|eN(x=π/2,t=π/2−1/N)| −5
10
|eN(x=π/N,t=1/N)| −6
10
0
10
1
10
2
N
10
3
10
Fig. 4.2. The convergence of the error in the Fourier method on the discrete trajectories x = π/2 and t = π/2 − 1/N , and x = π/N and t = 1/N .
of the discontinuities (x = 0, t = 0) and (x = π, t = 0). The difference, however, is that in the heat equation there is only a jump discontinuity in the second derivative of the solution at t = 0. For any t > 0, the discontinuity becomes smoothed out due to the diffusive operator. Thus, in the neighborhood of the singularities for 0 < t < , the original jump discontinuities become steep gradients that are analytically continuous.
1744
NATASHA FLYER AND PAUL N. SWARZTRAUBER 2
1.8
1.6
1.4
t
1.2
1
0.8
0.6
0.4
0.2
−0.8
−0.6
−0.4
−0.2
0
x
0.2
0.4
0.6
0.8
Fig. 4.3. The propagation of the error in the Chebyshev approximation for N = 128.
−1
10
−2
|eN(x=t)|
10
1/N2
−3
10
Chebyshev
−4
10
0
10
1
10
N
2
10
Fig. 4.4. The convergence of the error in the Chebyshev approximation along the characteristic t = x.
Spectral methods have difficulty resolving these gradients, but the Chebyshev method gives quartic convergence due to its quadratically clustered grid near the endpoints, as opposed to the uniform grid of the Fourier method. However, in the hyperbolic case, the discontinuities at t = 0 in the second derivative of the analytical solution are present for all time. Thus, the type of grid used is irrelevant, and any spectral method will converge quadratically.
1745
INITIAL-BOUNDARY VALUE PROBLEMS
4.3. Comparison with finite difference methods. The quadratic convergence of spectral methods naturally refocuses attention on second-order finite difference methods. The second-order centered finite difference method given by (3.30) yields the solution uN (xi , t) =
(4.6)
N
cos(λn t)
n=1,odd
where (4.7)
λn =
1 n − 4 Ψ(n/N ) sin nxi , n3 N
nπ 2N sin π 2N
is a second-order approximation to the eigenvalues n. The error is then given by (4.8)
eN (xi , t) = u(xi , t) − uN (xi , t) = TN (xi , t) + SN (xi , t) + EN (xi , t),
where N
TN (xi , t) =
(4.9)
n=1,odd
SN (xi , t) =
(4.10) (4.11)
cos nt − cos λn t sin nxi , n3
N
1 N4
nΨ(n/N ) cos λn t sin nxi ,
n=1,odd ∞
EN (xi , t) =
n=N +1,odd
cos nt sin nxi . n3
Using these formulas, we will analytically prove as well as computationally show that not only does a second-order finite difference method not give quadratic convergence but rather converges only slightly better than linearly. The dominant error term is TN (xi , t) which we demonstrate by first showing that |SN (xi , t)| < O(N −2 )
(4.12)
and
|EN (xi , t)| < O(N −2 ).
For t = x, we have |SN (t)| =
(4.13)
N 1 | nΨ(n/N ) cos λn t sin nt| 2N 4 n=1,odd
≤
4/3 , π 0 N where C is a constant. These bounds hold for any second-order finite difference operator in which there is a discontinuity in the second derivative of the solution. The fact that the convergence of the second-order finite difference scheme is not quadratic may not be surprising since the approximation is based on Taylor series expansions that assume continuity of the function at least up to the fourth derivative. However, the fact that it is 4/3 is surprising. For a variety of ICs that would correspond to violating different order compatibility conditions we propose the following conjecture for the upper bounds of the error term, TN (xi , t), due to a second-order finite difference approximation. Let (4.17)
TN (t) =
N n=1,odd
cos λn,N t − cos nt sin nt = n3
N
bn,N (t).
n=1,odd
When we enforce BCs, we violate some order compatibility condition and the coefficients of the error, |bn,N (t)|, will be algebraically decreasing due to the resulting
INITIAL-BOUNDARY VALUE PROBLEMS
1747
discontinuity in a derivative of the Taylor expansion of the IC about the boundary. Thus, we can always impose the bound |bn,N (t)| ≤
(4.18)
c , nβ
β ≥ 1,
where c is a constant. From (4.18) and generalizing the analysis in Appendix A we make the following conjecture whose derivation is relegated to Appendix B. Conjecture. The error term TN (xi , t) is bounded from above by N
n=1
|bn,N (t)|
≤ ≤ ≤
O(N (1−β)2/3 ), O(N −2 log N ), O(N −2 ),
β < 4, β = 4, β > 4.
For β < 4, the error from the IC dominates, while for β = 4 both the error from the IC and second-order difference approximation contribute equally. In the last case the error from the finite difference operator dominates. β > 4 implies that at least the fourth derivative of the IC is continuous. Thus, the best a second-order method can perform is O(N −2 ). β < 4 implies that the function is not C 4 , and therefore we would expect less than quadratic convergence since the error in a second-order method depends on the continuity of the fourth derivative of the function. However, the fact that β = 4 yields a convergence rate that is proportional to log(N ) is not only new but quite surprising. The above conjecture has been supported by the numerical evidence. For completeness, the convergence rate of the FD4 method that was used in section 3.4 is considered. However, we still see that convergence is less than quadratic and in fact appears to be only O(N −3/2 ), as demonstrated in Figure 4.5(b). The reduction in the convergence rates of the finite difference methods demonstrates the significant impact of incompatibilities which occur anytime the BCs are independent of the IC. In general, violation of the compatibility conditions will likely lead to algebraic convergence of spectral methods as well as reduce the convergence rate of finite difference schemes. 5. Conclusions. In this paper, we have applied the theory of compatibility conditions to the numerical solutions of PDEs and have demonstrated how the subtle singularities inherent in IBVPs impact the convergence properties of spectral and finite difference methods. The temporal derivatives of the solution as defined by the initial condition, differential operator, and Cauchy–Kowalesky theorem will not equal those determined by the independent BCs. The resulting singularities in the corners of the temporal-spatial domain, which are independent of the smoothness of the IC, disrupt the spectral convergence of the error normally associated with spectral methods, in a manner that differs significantly between parabolic and hyperbolic systems. For parabolic systems, the Chebyshev spectral method may be considered “selfhealing,” leading to a spectral rate of convergence for the error. However, this is not true for the Fourier method where aliasing is likely to induce an algebraic error in the initial Fourier representation. In the neighborhood of the singularities, all methods yield convergence rates that are algebraic. Thus, as measured in the maximum norm, the convergence of the approximate solution to the exact solution for spectral methods is both algebraic and nonuniform. Therefore, the usefulness of spectral methods may be diminished if we are interested in transient solutions for small t.
1748
NATASHA FLYER AND PAUL N. SWARZTRAUBER
For hyperbolic systems, the situation is more interesting but worse because the singularities at t = 0 on the boundary propagate throughout the temporal-spatial domain along the characteristic lines x ± ct. Thus, the algebraic rate of convergence for the error that existed in the neighborhood of the singularities for the parabolic case exists for all time in hyperbolic problems. As a result, the spectral method does not yield spectral convergence for hyperbolic IBVPs. This naturally draws attention to finite difference methods which require only O(N ) operations. However, it was shown that violation of the compatibility conditions can also reduce the order of convergence for finite difference schemes. The above study details the impact of the theory of compatibility conditions on numerical calculations for IBVPs. Application is broad because many scientific problems are calculated on finite domains where independent BCs are imposed. As a result, singularities in the corners of the temporal-spatial domain are likely to disrupt the accuracy of numerical methods. It is in this regard that the above study has analyzed the convergence rate of the error for spectral and finite difference methods and explored their application to solving IBVPs. Appendix A: Proof of 4/3 convergence. The term TN (xi , t), which accounts for the error in the eigenvalue approximation, dominates the error of the second-order finite difference scheme. Our goal is therefore to bound the error from above and below to prove the 4/3 convergence rate that was computationally observed in section 4.3. We will first prove the upper bound for the error term TN (x, t) along the characteristic t = x for x (0, π) |TN (t)| ≤ O(N −4/3 ).
(5.1) Proof. (5.2)
TN (t) =
N n=1,odd
cos λn,N t − cos nt sin nt = n3
N
bn,N (t),
n=1,odd
where (5.3)
λn,N = n
nπ ) sin( 2N nπ 2N
.
Since the lim supn→N cos(c(n)) sin(d(n)) = 1, we can immediately achieve an upper bound on the coefficients bn,N (5.4)
|bn,N | ≤
2 . n3
This bound is due to the discontinuity in the second derivative of the IC from violating the compatibility conditions. However, because the series is decreasing slowly as 1/n3 , this bound is only tight for the tail end of the series, when n is large. To describe the behavior for the initial terms (i.e., small n), we need to derive a secondary bound. This bound will come from the error in the second-order finite difference operator. Using trigonometric identities, we can rewrite (5.5)
cos λn,N t − cos nt = 2 sin
n − λn,N n + λn,N t sin t. 2 2
1749
INITIAL-BOUNDARY VALUE PROBLEMS
To find an upper bound on the product of the sine functions, we use the known result derived from the Taylor expansion of sin y/y for y [0, π] which gives (5.6)
| sin
n − λn,N π 2 n3 n − λn,N t| < t< t. 2 2 24 N 2
Therefore, |bn,N | ≤
(5.7)
|(cos λn,N t − cos nt) sin nt| π2 t ≤ . 3 n 24 N 2
The trick now is to divide the summation in (5.2) into two parts, γ
N
(5.8)
N
bn,N (t) =
n=1,odd
N
bn,N (t) +
bn,N (t),
n=N γ +1,odd
n=1,odd
where 0 < γ < 1. Equation (5.7) gives the tightest bound for the initial terms, while (5.4) gives a tighter bound for the tail end of the series. The bounds are minimized, giving equal contributions, for γ = 2/3. Therefore, we have N
(5.9) |TN (t)| = |
bn,N (t)| ≤
n=1,odd
2/3 N
n=1,odd
(5.10)
≤ N 2/3
(5.11)
≤
N
|bn,N (t)| +
|bn,N (t)|
n=N 2/3 +1,odd
π2 t + 24 N 2
N
2 3 n n=N 2/3 +1,odd 1 π2 t N −4/3 . + 4/3 = 1 + 24 N
π2 t 24 N 4/3
For the second part of the proof, we show that TN (x, t) for any average value of t along the characteristic t = x on x (0, π) is bounded from below by 1 π C (5.12) TN (t)dt > 4/3 . π 0 N Proof. Beginning with 1 π
(5.13)
0
π
TN (t)dt =
1 π
N n=1,odd
0
π
cos nt − cos λn t sin nt, n3
then by using trigonometric identities and integrating, we have (5.14)
(5.15)
1 π
0
π
TN (t)dt =
1 2
1 ≥ 2
N n=1,odd 2/3 N
n=1,odd
f ( n(1+λ2 n /n) π) + f ( n(1−λ2 n /n) π) n3 f ( n(1−λ2 n /n) π) , n3
where (5.16)
f (y) =
1 − cos(y) . y
1750
NATASHA FLYER AND PAUL N. SWARZTRAUBER
Next, we will bound the argument of f which is of the form 1 − siny y as shown below. By inspection, we can show that on y (0, π), y 2 /12 < 1 − siny y < y 2 /6, which yields the following bounds: nπ sin 2N π n(1 − λn /n) π = n 1 − nπ (5.17) 2 2 2N nπ π 3 n3 nπ 2N = 2 6 48 N 2 nπ π 3 n3 nπ 2N = (5.19) . > 2 12 96 N 2 In the range 0 < y < 1, f (y) can be bounded by its argument 1/4y < f (y) < 1/2y. Thus, for n < N 2/3 n(1 − λn /n) 1 n(1 − λn /n) (5.20) π 4.
(5.24)
|bn,N (t)| ≤
Similarly, for the tail end of the series as N becomes large, we have, using (5.23), (5.25)
∞ n=N γ +1,odd
∼
c N γ(1−β) . β−1
From these approximations, the conjecture in section 4.3 follows for the upper bound on the error term, TN (xi , t), due to the second-order finite difference scheme.
INITIAL-BOUNDARY VALUE PROBLEMS
1751
Acknowledgments. The authors gratefully thank Dr. Keith Lindsay at the National Center for Atmospheric Research for his discussions, help, and contributions, particularly with the convergence proofs in the appendices and the factorizations (5.8) and (5.9). REFERENCES [1] J.P. Boyd, Chebyshev and Fourier Spectral Methods, Springer-Verlag, New York, 1989. [2] J.P. Boyd and N. Flyer, Compatibility conditions for time-dependent partial differential equations and the rate of convergence of Chebyshev and Fourier spectral methods, Comput. Methods Appl. Mech. Engrg., 175 (1999), pp. 281–309. [3] A. Friedman, Partial Differential Equations of Parabolic Type, Prentice-Hall, Englewood Cliffs, NJ, 1964. [4] H.O. Kreiss, Initial boundary value problems for hyperbolic systems, Comm. Pure Appl. Math, 23 (1970), pp. 277–298. [5] H.O. Kreiss, V. Thomee, and O. Widlund, Smoothing of the initial data and rates of convergence for parabolic difference equations, Comm. Pure Appl. Math., 23 (1970), pp. 241–259. [6] H.O. Kreiss and J. Lorenz, Initial-Boundary Value Problems and the Navier-Stokes Equations, Academic Press, Boston, 1989. [7] R. Hersh, Mixed problems in several variables, J. Math. Mech., 12 (1963), pp. 317–334. [8] D. Gottlieb and S.A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, SIAM, Philadelphia, 1977. [9] B. Gustafsson, H.O. Kriess, and J. Oliger, Time-Dependent Problems and Difference Methods, Wiley-Interscience, New York, 1995. [10] O. Ladyzenskaja, On the convergence of Fourier series defining a solution of a mixed problem for hyperbolic equations, Dokl. Akad. Nauk SSSR, 85 (1952), pp. 481–484. [11] O. Ladyzenskaja, On the solvability of the fundamental boundary problems for equations of parabolic and hyperbolic type, Dokl. Akad. Nauk SSSR, 87 (1954), pp. 395–398. [12] O. Ladyzenskaja, V. Solonnikov, and N. Ural’ceva, Linear and Quasi-Linear Equations of Parabolic Type, Transl. Math. Monogr. 23, AMS, Providence, RI, 1968. [13] A. Majda, J. McDonough, and S. Osher, The Fourier method for nonsmooth initial data, Math. Comp., 32 (1978), pp. 1041–1081. [14] A. Majda and S. Osher, Propagation of error into regions of smoothness for accurate difference approximations to hyperbolic equations, Comm. Pure Appl. Math., 30 (1977), pp. 671–705. [15] S. Osher, Smoothing for spectral methods, in Spectral Methods for Partial Differential Equations, SIAM, Philadelphia, 1984, pp. 209–216. [16] J.B. Rauch and F.J. Massey, Differentiability of solutions to hyperbolic initial-boundary value problems, Trans. Amer. Math. Soc., 189 (1974), pp. 303–318. [17] R. Sakamoto, Hyperbolic Boundary Value Problems, Cambridge University Press, Cambridge, UK, 1982. [18] S. Smale, Smooth solutions of the heat and diffusion equations, Comment. Math. Helv., 55 (1980), pp. 1–12. [19] V.I. Smirnov, A Course of Higher Mathematics. Integral Equations and Partial Differential Equations, Vol. 4, Pergamon Press, London, 1964. [20] R. Temam, Behaviour at time t = 0 of the solutions of semi-linear evolution equations, J. Differential Equations, 43 (1982), pp. 73–92. [21] V. Thomee, Numerical Solution of Partial Differential Equations II, Academic Press, 1970, pp. 585–622. [22] V. Thomee and L. Wahlbin, Convergence rates of parabolic difference schemes for non-smooth data, Math. Comp., 28 (1974), pp. 1–13.
