MODIFIED NODAL CUBIC SPLINE COLLOCATION FOR POISSON'S ...
MODIFIED NODAL CUBIC SPLINE COLLOCATION FOR POISSON’S EQUATION ABEER ALI ABUSHAMA
† AND
BERNARD BIALECKI‡
Abstract. We present a new modified nodal cubic spline collocation scheme for solving the Dirichlet problem for Poisson’s equation on the unit square. We prove existence and uniqueness of a solution of the scheme and show how the solution can be computed on an (N + 1) × (N + 1) uniform partition of the square with cost O(N 2 logN ) using a direct fast Fourier transform method. Using two comparison functions, we derive an optimal fourth order error bound in the continuous maximum norm. We compare our scheme with other modified nodal cubic spline collocation schemes, in particular, the one proposed by Houstis et al. in [8]. We believe that our paper gives the first correct convergence analysis of a modified nodal cubic spline collocation for solving partial differential equations.
Key words. nodal collocation, cubic splines, convergence analysis, interpolants AMS subject classifications. 65N35, 65N12, 65N15, 65N22
1. Introduction. De Boor [7] proved that classical nodal cubic spline collocation for solving two-point boundary value problems is only second–order accurate and no better. For two-point boundary value problems, Archer [2] and independently Daniel and Swartz [6] developed a modified nodal cubic spline collocation (MNCSC) scheme which is fourth order accurate. The approximate solution in this scheme satisfies higher-order perturbations of the ordinary differential equation at the partition nodes. Based on the method of [2] and [6], Houstis et al. [8] derived a fourth order MNCSC scheme for solving elliptic boundary value problems on rectangles. For the Helmholtz equation, a direct fast Fourier transform (FFT) algorithm for solving this scheme was proposed recently in [3]. In this paper, we consider the Dirichlet boundary value problem for Poisson’s equation (1.1)
∆u = f in Ω,
u = 0 on ∂Ω,
where ∆ denotes the Laplacian, Ω = (0, 1) × (0, 1), and ∂Ω is the boundary of Ω. Let +1 ρx = {xi }N i=0 be a uniform partition of [0, 1] in the x-direction such that xi = ih, i = 0, . . . , N + 1, where h = 1/(N + 1). For the sake of simplicity, we assume that a +1 uniform partition ρy = {yj }N i=0 of [0, 1] in the y-direction is such that yj = xj . Let S3 be the space of cubic splines defined by S3 = {v ∈ C 2 [0, 1] : v|[xi−1 ,xi ] ∈ P3 , i = 1, . . . , N + 1}, where P3 denotes the set of all polynomials of degree ≤ 3, and let S D = {v ∈ S3 : v(0) = v(1) = 0}. Our MNCSC scheme for solving (1.1) is formulated as follows: Find uh ∈ S D ⊗ S D such that (1.2)
∆uh (xi , yj ) −
h2 2 2 h2 Dx Dy uh (xi , yj ) = f (xi , yj ) − ∆f (xi , yj ), 6 12 i, j = 0, . . . , N + 1.
† Department of Mathematical and Computer Sciences, Colorado School of Mines, Golden, Colorado 80401-1887, U.S.A. ([email protected]) ‡ Department of Mathematical and Computer Sciences, Colorado School of Mines, Golden, Colorado 80401-1887, U.S.A. ([email protected])
1
The scheme (1.2) is motivated by the fourth order finite difference method for (1.1), see, for example, equation (7) in section 4.5 of [9]. Using uh = u = 0 on ∂Ω and (1.1), we see that (1.2) is equivalent to: (1.3)
2Dx2 Dy2 uh (xi , yj ) = ∆f (xi , yj ),
i, j = 0, N + 1,
(1.4)
Dx2 uh (xi , yj ) −
h2 2 2 h2 Dx Dy uh (xi , yj ) = f (xi , yj ) − ∆f (xi , yj ), 6 12 i = 0, N + 1, j = 1, . . . , N,
(1.5)
Dy2 uh (xi , yj ) −
h2 2 2 h2 Dx Dy uh (xi , yj ) = f (xi , yj ) − ∆f (xi , yj ), 6 12 i = 1, . . . , N, j = 0, N + 1,
(1.6) ∆uh (xi , yj ) −
h2 h2 2 2 Dx Dy uh (xi , yj ) = f (xi , yj ) − ∆f (xi , yj ), i, j = 1, . . . , N. 6 12
The scheme (4.2)–(4.4) of [8] for (1.1) is: Find uh ∈ S D ⊗ S D satisfying (1.3) and (1.7)
Dx2 uh (xi , yj ) = f (xi , yj ),
i = 0, N + 1, j = 1, . . . , N,
(1.8)
Dy2 uh (xi , yj ) = f (xi , yj ),
i = 1, . . . , N, j = 0, N + 1,
(1.9)
(Lx + Ly )uh (xi , yj ) = f (xi , yj ), i, j = 1, . . . , N,
where, for i, j = 1, . . . , N , (1.10)
Lx v(xi , yj ) =
¤ 1 £ 2 Dx v(xi−1 , yj ) + 10Dx2 v(xi , yj ) + Dx2 v(xi+1 , yj ) , 12
Ly v(xi , yj ) =
¤ 1 £ 2 Dy v(xi , yj−1 ) + 10Dy2 v(xi , yj ) + Dy2 v(xi , yj+1 ) . 12
Our scheme and that of [8] are identical at the corners of Ω. However, they are different at the remaining partition nodes. While (1.4)–(1.6) involve perturbations of both the left- and right-hand sides, (1.9) involves a perturbation of the left-hand side only. Numerical results show that our scheme exhibits superconvergence phenomena while that of [8] does not. An outline of this paper is as follows. We give preliminaries in section 2. The matrix-vector form of our scheme, an existence and uniqueness proof of its solution, and a direct FFT algorithm for solving the scheme are presented in section 3. In section 4, using two comparison functions, we derive a fourth order error bound in the continuous maximum norm. In section 5, we give convergence analysis of the scheme in [4] that consists of (1.3)–(1.5) and (1.9). We also explain why convergence analysis of the scheme (1.3) and (1.7)–(1.9), given in [8], is incorrect. This is why, we believe, our paper gives the first correct convergence analysis of MNCSC for solving partial differential equations. Section 6 includes numerical results obtained using our scheme. 2
+1 2. Preliminaries. We extend the uniform partition ρx = {xi }N i=0 outside of [0, 1] using xi = ih, i = −3, −2, −1, N + 2, N + 3, N + 4, and introduce Ii = [xi−1 , xi ], N +2 i = −2, . . . , N + 4. Let {Bm }m=−1 be the basis for S3 defined by g1 [(x − xm−2 )/h], x ∈ Im−1 , g2 [(x − xm−1 )/h], x ∈ Im , g2 [(xm+1 − x)/h], x ∈ Im+1 , Bm (x) = (2.1) g1 [(xm+2 − x)/h], x ∈ Im+2 , 0, otherwise,
where (2.2)
g1 (x) = x3 ,
g2 (x) = 1 + 3x + 3x2 − 3x3 .
The basis functions are such that, for m = 0, . . . , N + 1, (2.3)
Bm−1 (xm ) = 1, Bm (xm ) = 4, Bm+1 (xm ) = 1, 00 00 00 (xm ) = −12/h2 , Bm+1 (xm ) = 6/h2 . Bm−1 (xm ) = 6/h2 , Bm
D N +1 Let {Bm }m=0 be the basis for S D defined by
(2.4)
B1D = B1 − B−1 , B0D = B0 − 4B−1 , D Bm = Bm , m = 2, . . . , N − 1, D D BN = BN − BN +2 , BN +1 = BN +1 − 4BN +2 .
It follows from (2.3) that (2.5)
B0D (x1 ) = 1, B1D (x1 ) = 4, D D BN (xN −1 ) = 1, BN (xN ) = 4,
£ D ¤00 (x ) = −36/h2 , B £ 0D ¤00 0 (x ) = 6/h2 , B (2.6)£ 0D ¤00 1 B (xN −1 ) = 6/h2 , ¤ £ N D 00 BN (xN +1 ) = 0,
B1D (x2 ) = 1, D BN +1 (xN ) = 1,
£ D ¤00 (x ) = 0, B £ 1D ¤00 0 (x ) = −12/h2 , B £ 1D ¤00 1 B (x ) = −12/h2 , ¤00 N £ N D BN +1 (xN +1 ) = −36/h2 .
£ D ¤00 (x ) = 6/h2 , B £ 1D ¤00 2 BN +1 (xN ) = 6/h2 ,
It also follows from (2.5), (2.6), (2.4), and (2.3) that, for i = 1, . . . , N, ½ h2 D 00 6, m = i, D (2.7) m = 0, . . . , N + 1. ] (xi ) = Bm (xi ) − [Bm 0, m 6= i, 6 Throughout the paper, C denotes a generic positive constant that is independent of u and h. D D N +1 Lemma 2.1. {Bm }m=0 of (2.4) satisfy max |Bm (x)| ≤ C, m = 0, . . . , N + 1. x∈[0,1]
Proof. For each fixed m = −1, . . . , N + 2, using Ii = [xi−1 , xi ] and xi = ih, we have (2.8)
0 ≤ (x − xm−2 )/h ≤ 1, x ∈ Im−1 , 0 ≤ (xm+1 − x)/h ≤ 1, x ∈ Im+1 ,
0 ≤ (x − xm−1 )/h ≤ 1, x ∈ Im , 0 ≤ (xm+2 − x)/h ≤ 1, x ∈ Im+2 .
Equations (2.2) and (2.8) give (2.9)
|g1 [(x − xm−2 )/h]| ≤ 1, x ∈ Im−1 , |g2 [(xm+1 − x)/h]| ≤ 7, x ∈ Im+1 ,
|g2 [(x − xm−1 )/h]| ≤ 7, x ∈ Im , |g1 [(xm+2 − x)/h]| ≤ 1, x ∈ Im+2 . 3
Using (2.1) and (2.9), we see that max |Bm (x)| ≤ C, m = −1, . . . , N + 2. Hence x∈[0,1]
D is a linear the required inequality follows from (2.4) which implies that each Bm N +2 combination of at most two of the functions {Bn }n=−1 . 2 N +1
D For {Bm }m=0 of (2.4), we introduce N × N matrices A and B defined by
(2.10)
N
N
D 00 ] (xi ), B = (bj,n )j,n=1 , bj,n = BnD (yj ). A = (ai,m )i,m=1 , ai,m = [Bm
It follows from (2.4), (2.3), (2.5), and (2.6) that A = 6h−2 T,
(2.11)
B = T + 6I,
where I is the identity matrix and the N × N matrix T −2 1 1 −2 1 . . .. .. ... (2.12) T = 1 −2 1 1 −2
is given by .
Lemma 2.2. If B[u1 , . . . , uN ]T = [v1 , . . . , vN ]T , where B is defined in (2.10), then max |ui | ≤ C max |vi |. 1≤i≤N 1≤i≤N X Proof. It follows from (2.10), (2.11), and (2.12) that |bi,i | − |bi,j | ≥ 2, i = 1, . . . , N . i6=j
Hence the required result follows, for example, from the discussion on page 21 in [1]. 2 In what follows, [u1,1 , . . . , uN,N ]T is the short notation for [u1,1 , . . . , u1,N , u2,1 , . . . , u2,N , . . . , uN,1 , . . . , uN,N ]T . Lemma 2.3. If u = [u1,1 , . . . , uN,N ]T and v = [v1,1 , . . . , vN,N ]T are such that (B ⊗ B)u = v, where B is defined in (2.10), then max |ui,j | ≤ C max |vi,j |. 1≤i,j≤N
1≤i,j≤N
Proof. Since B ⊗ B = (B ⊗ I)(I ⊗ B), we have (2.13)
v = (B ⊗ I)w,
w = (I ⊗ B)u.
Using ( 2.13) and Lemma 2.2, we obtain max |ui,j | ≤ C max |wi,j |,
1≤i,j≤N
max |wi,j | ≤ C max |vi,j |,
1≤i,j≤N
1≤i,j≤N
1≤i,j≤N
which imply the required inequality. 2 It is well known (see Theorem 4.5.2 of [10]) that for T of (2.12), we have (2.14)
QT Q = Λ,
QQ = I,
where the N × N matrices Λ and Q are given by (2.15)
(2.16)
N
Λ = diag(λi )i=1 ,
N
Q = (qi,j )i,j=1 ,
λi = −4sin2
qi,j =
µ
4
2 N +1
iπ , 2(N + 1)
¶1/2
sin
ijπ . N +1
Lemma 2.4. If v = [v1,1 , . . . , vN,N ]T and w = [w1,1 , . . . , wN,N ]T are such that ¶¸ · µ T h2 T T T (2.17) v = w, ⊗ I + I ⊗ + ⊗ h2 h2 6 h2 h2 where T is the matrix defined in (2.12), then
2 max vi,j ≤ Ch2
1≤i,j≤N
N N X X
2 wi,j .
i=1 j=1
Proof. The matrix in (2.17) arises in the fourth order finite difference method for (1.1). Hence the desired result follows, for example, from the last unnumbered equation on page 296 in [9]. 2 Finally, we observe that the matrix-vector form of (2.18)
φi,j =
N X
(1)
ci,m
m=1
N X
(2)
cj,n ψm,n ,
i, j = 1, . . . , N,
n=1
is φ = (C1 ⊗ C2 )ψ,
(2.19) ´N ³ (1) where C1 = ci,m
i,m=1
´N ³ (2) , C2 = cj,n
j,n=1
, and
φ = [φ1,1 , . . . , φN,N ]T ,
ψ = [ψ1,1 , . . . , , ψN,N ]T .
3. Matrix-Vector Form of Scheme. Since dim(S D ⊗ S D ) = (N + 2)2 , the scheme (1.3)–(1.6) involves (N + 2)2 equations in (N + 2)2 unknowns. Using the basis D N +1 }m=0 of (2.4) for the space S D , we have {Bm (3.1)
uh (x, y) =
+1 N +1 N X X
D (x)BnD (y). um,n Bm
m=0 n=0
Substituting (3.1) into (1.3), we obtain (3.2)
2
N +1 N +1 X X
D 00 ] (xi )[BnD ]00 (yj ) = ∆f (xi , yj ), um,n [Bm
i, j = 0, N + 1.
m=0 n=0
Using (2.6), we conclude that (3.2) gives ui,j =
(3.3)
h4 ∆f (xi , yj ), 2592
i, j = 0, N + 1.
Substituting (3.1) into (1.4), we obtain
(3.4)
+1 N +1 N X X
µ ¶ h2 D 00 um,n [Bm ] (xi ) BnD (yj ) − [BnD ]00 (yj ) 6 m=0 n=0 2 h i = 0, N + 1, j = 1, . . . , N. = f (xi , yj ) − ∆f (xi , yj ), 12
Using (2.6) and (2.7), we see that (3.4) gives (3.5)
ui,j = −
h2 h4 f (xi , yj ) + ∆f (xi , yj ), i = 0, N + 1, j = 1, . . . , N. 216 2592 5
Using (3.5) and symmetry with respect to x and y, we conclude that (1.5) gives (3.6)
ui,j = −
h2 h4 f (xi , yj ) + ∆f (xi , yj ), i = 1, . . . , N, j = 0, N + 1. 216 2592
Substituting (3.1) into (1.6), we obtain +1 N +1 N X X
· ¶ ¸ µ h2 D 00 D D 00 (xi ) − [Bm ] (xi )BnD (yj ) + Bm ] (xi ) [BnD ]00 (yj ) um,n [Bm 6 (3.7) m=0 n=0 h2 = f (xi , yj ) − ∆f (xi , yj ), i, j = 1, . . . , N. 12 +1 N Moving the terms involving {um,n }N n=0 , m = 0, N + 1, {um,n }m=1 , n = 0, N + 1, to the right-hand side of (3.7), we get
(3.8)
N N X X
um,n
m=1 n=1
µ
D 00 ] (xi )BnD (yj ) [Bm
· ¶ ¸ h2 D 00 D D 00 + Bm (xi ) − [Bm ] (xi ) [Bn ] (yj ) 6
= pi,j , i, j = 1, . . . , N, where pi,j = f (xi , yj ) − N +1 X
h2 ∆f (xi , yj ) 12 µ
· ¶ ¸ h2 D 00 D D 00 + Bm (xi ) − [Bm ] (xi ) [Bn ] (yj ) um,n − 6 m=0,N +1 n=0 · µ ¶ ¸ N X X h2 D 00 D D 00 (xi ) − [Bm ] (xi )BnD (yj ) + Bm um,n [Bm − ] (xi ) [BnD ]00 (yj ) . 6 m=1 X
D 00 ] (xi )BnD (yj ) [Bm
n=0,N +1
Using (2.18)–(2.19), we write (3.8) as ¶ ¸ · µ h2 (3.9) A ⊗ B + B − A ⊗ A u = p, 6 where u = [u1,1 , . . . , uN,N ]T , p = [p1,1 , . . . , pN,N ]T , and A, B are defined in (2.10). Using (2.11), we see that µ ¶ h2 6 A ⊗ B + B − A ⊗ A = 2 [6T ⊗ I + (6I + T ) ⊗ T ] , 6 h and hence the system (3.9) simplifies to (3.10)
6h−2 [6T ⊗ I + (6I + T ) ⊗ T ] u = p.
We are now ready to prove existence and uniqueness of uh in S D ⊗ S D that satisfies (1.3)–(1.6). Theorem 3.1. There exists unique uh in S D ⊗ S D satisfying (1.3)–(1.6). Proof. Since the number of equations in (1.3)–(1.6) is equal to the number of unknowns, we assume that the right-hand side in (1.3)–(1.6) is zero, and show that uh = 0 is the only solution of the resulting scheme. Using (3.1), (3.3), (3.5), and (3.6), we have (3.11) um,n = 0, m = 0, N + 1, n = 0, . . . , N + 1, m = 1, . . . , N, n = 0, N + 1. 6
Clearly (3.12) 6h
−2
T T h2 [6T ⊗ I + (6I + T ) ⊗ T ] = 36 2 ⊗ I + I ⊗ 2 + h h 6 ·
µ
T T ⊗ 2 h2 h
¶¸
.
Hence it follows from (3.10) with p replaced by 0, (3.12), and Lemma 2.4 that (3.13)
um,n = 0,
m, n = 1, . . . , N.
Equations (3.1), (3.11), and (3.13) give uh = 0. 2 Using Q of (2.16), we see that (3.10) is equivalent to (3.14) 6h−2 (Q ⊗ I) [6T ⊗ I + (6I + T ) ⊗ T ] (Q ⊗ I)(Q−1 ⊗ I)u = (Q ⊗ I)p. Introducing u0 = (Q−1 ⊗ I)u and p0 = (Q ⊗ I)p, and using (3.14) and (2.14), we obtain 6h−2 [6Λ ⊗ I + (6I + Λ) ⊗ T ]u0 = p0 ,
(3.15)
where Λ is defined in (2.15). The system (3.15) reduces to the N independent systems (3.16)
6h−2 [6λi I + (6 + λi )T ]u0i = p0i ,
i = 1, . . . , N,
where u0i = [u0i,1 , . . . , u0i,N ]T , p0i = [p0i,1 , . . . , p0i,N ]T , i = 1, . . . , N . We have the following algorithm for solving (3.10): Step 1. Compute p0 = (Q ⊗ I)p. Step 2. Solve the N systems in (3.16). Step 3. Compute u = (Q ⊗ I)u0 . Since the entries of Q in (2.16) are given in terms of sines, steps 1 and 3 are performed each using FFTs at a cost O(N 2 log N ). In step 2, the systems are tridiagonal, so this step is performed at a cost O(N 2 ). Thus the total cost of the algorithm is O(N 2 log N ). 4. Convergence Analysis. In what follows, C(u) denotes a generic positive constant that is independent of h, but depends on u. Our goal is to show that if u in C 6 (Ω) and uh in S D ⊗ S D are the solutions of (1.1) and (1.3)–(1.6), respectively, then ku − uh kC(Ω) ≤ C(u)h4 ,
(4.1)
where kgkC(Ω) = max |g(x)| for g in C(Ω). x∈Ω
To prove (4.1), for u in C 4 (Ω), we introduce two comparison functions, the spline interpolants S and Z in S D ⊗ S D of u defined respectively by (4.2)
(4.3)
Dx2 Dy2 S(xi , yj ) = Dx2 Dy2 u(xi , yj ),
i, j = 0, N + 1,
h2 2 2 h2 Dx Dy S(xi , yj ) = Dx2 u(xi , yj ) − Dx4 u(xi , yj ) 6 12 h2 2 2 − Dx Dy u(xi , yj ), i = 0, N + 1, j = 1, . . . , N, 6
Dx2 S(xi , yj ) −
7
(4.4)
(4.5)
h2 h2 Dy2 S(xi , yj ) − Dx2 Dy2 S(xi , yj ) = Dy2 u(xi , yj ) − Dy4 u(xi , yj ) 6 12 h2 i = 1, . . . , N, j = 0, N + 1, − Dx2 Dy2 u(xi , yj ), 6
S(xi , yj ) = u(xi , yj ),
i, j = 1, . . . , N,
and (4.6)
Dx2 Dy2 Z(xi , yj ) = Dx2 Dy2 u(xi , yj ),
i, j = 0, N + 1,
(4.7)
Dx2 Z(xi , yj ) = Dx2 u(xi , yj ),
i = 0, N + 1, j = 1, . . . , N,
(4.8)
Dy2 Z(xi , yj ) = Dy2 u(xi , yj ),
i = 1, . . . , N, j = 0, N + 1,
(4.9)
Z(xi , yj ) = u(xi , yj ),
i, j = 1, . . . , N.
It follows from (1.1) that (4.10)
f = Dx2 u + Dy2 u,
∆f = Dx4 u + Dy4 u + 2Dx2 Dy2 u.
Hence, using u = 0 on ∂Ω, we see that (1.3)–(1.5) reduce, respectively, to (4.11)
Dx2 Dy2 uh (xi , yj ) = Dx2 Dy2 u(xi , yj ),
i, j = 0, N + 1,
(4.12)
h2 h2 Dx2 uh (xi , yj ) − Dx2 Dy2 uh (xi , yj ) = Dx2 u(xi , yj ) − Dx4 u(xi , yj ) 6 12 h2 2 2 − Dx Dy u(xi , yj ), i = 0, N + 1, j = 1, . . . , N, 6
(4.13)
h2 h2 Dy2 uh (xi , yj ) − Dx2 Dy2 uh (xi , yj ) = Dy2 u(xi , yj ) − Dy4 u(xi , yj ) 6 12 h2 2 2 i = 1, . . . , N, j = 0, N + 1. − Dx Dy u(xi , yj ), 6
Comparing (4.11)–(4.13) and (4.2)–(4.4), we see that uh and S are defined in the same way for i = 0, N + 1, j = 0, . . . , N + 1, and i = 1, . . . , N , j = 0, N + 1. On the other hand, (4.6)–(4.8) are a simplified, tensor product, version of (4.2)–(4.4). The triangle inequality gives (4.14)
ku − uh kC(Ω) ≤ ku − ZkC(Ω) + kZ − SkC(Ω) + kS − uh kC(Ω) .
In what follows, we bound the three terms on the right-hand side of (4.14). 8
4.1. Bounding ku − ZkC(Ω) . We need the following results. Lemma 4.1. Let the interpolant Ix v in S3 of v in C 2 [0, 1] be defined by (4.15) (Ix v)00 (xi ) = v 00 (xi ), i = 0, N + 1,
Ix v(xi ) = v(xi ), i = 0, . . . , N + 1.
Then max |v(x) − Ix v(x)| ≤ C max |v 00 (x)|h2 .
(4.16)
x∈[0,1]
x∈[0,1]
If v ∈ C 4 [0, 1], then max |v(x) − Ix v(x)| ≤ C max |v (4) (x)|h4 .
(4.17)
x∈[0,1]
x∈[0,1]
Proof. First we prove (4.16). Using the discussion on page 404 in [5], we have (4.18) Ix v(x) = v(xi ) + Bi (x − xi ) + Ci (x − xi )2 + Di (x − xi )3 , x ∈ [xi , xi+1 ], for i = 0, . . . , N , where h 1 ri 1 h (4.19) Bi = − ri+1 − ri + [v(xi+1 ) − v(xi )] , Ci = , Di = (ri+1 − ri ), 6 3 h 2 6h and ri = (Ix v)00 (xi ). Equations (4.18) and (4.19) give (4.20)
h h ri Ix v(x) − v(x) = Ai (x)− ri+1 (x − xi ) − ri (x − xi ) + (x − xi )2 6 3 2 1 3 (ri+1 − ri )(x − xi ) , x ∈ [xi , xi+1 ], + 6h
where Ai (x) = v(xi ) − v(x) +
(4.21)
v(xi+1 ) − v(xi ) (x − xi ), x ∈ [xi , xi+1 ]. h
Using (4.20) and the triangle inequality, we obtain, for x ∈ [xi , xi+1 ], ¶ µ 4 1 (4.22) |Ix v(x) − v(x)| ≤ |Ai (x)| + h2 |ri | + |ri+1 | ≤ |Ai (x)| + h2 max |ri |. 3 3 0≤i≤N +1 We introduce
E=
N +1 (ei,j )i,j=0
1 1 =
4 .. .
1 .. . 1
..
. 4
, 1 1
r = [r0 , . . . , rN +1 ]T , p = [p0 , . . . , pN +1 ]T ,
where (4.23)
pi =
½
v 00 (xi ), i = 0, N + 1, h−2 [v(xi−1 ) − 2v(xi ) + v(xi+1 )], i = 1, . . . , N.
It follows X from the discussion on pages 400 and 401 in [5] that Er = p. Since |ei,i | − |ei,j | ≥ 1, i, j = 0, . . . , N + 1, the discussion on page 21 in [1] implies that i6=j
(4.24)
max
0≤i≤N +1
|ri | ≤ C 9
max
0≤i≤N +1
|pi |.
Using Taylor’s theorem, we obtain |v(xi−1 ) − 2v(xi ) + v(xi+1 )| ≤ Ch2 max |v 00 (x)|,
(4.25)
x∈[0,1]
i = 1, . . . , N.
It follows from (4.24), (4.23) and (4.25), that (4.26)
max
0≤i≤N +1
|ri | ≤ C max |v 00 (x)|. x∈[0,1]
Using Taylors’ theorem to expand v(x), x ∈ [xi , xi+1 ], around xi , we have (x − xi )2 00 v (ξi,x ), 2
v(x) = v(xi ) + (x − xi )v 0 (xi ) +
(4.27)
xi ≤ ξi,x ≤ x.
Using (4.21), (4.27), and the triangle inequality, we obtain, for x ∈ [xi , xi+1 ], ¯ ¯ ¯ (x − xi )2 00 ¯ h 00 ¯ (4.28) |Ai (x)| = ¯ v (ξi,x ) − (x − xi )v (ξi,xi+1 )¯¯ ≤ h2 max |v 00 (x)|. 2 2 x∈[0,1]
Inequality (4.16) follows from (4.22), (4.26), and (4.28). A proof of (4.17) is given in the proof of Theorem 2.3.4 in [1]. 2 Lemma 4.2. If u ∈ C 4 (Ω), Z in S D ⊗ S D is defined by (4.6)–(4.9), and Ix u and Iy u are defined in (4.15), then for (x, y) in Ω, we have (4.29)
Z(x, y) = Ix (Iy u)(x, y),
Dx2 (Iy u)(x, y) = Iy (Dx2 u)(x, y).
+3 Proof. Let {Ci }N i=0 be the basis for S3 such that
Ci (xj ) = δij , Ci00 (xj ) = 0,
i, j = 0, . . . , N + 1, i = 0, . . . , N + 1, j = 0, N + 1,
(4.30) CN +2 (xj ) = CN +3 (xj ) = 0,
j = 0, . . . , N + 1,
00 00 CN +2 (x0 ) = CN +3 (xN +1 ) = 1,
00 00 CN +2 (xN +1 ) = CN +3 (x0 ) = 0,
where δij is the Kronecker delta. Using (4.15) and (4.30), we have for (x, y) ∈ Ω, Ix (Iyu)(x, y) N +1 X = Ix u(x, yj )Cj (y) + Dy2 u(x, y0 )CN +2 (y) + Dy2 u(x, yN +1 )CN +3 (y) j=0
=
N +1 X
N +1 X
u(xi , yj )Cj (y) + Dy2 u(xi , y0 )CN +2 (y) + Dy2 u(xi , yN +1 )CN +3 (y) Ci (x)
i=0 j=0 N +1 X + Dx2 u(x0 , yj )Cj (y) + Dx2 Dy2 u(x0 , y0 )CN +2 (y) j=0
N +1 X ¤ Dx2 u(xN +1 , yj )Cj (y) +Dx2 Dy2 u(x0 , yN +1 )CN +3 (y) CN +2 (x) +
j=0 ¤ +Dx2 Dy2 u(xN +1 , y0 )CN +2 (y) + Dx2 Dy2 u(xN +1 , yN +1 )CN +3 (y) CN +3 (x).
10
Since u = 0 on ∂Ω, all terms involving C0 (x), CN +1 (x), C0 (y), CN +1 (y) drop out which implies that Ix (Iy u) ∈ S D ⊗ S D . Using (4.30), we verify that Ix (Iy u) satisfies (4.6)–(4.9), that is, (4.6)–(4.9) hold with Ix (Iy u) in place of Z. Hence, the uniqueness of the interpolant Z implies the first equation in (4.29). To prove the second equation in (4.29), we use (4.15) and (4.30) to see that for (x, y) ∈ Ω, Iy (Dx2 u)(x, y) N +1 X Dx2 u(x, yj )Cj (y) + Dx2 Dy2 u(x, y0 )CN +2 (y) + Dx2 Dy2 u(x, yN +1 )CN +3 (y) = j=0 N +1 X u(x, yj )Cj (y) + Dy2 u(x, y0 )CN +2 (y) + Dy2 u(x, yN +1 )CN +3 (y) = Dx2 j=0
= Dx2 (Iy u)(x, y).
2
Theorem 4.1. If u ∈ C 4 (Ω) and Z in S D ⊗ S D is defined by (4.6)–(4.9), then ku − ZkC(Ω) ≤ C(u)h4 . Proof. Using (4.29) and the triangle inequality, we have (4.31)
ku − ZkC(Ω)
≤ ku − Ix ukC(Ω) + kIx (u − Iy u) − (u − Iy u)kC(Ω) + ku − Iy ukC(Ω) .
For any fixed y in [0, 1], Ix u(·, y) is the cubic spline interpolant of u(·, y). Using this, symmetry with respect to x and y, and (4.17), we have (4.32)
ku − Ix ukC(Ω) ≤ C(u)h4 ,
ku − Iy ukC(Ω) ≤ C(u)h4 .
For any fixed y in [0, 1], Ix (u−Iy u)(·, y) is the cubic spline interpolant of (u−Iy u)(·, y). Hence it follows from (4.16) that (4.33)
kIx (u − Iy u) − (u − Iy u)kC(Ω) ≤ CkDx2 (u − Iy u)kC(Ω) h2 .
Using (4.29) and(4.16), we obtain (4.34)
kDx2 (u − Iy u)kC(Ω) = kDx2 u − Iy (Dx2 u)kC(Ω) ≤ CkDx2 Dy2 ukC(Ω) h2 .
Combining (4.33) and (4.34), we have (4.35)
kIx (u − Iy u) − (u − Iy u)kC(Ω) ≤ C(u)h4 .
The desired inequality now follows from (4.31), (4.32), and (4.35).
2
4.2. Bounding kZ − SkC(Ω) . We start by proving the following lemma. Lemma 4.3. If u ∈ C 4 (Ω) and (4.36) S(x, y) =
N +1 N +1 X X
D sm,n Bm (x)BnD (y), Z(x, y) =
m=0 n=0
N +1 N +1 X X
D zm,n Bm (x)BnD (y),
m=0 n=0
are defined by (4.2)–(4.5) and (4.6)–(4.9), respectively, then |sm,n − zm,n | ≤ C(u)h4 , 11
m, n = 0, . . . , N + 1.
Proof. Using (4.2), (4.6), and following the derivation of (3.3) from (1.3), we obtain (4.37)
sm,n = zm,n ,
m, n = 0, N + 1.
Next we prove the required inequality for m = 0, n = 1, . . . , N . Using (4.7), we have (4.38)
Dx2 (S − Z)(x0 , yj ) = Dx2 S(x0 , yj ) − Dx2 u(x0 , yj ),
j = 1, . . . , N.
It follows from (4.36), (4.37), and (2.6) that (4.39) Dx2 (S − Z)(x0 , yj ) = −36h−2
N X
(s0,n − z0,n )BnD (yj ),
j = 1, . . . , N.
n=1
Using (4.36), (2.6), (2.4), (2.3), and (2.5), we obtain, for j = 1, . . . , N , (4.40) Dx2 S(x0 , yj ) = −36h−2
N +1 X
s0,n BnD (yj ) = −36h−2 (s0,j−1 + 4s0,j + s0,j+1 ).
n=0
Substituting (4.39) and (4.40) into (4.38), and multiplying through by −h2 /36, we have (4.41)
N X
(s0,n − z0,n )BnD (yj ) = s0,j−1 + 4s0,j + s0,j+1 +
n=1
h2 2 D u(x0 , yj ) 36 x
for j = 1, . . . , N . Using (4.2), (4.3), and following the derivations of (3.3) from (1.3) and (3.5) from (1.4), we obtain s0,j =
(4.42)
h4 D2 D2 u(x0 , yj ), 1296 x y
j = 0, N + 1,
and (4.43)
s0,j
¸ · h2 2 2 h2 h2 4 2 =− Dx u(x0 , yj ) − Dx u(x0 , yj ) − Dx Dy u(x0 , yj ) 216 12 6
for j = 1, . . . , N . Since u = 0 on ∂Ω, (4.42) is the same as (4.43) with j = 0, N + 1. This observation and (4.43) imply that for j = 1, . . . , N , we have ¸ · h2 h2 h2 (4.44) s0,j±1 = − Dx2 u(x0 , yj±1 ) − Dx4 u(x0 , yj±1 ) − Dx2 Dy2 u(x0 , yj±1 ) . 216 12 6 Using Taylor’s theorem, we obtain (4.45) Dx2 u(x0 , yj±1 ) = Dx2 u(x0 , yj ) ± hDx2 Dy u(x0 , yj ) +
h2 2 2 D D u(x0 , ξj± ), 2 x y
where yj−1 ≤ ξj− ≤ yj , yj ≤ ξj+ ≤ yj+1 . Using (4.44), (4.43), and (4.45), we obtain ¯ ¯ ¯ ¯ h2 2 ¯ (4.46) ¯s0,j−1 + 4s0,j + s0,j+1 + Dx u(x0 , yj )¯¯ ≤ C(u)h4 , j = 1, . . . , N. 36 N
It follows from (4.46) that (4.41) is a system in {s0,n − z0,n }n=1 with the matrix B defined in(2.10) and with each entry on the right-hand side bounded in absolute value by C(u)h4 . Hence, Lemma 2.2 implies (4.47)
max |s0,n − z0,n | ≤ C(u)h4 .
1≤n≤N
12
Using (4.47) and symmetry with respect to x and y, we also have max |sN +1,n − zN +1,n | ≤ C(u)h4 ,
1≤n≤N
(4.48)
max |sm,n − zm,n | ≤ C(u)h4 ,
1≤m≤N
n = 0, N + 1.
Finally we prove the required inequality for m, n = 1, . . . , N. Using (4.5) and (4.9), we have (S − Z)(xi , yj ) = 0, i, j = 1, . . . , N , which, by (4.36) and (4.37), can be written as (4.49)
N X N X
D (sm,n − zm,n )Bm (xi )BnD (yj ) = di,j ,
i, j = 1, . . . N,
m=1 n=1
where
di,j =
X
N X
m=0,N +1 n=1
+
N X
X
m=1 n=0,N +1
D (zm,n − sm,n )Bm (xi )BnD (yj ).
Since for any fixed i, j, each of the above double sums reduces to at most three terms, using the triangle inequality, (4.47), (4.48), and Lemma 2.1, we obtain (4.50)
|di,j | ≤ C(u)h4 ,
i, j = 1, . . . , N. N
It follows from (2.18)–(2.19) that (4.49) is a system in {zm,n − sm,n }m,n=1 with the matrix B ⊗ B, where B is defined in (2.10). Hence, for m, n = 1, . . . , N, the required inequality follows from (4.50) and Lemma 2.3. 2 Theorem 4.2. If u ∈ C 4 (Ω) and S, Z in S D ⊗ S D are defined by (4.2)–(4.5) and (4.6)–(4.9), respectively, then kZ − SkC(Ω) ≤ C(u)h4 . Proof. Since Z − S is continuous on Ω, there is (x∗ , y∗ ) in Ω such that kZ − SkC(Ω) = |(Z − S)(x∗ , y∗ )| . Hence, (4.36) and the triangle inequality imply kZ − SkC(Ω) ≤
N +1 N +1 X X
D |sm,n − zm,n ||Bm (x∗ )||BnD (y∗ )|.
m=0 n=0
Since the above double sum reduces to at most nine terms, the required inequality follows from Lemmas 4.3 and 2.1. 2 4.3. Bounding kS − uh kC(Ω) and ku − uh kC(Ω) . We need the following results. Lemma 4.4. If u ∈ C 6 (Ω) and S in S D ⊗ S D is defined by (4.2)–(4.5), then for i = 0, N + 1, j = 1, . . . , N , ¯ 2 2 ¯ ¯Dx Dy S(xi , yj ) − Dx2 Dy2 u(xi , yj )¯ ≤ C(u)h2 , (4.51) (4.52)
¯ ¯ 2 ¯ ¯ 2 ¯Dx S(xi , yj ) − Dx2 u(xi , yj ) + h Dx4 u(xi , yj )¯ ≤ C(u)h4 . ¯ ¯ 12 13
Proof. We prove (4.51) for i = 0; for i = N + 1, (4.51) follows by symmetry with respect to x. Using (4.36), we obtain Dx2 Dy2 S(x0 , yj ) =
+1 N +1 N X X
m=0 n=0
¤00 £ £ D ¤00 (x0 ) BnD (yj ), sm,n Bm
j = 1, . . . , N,
and hence (2.4), (2.3), and (2.6) imply (4.53)
Dx2 Dy2 S(x0 , yj ) = −216h−4 (s0,j−1 − 2s0,j + s0,j+1 ), j = 1, . . . , N.
Equations (4.53), (4.43), and (4.44) give, for j = 1, . . . , N,
(4.54)
Dx2 Dy2£S(x0 , yj ) − Dx2 Dy2 u(x0 , yj ) = −Dx2 Dy2 u(x0 , yj ) ¤ +h−2 Dx2 u(x0 , yj−1 ) − 2Dx2 u(x0 , yj ) + Dx2 u(x0 , yj+1 ) ¤ 1 £ 4 − Dx u(x0 , yj−1 ) − 2Dx4 u(x0 , yj ) + Dx4 u(x0 , yj+1 ) 12 ¤ 1£ − Dx2 Dy2 u(x0 , yj−1 ) − 2Dx2 Dy2 u(x0 , yj ) + Dx2 Dy2 u(x0 , yj+1 ) . 6
Using Taylor’s theorem, we obtain
(4.55)
Dx2 u(x0 , yj±1 )
h2 2 2 D D u(x0 , yj ) 2 x y 3 4 h h ± Dx2 Dy3 u(x0 , yj ) + Dx2 Dy4 u(x0 , ξj± ), 3! 4!
= Dx2 u(x0 , yj ) ± hDx2 Dy u(x0 , yj ) +
(4.56) Dx4 u(x0 , yj±1 ) = Dx4 u(x0 , yj ) ± hDx4 Dy u(x0 , yj ) +
h2 4 2 D D u(x0 , ηj± ), 2 x y
(4.57) Dx2 Dy2 u(x0 , yj±1 ) = Dx2 Dy2 u(x0 , yj ) ± hDx2 Dy3 u(x0 , yj ) +
h2 2 4 D D u(x0 , κ± j ), 2 x y
+ + + where yj−1 ≤ ξj− , ηj− , κ− j ≤ yj , yj ≤ ξj , ηj , κj ≤ yj+1 . Equations (4.55)–(4.57) give ¯ −2 £ 2 ¯ ¤ ¯h Dx u(x¯0 , yj−1 ) − 2Dx2 u(x0 , yj ) + Dx2 u(x0 , yj+1 ) − Dx2¯Dy2 u(x0 , yj )¯ ≤ C(u)h2 , ¯ x4 u(x0 , yj−1 ) − 2Dx4 u(x0 , yj ) + Dx4 u(x0 , yj+1 )¯ ≤ C(u)h2 , ¯ 2 D ¯ ¯Dx Dy2 u(x0 , yj−1 ) − 2Dx2 Dy2 u(x0 , yj ) + Dx2 Dy2 u(x0 , yj+1 )¯ ≤ C(u)h2 ,
and hence (4.51) for i = 0 follows from (4.54) and the triangle inequality. Using (4.3) and (4.51), we obtain (4.52). 2 Lemma 4.5. If u ∈ C 6 (Ω) and S in S D ⊗ S D is defined by (4.2)–(4.5), then, for i, j = 1, . . . , N , we have ¯ ¯ ¯ ¯ 2 h2 4 4 2 ¯ ¯ (4.58) ¯Dx S(xi , yj ) − Dx u(xi , yj ) + 12 Dx u(xi , yj )¯ ≤ C(u)h , (4.59)
(4.60)
¯ ¯ 2 ¯ ¯ 2 ¯Dy S(xi , yj ) − Dy2 u(xi , yj ) + h Dy4 u(xi , yj )¯ ≤ C(u)h4 , ¯ ¯ 12 |Dx2 Dy2 S(xi , yj ) − Dx2 Dy2 u(xi , yj )| ≤ C(u)h2 . 14
Proof. First we prove (4.58). For i = 0, . . . , N + 1, j = 1, . . . , N , we introduce ¸ · h2 4 2 2 (4.61) di,j = Dx S(xi , yj ) − Dx u(xi , yj ) − Dx u(xi , yj ) . 12 Then (4.62)
di−1,j + 4di,j + di+1,j = φi,j − ψi,j ,
i, j = 1, . . . , N,
where φi,j (4.63)
ψi,j (4.64)
2 2 = Dx2 S(x · i−1 , yj ) + 4Dx2S(xi , yj ) + D ¸x S(xi+1 , yj ) h −6 Dx2 u(xi , yj ) + Dx4 u(xi , yj ) , 12
· ¸ h2 4 h2 4 2 = − Dx u(xi−1 , yj ) + 4 Dx u(xi , yj ) − Dx u(xi , yj ) 12 12 · ¸ h2 4 h2 4 2 2 +Dx u(xi+1 , yj ) − Dx u(xi+1 , yj ) − 6 Dx u(xi , yj ) + Dx u(xi , yj ) 12 12 = Dx2 u(xi−1 , yj ) − 2Dx2 u(xi , yj ) + Dx2 u(xi+1 , yj ) ¤ h2 £ 4 − Dx u(xi−1 , yj ) + 10Dx4 u(xi , yj ) + Dx4 u(xi+1 , yj ) . 12 Dx2 u(xi−1 , yj )
Since S(·, yj ) ∈ S3 , (2.1.7) in [1], (4.5), and S = u = 0 on ∂Ω, imply that (4.65)
Dx2 S(xi−1 , yj ) + 4Dx2 S(xi , yj ) + Dx2 S(xi+1 , yj ) = 6h−2 [u(xi−1 , yj ) − 2u(xi , yj ) + u(xi+1 , yj )] , i, j = 1, . . . , N.
Using Taylor’s theorem, we obtain h3 h2 u(xi±1 , yj ) = u(xi , yj ) ± hDx u(xi , yj ) + Dx2 u(xi , yj )± Dx3 u(xi , yj ) 2 3! h4 4 h5 5 h6 6 + Dx u(xi , yj )± Dx u(xi , yj ) + Dx u(ξi± , yj ), 4! 5! 6! where xi−1 ≤ ξi− ≤ xi , xi ≤ ξi+ ≤ xi+1 , and hence ¯ −2 ¯h [u(xi−1 , yj ) − 2u(xi , yj ) + u(xi+1 , yj )] ¸¯ · ¯ h2 (4.66) − Dx2 u(xi , yj ) + Dx4 u(xi , yj ) ¯¯ ≤ C(u)h4 , i, j = 1, . . . , N. 12 Using (4.63), (4.65), and (4.66), we obtain (4.67)
|φi,j | ≤ C(u)h4 , i, j = 1, . . . , N.
Using Taylor’s theorem, we obtain h2 Dx2 u(xi±1 , yj ) = Dx2 u(xi , yj ) ± hDx3 u(xi , yj ) + Dx4 u(xi , yj ) 2 h4 h3 ± Dx5 u(xi , yj ) + Dx6 u(ξi± , yj ), 3! 4! h2 4 4 5 Dx u(xi±1 , yj ) = Dx u(xi , yj ) ± hDx u(xi , yj ) + Dx6 u(ηi± , yj ), 2 15
where xi−1 ≤ ξi− , ηi− ≤ xi , xi ≤ ξi+ , ηi+ ≤ xi+1 , and hence (4.64) gives |ψi,j | ≤ C(u)h4 ,
(4.68)
i, j = 1, . . . , N.
Using (4.61) and (4.52), we have |di,j | ≤ C(u)h4 ,
(4.69)
i = 0, N + 1,
j = 1, . . . , N.
It follows from (4.67)–(4.69) that moving d0,j and dN +1,j to the right-hand side of N (4.62), we obtain, for each j = 1, . . . , N , a system in {di,j }i=1 with the matrix B of (2.10)–(2.12), and with each entry on the right-hand side bounded in absolute value by C(u)h4 . Hence (4.58) follows from (4.61) and Lemma 2.2, and (4.59) follows from (4.58) by symmetry with respect to x and y. Next we prove (4.60). Since S(x, ·) ∈ S3 for x ∈ [0, 1], (2.1.7) in [1] gives (4.70)
Dy2 S(x, yj−1 ) + 4Dy2 S(x, yj ) + Dy2 S(x, yj+1 ) = 6h−2 [S(x, yj−1 ) − 2S(x, yj ) + S(x, yj+1 )] , j = 1, . . . , N, x ∈ [0, 1].
Differentiating (4.70) twice with respect to x, we obtain, for j = 1, . . . , N , x ∈ [0, 1], (4.71)
2 2 yj ) + Dx2 Dy2 S(x, yj+1 )¤ Dx2 Dy2 S(x, £ 2yj−1 ) + 4Dx Dy S(x, 2 −2 Dx S(x, yj−1 ) − 2Dx S(x, yj ) + Dx2 S(x, yj+1 ) . = 6h
Using (4.71) with x = xi−1 , xi , xi+1 , we obtain, for i, j = 1, . . . , N , (4.72)
(4.73)
(4.74)
Dx2 Dy2 S(x , yj−1 ) + 4Dx2 Dy2 S(xi−1 , yj ) + Dx2 Dy2 S(xi−1 , yj+1 )¤ £ i−1 2 −2 Dx S(xi−1 , yj−1 ) − 2Dx2 S(xi−1 , yj ) + Dx2 S(xi−1 , yj+1 ) , = 6h 2 2 2 2 Dx2 Dy2 S(x i , yj ) + Dx Dy S(xi , yj+1 )¤ £ i ,2yj−1 ) + 4Dx Dy S(x −2 2 = 6h Dx S(xi , yj−1 ) − 2Dx S(xi , yj ) + Dx2 S(xi , yj+1 ) ,
, yj−1 ) + 4Dx2 Dy2 S(xi+1 , yj ) + Dx2 Dy2 S(xi+1 , yj+1 )¤ Dx2 Dy2 S(x £ i+1 2 −2 Dx S(xi+1 , yj−1 ) − 2Dx2 S(xi+1 , yj ) + Dx2 S(xi+1 , yj+1 ) . = 6h
Adding (4.72), (4.74) and (4.73) multiplied through by 4, and using (4.65) and S = u = 0 on ∂Ω, we obtain
(4.75)
Dx2 Dy2 S(xi−1 , yj−1 ) + 4Dx2 Dy2 S(xi−1 , yj ) + Dx2 Dy2 S(xi−1 , yj+1 ) +4Dx2 Dy2 S(xi , yj−1 ) + 16Dx2 Dy2 S(xi , yj ) + 4Dx2 Dy2 S(xi , yj+1 ) +Dx2 Dy2 S(xi+1 , yj−1 ) + 4Dx2 Dy2 S(xi+1 , yj ) + Dx2 Dy2 S(xi+1 , yj+1 ) = 36h−4 αi,j , i, j = 1, . . . , N,
where αi,j (4.76)
= u(xi−1 , yj−1 ) − 2u(xi , yj−1 ) + u(xi+1 , yj−1 ) −2u(xi−1 , yj ) + 4u(xi , yj ) − 2u(xi+1 , yj ) +u(xi−1 , yj+1 ) − 2u(xi , yj+1 ) + u(xi+1 , yj+1 ).
Using (4.76) and the discussion on pages 290–292 in [9], we have ¯ ¯ −4 ¯h αi,j − Dx2 Dy2 u(xi , yj )¯ ≤ C(u)h2 , (4.77) i, j = 1, . . . , N. 16
Equation (4.75) is equivalent to
(4.78)
Dx2 Dy2 (S − u)(xi−1 , yj−1 ) + 4Dx2 Dy2 (S − u)(xi−1 , yj ) +Dx2 Dy2 (S − u)(xi−1 , yj+1 ) + 4Dx2 Dy2 (S − u)(xi , yj−1 ) +16Dx2 Dy2 (S − u)(xi , yj ) + 4Dx2 Dy2 (S − u)(xi , yj+1 ) +Dx2 Dy2 (S − u)(xi+1 , yj−1 ) + 4Dx2 Dy2 (S − u)(xi+1 , yj ) +Dx2 Dy2 (S − u)(xi+1 , yj+1 ) = 36h−4 αi,j − βi,j , i, j = 1, . . . , N,
where βi,j (4.79)
= Dx2 Dy2 u(xi−1 , yj−1 ) + 4Dx2 Dy2 u(xi−1 , yj ) + Dx2 Dy2 u(xi−1 , yj+1 ) +4Dx2 Dy2 u(xi , yj−1 ) + 16Dx2 Dy2 u(xi , yj ) + 4Dx2 Dy2 u(xi , yj+1 ) +Dx2 Dy2 u(xi+1 , yj−1 ) + 4Dx2 Dy2 u(xi+1 , yj ) + Dx2 Dy2 u(xi+1 , yj+1 ).
Using Taylor’s theorem, we obtain Dx2 Dy2 u(xi−1 , yj±1 ) = Dx2 Dy2 u(xi , yj ) − hDx3 Dy2 u(xi , yj ) ± hDx2 Dy3 Dy u(xi , yj ) + ²± i,j , ± Dx2 Dy2 u(xi+1 , yj±1 ) = Dx2 Dy2 u(xi , yj ) + hDx3 Dy2 u(xi , yj ) ± hDx2 Dy3 Dy u(xi , yj ) + σi,j , ± 2 2 2 2 2 3 Dx Dy u(xi , yj±1 ) = Dx Dy u(xi , yj ) ± hDx Dy u(xi , yj ) + µi,j , ± Dx2 Dy2 u(xi±1 , yj ) = Dx2 Dy2 u(xi , yj ) ± hDx3 Dy2 u(xi , yj ) + νi,j . ¯ ¯ ¯ ±¯ ¯ ±¯ ¯ ±¯ 2 ¯ ¯ ¯ ¯ ¯ ¯ ¯ where ¯²± i,j , σi,j , µi,j , νi,j ≤ C(u)h , i, j = 1, . . . , N , and hence (4.79) gives ¯ ¯ ¯βi,j − 36Dx2 Dy2 u(xi , yj )¯ ≤ C(u)h2 , (4.80) i, j = 1, . . . , N.
It follows from (4.77) and (4.80) that the right-hand side of (4.78) is bounded in absolute value by C(u)h2 . Using (4.2) and moving terms involving Dx2 Dy2 (S − u)(xi , yj ), i = 0, N + 1, j = 1, . . . , N, i = 1, . . . , N, j = 0, N + 1, N
to the right-hand side of (4.78), we obtain a system in {Dx2 Dy2 (S − u)(xi , yj )}i,j=1 with the matrix B ⊗ B, where B is given in (2.10)–(2.12). By (4.51) and symmetry with respect to x and y, each entry on the right-hand side in this system is bounded in absolute value by C(u)h2 . Therefore, (4.60) follows from Lemma 2.3. 2 Lemma 4.6. If u ∈ C 6 (Ω) and (4.81) uh (x, y) =
+1 N +1 N X X
D (x)BnD (y), um,n Bm
S(x, y) =
+1 N +1 N X X
D (x)BnD (y), sm,n Bm
m=0 n=0
m=0 n=0
are defined by (1.3)–(1.6) and (4.2)–(4.5), respectively, then (4.82)
max
1≤m,n≤N
|sm,n − um,n | ≤ C(u)h4 , m, n = 1, . . . , N.
Proof. Using (4.2)–(4.4), (4.11)–(4.13), and following the derivations of (3.3) from (1.3), (3.5) from (1.4), and (3.6) from (1.5), we conclude that (4.83) sm,n = um,n , m = 0, N + 1, n = 0, . . . , N + 1, m = 1, . . . , N, n = 0, N + 1. We define {wi,j }N i,j=1 by (4.84) ∆(S − uh )(xi , yj ) −
h2 2 2 D D (S − uh )(xi , yj ) = wi,j , i, j = 1, . . . , N. 6 x y 17
Using (4.84), (1.6), and (4.10), we obtain h2 wi,j = Dx2 S(xi , yj ) + Dy2 S(xi , yj ) − Dx2 Dy2 S(xi , yj ) 6 ¤ h2 £ 4 −Dx2 u(xi , yj ) − Dy2 u(xi , yj ) + Dx u(xi , yj ) + Dy4 u(xi , yj ) + 2Dx2 Dy2 u(xi , yj ) , 12 and hence (4.58)–(4.60) and the triangle inequality imply that (4.85)
|wi,j | ≤ C(u)h4 , i, j = 1, . . . , N.
Introducing v = [s1,1 − u1,1 , . . . , sN,N − uN,N ]T , w = [w1,1 , . . . , wN,N ]T , using (4.84), (4.81), (4.83), and following the derivation of (3.10) from (1.6), we obtain (4.86)
6h−2 [6T ⊗ I + (6I + T ) ⊗ T ] v = w.
Since h = 1/(N + 1), (4.85) gives h2
N X
2 wi,j ≤ C 2 (u)h8 and hence (4.82) follows from
i,j=1
(4.86), (3.12), and Lemma 2.4. 2 6 Theorem 4.3. If u ∈ C (Ω) and uh and S in S D ⊗S D are defined by (1.3)–(1.6) and (4.2)–(4.5), respectively, then kS − uh kC(Ω) ≤ C(u)h4 . Proof. Since uh − S is continuous on Ω, there is (x∗ , y∗ ) in Ω such that kuh − SkC(Ω) = |(S − uh )(x∗ , y∗ )|. Hence, (3.1), (4.36), (4.83), and the triangle inequality give kuh − SkC(Ω) ≤
N N X X
D |sm,n − um,n ||Bm (x∗ )||BnD (y∗ )|.
m=1 n=1
Since the above double sum reduces to at the most nine terms, the desired result follows from Lemmas 4.6 and 2.1. 2 Theorem 4.4. If u in C 6 (Ω) and uh in S D ⊗ S D are the solutions of (1.1) and (1.3)–(1.6), respectively, then ku − uh kC(Ω) ≤ C(u)h4 . Proof. The required inequality follows from (4.14) and Theorems 4.1, 4.2, 4.3. 2 5. Other Schemes. Consider the scheme for solving (1.1) formulated as follows: Find uh ∈ S D ⊗ S D satisfying (1.3)–(1.5) and (1.9). This scheme is essentially the same as the scheme (4.1)–(4.3) in [4], except that (1.5) is replaced in [4] with 1 [13Dy2 uh (xi , yj ) − 2Dy2 uh (xi , yj+1 ) + Dy2 uh (xi , yj+2 )] = f (xi , yj ), j = 0, 12 1 [D2 uh (xi , yj−2 ) − 2Dy2 uh (xi , yj−1 ) + 13Dy2 uh (xi , yj )] = f (xi , yj ), j = N + 1, 12 y where i = 1, . . . , N . It follows from (3.1), the discussion in section 3, and (2.9) of [4] that the the matrix-vector form of (1.3)–(1.5) and (1.9) is (5.1)
(A ⊗ B + B ⊗ A)u = p, 18
where u = [u1,1 , . . . , uN,N ]T , p = [p1,1 , . . . , pN,N ]T , pi,j = f (xi , yj )−
N +1 X
X
m=0,N +1 n=0
−
N X
X
m=1 n=0,N +1
¤ £ D D D (yj ) (xi )Ly Bm (xi )BnD (yj ) + Bm um,n Lx Bm
¤ £ D D (xi )Ly BnD (yj ) , (xi )BnD (yj ) + Bm um,n Lx Bm
+1 N {ui,j }N j=0 , i = 0, N + 1, {ui,j }i=1 , j = 0, N + 1, are given in (3.3), (3.5), (3.6),
A=
(5.2)
1 (T 2 + 12T ), B = T + 6I, 2h2
and T is defined in (2.12). Lemma 5.1. Assume A, B are as in (5.2) and v = [v1,1 , . . . , vN,N ]T , w = [w1,1 , . . . , wN,N ]T are such that (A ⊗ B + B ⊗ A)v = w. Then 2 max vi,j ≤ Ch2
1≤i,j≤N
N N X X 2 wi,j . i=1 j=1
Proof. It follows from (5.2) that 1 3 6 36 A ⊗ B + B ⊗ A = 2 (T 2 ⊗ T ) + 2 (T 2 ⊗ I) + 2 (T ⊗ T ) + 2 (T ⊗ I) 2h h h h (5.3) 3 6 36 1 + 2 (T ⊗ T 2 ) + 2 (I ⊗ T 2 ) + 2 (T ⊗ T ) + 2 (I ⊗ T ) = 36[r(T ) + s(T )], 2h h h h where for an N × N matrix P , r(P ) = h−2 (P ⊗ I + I ⊗ P ),
(5.4) (5.5)
s(P ) =
¢ ¢ 1 1 ¡ 2 1 ¡ 2 (P ⊗ P ) + P ⊗ P + P ⊗ P2 + P ⊗ I + I ⊗ P2 . 2 2 2 3h 72h 12h
First, we will show that (5.6)
([r(T ) + s(T )]z, z) ≤
2 (r(T )z, z), 9
2
z ∈ RN ,
2
where (·, ·) is the standard inner product in RN . It follows from (2.14) and QT = Q for Q of (2.16) that (r(T )z, z) = ([Q ⊗ Q]r(Λ)[Q ⊗ Q]z, z) = (r(Λ)[Q ⊗ Q]z, [Q ⊗ Q]z), (s(T )z, z) = ([Q ⊗ Q]s(Λ)[Q ⊗ Q]z, z) = (s(Λ)[Q ⊗ Q]z, [Q ⊗ Q]z), where Λ is given in (2.15). Hence (5.6) is equivalent to ([r(Λ) + s(Λ)]z, z) ≤
2 (r(Λ)z, z), 9
which, by (5.4), (5.5), and (2.15), is in turn equivalent to (5.7)
g(λi , λj ) ≤ 0,
i, j = 1, . . . , N, 19
2
z ∈ RN ,
where g(x, y) =
1 1 1 7 (x + y) + xy + (x2 y + xy 2 ) + (x2 + y 2 ). 9 3 72 12
It follows from (2.15) that −4 ≤ λi ≤ 0, i = 1, . . . , N . Hence, (5.7) follows from g(x, y) ≤ 0,
x, y ∈ [−4, 0],
which is established using elementary calculus. The matrices r(T ) and s(T ) are symmetric, r(T )s(T ) = s(T )r(T ), and −r(T ) is positive definite. Hence, (5.6) and 6) on page 135 in [9] imply that (5.8)
kr(T )zk2 ≤
9 k[r(T ) + s(T )]zk2 , 2
2
z ∈ RN .
where k · k2 is the two vector norm. It is known (see, for example, the embedding theorem on page 281 in [9]) that (5.9)
2 ≤ max zi,j
1≤i,j≤N
1 2 h kr(T )zk22 , 4
2
z = [z1,1 , . . . , zN,N ]T ∈ RN .
Hence the desired result follows from (5.9), (5.8), and (5.3). 2 Theorem 5.1. If u ∈ C 6 (Ω) and uh and S are defined by (1.3)–(1.5) and (1.9), and (4.2)–(4.5), respectively, then kS − uh kC(Ω) ≤ C(u)h4 . Proof. Following the proof of Lemma 4.6, we define {wi,j }N i,j=1 by (5.10)
(Lx + Ly )(S − uh )(xi , yj ) = wi,j , i, j = 1, . . . , N.
Using (5.10), (1.9), (1.1), we obtain wi,j = Lx S(xi , yj ) − Dx2 u(xi , yj ) + Ly S(xi , yj ) − Dy2 u(xi , yj ). Equations (1.10), (4.58), and (4.52) give, for i, j = 1, . . . , N , Lx S(xi , yj ) − Dx2 u(xi , yj ) = Dx2 S(xi , yj ) − Dx2 u(xi , yj ) 1 + [Dx2 S(xi−1 , yj ) − 2Dx2 S(xi , yj ) + Dx2 S(xi+1 , yj )] 12 2 h 1 = − Dx4 u(xi , yj ) + [Dx2 u(xi−1 , yj ) − 2Dx2 u(xi , yj ) + Dx2 u(xi+1 , yj )] 12 12 h2 − [D4 u(xi−1 , yj ) − 2Dx4 u(xi , yj ) + Dx4 u(xi+1 , yj )] + ²i,j , 144 x where |²i,j | ≤ C(u)h4 , i, j = 1, . . . , N . Hence Taylor’s theorem and similar considerations for Ly S(xi , yj ) − Dy2 u(xi , yj ) show that (4.85) holds. It follows from (4.81) and (4.83) that the matrix-vector form of (5.10) is (A ⊗ B + B ⊗ A)v = w, where A, B are as in (5.2), v = [s1,1 −u1,1 , . . . , sN,N −uN,N ]T , w = [w1,1 , . . . , wN,N ]T . Hence Lemma 5.1 implies (4.82) and the desired result follows from the proof of Theorem 4.3. 2 Theorem 5.2. If u in C 6 (Ω) and uh in S D ⊗ S D are the solutions of (1.1), and (1.3)–(1.5) and (1.9), respectively, then (5.11)
ku − uh kC(Ω) ≤ C(u)h4 . 20
Proof. The required inequality follows from (4.14) and Theorems 4.1, 4.2, 5.1. 2 It is claimed in Theorem 4.1 of [8] that for the scheme (1.3) and (1.7)–(1.9), one has (5.11) provided that u ∈ C 6 (Ω). The proof of this claim in [8] is based on using Z defined in (4.6)–(4.9) as a comparison function. It is claimed, for example, in Lemma 2.1 of [8] that Z has properties (4.58) and (4.59), that is, (4.58) and (4.59) hold with Z in place of S. Unfortunately, numerical examples indicate that such property does not hold even in one dimensional case. Specifically, for u(x) = x(x − 1)ex and Z ∈ S D such that Z(xi ) = u(xi ), i = 1, . . . , N,
Z 00 (xi ) = u00 (xi ), i = 0, N + 1,
we only have ¯ ¯ ¯ ¯ h2 max ¯¯Z 00 (xi ) − u00 (xi ) + u(4) (xi )¯¯ = Ch2 1≤i≤N 12
and not better. It should be noted that the convergence analysis of [6] for two-point boundary value problems involves the comparison function S ∈ S D defined by S(xi ) = u(xi ), i = 1, . . . , N,
S 00 (xi ) = u00 (xi ) −
h2 (4) u (xi ), i = 0, N + 1, 12
which, in part, was motivation for the definition (4.2)–(4.5). The convergence analysis of the scheme (4.2)–(4.4) in [8] remains an open problem. We believe that such analysis may require proving stability not only with respect to the right-hand side but also with respect to the boundary conditions. 6. Numerical Results. We used scheme (1.3)–(1.6) and algorithm of section 3 to solve a test problem (1.1). The computations were carried out in double precision. We determined the nodal and global errors using the formulas kwkh =
max
0≤i,j≤N +1
kwkC(Ω) ≈
|w(xi , yj )|,
max
0≤i,j≤501
|w(ti , tj )|,
where ti = i/501, i = 1, . . . , 501. Convergence rates were determined using the formula rate =
log(eN/2 /eN ) , log[(N + 1)/(N/2 + 1)]
where eN is the error corresponding to the partition ρx × ρy . We took f in (1.1) corresponding to the exact solution u(x, y) = 3exy (x2 − x)(y 2 − y). We see from the results in Tables 1 and 2 that the scheme (1.3)–(1.6) produces fourth order accuracy for u in both the discrete and the continuous maximum norms. We also observe superconvergence phenomena since the derivative approximations at the partition nodes are of order four. REFERENCES [1] J. H. Ahlberg, E. N. Nilson, and J. L. Walsh, The Theory of Splines and Their Applications, Academic Press, New York, 1967. 21
Table 1 Nodal errors and convergence rates for u, ux , uy , and uxy
N 4 8 16 32 64 128
ku − uh kh Error Rate 7.305–05 6.574–06 4.097 5.036–07 4.040 3.585–08 3.983 2.380–09 4.001 1.534–10 4.000
k(u − uh )x kh Error Rate 9.219–04 8.789–05 3.998 6.715–06 4.044 4.712–07 4.005 3.129–08 4.001 2.017–09 4.000
k(u − uh )y kh Error Rate 9.219–04 8.789–05 3.998 6.715–06 4.044 4.712–07 4.005 3.129–08 4.001 2.017–09 4.000
k(u − uh )xy kh Error Rate 1.673–02 2.134–03 3.504 2.158–04 3.603 1.881–05 3.678 1.497–06 3.734 1.127–07 3.774
Table 2 Global errors and convergence rates for u, ux , uy , and uxy
N 4 8 16 32 64 128
ku − uh kC(Ω) Error Rate 8.630–05 8.606–06 3.922 6.776–07 3.996 4.763–08 4.003 3.157–09 4.003 2.038–10 3.998
k(u − uh )x kC(Ω) Error Rate 1.125–03 1.186–04 3.827 1.232–05 3.561 1.267–06 3.429 1.308–07 3.350 1.467–08 3.192
k(u − uh )y kC(Ω) Error Rate 1.125–03 1.186–0 3.827 1.232–05 3.561 1.267–06 3.429 1.308–07 3.350 1.467–08 3.192
k(u − uh )xy kC(Ω) Error Rate 1.654–02 2.169–03 3.456 2.593–04 3.339 2.997–05 3.253 3.308–06 3.251 3.823–07 3.148
[2] D. Archer, An O(h4 ) cubic spline collocation method for quasilinear parabolic equations, SIAM J. Number. Anal., 14 (1977), 620–637. [3] B. Bialecki, G. Fairweather, and A. Karageorghis, Matrix decomposition algorithms for modified spline collocation for Helmholtz problems, SIAM J. Sci. Comput., 24 (2003), 1733–1753. [4] B. Bialecki, G. Fairweather, and A. Karageorghis, Optimal superconvergent one step nodal cubic spline collocation methods, SIAM J. Sci. Comput., 27 (2005), 575–598. [5] W. Cheney, and D. Kincaid, Numerical Mathematics and Computing, Brooks Cole, California, 1999. [6] J. W. Daniel and B. K. Swartz, Extrapolated collocation for two–point boundary–value problems using cubic splines, J. Inst. Math Appl., 16 (1975), 161–174. [7] C. de Boor, The Method of Projections as Applied to the Numerical Solution of Two Point Boundary Value Problems Using Cubic Splines, Ph.D. thesis, University of Michigan, Ann Arbor, Michigan, 1966. [8] E. N. Houstis, E. A. Vavalis, and J. R. Rice, Convergence of O(h4 ) cubic spline collocation methods for elliptic partial differential equations, SIAM J. Numer. Anal., 25 (1988), 54–74. [9] A. A. Samarski, The Theory of Difference Schemes, Marcel Dekker, Inc., New York, Basel, 2001. [10] C. Van Loan, Computational Frameworks for the Fast Fourier Transform, SIAM, Philadelphia, 1992.
22
† AND
BERNARD BIALECKI‡
Abstract. We present a new modified nodal cubic spline collocation scheme for solving the Dirichlet problem for Poisson’s equation on the unit square. We prove existence and uniqueness of a solution of the scheme and show how the solution can be computed on an (N + 1) × (N + 1) uniform partition of the square with cost O(N 2 logN ) using a direct fast Fourier transform method. Using two comparison functions, we derive an optimal fourth order error bound in the continuous maximum norm. We compare our scheme with other modified nodal cubic spline collocation schemes, in particular, the one proposed by Houstis et al. in [8]. We believe that our paper gives the first correct convergence analysis of a modified nodal cubic spline collocation for solving partial differential equations.
Key words. nodal collocation, cubic splines, convergence analysis, interpolants AMS subject classifications. 65N35, 65N12, 65N15, 65N22
1. Introduction. De Boor [7] proved that classical nodal cubic spline collocation for solving two-point boundary value problems is only second–order accurate and no better. For two-point boundary value problems, Archer [2] and independently Daniel and Swartz [6] developed a modified nodal cubic spline collocation (MNCSC) scheme which is fourth order accurate. The approximate solution in this scheme satisfies higher-order perturbations of the ordinary differential equation at the partition nodes. Based on the method of [2] and [6], Houstis et al. [8] derived a fourth order MNCSC scheme for solving elliptic boundary value problems on rectangles. For the Helmholtz equation, a direct fast Fourier transform (FFT) algorithm for solving this scheme was proposed recently in [3]. In this paper, we consider the Dirichlet boundary value problem for Poisson’s equation (1.1)
∆u = f in Ω,
u = 0 on ∂Ω,
where ∆ denotes the Laplacian, Ω = (0, 1) × (0, 1), and ∂Ω is the boundary of Ω. Let +1 ρx = {xi }N i=0 be a uniform partition of [0, 1] in the x-direction such that xi = ih, i = 0, . . . , N + 1, where h = 1/(N + 1). For the sake of simplicity, we assume that a +1 uniform partition ρy = {yj }N i=0 of [0, 1] in the y-direction is such that yj = xj . Let S3 be the space of cubic splines defined by S3 = {v ∈ C 2 [0, 1] : v|[xi−1 ,xi ] ∈ P3 , i = 1, . . . , N + 1}, where P3 denotes the set of all polynomials of degree ≤ 3, and let S D = {v ∈ S3 : v(0) = v(1) = 0}. Our MNCSC scheme for solving (1.1) is formulated as follows: Find uh ∈ S D ⊗ S D such that (1.2)
∆uh (xi , yj ) −
h2 2 2 h2 Dx Dy uh (xi , yj ) = f (xi , yj ) − ∆f (xi , yj ), 6 12 i, j = 0, . . . , N + 1.
† Department of Mathematical and Computer Sciences, Colorado School of Mines, Golden, Colorado 80401-1887, U.S.A. ([email protected]) ‡ Department of Mathematical and Computer Sciences, Colorado School of Mines, Golden, Colorado 80401-1887, U.S.A. ([email protected])
1
The scheme (1.2) is motivated by the fourth order finite difference method for (1.1), see, for example, equation (7) in section 4.5 of [9]. Using uh = u = 0 on ∂Ω and (1.1), we see that (1.2) is equivalent to: (1.3)
2Dx2 Dy2 uh (xi , yj ) = ∆f (xi , yj ),
i, j = 0, N + 1,
(1.4)
Dx2 uh (xi , yj ) −
h2 2 2 h2 Dx Dy uh (xi , yj ) = f (xi , yj ) − ∆f (xi , yj ), 6 12 i = 0, N + 1, j = 1, . . . , N,
(1.5)
Dy2 uh (xi , yj ) −
h2 2 2 h2 Dx Dy uh (xi , yj ) = f (xi , yj ) − ∆f (xi , yj ), 6 12 i = 1, . . . , N, j = 0, N + 1,
(1.6) ∆uh (xi , yj ) −
h2 h2 2 2 Dx Dy uh (xi , yj ) = f (xi , yj ) − ∆f (xi , yj ), i, j = 1, . . . , N. 6 12
The scheme (4.2)–(4.4) of [8] for (1.1) is: Find uh ∈ S D ⊗ S D satisfying (1.3) and (1.7)
Dx2 uh (xi , yj ) = f (xi , yj ),
i = 0, N + 1, j = 1, . . . , N,
(1.8)
Dy2 uh (xi , yj ) = f (xi , yj ),
i = 1, . . . , N, j = 0, N + 1,
(1.9)
(Lx + Ly )uh (xi , yj ) = f (xi , yj ), i, j = 1, . . . , N,
where, for i, j = 1, . . . , N , (1.10)
Lx v(xi , yj ) =
¤ 1 £ 2 Dx v(xi−1 , yj ) + 10Dx2 v(xi , yj ) + Dx2 v(xi+1 , yj ) , 12
Ly v(xi , yj ) =
¤ 1 £ 2 Dy v(xi , yj−1 ) + 10Dy2 v(xi , yj ) + Dy2 v(xi , yj+1 ) . 12
Our scheme and that of [8] are identical at the corners of Ω. However, they are different at the remaining partition nodes. While (1.4)–(1.6) involve perturbations of both the left- and right-hand sides, (1.9) involves a perturbation of the left-hand side only. Numerical results show that our scheme exhibits superconvergence phenomena while that of [8] does not. An outline of this paper is as follows. We give preliminaries in section 2. The matrix-vector form of our scheme, an existence and uniqueness proof of its solution, and a direct FFT algorithm for solving the scheme are presented in section 3. In section 4, using two comparison functions, we derive a fourth order error bound in the continuous maximum norm. In section 5, we give convergence analysis of the scheme in [4] that consists of (1.3)–(1.5) and (1.9). We also explain why convergence analysis of the scheme (1.3) and (1.7)–(1.9), given in [8], is incorrect. This is why, we believe, our paper gives the first correct convergence analysis of MNCSC for solving partial differential equations. Section 6 includes numerical results obtained using our scheme. 2
+1 2. Preliminaries. We extend the uniform partition ρx = {xi }N i=0 outside of [0, 1] using xi = ih, i = −3, −2, −1, N + 2, N + 3, N + 4, and introduce Ii = [xi−1 , xi ], N +2 i = −2, . . . , N + 4. Let {Bm }m=−1 be the basis for S3 defined by g1 [(x − xm−2 )/h], x ∈ Im−1 , g2 [(x − xm−1 )/h], x ∈ Im , g2 [(xm+1 − x)/h], x ∈ Im+1 , Bm (x) = (2.1) g1 [(xm+2 − x)/h], x ∈ Im+2 , 0, otherwise,
where (2.2)
g1 (x) = x3 ,
g2 (x) = 1 + 3x + 3x2 − 3x3 .
The basis functions are such that, for m = 0, . . . , N + 1, (2.3)
Bm−1 (xm ) = 1, Bm (xm ) = 4, Bm+1 (xm ) = 1, 00 00 00 (xm ) = −12/h2 , Bm+1 (xm ) = 6/h2 . Bm−1 (xm ) = 6/h2 , Bm
D N +1 Let {Bm }m=0 be the basis for S D defined by
(2.4)
B1D = B1 − B−1 , B0D = B0 − 4B−1 , D Bm = Bm , m = 2, . . . , N − 1, D D BN = BN − BN +2 , BN +1 = BN +1 − 4BN +2 .
It follows from (2.3) that (2.5)
B0D (x1 ) = 1, B1D (x1 ) = 4, D D BN (xN −1 ) = 1, BN (xN ) = 4,
£ D ¤00 (x ) = −36/h2 , B £ 0D ¤00 0 (x ) = 6/h2 , B (2.6)£ 0D ¤00 1 B (xN −1 ) = 6/h2 , ¤ £ N D 00 BN (xN +1 ) = 0,
B1D (x2 ) = 1, D BN +1 (xN ) = 1,
£ D ¤00 (x ) = 0, B £ 1D ¤00 0 (x ) = −12/h2 , B £ 1D ¤00 1 B (x ) = −12/h2 , ¤00 N £ N D BN +1 (xN +1 ) = −36/h2 .
£ D ¤00 (x ) = 6/h2 , B £ 1D ¤00 2 BN +1 (xN ) = 6/h2 ,
It also follows from (2.5), (2.6), (2.4), and (2.3) that, for i = 1, . . . , N, ½ h2 D 00 6, m = i, D (2.7) m = 0, . . . , N + 1. ] (xi ) = Bm (xi ) − [Bm 0, m 6= i, 6 Throughout the paper, C denotes a generic positive constant that is independent of u and h. D D N +1 Lemma 2.1. {Bm }m=0 of (2.4) satisfy max |Bm (x)| ≤ C, m = 0, . . . , N + 1. x∈[0,1]
Proof. For each fixed m = −1, . . . , N + 2, using Ii = [xi−1 , xi ] and xi = ih, we have (2.8)
0 ≤ (x − xm−2 )/h ≤ 1, x ∈ Im−1 , 0 ≤ (xm+1 − x)/h ≤ 1, x ∈ Im+1 ,
0 ≤ (x − xm−1 )/h ≤ 1, x ∈ Im , 0 ≤ (xm+2 − x)/h ≤ 1, x ∈ Im+2 .
Equations (2.2) and (2.8) give (2.9)
|g1 [(x − xm−2 )/h]| ≤ 1, x ∈ Im−1 , |g2 [(xm+1 − x)/h]| ≤ 7, x ∈ Im+1 ,
|g2 [(x − xm−1 )/h]| ≤ 7, x ∈ Im , |g1 [(xm+2 − x)/h]| ≤ 1, x ∈ Im+2 . 3
Using (2.1) and (2.9), we see that max |Bm (x)| ≤ C, m = −1, . . . , N + 2. Hence x∈[0,1]
D is a linear the required inequality follows from (2.4) which implies that each Bm N +2 combination of at most two of the functions {Bn }n=−1 . 2 N +1
D For {Bm }m=0 of (2.4), we introduce N × N matrices A and B defined by
(2.10)
N
N
D 00 ] (xi ), B = (bj,n )j,n=1 , bj,n = BnD (yj ). A = (ai,m )i,m=1 , ai,m = [Bm
It follows from (2.4), (2.3), (2.5), and (2.6) that A = 6h−2 T,
(2.11)
B = T + 6I,
where I is the identity matrix and the N × N matrix T −2 1 1 −2 1 . . .. .. ... (2.12) T = 1 −2 1 1 −2
is given by .
Lemma 2.2. If B[u1 , . . . , uN ]T = [v1 , . . . , vN ]T , where B is defined in (2.10), then max |ui | ≤ C max |vi |. 1≤i≤N 1≤i≤N X Proof. It follows from (2.10), (2.11), and (2.12) that |bi,i | − |bi,j | ≥ 2, i = 1, . . . , N . i6=j
Hence the required result follows, for example, from the discussion on page 21 in [1]. 2 In what follows, [u1,1 , . . . , uN,N ]T is the short notation for [u1,1 , . . . , u1,N , u2,1 , . . . , u2,N , . . . , uN,1 , . . . , uN,N ]T . Lemma 2.3. If u = [u1,1 , . . . , uN,N ]T and v = [v1,1 , . . . , vN,N ]T are such that (B ⊗ B)u = v, where B is defined in (2.10), then max |ui,j | ≤ C max |vi,j |. 1≤i,j≤N
1≤i,j≤N
Proof. Since B ⊗ B = (B ⊗ I)(I ⊗ B), we have (2.13)
v = (B ⊗ I)w,
w = (I ⊗ B)u.
Using ( 2.13) and Lemma 2.2, we obtain max |ui,j | ≤ C max |wi,j |,
1≤i,j≤N
max |wi,j | ≤ C max |vi,j |,
1≤i,j≤N
1≤i,j≤N
1≤i,j≤N
which imply the required inequality. 2 It is well known (see Theorem 4.5.2 of [10]) that for T of (2.12), we have (2.14)
QT Q = Λ,
QQ = I,
where the N × N matrices Λ and Q are given by (2.15)
(2.16)
N
Λ = diag(λi )i=1 ,
N
Q = (qi,j )i,j=1 ,
λi = −4sin2
qi,j =
µ
4
2 N +1
iπ , 2(N + 1)
¶1/2
sin
ijπ . N +1
Lemma 2.4. If v = [v1,1 , . . . , vN,N ]T and w = [w1,1 , . . . , wN,N ]T are such that ¶¸ · µ T h2 T T T (2.17) v = w, ⊗ I + I ⊗ + ⊗ h2 h2 6 h2 h2 where T is the matrix defined in (2.12), then
2 max vi,j ≤ Ch2
1≤i,j≤N
N N X X
2 wi,j .
i=1 j=1
Proof. The matrix in (2.17) arises in the fourth order finite difference method for (1.1). Hence the desired result follows, for example, from the last unnumbered equation on page 296 in [9]. 2 Finally, we observe that the matrix-vector form of (2.18)
φi,j =
N X
(1)
ci,m
m=1
N X
(2)
cj,n ψm,n ,
i, j = 1, . . . , N,
n=1
is φ = (C1 ⊗ C2 )ψ,
(2.19) ´N ³ (1) where C1 = ci,m
i,m=1
´N ³ (2) , C2 = cj,n
j,n=1
, and
φ = [φ1,1 , . . . , φN,N ]T ,
ψ = [ψ1,1 , . . . , , ψN,N ]T .
3. Matrix-Vector Form of Scheme. Since dim(S D ⊗ S D ) = (N + 2)2 , the scheme (1.3)–(1.6) involves (N + 2)2 equations in (N + 2)2 unknowns. Using the basis D N +1 }m=0 of (2.4) for the space S D , we have {Bm (3.1)
uh (x, y) =
+1 N +1 N X X
D (x)BnD (y). um,n Bm
m=0 n=0
Substituting (3.1) into (1.3), we obtain (3.2)
2
N +1 N +1 X X
D 00 ] (xi )[BnD ]00 (yj ) = ∆f (xi , yj ), um,n [Bm
i, j = 0, N + 1.
m=0 n=0
Using (2.6), we conclude that (3.2) gives ui,j =
(3.3)
h4 ∆f (xi , yj ), 2592
i, j = 0, N + 1.
Substituting (3.1) into (1.4), we obtain
(3.4)
+1 N +1 N X X
µ ¶ h2 D 00 um,n [Bm ] (xi ) BnD (yj ) − [BnD ]00 (yj ) 6 m=0 n=0 2 h i = 0, N + 1, j = 1, . . . , N. = f (xi , yj ) − ∆f (xi , yj ), 12
Using (2.6) and (2.7), we see that (3.4) gives (3.5)
ui,j = −
h2 h4 f (xi , yj ) + ∆f (xi , yj ), i = 0, N + 1, j = 1, . . . , N. 216 2592 5
Using (3.5) and symmetry with respect to x and y, we conclude that (1.5) gives (3.6)
ui,j = −
h2 h4 f (xi , yj ) + ∆f (xi , yj ), i = 1, . . . , N, j = 0, N + 1. 216 2592
Substituting (3.1) into (1.6), we obtain +1 N +1 N X X
· ¶ ¸ µ h2 D 00 D D 00 (xi ) − [Bm ] (xi )BnD (yj ) + Bm ] (xi ) [BnD ]00 (yj ) um,n [Bm 6 (3.7) m=0 n=0 h2 = f (xi , yj ) − ∆f (xi , yj ), i, j = 1, . . . , N. 12 +1 N Moving the terms involving {um,n }N n=0 , m = 0, N + 1, {um,n }m=1 , n = 0, N + 1, to the right-hand side of (3.7), we get
(3.8)
N N X X
um,n
m=1 n=1
µ
D 00 ] (xi )BnD (yj ) [Bm
· ¶ ¸ h2 D 00 D D 00 + Bm (xi ) − [Bm ] (xi ) [Bn ] (yj ) 6
= pi,j , i, j = 1, . . . , N, where pi,j = f (xi , yj ) − N +1 X
h2 ∆f (xi , yj ) 12 µ
· ¶ ¸ h2 D 00 D D 00 + Bm (xi ) − [Bm ] (xi ) [Bn ] (yj ) um,n − 6 m=0,N +1 n=0 · µ ¶ ¸ N X X h2 D 00 D D 00 (xi ) − [Bm ] (xi )BnD (yj ) + Bm um,n [Bm − ] (xi ) [BnD ]00 (yj ) . 6 m=1 X
D 00 ] (xi )BnD (yj ) [Bm
n=0,N +1
Using (2.18)–(2.19), we write (3.8) as ¶ ¸ · µ h2 (3.9) A ⊗ B + B − A ⊗ A u = p, 6 where u = [u1,1 , . . . , uN,N ]T , p = [p1,1 , . . . , pN,N ]T , and A, B are defined in (2.10). Using (2.11), we see that µ ¶ h2 6 A ⊗ B + B − A ⊗ A = 2 [6T ⊗ I + (6I + T ) ⊗ T ] , 6 h and hence the system (3.9) simplifies to (3.10)
6h−2 [6T ⊗ I + (6I + T ) ⊗ T ] u = p.
We are now ready to prove existence and uniqueness of uh in S D ⊗ S D that satisfies (1.3)–(1.6). Theorem 3.1. There exists unique uh in S D ⊗ S D satisfying (1.3)–(1.6). Proof. Since the number of equations in (1.3)–(1.6) is equal to the number of unknowns, we assume that the right-hand side in (1.3)–(1.6) is zero, and show that uh = 0 is the only solution of the resulting scheme. Using (3.1), (3.3), (3.5), and (3.6), we have (3.11) um,n = 0, m = 0, N + 1, n = 0, . . . , N + 1, m = 1, . . . , N, n = 0, N + 1. 6
Clearly (3.12) 6h
−2
T T h2 [6T ⊗ I + (6I + T ) ⊗ T ] = 36 2 ⊗ I + I ⊗ 2 + h h 6 ·
µ
T T ⊗ 2 h2 h
¶¸
.
Hence it follows from (3.10) with p replaced by 0, (3.12), and Lemma 2.4 that (3.13)
um,n = 0,
m, n = 1, . . . , N.
Equations (3.1), (3.11), and (3.13) give uh = 0. 2 Using Q of (2.16), we see that (3.10) is equivalent to (3.14) 6h−2 (Q ⊗ I) [6T ⊗ I + (6I + T ) ⊗ T ] (Q ⊗ I)(Q−1 ⊗ I)u = (Q ⊗ I)p. Introducing u0 = (Q−1 ⊗ I)u and p0 = (Q ⊗ I)p, and using (3.14) and (2.14), we obtain 6h−2 [6Λ ⊗ I + (6I + Λ) ⊗ T ]u0 = p0 ,
(3.15)
where Λ is defined in (2.15). The system (3.15) reduces to the N independent systems (3.16)
6h−2 [6λi I + (6 + λi )T ]u0i = p0i ,
i = 1, . . . , N,
where u0i = [u0i,1 , . . . , u0i,N ]T , p0i = [p0i,1 , . . . , p0i,N ]T , i = 1, . . . , N . We have the following algorithm for solving (3.10): Step 1. Compute p0 = (Q ⊗ I)p. Step 2. Solve the N systems in (3.16). Step 3. Compute u = (Q ⊗ I)u0 . Since the entries of Q in (2.16) are given in terms of sines, steps 1 and 3 are performed each using FFTs at a cost O(N 2 log N ). In step 2, the systems are tridiagonal, so this step is performed at a cost O(N 2 ). Thus the total cost of the algorithm is O(N 2 log N ). 4. Convergence Analysis. In what follows, C(u) denotes a generic positive constant that is independent of h, but depends on u. Our goal is to show that if u in C 6 (Ω) and uh in S D ⊗ S D are the solutions of (1.1) and (1.3)–(1.6), respectively, then ku − uh kC(Ω) ≤ C(u)h4 ,
(4.1)
where kgkC(Ω) = max |g(x)| for g in C(Ω). x∈Ω
To prove (4.1), for u in C 4 (Ω), we introduce two comparison functions, the spline interpolants S and Z in S D ⊗ S D of u defined respectively by (4.2)
(4.3)
Dx2 Dy2 S(xi , yj ) = Dx2 Dy2 u(xi , yj ),
i, j = 0, N + 1,
h2 2 2 h2 Dx Dy S(xi , yj ) = Dx2 u(xi , yj ) − Dx4 u(xi , yj ) 6 12 h2 2 2 − Dx Dy u(xi , yj ), i = 0, N + 1, j = 1, . . . , N, 6
Dx2 S(xi , yj ) −
7
(4.4)
(4.5)
h2 h2 Dy2 S(xi , yj ) − Dx2 Dy2 S(xi , yj ) = Dy2 u(xi , yj ) − Dy4 u(xi , yj ) 6 12 h2 i = 1, . . . , N, j = 0, N + 1, − Dx2 Dy2 u(xi , yj ), 6
S(xi , yj ) = u(xi , yj ),
i, j = 1, . . . , N,
and (4.6)
Dx2 Dy2 Z(xi , yj ) = Dx2 Dy2 u(xi , yj ),
i, j = 0, N + 1,
(4.7)
Dx2 Z(xi , yj ) = Dx2 u(xi , yj ),
i = 0, N + 1, j = 1, . . . , N,
(4.8)
Dy2 Z(xi , yj ) = Dy2 u(xi , yj ),
i = 1, . . . , N, j = 0, N + 1,
(4.9)
Z(xi , yj ) = u(xi , yj ),
i, j = 1, . . . , N.
It follows from (1.1) that (4.10)
f = Dx2 u + Dy2 u,
∆f = Dx4 u + Dy4 u + 2Dx2 Dy2 u.
Hence, using u = 0 on ∂Ω, we see that (1.3)–(1.5) reduce, respectively, to (4.11)
Dx2 Dy2 uh (xi , yj ) = Dx2 Dy2 u(xi , yj ),
i, j = 0, N + 1,
(4.12)
h2 h2 Dx2 uh (xi , yj ) − Dx2 Dy2 uh (xi , yj ) = Dx2 u(xi , yj ) − Dx4 u(xi , yj ) 6 12 h2 2 2 − Dx Dy u(xi , yj ), i = 0, N + 1, j = 1, . . . , N, 6
(4.13)
h2 h2 Dy2 uh (xi , yj ) − Dx2 Dy2 uh (xi , yj ) = Dy2 u(xi , yj ) − Dy4 u(xi , yj ) 6 12 h2 2 2 i = 1, . . . , N, j = 0, N + 1. − Dx Dy u(xi , yj ), 6
Comparing (4.11)–(4.13) and (4.2)–(4.4), we see that uh and S are defined in the same way for i = 0, N + 1, j = 0, . . . , N + 1, and i = 1, . . . , N , j = 0, N + 1. On the other hand, (4.6)–(4.8) are a simplified, tensor product, version of (4.2)–(4.4). The triangle inequality gives (4.14)
ku − uh kC(Ω) ≤ ku − ZkC(Ω) + kZ − SkC(Ω) + kS − uh kC(Ω) .
In what follows, we bound the three terms on the right-hand side of (4.14). 8
4.1. Bounding ku − ZkC(Ω) . We need the following results. Lemma 4.1. Let the interpolant Ix v in S3 of v in C 2 [0, 1] be defined by (4.15) (Ix v)00 (xi ) = v 00 (xi ), i = 0, N + 1,
Ix v(xi ) = v(xi ), i = 0, . . . , N + 1.
Then max |v(x) − Ix v(x)| ≤ C max |v 00 (x)|h2 .
(4.16)
x∈[0,1]
x∈[0,1]
If v ∈ C 4 [0, 1], then max |v(x) − Ix v(x)| ≤ C max |v (4) (x)|h4 .
(4.17)
x∈[0,1]
x∈[0,1]
Proof. First we prove (4.16). Using the discussion on page 404 in [5], we have (4.18) Ix v(x) = v(xi ) + Bi (x − xi ) + Ci (x − xi )2 + Di (x − xi )3 , x ∈ [xi , xi+1 ], for i = 0, . . . , N , where h 1 ri 1 h (4.19) Bi = − ri+1 − ri + [v(xi+1 ) − v(xi )] , Ci = , Di = (ri+1 − ri ), 6 3 h 2 6h and ri = (Ix v)00 (xi ). Equations (4.18) and (4.19) give (4.20)
h h ri Ix v(x) − v(x) = Ai (x)− ri+1 (x − xi ) − ri (x − xi ) + (x − xi )2 6 3 2 1 3 (ri+1 − ri )(x − xi ) , x ∈ [xi , xi+1 ], + 6h
where Ai (x) = v(xi ) − v(x) +
(4.21)
v(xi+1 ) − v(xi ) (x − xi ), x ∈ [xi , xi+1 ]. h
Using (4.20) and the triangle inequality, we obtain, for x ∈ [xi , xi+1 ], ¶ µ 4 1 (4.22) |Ix v(x) − v(x)| ≤ |Ai (x)| + h2 |ri | + |ri+1 | ≤ |Ai (x)| + h2 max |ri |. 3 3 0≤i≤N +1 We introduce
E=
N +1 (ei,j )i,j=0
1 1 =
4 .. .
1 .. . 1
..
. 4
, 1 1
r = [r0 , . . . , rN +1 ]T , p = [p0 , . . . , pN +1 ]T ,
where (4.23)
pi =
½
v 00 (xi ), i = 0, N + 1, h−2 [v(xi−1 ) − 2v(xi ) + v(xi+1 )], i = 1, . . . , N.
It follows X from the discussion on pages 400 and 401 in [5] that Er = p. Since |ei,i | − |ei,j | ≥ 1, i, j = 0, . . . , N + 1, the discussion on page 21 in [1] implies that i6=j
(4.24)
max
0≤i≤N +1
|ri | ≤ C 9
max
0≤i≤N +1
|pi |.
Using Taylor’s theorem, we obtain |v(xi−1 ) − 2v(xi ) + v(xi+1 )| ≤ Ch2 max |v 00 (x)|,
(4.25)
x∈[0,1]
i = 1, . . . , N.
It follows from (4.24), (4.23) and (4.25), that (4.26)
max
0≤i≤N +1
|ri | ≤ C max |v 00 (x)|. x∈[0,1]
Using Taylors’ theorem to expand v(x), x ∈ [xi , xi+1 ], around xi , we have (x − xi )2 00 v (ξi,x ), 2
v(x) = v(xi ) + (x − xi )v 0 (xi ) +
(4.27)
xi ≤ ξi,x ≤ x.
Using (4.21), (4.27), and the triangle inequality, we obtain, for x ∈ [xi , xi+1 ], ¯ ¯ ¯ (x − xi )2 00 ¯ h 00 ¯ (4.28) |Ai (x)| = ¯ v (ξi,x ) − (x − xi )v (ξi,xi+1 )¯¯ ≤ h2 max |v 00 (x)|. 2 2 x∈[0,1]
Inequality (4.16) follows from (4.22), (4.26), and (4.28). A proof of (4.17) is given in the proof of Theorem 2.3.4 in [1]. 2 Lemma 4.2. If u ∈ C 4 (Ω), Z in S D ⊗ S D is defined by (4.6)–(4.9), and Ix u and Iy u are defined in (4.15), then for (x, y) in Ω, we have (4.29)
Z(x, y) = Ix (Iy u)(x, y),
Dx2 (Iy u)(x, y) = Iy (Dx2 u)(x, y).
+3 Proof. Let {Ci }N i=0 be the basis for S3 such that
Ci (xj ) = δij , Ci00 (xj ) = 0,
i, j = 0, . . . , N + 1, i = 0, . . . , N + 1, j = 0, N + 1,
(4.30) CN +2 (xj ) = CN +3 (xj ) = 0,
j = 0, . . . , N + 1,
00 00 CN +2 (x0 ) = CN +3 (xN +1 ) = 1,
00 00 CN +2 (xN +1 ) = CN +3 (x0 ) = 0,
where δij is the Kronecker delta. Using (4.15) and (4.30), we have for (x, y) ∈ Ω, Ix (Iyu)(x, y) N +1 X = Ix u(x, yj )Cj (y) + Dy2 u(x, y0 )CN +2 (y) + Dy2 u(x, yN +1 )CN +3 (y) j=0
=
N +1 X
N +1 X
u(xi , yj )Cj (y) + Dy2 u(xi , y0 )CN +2 (y) + Dy2 u(xi , yN +1 )CN +3 (y) Ci (x)
i=0 j=0 N +1 X + Dx2 u(x0 , yj )Cj (y) + Dx2 Dy2 u(x0 , y0 )CN +2 (y) j=0
N +1 X ¤ Dx2 u(xN +1 , yj )Cj (y) +Dx2 Dy2 u(x0 , yN +1 )CN +3 (y) CN +2 (x) +
j=0 ¤ +Dx2 Dy2 u(xN +1 , y0 )CN +2 (y) + Dx2 Dy2 u(xN +1 , yN +1 )CN +3 (y) CN +3 (x).
10
Since u = 0 on ∂Ω, all terms involving C0 (x), CN +1 (x), C0 (y), CN +1 (y) drop out which implies that Ix (Iy u) ∈ S D ⊗ S D . Using (4.30), we verify that Ix (Iy u) satisfies (4.6)–(4.9), that is, (4.6)–(4.9) hold with Ix (Iy u) in place of Z. Hence, the uniqueness of the interpolant Z implies the first equation in (4.29). To prove the second equation in (4.29), we use (4.15) and (4.30) to see that for (x, y) ∈ Ω, Iy (Dx2 u)(x, y) N +1 X Dx2 u(x, yj )Cj (y) + Dx2 Dy2 u(x, y0 )CN +2 (y) + Dx2 Dy2 u(x, yN +1 )CN +3 (y) = j=0 N +1 X u(x, yj )Cj (y) + Dy2 u(x, y0 )CN +2 (y) + Dy2 u(x, yN +1 )CN +3 (y) = Dx2 j=0
= Dx2 (Iy u)(x, y).
2
Theorem 4.1. If u ∈ C 4 (Ω) and Z in S D ⊗ S D is defined by (4.6)–(4.9), then ku − ZkC(Ω) ≤ C(u)h4 . Proof. Using (4.29) and the triangle inequality, we have (4.31)
ku − ZkC(Ω)
≤ ku − Ix ukC(Ω) + kIx (u − Iy u) − (u − Iy u)kC(Ω) + ku − Iy ukC(Ω) .
For any fixed y in [0, 1], Ix u(·, y) is the cubic spline interpolant of u(·, y). Using this, symmetry with respect to x and y, and (4.17), we have (4.32)
ku − Ix ukC(Ω) ≤ C(u)h4 ,
ku − Iy ukC(Ω) ≤ C(u)h4 .
For any fixed y in [0, 1], Ix (u−Iy u)(·, y) is the cubic spline interpolant of (u−Iy u)(·, y). Hence it follows from (4.16) that (4.33)
kIx (u − Iy u) − (u − Iy u)kC(Ω) ≤ CkDx2 (u − Iy u)kC(Ω) h2 .
Using (4.29) and(4.16), we obtain (4.34)
kDx2 (u − Iy u)kC(Ω) = kDx2 u − Iy (Dx2 u)kC(Ω) ≤ CkDx2 Dy2 ukC(Ω) h2 .
Combining (4.33) and (4.34), we have (4.35)
kIx (u − Iy u) − (u − Iy u)kC(Ω) ≤ C(u)h4 .
The desired inequality now follows from (4.31), (4.32), and (4.35).
2
4.2. Bounding kZ − SkC(Ω) . We start by proving the following lemma. Lemma 4.3. If u ∈ C 4 (Ω) and (4.36) S(x, y) =
N +1 N +1 X X
D sm,n Bm (x)BnD (y), Z(x, y) =
m=0 n=0
N +1 N +1 X X
D zm,n Bm (x)BnD (y),
m=0 n=0
are defined by (4.2)–(4.5) and (4.6)–(4.9), respectively, then |sm,n − zm,n | ≤ C(u)h4 , 11
m, n = 0, . . . , N + 1.
Proof. Using (4.2), (4.6), and following the derivation of (3.3) from (1.3), we obtain (4.37)
sm,n = zm,n ,
m, n = 0, N + 1.
Next we prove the required inequality for m = 0, n = 1, . . . , N . Using (4.7), we have (4.38)
Dx2 (S − Z)(x0 , yj ) = Dx2 S(x0 , yj ) − Dx2 u(x0 , yj ),
j = 1, . . . , N.
It follows from (4.36), (4.37), and (2.6) that (4.39) Dx2 (S − Z)(x0 , yj ) = −36h−2
N X
(s0,n − z0,n )BnD (yj ),
j = 1, . . . , N.
n=1
Using (4.36), (2.6), (2.4), (2.3), and (2.5), we obtain, for j = 1, . . . , N , (4.40) Dx2 S(x0 , yj ) = −36h−2
N +1 X
s0,n BnD (yj ) = −36h−2 (s0,j−1 + 4s0,j + s0,j+1 ).
n=0
Substituting (4.39) and (4.40) into (4.38), and multiplying through by −h2 /36, we have (4.41)
N X
(s0,n − z0,n )BnD (yj ) = s0,j−1 + 4s0,j + s0,j+1 +
n=1
h2 2 D u(x0 , yj ) 36 x
for j = 1, . . . , N . Using (4.2), (4.3), and following the derivations of (3.3) from (1.3) and (3.5) from (1.4), we obtain s0,j =
(4.42)
h4 D2 D2 u(x0 , yj ), 1296 x y
j = 0, N + 1,
and (4.43)
s0,j
¸ · h2 2 2 h2 h2 4 2 =− Dx u(x0 , yj ) − Dx u(x0 , yj ) − Dx Dy u(x0 , yj ) 216 12 6
for j = 1, . . . , N . Since u = 0 on ∂Ω, (4.42) is the same as (4.43) with j = 0, N + 1. This observation and (4.43) imply that for j = 1, . . . , N , we have ¸ · h2 h2 h2 (4.44) s0,j±1 = − Dx2 u(x0 , yj±1 ) − Dx4 u(x0 , yj±1 ) − Dx2 Dy2 u(x0 , yj±1 ) . 216 12 6 Using Taylor’s theorem, we obtain (4.45) Dx2 u(x0 , yj±1 ) = Dx2 u(x0 , yj ) ± hDx2 Dy u(x0 , yj ) +
h2 2 2 D D u(x0 , ξj± ), 2 x y
where yj−1 ≤ ξj− ≤ yj , yj ≤ ξj+ ≤ yj+1 . Using (4.44), (4.43), and (4.45), we obtain ¯ ¯ ¯ ¯ h2 2 ¯ (4.46) ¯s0,j−1 + 4s0,j + s0,j+1 + Dx u(x0 , yj )¯¯ ≤ C(u)h4 , j = 1, . . . , N. 36 N
It follows from (4.46) that (4.41) is a system in {s0,n − z0,n }n=1 with the matrix B defined in(2.10) and with each entry on the right-hand side bounded in absolute value by C(u)h4 . Hence, Lemma 2.2 implies (4.47)
max |s0,n − z0,n | ≤ C(u)h4 .
1≤n≤N
12
Using (4.47) and symmetry with respect to x and y, we also have max |sN +1,n − zN +1,n | ≤ C(u)h4 ,
1≤n≤N
(4.48)
max |sm,n − zm,n | ≤ C(u)h4 ,
1≤m≤N
n = 0, N + 1.
Finally we prove the required inequality for m, n = 1, . . . , N. Using (4.5) and (4.9), we have (S − Z)(xi , yj ) = 0, i, j = 1, . . . , N , which, by (4.36) and (4.37), can be written as (4.49)
N X N X
D (sm,n − zm,n )Bm (xi )BnD (yj ) = di,j ,
i, j = 1, . . . N,
m=1 n=1
where
di,j =
X
N X
m=0,N +1 n=1
+
N X
X
m=1 n=0,N +1
D (zm,n − sm,n )Bm (xi )BnD (yj ).
Since for any fixed i, j, each of the above double sums reduces to at most three terms, using the triangle inequality, (4.47), (4.48), and Lemma 2.1, we obtain (4.50)
|di,j | ≤ C(u)h4 ,
i, j = 1, . . . , N. N
It follows from (2.18)–(2.19) that (4.49) is a system in {zm,n − sm,n }m,n=1 with the matrix B ⊗ B, where B is defined in (2.10). Hence, for m, n = 1, . . . , N, the required inequality follows from (4.50) and Lemma 2.3. 2 Theorem 4.2. If u ∈ C 4 (Ω) and S, Z in S D ⊗ S D are defined by (4.2)–(4.5) and (4.6)–(4.9), respectively, then kZ − SkC(Ω) ≤ C(u)h4 . Proof. Since Z − S is continuous on Ω, there is (x∗ , y∗ ) in Ω such that kZ − SkC(Ω) = |(Z − S)(x∗ , y∗ )| . Hence, (4.36) and the triangle inequality imply kZ − SkC(Ω) ≤
N +1 N +1 X X
D |sm,n − zm,n ||Bm (x∗ )||BnD (y∗ )|.
m=0 n=0
Since the above double sum reduces to at most nine terms, the required inequality follows from Lemmas 4.3 and 2.1. 2 4.3. Bounding kS − uh kC(Ω) and ku − uh kC(Ω) . We need the following results. Lemma 4.4. If u ∈ C 6 (Ω) and S in S D ⊗ S D is defined by (4.2)–(4.5), then for i = 0, N + 1, j = 1, . . . , N , ¯ 2 2 ¯ ¯Dx Dy S(xi , yj ) − Dx2 Dy2 u(xi , yj )¯ ≤ C(u)h2 , (4.51) (4.52)
¯ ¯ 2 ¯ ¯ 2 ¯Dx S(xi , yj ) − Dx2 u(xi , yj ) + h Dx4 u(xi , yj )¯ ≤ C(u)h4 . ¯ ¯ 12 13
Proof. We prove (4.51) for i = 0; for i = N + 1, (4.51) follows by symmetry with respect to x. Using (4.36), we obtain Dx2 Dy2 S(x0 , yj ) =
+1 N +1 N X X
m=0 n=0
¤00 £ £ D ¤00 (x0 ) BnD (yj ), sm,n Bm
j = 1, . . . , N,
and hence (2.4), (2.3), and (2.6) imply (4.53)
Dx2 Dy2 S(x0 , yj ) = −216h−4 (s0,j−1 − 2s0,j + s0,j+1 ), j = 1, . . . , N.
Equations (4.53), (4.43), and (4.44) give, for j = 1, . . . , N,
(4.54)
Dx2 Dy2£S(x0 , yj ) − Dx2 Dy2 u(x0 , yj ) = −Dx2 Dy2 u(x0 , yj ) ¤ +h−2 Dx2 u(x0 , yj−1 ) − 2Dx2 u(x0 , yj ) + Dx2 u(x0 , yj+1 ) ¤ 1 £ 4 − Dx u(x0 , yj−1 ) − 2Dx4 u(x0 , yj ) + Dx4 u(x0 , yj+1 ) 12 ¤ 1£ − Dx2 Dy2 u(x0 , yj−1 ) − 2Dx2 Dy2 u(x0 , yj ) + Dx2 Dy2 u(x0 , yj+1 ) . 6
Using Taylor’s theorem, we obtain
(4.55)
Dx2 u(x0 , yj±1 )
h2 2 2 D D u(x0 , yj ) 2 x y 3 4 h h ± Dx2 Dy3 u(x0 , yj ) + Dx2 Dy4 u(x0 , ξj± ), 3! 4!
= Dx2 u(x0 , yj ) ± hDx2 Dy u(x0 , yj ) +
(4.56) Dx4 u(x0 , yj±1 ) = Dx4 u(x0 , yj ) ± hDx4 Dy u(x0 , yj ) +
h2 4 2 D D u(x0 , ηj± ), 2 x y
(4.57) Dx2 Dy2 u(x0 , yj±1 ) = Dx2 Dy2 u(x0 , yj ) ± hDx2 Dy3 u(x0 , yj ) +
h2 2 4 D D u(x0 , κ± j ), 2 x y
+ + + where yj−1 ≤ ξj− , ηj− , κ− j ≤ yj , yj ≤ ξj , ηj , κj ≤ yj+1 . Equations (4.55)–(4.57) give ¯ −2 £ 2 ¯ ¤ ¯h Dx u(x¯0 , yj−1 ) − 2Dx2 u(x0 , yj ) + Dx2 u(x0 , yj+1 ) − Dx2¯Dy2 u(x0 , yj )¯ ≤ C(u)h2 , ¯ x4 u(x0 , yj−1 ) − 2Dx4 u(x0 , yj ) + Dx4 u(x0 , yj+1 )¯ ≤ C(u)h2 , ¯ 2 D ¯ ¯Dx Dy2 u(x0 , yj−1 ) − 2Dx2 Dy2 u(x0 , yj ) + Dx2 Dy2 u(x0 , yj+1 )¯ ≤ C(u)h2 ,
and hence (4.51) for i = 0 follows from (4.54) and the triangle inequality. Using (4.3) and (4.51), we obtain (4.52). 2 Lemma 4.5. If u ∈ C 6 (Ω) and S in S D ⊗ S D is defined by (4.2)–(4.5), then, for i, j = 1, . . . , N , we have ¯ ¯ ¯ ¯ 2 h2 4 4 2 ¯ ¯ (4.58) ¯Dx S(xi , yj ) − Dx u(xi , yj ) + 12 Dx u(xi , yj )¯ ≤ C(u)h , (4.59)
(4.60)
¯ ¯ 2 ¯ ¯ 2 ¯Dy S(xi , yj ) − Dy2 u(xi , yj ) + h Dy4 u(xi , yj )¯ ≤ C(u)h4 , ¯ ¯ 12 |Dx2 Dy2 S(xi , yj ) − Dx2 Dy2 u(xi , yj )| ≤ C(u)h2 . 14
Proof. First we prove (4.58). For i = 0, . . . , N + 1, j = 1, . . . , N , we introduce ¸ · h2 4 2 2 (4.61) di,j = Dx S(xi , yj ) − Dx u(xi , yj ) − Dx u(xi , yj ) . 12 Then (4.62)
di−1,j + 4di,j + di+1,j = φi,j − ψi,j ,
i, j = 1, . . . , N,
where φi,j (4.63)
ψi,j (4.64)
2 2 = Dx2 S(x · i−1 , yj ) + 4Dx2S(xi , yj ) + D ¸x S(xi+1 , yj ) h −6 Dx2 u(xi , yj ) + Dx4 u(xi , yj ) , 12
· ¸ h2 4 h2 4 2 = − Dx u(xi−1 , yj ) + 4 Dx u(xi , yj ) − Dx u(xi , yj ) 12 12 · ¸ h2 4 h2 4 2 2 +Dx u(xi+1 , yj ) − Dx u(xi+1 , yj ) − 6 Dx u(xi , yj ) + Dx u(xi , yj ) 12 12 = Dx2 u(xi−1 , yj ) − 2Dx2 u(xi , yj ) + Dx2 u(xi+1 , yj ) ¤ h2 £ 4 − Dx u(xi−1 , yj ) + 10Dx4 u(xi , yj ) + Dx4 u(xi+1 , yj ) . 12 Dx2 u(xi−1 , yj )
Since S(·, yj ) ∈ S3 , (2.1.7) in [1], (4.5), and S = u = 0 on ∂Ω, imply that (4.65)
Dx2 S(xi−1 , yj ) + 4Dx2 S(xi , yj ) + Dx2 S(xi+1 , yj ) = 6h−2 [u(xi−1 , yj ) − 2u(xi , yj ) + u(xi+1 , yj )] , i, j = 1, . . . , N.
Using Taylor’s theorem, we obtain h3 h2 u(xi±1 , yj ) = u(xi , yj ) ± hDx u(xi , yj ) + Dx2 u(xi , yj )± Dx3 u(xi , yj ) 2 3! h4 4 h5 5 h6 6 + Dx u(xi , yj )± Dx u(xi , yj ) + Dx u(ξi± , yj ), 4! 5! 6! where xi−1 ≤ ξi− ≤ xi , xi ≤ ξi+ ≤ xi+1 , and hence ¯ −2 ¯h [u(xi−1 , yj ) − 2u(xi , yj ) + u(xi+1 , yj )] ¸¯ · ¯ h2 (4.66) − Dx2 u(xi , yj ) + Dx4 u(xi , yj ) ¯¯ ≤ C(u)h4 , i, j = 1, . . . , N. 12 Using (4.63), (4.65), and (4.66), we obtain (4.67)
|φi,j | ≤ C(u)h4 , i, j = 1, . . . , N.
Using Taylor’s theorem, we obtain h2 Dx2 u(xi±1 , yj ) = Dx2 u(xi , yj ) ± hDx3 u(xi , yj ) + Dx4 u(xi , yj ) 2 h4 h3 ± Dx5 u(xi , yj ) + Dx6 u(ξi± , yj ), 3! 4! h2 4 4 5 Dx u(xi±1 , yj ) = Dx u(xi , yj ) ± hDx u(xi , yj ) + Dx6 u(ηi± , yj ), 2 15
where xi−1 ≤ ξi− , ηi− ≤ xi , xi ≤ ξi+ , ηi+ ≤ xi+1 , and hence (4.64) gives |ψi,j | ≤ C(u)h4 ,
(4.68)
i, j = 1, . . . , N.
Using (4.61) and (4.52), we have |di,j | ≤ C(u)h4 ,
(4.69)
i = 0, N + 1,
j = 1, . . . , N.
It follows from (4.67)–(4.69) that moving d0,j and dN +1,j to the right-hand side of N (4.62), we obtain, for each j = 1, . . . , N , a system in {di,j }i=1 with the matrix B of (2.10)–(2.12), and with each entry on the right-hand side bounded in absolute value by C(u)h4 . Hence (4.58) follows from (4.61) and Lemma 2.2, and (4.59) follows from (4.58) by symmetry with respect to x and y. Next we prove (4.60). Since S(x, ·) ∈ S3 for x ∈ [0, 1], (2.1.7) in [1] gives (4.70)
Dy2 S(x, yj−1 ) + 4Dy2 S(x, yj ) + Dy2 S(x, yj+1 ) = 6h−2 [S(x, yj−1 ) − 2S(x, yj ) + S(x, yj+1 )] , j = 1, . . . , N, x ∈ [0, 1].
Differentiating (4.70) twice with respect to x, we obtain, for j = 1, . . . , N , x ∈ [0, 1], (4.71)
2 2 yj ) + Dx2 Dy2 S(x, yj+1 )¤ Dx2 Dy2 S(x, £ 2yj−1 ) + 4Dx Dy S(x, 2 −2 Dx S(x, yj−1 ) − 2Dx S(x, yj ) + Dx2 S(x, yj+1 ) . = 6h
Using (4.71) with x = xi−1 , xi , xi+1 , we obtain, for i, j = 1, . . . , N , (4.72)
(4.73)
(4.74)
Dx2 Dy2 S(x , yj−1 ) + 4Dx2 Dy2 S(xi−1 , yj ) + Dx2 Dy2 S(xi−1 , yj+1 )¤ £ i−1 2 −2 Dx S(xi−1 , yj−1 ) − 2Dx2 S(xi−1 , yj ) + Dx2 S(xi−1 , yj+1 ) , = 6h 2 2 2 2 Dx2 Dy2 S(x i , yj ) + Dx Dy S(xi , yj+1 )¤ £ i ,2yj−1 ) + 4Dx Dy S(x −2 2 = 6h Dx S(xi , yj−1 ) − 2Dx S(xi , yj ) + Dx2 S(xi , yj+1 ) ,
, yj−1 ) + 4Dx2 Dy2 S(xi+1 , yj ) + Dx2 Dy2 S(xi+1 , yj+1 )¤ Dx2 Dy2 S(x £ i+1 2 −2 Dx S(xi+1 , yj−1 ) − 2Dx2 S(xi+1 , yj ) + Dx2 S(xi+1 , yj+1 ) . = 6h
Adding (4.72), (4.74) and (4.73) multiplied through by 4, and using (4.65) and S = u = 0 on ∂Ω, we obtain
(4.75)
Dx2 Dy2 S(xi−1 , yj−1 ) + 4Dx2 Dy2 S(xi−1 , yj ) + Dx2 Dy2 S(xi−1 , yj+1 ) +4Dx2 Dy2 S(xi , yj−1 ) + 16Dx2 Dy2 S(xi , yj ) + 4Dx2 Dy2 S(xi , yj+1 ) +Dx2 Dy2 S(xi+1 , yj−1 ) + 4Dx2 Dy2 S(xi+1 , yj ) + Dx2 Dy2 S(xi+1 , yj+1 ) = 36h−4 αi,j , i, j = 1, . . . , N,
where αi,j (4.76)
= u(xi−1 , yj−1 ) − 2u(xi , yj−1 ) + u(xi+1 , yj−1 ) −2u(xi−1 , yj ) + 4u(xi , yj ) − 2u(xi+1 , yj ) +u(xi−1 , yj+1 ) − 2u(xi , yj+1 ) + u(xi+1 , yj+1 ).
Using (4.76) and the discussion on pages 290–292 in [9], we have ¯ ¯ −4 ¯h αi,j − Dx2 Dy2 u(xi , yj )¯ ≤ C(u)h2 , (4.77) i, j = 1, . . . , N. 16
Equation (4.75) is equivalent to
(4.78)
Dx2 Dy2 (S − u)(xi−1 , yj−1 ) + 4Dx2 Dy2 (S − u)(xi−1 , yj ) +Dx2 Dy2 (S − u)(xi−1 , yj+1 ) + 4Dx2 Dy2 (S − u)(xi , yj−1 ) +16Dx2 Dy2 (S − u)(xi , yj ) + 4Dx2 Dy2 (S − u)(xi , yj+1 ) +Dx2 Dy2 (S − u)(xi+1 , yj−1 ) + 4Dx2 Dy2 (S − u)(xi+1 , yj ) +Dx2 Dy2 (S − u)(xi+1 , yj+1 ) = 36h−4 αi,j − βi,j , i, j = 1, . . . , N,
where βi,j (4.79)
= Dx2 Dy2 u(xi−1 , yj−1 ) + 4Dx2 Dy2 u(xi−1 , yj ) + Dx2 Dy2 u(xi−1 , yj+1 ) +4Dx2 Dy2 u(xi , yj−1 ) + 16Dx2 Dy2 u(xi , yj ) + 4Dx2 Dy2 u(xi , yj+1 ) +Dx2 Dy2 u(xi+1 , yj−1 ) + 4Dx2 Dy2 u(xi+1 , yj ) + Dx2 Dy2 u(xi+1 , yj+1 ).
Using Taylor’s theorem, we obtain Dx2 Dy2 u(xi−1 , yj±1 ) = Dx2 Dy2 u(xi , yj ) − hDx3 Dy2 u(xi , yj ) ± hDx2 Dy3 Dy u(xi , yj ) + ²± i,j , ± Dx2 Dy2 u(xi+1 , yj±1 ) = Dx2 Dy2 u(xi , yj ) + hDx3 Dy2 u(xi , yj ) ± hDx2 Dy3 Dy u(xi , yj ) + σi,j , ± 2 2 2 2 2 3 Dx Dy u(xi , yj±1 ) = Dx Dy u(xi , yj ) ± hDx Dy u(xi , yj ) + µi,j , ± Dx2 Dy2 u(xi±1 , yj ) = Dx2 Dy2 u(xi , yj ) ± hDx3 Dy2 u(xi , yj ) + νi,j . ¯ ¯ ¯ ±¯ ¯ ±¯ ¯ ±¯ 2 ¯ ¯ ¯ ¯ ¯ ¯ ¯ where ¯²± i,j , σi,j , µi,j , νi,j ≤ C(u)h , i, j = 1, . . . , N , and hence (4.79) gives ¯ ¯ ¯βi,j − 36Dx2 Dy2 u(xi , yj )¯ ≤ C(u)h2 , (4.80) i, j = 1, . . . , N.
It follows from (4.77) and (4.80) that the right-hand side of (4.78) is bounded in absolute value by C(u)h2 . Using (4.2) and moving terms involving Dx2 Dy2 (S − u)(xi , yj ), i = 0, N + 1, j = 1, . . . , N, i = 1, . . . , N, j = 0, N + 1, N
to the right-hand side of (4.78), we obtain a system in {Dx2 Dy2 (S − u)(xi , yj )}i,j=1 with the matrix B ⊗ B, where B is given in (2.10)–(2.12). By (4.51) and symmetry with respect to x and y, each entry on the right-hand side in this system is bounded in absolute value by C(u)h2 . Therefore, (4.60) follows from Lemma 2.3. 2 Lemma 4.6. If u ∈ C 6 (Ω) and (4.81) uh (x, y) =
+1 N +1 N X X
D (x)BnD (y), um,n Bm
S(x, y) =
+1 N +1 N X X
D (x)BnD (y), sm,n Bm
m=0 n=0
m=0 n=0
are defined by (1.3)–(1.6) and (4.2)–(4.5), respectively, then (4.82)
max
1≤m,n≤N
|sm,n − um,n | ≤ C(u)h4 , m, n = 1, . . . , N.
Proof. Using (4.2)–(4.4), (4.11)–(4.13), and following the derivations of (3.3) from (1.3), (3.5) from (1.4), and (3.6) from (1.5), we conclude that (4.83) sm,n = um,n , m = 0, N + 1, n = 0, . . . , N + 1, m = 1, . . . , N, n = 0, N + 1. We define {wi,j }N i,j=1 by (4.84) ∆(S − uh )(xi , yj ) −
h2 2 2 D D (S − uh )(xi , yj ) = wi,j , i, j = 1, . . . , N. 6 x y 17
Using (4.84), (1.6), and (4.10), we obtain h2 wi,j = Dx2 S(xi , yj ) + Dy2 S(xi , yj ) − Dx2 Dy2 S(xi , yj ) 6 ¤ h2 £ 4 −Dx2 u(xi , yj ) − Dy2 u(xi , yj ) + Dx u(xi , yj ) + Dy4 u(xi , yj ) + 2Dx2 Dy2 u(xi , yj ) , 12 and hence (4.58)–(4.60) and the triangle inequality imply that (4.85)
|wi,j | ≤ C(u)h4 , i, j = 1, . . . , N.
Introducing v = [s1,1 − u1,1 , . . . , sN,N − uN,N ]T , w = [w1,1 , . . . , wN,N ]T , using (4.84), (4.81), (4.83), and following the derivation of (3.10) from (1.6), we obtain (4.86)
6h−2 [6T ⊗ I + (6I + T ) ⊗ T ] v = w.
Since h = 1/(N + 1), (4.85) gives h2
N X
2 wi,j ≤ C 2 (u)h8 and hence (4.82) follows from
i,j=1
(4.86), (3.12), and Lemma 2.4. 2 6 Theorem 4.3. If u ∈ C (Ω) and uh and S in S D ⊗S D are defined by (1.3)–(1.6) and (4.2)–(4.5), respectively, then kS − uh kC(Ω) ≤ C(u)h4 . Proof. Since uh − S is continuous on Ω, there is (x∗ , y∗ ) in Ω such that kuh − SkC(Ω) = |(S − uh )(x∗ , y∗ )|. Hence, (3.1), (4.36), (4.83), and the triangle inequality give kuh − SkC(Ω) ≤
N N X X
D |sm,n − um,n ||Bm (x∗ )||BnD (y∗ )|.
m=1 n=1
Since the above double sum reduces to at the most nine terms, the desired result follows from Lemmas 4.6 and 2.1. 2 Theorem 4.4. If u in C 6 (Ω) and uh in S D ⊗ S D are the solutions of (1.1) and (1.3)–(1.6), respectively, then ku − uh kC(Ω) ≤ C(u)h4 . Proof. The required inequality follows from (4.14) and Theorems 4.1, 4.2, 4.3. 2 5. Other Schemes. Consider the scheme for solving (1.1) formulated as follows: Find uh ∈ S D ⊗ S D satisfying (1.3)–(1.5) and (1.9). This scheme is essentially the same as the scheme (4.1)–(4.3) in [4], except that (1.5) is replaced in [4] with 1 [13Dy2 uh (xi , yj ) − 2Dy2 uh (xi , yj+1 ) + Dy2 uh (xi , yj+2 )] = f (xi , yj ), j = 0, 12 1 [D2 uh (xi , yj−2 ) − 2Dy2 uh (xi , yj−1 ) + 13Dy2 uh (xi , yj )] = f (xi , yj ), j = N + 1, 12 y where i = 1, . . . , N . It follows from (3.1), the discussion in section 3, and (2.9) of [4] that the the matrix-vector form of (1.3)–(1.5) and (1.9) is (5.1)
(A ⊗ B + B ⊗ A)u = p, 18
where u = [u1,1 , . . . , uN,N ]T , p = [p1,1 , . . . , pN,N ]T , pi,j = f (xi , yj )−
N +1 X
X
m=0,N +1 n=0
−
N X
X
m=1 n=0,N +1
¤ £ D D D (yj ) (xi )Ly Bm (xi )BnD (yj ) + Bm um,n Lx Bm
¤ £ D D (xi )Ly BnD (yj ) , (xi )BnD (yj ) + Bm um,n Lx Bm
+1 N {ui,j }N j=0 , i = 0, N + 1, {ui,j }i=1 , j = 0, N + 1, are given in (3.3), (3.5), (3.6),
A=
(5.2)
1 (T 2 + 12T ), B = T + 6I, 2h2
and T is defined in (2.12). Lemma 5.1. Assume A, B are as in (5.2) and v = [v1,1 , . . . , vN,N ]T , w = [w1,1 , . . . , wN,N ]T are such that (A ⊗ B + B ⊗ A)v = w. Then 2 max vi,j ≤ Ch2
1≤i,j≤N
N N X X 2 wi,j . i=1 j=1
Proof. It follows from (5.2) that 1 3 6 36 A ⊗ B + B ⊗ A = 2 (T 2 ⊗ T ) + 2 (T 2 ⊗ I) + 2 (T ⊗ T ) + 2 (T ⊗ I) 2h h h h (5.3) 3 6 36 1 + 2 (T ⊗ T 2 ) + 2 (I ⊗ T 2 ) + 2 (T ⊗ T ) + 2 (I ⊗ T ) = 36[r(T ) + s(T )], 2h h h h where for an N × N matrix P , r(P ) = h−2 (P ⊗ I + I ⊗ P ),
(5.4) (5.5)
s(P ) =
¢ ¢ 1 1 ¡ 2 1 ¡ 2 (P ⊗ P ) + P ⊗ P + P ⊗ P2 + P ⊗ I + I ⊗ P2 . 2 2 2 3h 72h 12h
First, we will show that (5.6)
([r(T ) + s(T )]z, z) ≤
2 (r(T )z, z), 9
2
z ∈ RN ,
2
where (·, ·) is the standard inner product in RN . It follows from (2.14) and QT = Q for Q of (2.16) that (r(T )z, z) = ([Q ⊗ Q]r(Λ)[Q ⊗ Q]z, z) = (r(Λ)[Q ⊗ Q]z, [Q ⊗ Q]z), (s(T )z, z) = ([Q ⊗ Q]s(Λ)[Q ⊗ Q]z, z) = (s(Λ)[Q ⊗ Q]z, [Q ⊗ Q]z), where Λ is given in (2.15). Hence (5.6) is equivalent to ([r(Λ) + s(Λ)]z, z) ≤
2 (r(Λ)z, z), 9
which, by (5.4), (5.5), and (2.15), is in turn equivalent to (5.7)
g(λi , λj ) ≤ 0,
i, j = 1, . . . , N, 19
2
z ∈ RN ,
where g(x, y) =
1 1 1 7 (x + y) + xy + (x2 y + xy 2 ) + (x2 + y 2 ). 9 3 72 12
It follows from (2.15) that −4 ≤ λi ≤ 0, i = 1, . . . , N . Hence, (5.7) follows from g(x, y) ≤ 0,
x, y ∈ [−4, 0],
which is established using elementary calculus. The matrices r(T ) and s(T ) are symmetric, r(T )s(T ) = s(T )r(T ), and −r(T ) is positive definite. Hence, (5.6) and 6) on page 135 in [9] imply that (5.8)
kr(T )zk2 ≤
9 k[r(T ) + s(T )]zk2 , 2
2
z ∈ RN .
where k · k2 is the two vector norm. It is known (see, for example, the embedding theorem on page 281 in [9]) that (5.9)
2 ≤ max zi,j
1≤i,j≤N
1 2 h kr(T )zk22 , 4
2
z = [z1,1 , . . . , zN,N ]T ∈ RN .
Hence the desired result follows from (5.9), (5.8), and (5.3). 2 Theorem 5.1. If u ∈ C 6 (Ω) and uh and S are defined by (1.3)–(1.5) and (1.9), and (4.2)–(4.5), respectively, then kS − uh kC(Ω) ≤ C(u)h4 . Proof. Following the proof of Lemma 4.6, we define {wi,j }N i,j=1 by (5.10)
(Lx + Ly )(S − uh )(xi , yj ) = wi,j , i, j = 1, . . . , N.
Using (5.10), (1.9), (1.1), we obtain wi,j = Lx S(xi , yj ) − Dx2 u(xi , yj ) + Ly S(xi , yj ) − Dy2 u(xi , yj ). Equations (1.10), (4.58), and (4.52) give, for i, j = 1, . . . , N , Lx S(xi , yj ) − Dx2 u(xi , yj ) = Dx2 S(xi , yj ) − Dx2 u(xi , yj ) 1 + [Dx2 S(xi−1 , yj ) − 2Dx2 S(xi , yj ) + Dx2 S(xi+1 , yj )] 12 2 h 1 = − Dx4 u(xi , yj ) + [Dx2 u(xi−1 , yj ) − 2Dx2 u(xi , yj ) + Dx2 u(xi+1 , yj )] 12 12 h2 − [D4 u(xi−1 , yj ) − 2Dx4 u(xi , yj ) + Dx4 u(xi+1 , yj )] + ²i,j , 144 x where |²i,j | ≤ C(u)h4 , i, j = 1, . . . , N . Hence Taylor’s theorem and similar considerations for Ly S(xi , yj ) − Dy2 u(xi , yj ) show that (4.85) holds. It follows from (4.81) and (4.83) that the matrix-vector form of (5.10) is (A ⊗ B + B ⊗ A)v = w, where A, B are as in (5.2), v = [s1,1 −u1,1 , . . . , sN,N −uN,N ]T , w = [w1,1 , . . . , wN,N ]T . Hence Lemma 5.1 implies (4.82) and the desired result follows from the proof of Theorem 4.3. 2 Theorem 5.2. If u in C 6 (Ω) and uh in S D ⊗ S D are the solutions of (1.1), and (1.3)–(1.5) and (1.9), respectively, then (5.11)
ku − uh kC(Ω) ≤ C(u)h4 . 20
Proof. The required inequality follows from (4.14) and Theorems 4.1, 4.2, 5.1. 2 It is claimed in Theorem 4.1 of [8] that for the scheme (1.3) and (1.7)–(1.9), one has (5.11) provided that u ∈ C 6 (Ω). The proof of this claim in [8] is based on using Z defined in (4.6)–(4.9) as a comparison function. It is claimed, for example, in Lemma 2.1 of [8] that Z has properties (4.58) and (4.59), that is, (4.58) and (4.59) hold with Z in place of S. Unfortunately, numerical examples indicate that such property does not hold even in one dimensional case. Specifically, for u(x) = x(x − 1)ex and Z ∈ S D such that Z(xi ) = u(xi ), i = 1, . . . , N,
Z 00 (xi ) = u00 (xi ), i = 0, N + 1,
we only have ¯ ¯ ¯ ¯ h2 max ¯¯Z 00 (xi ) − u00 (xi ) + u(4) (xi )¯¯ = Ch2 1≤i≤N 12
and not better. It should be noted that the convergence analysis of [6] for two-point boundary value problems involves the comparison function S ∈ S D defined by S(xi ) = u(xi ), i = 1, . . . , N,
S 00 (xi ) = u00 (xi ) −
h2 (4) u (xi ), i = 0, N + 1, 12
which, in part, was motivation for the definition (4.2)–(4.5). The convergence analysis of the scheme (4.2)–(4.4) in [8] remains an open problem. We believe that such analysis may require proving stability not only with respect to the right-hand side but also with respect to the boundary conditions. 6. Numerical Results. We used scheme (1.3)–(1.6) and algorithm of section 3 to solve a test problem (1.1). The computations were carried out in double precision. We determined the nodal and global errors using the formulas kwkh =
max
0≤i,j≤N +1
kwkC(Ω) ≈
|w(xi , yj )|,
max
0≤i,j≤501
|w(ti , tj )|,
where ti = i/501, i = 1, . . . , 501. Convergence rates were determined using the formula rate =
log(eN/2 /eN ) , log[(N + 1)/(N/2 + 1)]
where eN is the error corresponding to the partition ρx × ρy . We took f in (1.1) corresponding to the exact solution u(x, y) = 3exy (x2 − x)(y 2 − y). We see from the results in Tables 1 and 2 that the scheme (1.3)–(1.6) produces fourth order accuracy for u in both the discrete and the continuous maximum norms. We also observe superconvergence phenomena since the derivative approximations at the partition nodes are of order four. REFERENCES [1] J. H. Ahlberg, E. N. Nilson, and J. L. Walsh, The Theory of Splines and Their Applications, Academic Press, New York, 1967. 21
Table 1 Nodal errors and convergence rates for u, ux , uy , and uxy
N 4 8 16 32 64 128
ku − uh kh Error Rate 7.305–05 6.574–06 4.097 5.036–07 4.040 3.585–08 3.983 2.380–09 4.001 1.534–10 4.000
k(u − uh )x kh Error Rate 9.219–04 8.789–05 3.998 6.715–06 4.044 4.712–07 4.005 3.129–08 4.001 2.017–09 4.000
k(u − uh )y kh Error Rate 9.219–04 8.789–05 3.998 6.715–06 4.044 4.712–07 4.005 3.129–08 4.001 2.017–09 4.000
k(u − uh )xy kh Error Rate 1.673–02 2.134–03 3.504 2.158–04 3.603 1.881–05 3.678 1.497–06 3.734 1.127–07 3.774
Table 2 Global errors and convergence rates for u, ux , uy , and uxy
N 4 8 16 32 64 128
ku − uh kC(Ω) Error Rate 8.630–05 8.606–06 3.922 6.776–07 3.996 4.763–08 4.003 3.157–09 4.003 2.038–10 3.998
k(u − uh )x kC(Ω) Error Rate 1.125–03 1.186–04 3.827 1.232–05 3.561 1.267–06 3.429 1.308–07 3.350 1.467–08 3.192
k(u − uh )y kC(Ω) Error Rate 1.125–03 1.186–0 3.827 1.232–05 3.561 1.267–06 3.429 1.308–07 3.350 1.467–08 3.192
k(u − uh )xy kC(Ω) Error Rate 1.654–02 2.169–03 3.456 2.593–04 3.339 2.997–05 3.253 3.308–06 3.251 3.823–07 3.148
[2] D. Archer, An O(h4 ) cubic spline collocation method for quasilinear parabolic equations, SIAM J. Number. Anal., 14 (1977), 620–637. [3] B. Bialecki, G. Fairweather, and A. Karageorghis, Matrix decomposition algorithms for modified spline collocation for Helmholtz problems, SIAM J. Sci. Comput., 24 (2003), 1733–1753. [4] B. Bialecki, G. Fairweather, and A. Karageorghis, Optimal superconvergent one step nodal cubic spline collocation methods, SIAM J. Sci. Comput., 27 (2005), 575–598. [5] W. Cheney, and D. Kincaid, Numerical Mathematics and Computing, Brooks Cole, California, 1999. [6] J. W. Daniel and B. K. Swartz, Extrapolated collocation for two–point boundary–value problems using cubic splines, J. Inst. Math Appl., 16 (1975), 161–174. [7] C. de Boor, The Method of Projections as Applied to the Numerical Solution of Two Point Boundary Value Problems Using Cubic Splines, Ph.D. thesis, University of Michigan, Ann Arbor, Michigan, 1966. [8] E. N. Houstis, E. A. Vavalis, and J. R. Rice, Convergence of O(h4 ) cubic spline collocation methods for elliptic partial differential equations, SIAM J. Numer. Anal., 25 (1988), 54–74. [9] A. A. Samarski, The Theory of Difference Schemes, Marcel Dekker, Inc., New York, Basel, 2001. [10] C. Van Loan, Computational Frameworks for the Fast Fourier Transform, SIAM, Philadelphia, 1992.
22
