Domain Decomposition Spectral Method for Mixed ... - Semantic Scholar
J Sci Comput (2012) 53:451–480 DOI 10.1007/s10915-012-9581-z
Domain Decomposition Spectral Method for Mixed Inhomogeneous Boundary Value Problems of High Order Differential Equations on Unbounded Domains Chao Zhang · Ben-yu Guo
Received: 30 December 2010 / Revised: 6 February 2012 / Accepted: 12 February 2012 / Published online: 28 February 2012 © Springer Science+Business Media, LLC 2012
Abstract In this paper, we develop domain decomposition spectral method for mixed inhomogeneous boundary value problems of high order differential equations defined on unbounded domains. We introduce an orthogonal family of new generalized Laguerre functions, with the weight function x α , α being any real number. The corresponding quasiorthogonal approximation and Gauss-Radau type interpolation are investigated, which play important roles in the related spectral and collocation methods. As examples of applications, we propose the domain decomposition spectral methods for two fourth order problems, and the spectral method with essential imposition of boundary conditions. The spectral accuracy is proved. Numerical results demonstrate the effectiveness of suggested algorithms. Keywords Laguerre quasi-orthogonal approximation · Laguerre-Gauss-Radau type interpolation · Domain decomposition spectral methods · High order problems with mixed inhomogeneous boundary conditions
The work of C. Zhang is supported in part by NSF of China No. 11171227 and Research Fund for young teachers of Jiangsu Normal University No. 11XLR27. The work of B.-y. Guo is supported in part by NSF of China No. 11171227, Fund for Doctor Authority of Chinese Educational Ministry No. 20080270001, Shanghai Leading Academic Discipline Project No. S30405, and Fund for E-institute of Shanghai Universities No. E03004. C. Zhang Jiangsu Normal University, Xuzhou 221116, China e-mail: [email protected] B.-y. Guo () Shanghai Normal University, Shanghai 200234, China e-mail: [email protected] B.-y. Guo Scientific Computing Key Laboratory, Division of Computational Science of E-institute, Shanghai Universities, Shanghai, China
452
J Sci Comput (2012) 53:451–480
1 Introduction Spectral and pseudospectral methods have been used successfully for scientific computation, see [3, 5–7, 9, 10, 20] and the references therein. During the past decade, more and more attentions were paid to spectral and pseudospectral methods for unbounded domains. Among those approaches, four direct and commonly used methods are the Hermite spectral and pseudospectral methods, and the Laguerre spectral and pseudospectral methods, see [8, 11, 13, 15, 21, 23] and the references therein. The existing results on the Laguerre approximation are mostly suitable for second order problems. However, it is also interesting to consider high order problems, cf. [4, 15]. Recently, Guo, Sun and Zhang [19] proposed the generalized Laguerre quasi-orthogonal approximation with the weight function x α e−βx , where α is any real number and β > 0. By adjusting the parameters α and β, it matches considered functions at the origin properly, and simulates the asymptotic behaviors more reasonably. This approach provides a tool for treating high order problems with solutions growing up at infinity. Whereas, in many practical cases, such as some problems arising in fluid dynamics and astrophysics, the solutions do not grow up at infinity. In those cases, the above framework is no longer the most appreciate. On the other hand, the solutions may vary rapidly near the fixed boundaries. For raising numerical accuracy and saving work, we might use the Legendre approximation near the boundaries and the Laguerre approximation on the remaining domains. It leads to domain decomposition spectral method, cf. [2, 12]. Furthermore, for solving certain problems in statistical physics, such as the Fokker-Planck equation in an infinite channel, we also need domain decomposition spectral method, cf. [14]. But the presence of weight x α e−βx brings difficulties in domain decomposition spectral method for high order problems. In this paper, we consider another spectral method on unbounded domains. We shall introduce a new family of generalized Laguerre functions, which are mutually orthogonal with the weight function x α , α being any real number. Then we develop the related orthogonal approximation, which is more suitable for approximated functions decaying or growing algebraically at infinity. We next propose the corresponding quasi-orthogonal approximation, with which we could deal with various mixed inhomogeneous boundary value problems easily. In particular, the Laguerre quasi-orthogonal approximation with negative integer α, serves as a basic tool in domain decomposition spectral method of high order problems, since it keeps certain continuity on all common boundaries of adjacent subdomains. Moreover, such new approximation is applicable to numerical solutions of some exterior problems. We also study the related Laguerre-Gauss-Radau type interpolation, which is the mathematical foundation of the new Laguerre pseudospectral and collocation methods. As examples of applications, we provide the domain decomposition spectral methods for two high order problems, and prove their spectral accuracy. We also propose spectral method with essential imposition of boundary conditions for a model problem, with the error analysis. Numerical results demonstrate the high accuracy of suggested algorithms. It is noted that Shen [21], and Guo and Ma [12] considered the case with α = 0, β = 1, and Guo and Zhang [17] studied the case with α > −1, β > 0. The related results are available for second order problems. Our new results improve and generalize the existing results, and so enlarge the applications of Laguerre approximation essentially. In particular, they are very useful for domain decomposition spectral method of high order problems with inhomogeneous boundary conditions, as well as exterior problems. This paper is organized as follows. The next section is for preliminaries. In Sect. 3, we develop the new Laguerre quasi-orthogonal approximation. In Sect. 4, we study the new Laguerre-Gauss-Radau type interpolation. Section 5 is devoted to domain decomposition
J Sci Comput (2012) 53:451–480
453
spectral methods and spectral method with essential imposition of boundary conditions. In Sect. 6, we present some numerical results. The final section is for concluding remarks.
2 Preliminaries In this section, we recall the recent results on the generalized Laguerre orthogonal approximation and the Laguerre-Gauss-Radau type interpolation. Let = {x | 0 < x < ∞} and χ (x) be certain a weight function. For integer r ≥ 0, we define the weighted Sobolev space Hχr () in the usual way, with the inner product (· , ·)r,χ , , the semi-norm | · |r,χ , and the norm · r,χ , . In particular, the inner product and the norm dk v of L2χ () are denoted by (· , ·)χ , and · χ , , respectively. For simplicity, we denote dx k by ∂xk v. For integer r ≥ 1, r 0 Hχ ()
= v ∈ Hχr () | ∂xk v(0) = 0, 0 ≤ k ≤ r − 1 .
We omit the subscript χ in notations whenever χ (x) ≡ 1. The scaled generalized Laguerre polynomials of degree l ≥ 0 were given in [16], as (α,β)
Ll
(x) =
1 −α βx l l+α −βx x e ∂x x e , l!
α > −1, β > 0, l ≥ 0.
Now, let α be any real number. Denote by [α] the largest integer ≤ α. Let l¯α = [−α] for α ≤ −1, and l¯α = 0 for α > −1. Meanwhile, lα = l − [−α] for α ≤ −1, and lα = l for α > −1. The new generalized Laguerre functions with real parameter α of degree l ≥ 0 were introduced in [19], namely, (α,β) Ll (x)
=
x −α Llα (x), (α,β) Ll (x), (−α,β)
α ≤ −1, l ≥ l¯α = [−α], α > −1, l ≥ l¯α = 0.
(2.1)
They satisfy the following Sturm-Liouville problem (see (4.6) of [19]), (α,β) (α,β) (α,β) ∂x x α+1 e−βx ∂x Ll (x) + λl x α e−βx Ll (x) = 0,
l ≥ l¯α ,
(2.2)
with the eigenvalues (α,β) λl
=
β(lα − α) = β(l − [−α] − α), βlα = βl,
for α ≤ −1, for α > −1.
(2.3)
Next, let ωα,β (x) = x α e−βx , and (α,β)
γl
=
(l + α + 1) , β α+1 l!
α > −1, β > 0, l ≥ 0.
According to (4.8) of [19], we have (α,β)
Ll
(α,β)
(x)Ll
(α,β)
(x)ωα,β (x)dx = ηl
δl,l
(2.4)
454
J Sci Comput (2012) 53:451–480
where δl,l is the Kronecker symbol, and (−α,β) (−α,β) = γl−[−α] , for α ≤ −1, l ≥ l α , γα (α,β) ηl = l(α,β) γl , for α > −1, l ≥ 0. (α,β)
The set of all Ll
(2.5)
(x) is complete in L2ωα,β (). Hence, for any v ∈ L2ωα,β (), we have v(x) =
∞
(α,β)
vˆl
(α,β)
Ll
(2.6)
(x)
l=l¯α
with (α,β) vˆl
=
1
(α,β)
v(x)Ll
(α,β) ηl
(2.7)
(x)ωα,β (x)dx.
Now, let N ≥ l¯α be non-negative integer, and (α,β) (α,β) (α,β) (α,β) QN () = span Ll¯ (x), Ll¯ +1 (x), . . . , LN (x) . α
α
(α,β)
The orthogonal projection PN,α,β, : L2ωα,β () → QN (PN,α,β, v − v, φ)ωα,β , = 0,
() is defined by (α,β)
∀φ ∈ QN
().
Due to Theorem 4.2 of [19], we have the following basic result. Lemma 2.1 If ∂xν v ∈ L2ω−m+ν,β () for ν = 0, k, r, and integers 1 ≤ m ≤ min(r, N ), then for 0 ≤ k ≤ r ≤ N + 1, k k−r ∂ (PN,−m,β, v − v) ≤ c(βN ) 2 ∂xr v ω . (2.8) x ω , , −m+r,β
−m+k,β
Hereafter, we denote by c a generic positive constant independent of any function and N . Next, we introduce the following space with integer r ≥ 0, Hωr −m,β ,A () = v | v is measurable on and vHωr
−1 and β > 0, we (α,β) (α,β) denote by ξG,N,j, (0 ≤ j ≤ N ) the zeros of the polynomial LN+1 (x), which are arranged in (α,β) ascending order. Meanwhile, ωG,N,j, (0 ≤ j ≤ N ) stand for the corresponding Christoffel numbers such that φ(x)ωα,β (x)dx =
N (α,β) (α,β) φ ξG,N,j, ωG,N,j, ,
∀φ ∈ P2N+1 ().
(2.11)
j =0
We know from (2.10) of [18] that (α,β)
ωG,N,j, =
1
(N + α + 2) , (α,β) (α,β) (α,β) α β (N + 2) ξG,N,j, [∂x LN+1 (ξG,N,j, )]2
0 ≤ j ≤ N.
Next, let m ≥ 0, and (m,β)
(m,β)
ξN,j, = ξG,N−m,j, , (m,β) −2m (m,β) (m,β) ωN,j, = ωG,N−m,j, ξG,N−m,j, ,
(2.12)
0 ≤ j ≤ N − m.
Thanks to (5.9) of [19], we have N−m
(m,β) (m,β) (m,β) φ ξN,j, ψ ξN,j, ωN,j, = (φ, ψ)ω−m,β , ,
(−m,β)
∀φ ∈ QN
(−m,β)
(), ψ ∈ QN+1 ().
j =0
(2.13) For any integer r ≥ 0, we denote by C r () the space consisting of all r-times differential functions. Further, for integer r ≥ m − 1, r () = v ∈ C r () | ∂xk v(0) = 0, 0 ≤ k ≤ m − 1 . C0,m,β m−1 () and m ≥ 1, the new Laguerre-Gauss-Radau type interpolation For any v ∈ C0,m,β (−m,β) IN,−m,β, v ∈ QN () is determined uniquely by
(m,β)
(m,β)
IN,−m,β, v ξN,j, = v ξN,j, ,
0 ≤ j ≤ N − m.
(2.14)
Thanks to Lemma 5.2 of [19], we have the following result. Lemma 2.3 If v ∈ 0 Hωm−m,β ,A () ∩ Hωk−m+k ,β (), ∂xr v ∈ L2ω−m+r,β () and integers 1 ≤ m ≤ min(r, N ), then for 0 ≤ k ≤ r ≤ N + 1, k ∂ (IN,−m,β, v − v) x
ω−m+k,β ,
1 1 k+1−r ≤ c β − 2 + 1 (ln N ) 2 (βN ) 2 ∂xr v ω
−m+r,β ,
In the end of this section, we derive an equality, which will be used in Sect. 5.
.
(2.15)
456
J Sci Comput (2012) 53:451–480
Proposition 2.1 If u ∈ 0 Hω20,β (), then 2 2 ∂ u x
2 4 2 − β x 2 ∂ e 2 u + β ∂x e− β2 x u 2 + β e− β2 x u 2 . = x ω0,β , 2 16
(2.16)
β
Proof Let v(x) = e− 2 x u(x). Then
β β2 ∂x2 u(x) = ∂x2 v(x) + β∂x v(x) + v(x) e 2 x . 4 Since u ∈ 0 Hω20,β (), we have v ∈ 0 H 2 (). Therefore,
v(x)∂x v(x)dx =
∂x v(x)∂x2 v(x)dx = 0,
v(x)∂x2 v(x)dx = −∂x v2 .
Consequently, 2 2 ∂ u x ω
2 β2 v(x) dx 4 4 2 β β2 = ∂x2 v + β 2 ∂x v2 + v2 + v(x)∂x2 v(x)dx 16 2
0,β ,
=
∂x2 v(x) + β∂x v(x) +
2 β2 β4 = ∂x2 v + ∂x v2 + v2 . 2 16
The proof is completed.
3 New Laguerre Quasi-Orthogonal Approximation In this section, we develop the new Laguerre quasi-orthogonal approximation. 3.1 New Laguerre Orthogonal Approximation For any real number α and β > 0, the new generalized Laguerre functions are given by L˜ l
(α,β)
β
(x) = e− 2 x Ll
(α,β)
(x),
l ≥ l¯α .
(α,β) Due to (2.2), L˜ l (x) is the lth eigenfunction of the following Sturm-Liouville problem,
β β (α,β) ∂x x α+1 e−βx ∂x e 2 x v(x) + λl x α e− 2 x v(x) = 0,
l ≥ l¯α .
(3.1)
(α,β) Let ω˜ α (x) = x α . By virtue of (2.4), the set of all L˜ l (x) conforms a complete L2ω˜ α ()orthogonal system, namely, (α,β) (α,β) (α,β) L˜ l (x)L˜ l (x)ω˜ α (x)dx = ηl δl,l . (3.2)
J Sci Comput (2012) 53:451–480
457
Thus, for any v ∈ L2ω˜ α (), we have v(x) =
∞
(α,β)
v˜l
L˜ l
(α,β)
(3.3)
(x)
l=l¯α
with (α,β)
v˜l
=
1
v(x)L˜ l
(α,β)
(α,β)
ηl
(x)ω˜ α (x)dx.
(3.4)
Next, let β (α,β) (α,β) Q˜ N () = e− 2 x ψ | ψ ∈ QN () . ˜ (α,β) The L2ω˜ α ()-orthogonal projection P˜N,α,β, : L2ω˜ α () → Q () is defined by N (P˜N,α,β, v − v, φ)ω˜ α , = 0,
(α,β) ∀φ ∈ Q˜ N (),
or equivalently, P˜N,α,β, v(x) =
N
(α,β)
v˜l
L˜ l
(α,β)
(3.5)
(x).
l=l¯α
We now estimate the approximation error. To do this, we introduce the Sturm-Liouville operator: β β A˜ α,β v(x) = −x −α e 2 x ∂x x α+1 e−βx ∂x e 2 x v(x) . Evidently, A˜ α,β L˜ l
(α,β)
(α,β)
(x) = λl
L˜ l
(α,β)
(x),
l ≥ l¯α .
(3.6)
By integration by parts, we obtain
(α,β)
v, L˜ l
ω˜ α ,
(α,β) −1 (α,β) −1 (α,β) (α,β) = λl = λl A˜ α,β v, L˜ l . v, A˜ α,β L˜ l ω˜ α , ω˜ α ,
(3.7)
Therefore, if u, v are in the domain of the operator A˜ α,β , then (A˜ α,β u, v)ω˜ α , = (u, A˜ α,β v)ω˜ α , . Accordingly, A˜ α,β is a positive definite and self-conjugate operator. Thus, we could define the following Sobolev-type spaces with integer r ≥ 0, D A˜ rα,β = v | A˜ kα,β v ∈ L2ω˜ α (), 0 ≤ k ≤ r , equipped with the following semi-norm and norm, |v|D(A˜ r
α,β
r ˜ ) = Aα,β v
ω˜ α ,
vD(A˜ r
,
α,β )
=
r k=0
12 |v|2D(A˜ k ) α,β
.
Theorem 3.1 For any v ∈ D(A˜ rα,β ) and integers 0 ≤ μ ≤ r ≤ N + 1, |P˜N,α,β, v − v|D(A˜ μ
α,β )
≤ c(βN )μ−r |v|D(A˜ rα,β ) .
(3.8)
458
J Sci Comput (2012) 53:451–480
Proof We use (3.5) and (3.6) to deduce that |P˜N,α,β, v − v|2D(A˜ μ
α,β )
=
∞
(α,β) 2μ (α,β) (α,β) 2 ηl v˜l
λl
l=N+1
=
∞
(α,β) 2μ (α,β) −1 (α,β) 2 . ηl v, L˜ l ω˜ α ,
λl
l=N+1
Further, using (3.7) repeatedly, yields
(α,β) v, L˜ l
ω˜ α ,
(α,β) −r r (α,β) = λl A˜ α,β v, L˜ l . ω˜ α ,
The above two equalities, together with (2.3), lead to that |P˜N,α,β, v − v|2D(A˜ μ
α,β )
∞
=
(α,β) 2μ−2r (α,β) −1 ˜ r (α,β) 2 Aα,β v, L˜ l ηl ω˜ α ,
λl
l=N+1 ∞ (α,β) 2μ−2r (α,β) −1 r (α,β) 2 ≤ c λN+1 A˜ α,β v, L˜ l ηl ω˜
α ,
l=N+1
≤ c(βN )2μ−2r
∞
(α,β) −1
ηl
(α,β) 2 . ω˜ α ,
A˜ rα,β v, L˜ l
l=l¯α
On the other hand, thanks to (3.2)–(3.4), we have ∞ r 2 (α,β) −1 r (α,β) 2 A˜ v = A˜ α,β v, L˜ l ηl α,β ω˜ α , ω˜ l=l¯α
α ,
.
Consequently, |P˜N,α,β, v − v|D(A˜ μ
α,β )
≤ c(βN )μ−r A˜ rα,β v ω˜ α , = c(βN )μ−r |v|D(A˜ rα,β ) .
The proof is completed. 3.2 Laguerre Orthogonal Approximation with Negative Integer α
We now focus on a special Laguerre orthogonal approximation, which is the mostly appreciated for numerical solutions of high order differential equations defined on unbounded domains. In the forthcoming discussions, we assume integer m ≥ 1. Clearly, ∂xk L˜ l
(−m,β)
(0) = 0,
0 ≤ k ≤ m − 1, l ≥ m.
Proposition 3.1 For any v ∈ 0 Hω1α,β () and α ≤ 0, vωα,β , ≤
2 ∂x vωα,β , . β
(3.9)
J Sci Comput (2012) 53:451–480
459
Proof For any x ∈ , we have v 2 (x)x α e−βx =
x
∂ξ v 2 (ξ )ξ α e−βξ dξ
0
x
=2
α −βξ
v(ξ )∂ξ v(ξ )ξ e 0
x
dξ + α
v 2 (ξ )ξ α−1 e−βξ dξ
0
x
−β
v 2 (ξ )ξ α e−βξ dξ.
(3.10)
0
Since v ∈ 0 Hω1α,β (), we have v 2 (x)x α e−βx → 0 as x → ∞. Letting x → ∞ in (3.10), and using the Cauchy-Schwarz inequality, we derive that v2ωα,β , ≤
4 2α 4 v2ωα−1,β , ≤ 2 ∂x v2ωα,β , . ∂x v2ωα,β , + 2 β β β
This is the desired result. We next introduce the following space with integer r ≥ 0, Hωr˜ −m ,A () = v | v is measurable on and vHωr˜
−m ,A
max(m, n), we introduce the spaces r r (I ) = v | v is measurable on I and vHm,n,A Hm,n,A (I ) < ∞ , r r H0,m,n,A (I ) = v ∈ Hm,n,A (I ) | ∂xk v(−2) = 0 for 0 ≤ k ≤ n − 1, and ∂xk v(0) = 0 for 0 ≤ k ≤ m − 1 ,
466
J Sci Comput (2012) 53:451–480
equipped with the following semi-norm and norm, r r |v|Hm,n,A (I ) = ∂x v χ (−m+r,−n+r) ,I ,
r vHm,n,A (I ) =
r k=0
12 |v|2H k (I ) m,n,A
.
Next, Jl(m,n) (x) stands for the standard Jacobi polynomial of degree l, and (m,n) Yl(m,n) (x) = χ (m,n) (x)Jl−m−n (x + 1), x ∈ I, l ≥ m + n, (I ) = span Yl(m,n) (x), m + n ≤ l ≤ N . Q(m,n) N μ,0 μ For integers max(m, n) ≤ μ ≤ m + n, the projection PN,m,n,I : H0,m,n,A (I ) → Q(m,n) (I ) is N defined by
μ,0 v , ∂xμ φ χ (−m+μ,−n+μ) = 0, ∂xμ v − PN,m,n,I
∀φ ∈ Q(m,n) (I ). N
− + We also introduce the polynomials qm,n,j,I (x), qm,n,j,I (x) ∈ Pm+n−1 (I ) as follows,
− (x) = qm,n,j,I
n−1−j (m + l − 1)! 1 m (−x) (2 + x)l+j , m l l!(m − 1)! 2 j! 2 l=0
m−1−j (n + l − 1)! (−1)j + (−x)l+j . (x) = n (2 + x)n qm,n,j,I l l!(n − 1)! 2 j! 2 l=0
(5.1)
It can be checked that − (−2) = δk,j , ∂xk qm,n,j,I + (−2) ∂xl qm,n,j,I
− ∂xl qm,n,j,I (0) = 0, + ∂xk qm,n,j,I (0)
= 0,
= δk,j ,
0 ≤ j, k ≤ n − 1, 0 ≤ l ≤ m − 1, 0 ≤ j, k ≤ m − 1, 0 ≤ l ≤ n − 1.
(5.2)
r For any v ∈ Hm,n,A (I ) and r ≥ max(m, n), we set
vm,n,b,I (x) =
n−1 j =0
− ∂xj v(−2)qm,n,j,I (x) +
m−1
+ ∂xj v(0)qm,n,j,I (x).
j =0
Let v(x) ˜ = v(x) − vm,n,b,I (x). The Jacobi quasi-orthogonal projection is defined by (cf. (2.18) of [19]) μ μ,0 PN,m,n,I v(x) = PN,m,n,I v(x) ˜ + vm,n,b,I (x).
(5.3)
By using the results (2.21) and (2.24) of [19], we reach the following result. μ Proposition 5.1 If v ∈ Hm,n,A (I ), ∂xr v ∈ L2χ (−m+r,−n+r) (I ), integers m, n, r ≥ 1, N ≥ m + n, 0 ≤ k ≤ r ≤ N + 1, max(m, n, k) ≤ μ ≤ m + n, and r ≥ m + n (or max(m, n) ≤ r ≤ m + n − 1, 1 ≤ m, n ≤ 4), then
k μ ∂ P x
N,m,n,I v
− v χ (−m+k,−n+k) ,I ≤ cN k−r ∂xr v χ (−m+r,−n+r) ,I .
(5.4)
J Sci Comput (2012) 53:451–480
467
Now, let = {x | −2 < x < ∞}. We denote the inner product and norm of L2χ () by (· , ·)χ and · χ , respectively. We omit χ in notations, whenever χ ≡ 1. Let W (), W () H 2 () and λ ≥ 0. We introduce the bilinear form Aλ (u, v) = ∂x2 u, ∂x2 v + λ(u, v), ∀u ∈ W (), v ∈ W (). Next, we set N = (N1 , N2 ), and ˜ (0,β) WN () = φ ∈ W ()|φ|I ∈ PN1 (I ), φ| ∈ Q N2 () , ˜ (0,β) W N () = φ ∈ W ()|φ|I ∈ PN1 (I ), φ| ∈ Q N2 () . 2 We define the operator ∗ PN,λ : W () → WN (), by
2 Aλ ∗ PN,λ v − v, φ = 0,
∀φ ∈ W N ().
(5.5)
2 Proposition 5.2 If v ∈ W (), w ∈ WN () and ∗ PN,λ v − w ∈ W N (), then
2 2 Aλ v − ∗ PN,λ v, v − ∗ PN,λ v ≤ Aλ (v − w, v − w).
(5.6)
Proof A direct calculation shows 2 2 2 2 Aλ (v − w, v − w) = Aλ v − ∗ PN,λ v, v − ∗ PN,λ v + Aλ ∗ PN,λ v − w, ∗ PN,λ v−w 2 2 + 2Aλ v − ∗ PN,λ v, ∗ PN,λ v−w . Thanks to (5.5), we have
2 2 Aλ v − ∗ PN,λ v, ∗ PN,λ v − w = 0.
Since λ ≥ 0, we have
2 2 Aλ ∗ PN,λ v − w, ∗ PN,λ v − w ≥ 0.
Then, the desired result (5.6) follows from the previous statements immediately.
5.2 Dirichlet Problem of High Order Equation We first consider the following simple model problem, ∂x4 U (x) + λU (x) = f (x), x ∈ , ∂x U (−2) = b, U (−2) = a,
(5.7)
where f ∈ L2 (), a, b and λ > 0 are given constants. In addition, we require that 1 1 x 2 ∂xk U (x) → 0, as x → ∞, 0 ≤ k ≤ 2, and x − 2 ∂x3 U (x) → 0, as x → ∞. Let 1 H˜ 2 () = H 2 () ∩ v | x − 2 ∂x3 v(x) → 0, as x → ∞ , V () = 0 H 2 (). V () = v | v ∈ H˜ 2 () and ∂x v(−2) = b, v(−2) = a , Multiplying (5.7) by v ∈ V () and integrating the resulting equation over , we obtain the weak formulation of (5.7). It is to look for U ∈ V () such that
468
J Sci Comput (2012) 53:451–480
Aλ (U, v) = (f, v),
∀v ∈ V ().
(5.8)
It can be shown that problem (5.8) admits a unique solution. We now construct the domain decomposition spectral scheme for (5.8). Let ˜ (0,β) VN () = φ ∈ V ()|φ|I ∈ PN1 (I ), φ| ∈ Q N2 () , ˜ (0,β) V N () = φ ∈ V ()|φ|I ∈ PN1 (I ), φ| ∈ Q N2 () . The domain decomposition spectral method for (5.8), is to seek uN ∈ VN () such that Aλ (uN , φ) = (f, φ),
∀φ ∈ V N ().
(5.9)
For the existence of solutions of (5.9), it suffices to check the uniqueness of solution. Assume (2) (2) that u(1) ˜ N (x) = u(1) N (x) and uN (x) are solutions of (5.9), and u N (x) − uN (x) ∈ V N (). Then Aλ (u˜ N , φ) = 0,
∀φ ∈ V N ().
(5.10)
Putting φ = u˜ N in (5.10), we obtain 2 2 ∂ u˜ N + λu˜ N 2 = Aλ (u˜ N , u˜ N ) = 0. x Thus, u˜ N (x) ≡ 0. This implies the existence and the uniqueness of solution of (5.9). For numerical analysis of the proposed method, we introduce the composite quasiorthogonal projection PN v(x), defined by PN21 ,2,2,I v(x), x ∈ I, PN v(x) = (5.11) P˜N22 ,−2,β, v(x), x ∈ . Obviously, PN U ∈ VN (). By using (5.4), we deduce that for integer r1 ≥ 2, k 2 2 ∂ U − P 2 ≤ c ∂ k U − P 2 (−2+k,−2+k) x N1 ,2,2,I U x N1 ,2,2,I U I χ ,I 2 ≤ cN1 2(k−r1 ) ∂xr1 U χ (−2+r1 ,−2+r1 ) ,I , k = 0, 1, 2.
(5.12)
On the other hand, according to (3.14) and (3.16), we set
β β U˜ (x) = U (x) − Ub,2, (x), Ub,2, (x) = e− 2 x U (0) + ∂x U (0) + U (0) x , 2 β β P˜N22 ,−2,β, U (x) = e− 2 x 0 PN22 ,−2,β, e 2 x U˜ (x) + Ub,2, (x). We use (2.16) and (2.10) to derive that for integer r2 ≥ 2, k 2 β x 2 ∂ U − P˜ 2 = ∂ k e− β2 x e β2 x U˜ − 0 P 2 ˜ 2 U x N2 ,−2,β, U x N2 ,−2,β, e β β 2 ≤ cβ 2k−4 ∂x2 e 2 x U˜ − 0 PN22 ,−2,β, e 2 x U˜ ω
0,β ,
β 2 ≤ cβ 2k−4 (βN2 )2−r2 ∂xr2 e 2 x U˜ ω β 2 ≤ cβ 2k−4 (βN2 )2−r2 ∂xr2 e 2 x U ω
−2+r2 ,β ,
−2+r2 ,β ,
,
k = 0, 1, 2. (5.13)
J Sci Comput (2012) 53:451–480
469
We now estimate the error of numerical solution. For this purpose, we introduce the 2 auxiliary projection P N,λ : V () → VN (), defined by
2
Aλ P N,λ v − v, φ = 0,
∀φ ∈ V N ().
(5.14)
∀φ ∈ V N ().
(5.15)
∀φ ∈ V N ().
(5.16)
We have from (5.8) and (5.14) that
2
Aλ P N,λ U, φ = (f, φ),
Subtracting (5.15) from (5.9), yields
2
Aλ uN − P N,λ U, φ = 0, 2
Taking φ = uN − P N,λ U ∈ V N () in (5.16), we obtain
2
2
Aλ uN − P N,λ U, uN − P N,λ U = 0. 2
This implies uN = P N,λ U . Let r = (r1 , r2 ). For simplicity, we shall use the following notations, 2 1 2 A(v, k) = ∂xk v χ (−2+k,−2+k) ,I + ∂xk v 2 , 2 B(v, N, r) = N1 4−2r1 ∂ r1 U (−2+r ,−2+r ) x
+ 1+β
−4
1
χ
1 ,I
β 2 (βN2 )2−r2 ∂xr2 e 2 x U ω
12 −2+r2 ,β ,
.
Furthermore, we put in Proposition 5.2 that W () = V (), v = U,
W () = V (), w = PN U,
2 ∗ PN,λ U
WN () = VN (),
W N () = V N (),
2
= P N,λ U.
Then, we use the Gagliardo-Nirenberg inequality, projection theorem with (5.14), (5.12) and (5.13) successively, to obtain that for integers r1 , r2 ≥ 2, U − uN H 2 () ≤ c ∂x2 (U − uN ) + λU − uN 2 2 = c ∂x2 U − P N,λ U + λ U − P N,λ U ≤ c ∂x2 (U − PN U ) + λU − PN U ≤ c A(U − PN U, 2) + λA(U − PN U, 0) ≤ cB(U, N, r).
(5.17)
5.3 Mixed Inhomogeneous Boundary Value Problem We now consider the following mixed inhomogeneous boundary value problem, ∂x4 U (x) + λU (x) = f (x), x ∈ , U (−2) = a, ∂x2 U (−2) = b,
(5.18)
470
J Sci Comput (2012) 53:451–480 1
where f ∈ L2 (), a, b and λ > 0 are given constants. We also require that x 2 ∂xk U (x) → 0, 1 as x → ∞, 0 ≤ k ≤ 2, and x − 2 ∂x3 U (x) → 0, as x → ∞. Let V () = v | v ∈ H˜ 2 () and v(−2) = a , V () = v | v ∈ H 2 () and v(−2) = 0 . The weak formulation of (5.18) is to look for U ∈ V () such that Aλ (U, v) + b∂x v(−2) = (f, v),
∀v ∈ V ().
(5.19)
It can be shown that problem (5.19) admits a unique solution. Next, let ˜ (0,β) VN () = φ ∈ V () | φ|I ∈ PN (I ), φ| ∈ Q () , N 1
2
˜ (0,β) V N () = φ ∈ V () | φ|I ∈ PN1 (I ), φ| ∈ Q N2 () .
The domain decomposition spectral method for (5.19), is to seek uN ∈ VN () such that Aλ (uN , φ) + b∂x φ(−2) = (f, φ),
∀φ ∈ V N ().
(5.20)
Like (5.9), the above problem has a unique solution. We now estimate the error of numerical solution. Let PN v(x) be the same as in (5.11). 2 We also introduce the auxiliary operator P N,λ : V () → VN (), defined by
2
Aλ P N,λ v − v, φ = 0,
∀φ ∈ V N ().
(5.21)
We have from (5.19) and (5.21) that
2
Aλ P N,λ U, φ + b∂x φ(−2) = (f, φ),
∀φ ∈ V N ().
(5.22)
Subtracting (5.22) from (5.20), yields
2
Aλ uN − P N,λ U, φ = 0,
∀φ ∈ V N ().
(5.23)
2
Taking φ = uN − P N,λ U ∈ V N () in (5.23), we reach
2
2
Aλ uN − P N,λ U, uN − P N,λ U = 0. 2
Thereby, uN = P N,λ U . Finally, by virtue of the Gagliardo-Nirenberg inequality, projection theorem with (5.21), (5.12), (5.13) and Proposition 5.2 with W () = V (), v = U,
W () = V (), w = PN U,
2 ∗ PN,λ U
we obtain that for integers r1 , r2 ≥ 2,
WN () = VN (), =
2 P N,λ U,
W N () = V N (),
J Sci Comput (2012) 53:451–480
471
U − uN H 2 () ≤ c ∂x2 (U − uN ) + λU − uN = c ∂x2 (U − P N,λ U ) + λU − P N,λ U ≤ c ∂x2 (U − PN U ) + λU − PN U ≤ c A(U − PN U, 2) + λA(U − PN U, 0) ≤ cB(U, N, r).
(5.24)
Remark 5.1 We may use spectral method without domain decomposition to solve (5.18). In this case, we let (α,β) () = span Ll (x + 2), l¯α ≤ l ≤ N , β (α,β) (α,β) Q˜ N () = span e− 2 (x+2) φ, φ ∈ QN () , (α,β)
QN
VN () = V () ∩ Q˜ N
(0,β)
(),
V N () = V () ∩ Q˜ N
(0,β)
().
The spectral method for (5.18) is to find uN ∈ VN () such that Aλ (uN , φ) + b∂x φ(−2) = (f, φ),
∀φ ∈ V N ().
(5.25)
∗ (x) = ωα,β (x + 2). By an argument similar to the derivation of (5.24), we obtain Let ωα,β
β 2 U − uN 2H 2 () ≤ c 1 + β −4 (βN )2−r ∂xr e 2 x U ω∗
−2+r,β
.
5.4 Spectral Method with Essential Imposition of Boundary Conditions The Laguerre quasi-orthogonal approximation also plays an important role in spectral method with essential imposition of mixed inhomogeneous boundary conditions. To do this, we set V () = v | v ∈ H˜ 2 () and v(−2) = a ,
V () = v | v ∈ H 2 () and v(−2) = 0 .
The weak formulation of (5.18) is to look for U ∈ V () such that Aλ (U, v) + b∂x v(−2) = (f, v),
∀v ∈ V ().
(5.26)
This problem admits a unique solution. Let () ∩ v | ∂x2 v(−2) = b , (0,β) V N () = V () ∩ Q˜ N () ∩ v | ∂x2 v(−2) = 0 .
˜N VN () = V () ∩ Q
(0,β)
The spectral method for (5.26), is to seek uN ∈ VN () such that Aλ (uN , φ) + b∂x φ(−2) = (f, φ),
∀φ ∈ V N ().
Like (5.9) and (5.20), problem (5.27) possess a unique solution.
(5.27)
472
J Sci Comput (2012) 53:451–480 2
We now deal with the convergence of (5.27). We introduce the projection P N,λ : V () → VN (), defined by 2 Aλ P N,λ v − v, φ = 0, ∀φ ∈ V N (). (5.28) We have from (5.26) and (5.28) that
2
Aλ P N,λ U, φ + b∂x φ(−2) = (f, φ),
∀φ ∈ V N ().
(5.29)
Subtracting (5.29) from (5.27), yields
2
Aλ uN − P N,λ U, φ = 0,
∀φ ∈ V N ().
(5.30)
2
Taking φ = uN − P N,λ U ∈ V N () in (5.30), we obtain
2
2
Aλ uN − P N,λ U, uN − P N,λ U = 0. 2
This implies uN = P N,λ U . For error estimation of numerical solution, we need the Laguerre quasi-orthogonal projection defined on . According to (3.14) and (3.16), we introduce the function U˜ (x) = U (x) − Ub,3, (x), with
β β Ub,3, (x) = e− 2 (x+2) U (−2) + ∂x U (−2) + U (−2) (x + 2) 2
1 2 β β2 2 + ∂ U (−2) + ∂x U (−2) + U (−2) (x + 2) . 2 x 2 8 ∗ (x) be the same as in Remark 5.1. We define the orthogonal projecLet QN () and ωα,β 3 tion 0 PN,−3,β, by (α,β)
3 3 ∂x 0 PN,−3,β, v − v , ∂x3 v ω∗ = 0, 0,β
(−3,β)
∀φ ∈ QN
().
Like (3.16), the Laguerre quasi-orthogonal projection defined on , is defined by β (x+2) β 3 3 U (x) = e− 2 (x+2) 0 PN,−3,β, U˜ (x) + Ub,3, (x). e2 P˜N,−3,β, Thanks to Lemma 2.2 of [16], for any v ∈ 0 Hω1∗ (), 0,β
v2ω∗ ≤ 0,β
4 ∂x v2ω∗ . 0,β β2
(5.31)
Accordingly, we use an equality like (2.16), (5.31) and an inequality like (2.10) successively to derive that for integer r ≥ 3, k 2 k − β (x+2) β (x+2) β (x+2) 2 3 ∂ U − P˜ 3 = ∂ e 2 U˜ − 0 PN,−3,β, U˜ e2 e2 x N,−3,β, U x β β (x+2) 2 3 ≤ cβ 2k−4 ∂x2 e 2 (x+2) U˜ − 0 PN,−3,β, U˜ ω∗ e2 0,β
≤
cβ 2k−6 ∂x3
e
β 2 (x+2)
U˜
3 − 0 PN,−3,β,
e
β 2 (x+2)
2 U˜ ∗
ω0,β
J Sci Comput (2012) 53:451–480
473
β 2 ≤ cβ 2k−6 (βN )3−r ∂xr e 2 (x+2) U˜ ω∗
−3+r,β
β 2 ≤ cβ 2k−6 (βN )3−r ∂xr e 2 (x+2) U ω∗
−3+r,β
,
k = 0, 1, 2. (5.32)
On the other hand, like Proposition 5.2, we can prove that
3 3 Aλ U − P N,λ U, U − P N,λ U ≤ Aλ U − P˜N,−3,β, U, U − P˜N,−3,β, U . 2
2
Finally, we use the Gagliardo-Nirenberg inequality, projection theorem with (5.28), and (5.32) successively, to obtain that for integer r ≥ 3, 2 U − uN 2H 2 () ≤ c ∂x2 (U − uN ) + λU − uN 2 2 = c ∂x2 (U − P N,λ U ) + λU − P N,λ U 2 2 2 3 3 ≤ c ∂x2 U − P˜N,−3,β, U + λ U − P˜N,−3,β, U
β 2 c λ . (5.33) ≤ 2 1 + 4 (βN )3−r ∂xr e 2 x U ω∗ −3+r,β β β For other spectral and collocation methods with essential imposition of boundary conditions, we refer to the work of [1, 22]. 6 Numerical Results In this section, we present some numerical results to illustrate the efficiency of proposed methods. 6.1 Spectral Scheme (5.9) We need six kinds of basis functions. The first kind of basis functions are local, and correspond to subdomains I and respectively, namely, 1 Y (2,2) (x), x ∈ I, (6.1) ψ1,l (x) = 4(l−2)(l−3) l 0, otherwise, 2 (−2,β) β L˜ (x), x ∈ , ψ2,l (x) = l(l−1) l (6.2) 0, otherwise. Next, we define the basis function for matching numerical solution at x = 0, as (x), x ∈ I, q+ σ (x) = 2,2,0,Iβ β (1 + 2 x)e− 2 x , x ∈ .
(6.3)
Similarly, we define the basis function for matching the first derivative of numerical solution at x = 0, as (x), x ∈ I , q+ (6.4) ρ(x) = 2,2,1,I β x ∈ . xe− 2 x ,
474
J Sci Comput (2012) 53:451–480
Obviously, σ (0) = ∂x ρ(0) = 1 and ∂x σ (0) = ρ(0) = 0. Furthermore, we take the basis function for matching numerical solution at x = −2, as (x), x ∈ I , q− (6.5) ζ (x) = 2,2,0,I 0, otherwise. We also take the basis function for matching the first order derivative of numerical solution at x = −2, as (x), x ∈ I , q− (6.6) ξ(x) = 2,2,1,I 0, otherwise. Clearly, ζ (−2) = ∂x ξ(−2) = 1 and ∂x ζ (−2) = ξ(−2) = 0. In actual computation, we expand the numerical solution as uN (x) =
N1
uˆ 1,l ψ1,l (x) +
l=4
N2
uˆ 2,l ψ2,l (x) + uN (0)σ (x) + ∂x uN (0)ρ(x)
l=2
+ aζ (x) + bξ(x) ∈ VN ().
(6.7)
Let g(φ) = (f, φ). Inserting (6.7) into the domain decomposition spectral scheme (5.9), we obtain Aλ (uN , φ) = g(φ),
∀φ ∈ V N ().
(6.8)
where the test functions φ stand for the basis functions given by (6.1)–(6.5). We now use algorithm (6.8) to solve problem (5.7) with λ = 1. Let xLe,N1 ,j and ωLe,N1 ,j be the nodes and weights of the Legendre-Gauss-Lobatto numerical quadrature on the interval (−1, 1), respectively. Meanwhile, xLa,N2 ,j and ωLa,N2 ,j stand for the nodes and weights of the Leguerre-Gauss-Radau numerical quadrature on the interval (0, ∞), respectively. We measure the numerical errors by the following quantities:
N 1 2 U (xLe,N1 ,j − 1) − uN (xLe,N1 ,j − 1) ωLe,N1 ,j U − uN N,ave = j =0
+
N2
2 U (xLa,N2 ,j ) − uN (xLa,N2 ,j ) ωLa,N2 ,j
12 ,
j =1
U − uN N,max = max max U (xLe,N1 ,j − 1) − uN (xLe,N1 ,j − 1), 0≤j ≤N1
max U (xLa,N2 ,j ) − uN (xLa,N2 ,j ) .
1≤j ≤N2
We first take the test function U (x) = e−θx sin γ x,
θ, γ > 0,
(6.9)
which decays exponentially as x goes to infinity. In Tables 1 and 2, we list the values of U − uN N,ave and U − uN N,max , with θ = γ = 1, N2 = 4N1 and β = 1, 2, v.s. the modes N1 . The numerical results demonstrate that the numerical errors decay fast as N1 increases. This confirms the theoretical analysis in Sect. 5.2.
J Sci Comput (2012) 53:451–480 Table 1 The values U − uN N,ave of algorithm (6.8)
Table 2 The values U − uN N,max of algorithm (6.8)
475 N1 = 5
N1 = 10
N1 = 15
β =1
3.31E–3
4.85E–9
1.81E–13
β =2
2.86E–3
4.84E–9
6.88E–14
N1 = 5
N1 = 10
N1 = 15
β =1
1.81E–3
2.62E–9
3.51E–14
β =2
1.95E–3
2.62E–9
2.35E–14
We also see that a suitable choice of parameter β could raise the numerical accuracy. However, so far, there is no theoretical result on the best choice of β. Generally speaking, it seems reasonable to take certain β so that the asymptotic behaviors of numerical solutions are more similar to those of exact solutions. Clearly, the test function (6.9) with θ = 1 decays like e−x , while the numerical solution with β = 2 does so also. Thus it is more accurate. For a given problem, the better choice of β depends on some parameters involved in the underlying differential equation, which affect the asymptotic behavior of solution. For instance, we consider the following simple equation, ∂x4 U (x) + λU (x) = f (x),
λ > 0, 0 < x < ∞.
In addition, ∂xm U (x) → 0 as x → ∞, m = 1, 2. If f (x) fulfills certain condition, then the √ √ 4 solution U (x) decays like e− λx as x → ∞. In this case, we may take β = 2 4 λ. On the other hand, if the solution varies faster near the fixed boundary, then it is better to take suitably bigger β, since the distances between two adjacent roots of the scaled Laguerre functions with bigger β are small than those with small β. We refer to [17, 18] for more discussions. We next use (6.8) to solve (5.7) with the test function U (x) =
sin x + 1 , (x + 3)k
k ≥ 1,
(6.10)
which decays algebraically as x goes to infinity. In Tables 3 and 4, we list the values of U − uN N,ave and U − uN N,max , with k = 3, N2 = 4N1 and β = 1, 2, v.s. the mode N1 . The numerical results show again the high accuracy of algorithm (6.8), and the affecting of parameter β. In Table 5, we present the values of U − uN N,ave and U − uN N,max , with k = 3, 4 and β = 2, v.s. the mode N1 . The numerical results indicate that the smoother the exact solution, the better the numerical results. This coincides with the theoretical analysis. In fact, for the test function U (x) given by (6.10), β the norm ∂xr (e 2 x U )2ω−2+r,β , is finite, as long as r < 2k + 1. Thereby, the estimate (5.17) implies U − uN H 2 () β ≤ c N1 2−2k ∂x2k U χ (−2+2k,−2+2k) ,I + 1 + β −2 (βN2 )1−k ∂x2k e 2 x U ω
−2+2k,β ,
.
476
J Sci Comput (2012) 53:451–480
Table 3 The values U − uN N,ave of algorithm (6.8)
Table 4 The values U − uN N,max of algorithm (6.8)
N1 = 5
N1 = 10
N1 = 15
N1 = 20
N1 = 25
β =1
1.74E–2
6.70E–5
5.89E–6
2.74E–6
1.48E–6
β =2
1.80E–2
6.55E–5
2.90E–6
1.50E–6
8.42E–7
N1 = 5
N1 = 10
N1 = 15
N1 = 20
N1 = 25
β =1
1.12E–2
3.63E–5
2.57E–6
1.20E–6
6.34E–7
β =2
1.12E–2
3.63E–5
1.53E–6
6.50E–7
3.48E–7
Table 5 The errors of algorithm (6.8) k=3
k=4
U − uN N,ave
U − uN N,max
U − uN N,ave
U − uN N,max
N1 = 15
2.90E–6
1.53E–6
7.41E–7
3.32E–7
N1 = 20
1.50E–6
6.50E–7
1.40E–8
6.35E–9
N1 = 25
8.42E–7
3.48E–7
6.36E–9
2.65E–9
6.2 Spectral scheme (5.20) We need the same six kinds of basis functions as in (6.1)–(6.6). In actual computation, we expand the numerical solution as uN (x) =
N1 l=4
uˆ 1,l ψ1,l (x) +
N2
uˆ 2,l ψ2,l (x) + uN (0)σ (x) + ∂x uN (0)ρ(x)
l=2
+ aζ (x) + ∂x uN (−2)ξ(x) ∈ VN ().
(6.11)
Let g(φ) = (f, φ) − b∂x φ(−2). Inserting (6.11) into (5.20), we obtain Aλ (uN , φ) = g(φ),
∀φ ∈ V N ().
(6.12)
where the test functions φ stand for the basis functions given by (6.1)–(6.6). We now use algorithm (6.12) to solve the mixed inhomogeneous boundary value problem (5.18) with λ = 1. We first take the test function (6.9). In Tables 6 and 7, we list the values of U − uN N,ave and U − uN N,max , with σ = γ = 1, N2 = 4N1 and β = 1, 2, v.s. the modes N1 . Evidently, the errors decay very fast as N1 increases. Also, a suitable choice of parameter β leads to more accurate numerical results. We next use (6.12) to solve (5.18) with the test function (6.10). In Tables 8 and 9, we list the values of U − uN N,ave and U − uN N,max , with k = 3, N2 = 4N1 and β = 1, 2, v.s. the mode N1 . In Table 10, we present the values of U − uN N,ave and U − uN N,max , with k = 3, 4 and β = 2, v.s. the mode N1 . This coincides with the theoretical analysis again. For comparison, we consider problem (5.18) with the test function (6.9) and θ = 4, γ = 2. Clearly, the solution oscillates seriously near the endpoint x = −2. We use (6.12) with N2 = 4N1 and (5.25) with N = 5N1 , to solve the same problem. The numerical errors of
J Sci Comput (2012) 53:451–480 Table 6 The values U − uN N,ave of algorithm (6.12)
Table 7 The values U − uN N,max of algorithm (6.12)
Table 8 The values U − uN N,ave of algorithm (6.8)
Table 9 The values U − uN N,max of algorithm (6.12)
477 N1 = 5
N1 = 10
N1 = 15
β =1
3.65E–3
4.85E–9
1.86E–13
β =2
2.86E–3
4.84E–9
7.55E–14
N1 = 5
N1 = 10
N1 = 15
β =1
1.84E–3
2.62E–9
3.51E–14
β =2
1.82E–3
2.62E–9
2.43E–14
N1 = 5
N1 = 10
N1 = 15
N1 = 20
N1 = 25
β =1
6.45E–2
6.74E–5
5.89E–6
2.74E–6
1.48E–6
β =2
6.76E–2
6.58E–5
2.90E–6
1.49E–6
8.43E–7
N1 = 5
N1 = 10
N1 = 15
N1 = 20
N1 = 25
β =1
4.00E–2
3.72E–5
2.57E–6
1.20E–6
6.34E–7
β =2
4.00E–2
3.72E–5
1.53E–6
6.50E–7
3.49E–7
Table 10 The errors of algorithm (6.12) k=3
k=4
U − uN N,ave
U − uN N,max
U − uN N,ave
U − uN N,max
N1 = 15
2.90E–6
1.53E–6
7.41E–7
3.31E–7
N1 = 20
1.49E–6
6.50E–7
1.40E–8
6.35E–9
N1 = 25
8.43E–7
3.49E–7
6.36E–9
2.65E–9
(5.25) are measured by
U − uN N,ave,∗ =
N
U (xLa,N,j − 2) − uN (xLa,N,j
2 − 2) ωLa,N,j
12 ,
j =0
U − uN N,max,∗ = max U (xLa,N,j − 2) − uN (xLa,N,j − 2). 0≤j ≤N
In Table 11, we list the values of U − uN N,ave and U − uN N,max of algorithm (6.12), and the values of U − uN N,ave,∗ and U − uN N,max,∗ of algorithm (5.25). The parameter β = 2. Obviously, the spectral method with domain decomposition provides better numerical results.
478
J Sci Comput (2012) 53:451–480
Table 11 The errors of (6.12) with N2 = 4N1 and (5.25) with N = 5N1 Algorithm (6.12)
Algorithm (5.25)
U − uN N,ave
U − uN N,max
U − uN N,ave,∗
U − uN N,max,∗
N1 = 20
3.88E–10
1.47E–10
6.33E–9
2.18E–9
N1 = 25
7.64E–11
1.65E–11
1.53E–8
5.30E–9
6.3 Spectral scheme (5.27) We define the basis functions as follows,
β β2 β ψ0 (x) = 1 + (x + 2) + (x + 2)2 e− 2 (x+2) , 2 8
β β ψ1 (x) = x + 2 + (x + 2)2 e− 2 (x+2) , 2
(6.13)
β 1 ψ2 (x) = (x + 2)2 e− 2 (x+2) , 2
ψl (x) =
β3 (−3,β) L˜ (x + 2), l(l − 1)(l − 2) l
3 ≤ l ≤ N.
Evidently, ∂xk ψl (−2) = δk,l for 0 ≤ k, l ≤ 2, and ∂xk ψl (−2) = 0 for 0 ≤ k ≤ 2, 3 ≤ l ≤ N . In actual computation, we expand the numerical solution as uN (x) =
N
uˆ N,l ψl (x) + bψ2 (x) + ∂x uN (−2)ψ1 (x) + aψ0 (x) ∈ VN ().
l=3
Let g(φ) be the same as in (6.12). Then, the scheme (5.27) with essential imposition of boundary condition, becomes Aλ (uN , φ) = g(φ),
∀φ ∈ V N ().
(6.14)
Here we take the test functions φ, as presented in (6.13). We now use (6.14) to solve (5.18) with λ = 1 and the test function (6.9) with θ = γ = 1. In Table 12, we present the values of U − uN N,ave,∗ and U − uN N,max,∗ of algorithms (5.25) and (6.14), with β = 2, v.s. the mode N . We find that like [1], the algorithm (6.14) with essential imposition of boundary condition, provides better numerical results, than scheme (5.25) without domain decomposition and essential imposition of boundary condition. But, the error estimate (5.33) is not as good as the estimate for (5.25), see Remark 5.1. How to improve the estimate (5.33) is still an open problem.
7 Concluding Remarks In this paper, we developed the generalized Laguerre quasi-orthogonal approximation by using the new basis functions, which are mutually orthogonal with the weight function x α , α being any real number. This new approximation matches the considered function at the origin properly, and fits the asymptotic behaviors (growing or decaying algebraically at the
J Sci Comput (2012) 53:451–480
479
Table 12 The errors of algorithms (5.25) and (6.14) Algorithm (5.25)
Algorithm (6.14)
U − uN N,ave,∗
U − uN N,max,∗
U − uN N,ave,∗
U − uN N,max,∗
N = 100
1.44E–12
5.80E–13
1.36E–12
5.07E–13
N = 125
2.65E–11
9.08E–12
1.15E–12
3.89E–13
infinity) reasonably. By using this approximation, it is more convenient to deal with various mixed inhomogeneous boundary value problems of high order differential equations. Moreover, it plays an important role in the domain decomposition spectral method of high order problems defined on unbounded domains. We also studied the corresponding LaguerreGauss-Radau type interpolation, which serves as an important tool for the new Laguerre pseudospectral and collocation methods. As examples of applications, we provided the domain decomposition spectral methods for two high order problems with mixed inhomogeneous boundary conditions. The numerical results demonstrated their high accuracy. We also developed the spectral method with essential imposition of boundary condition, which often leads to better numerical results. Although we only considered several model problems, we could use the techniques developed in this work, coupled with other tricks, to treat other problems, such as the stream function form of the Navier Stokes equations in an infinite strap or outside certain obstacles, the infinite beam equations and so on.
References 1. Auteri, F., Parolini, N., Quartapelle, L.: Essential imposition of Neumann condition in Galerkin-Legendre elliptic solver. J. Comput. Phys. 185, 427–444 (2003) 2. Azaiez, M., Shen, J., Xu, C., Zhuang, Q.: A Laguerre-Legendre spectral method for the Stokes problems in a semi-infinite channel. SIAM J. Numer. Anal. 47, 271–292 (2008) 3. Bernardi, C., Maday, Y.: Spectral methods. In: Ciarlet, P.G., Lions, J.L. (eds.) Handbook of Numerical Analysis. Techniques of Scientific Computing, vol. 5 pp. 209–486. Amsterdam, Elsevier (1997) 4. Bernardi, C., Coppoletta, G., Maday, Y.: Some spectral approximations of multi-dimensional fourthorder problems. Internal report 90021, Laboratoire d’Analyse Numérique, Université Pierre et Marie Curie, Paris (1990) 5. Boyd, J.P.: Chebyshev and Fourier Spectral Methods, 2nd edn. Dover, New York (2001) 6. Canuto, C., Hussaini, M.Y., Quarteroni, A., Zang, T.A.: Spectral Methods, Fundamentals in Single Domains. Springer, Berlin (2006) 7. Funaro, D.: Polynomial Approximations of Differential Equations. Springer, Berlin (1992) 8. Funaro, D., Kavian, O.: Approximation of some diffusion evolution equations in unbounded domains by Hermite function. Math. Comput. 57, 597–619 (1999) 9. Gottlieb, D., Orszag, S.A.: Numerical Analysis of Spectral Methods: Theory and Applications. SIAM/CBMS, Philadelphia (1977) 10. Guo, B.-y.: Spectral Methods and Their Applications. World Scientific, Singapore (1998) 11. Guo, B.-y.: Error estimation of Hermite spectral method for nonlinear partial differential equations. Math. Comput. 68, 1067–1078 (1999) 12. Guo, B.-y., Ma, H.-p.: Composite Legendre-Laguerre approximation in unbounded domains. J. Comput. Math. 19, 101–112 (2001) 13. Guo, B.-y., Shen, J.: Laguerre-Galerkin method for nonlinear partial differential equations on a semiinfinite interval. Numer. Math. 86, 635–654 (2000) 14. Guo, B.-y., Wang, T.-j.: Composite generalized Laguerre-Legendre spectral method with domain decomposition and its application to Fokker-Planck equation in an infinite channel. Math. Comput. 78, 129–151 (2009) 15. Guo, B.-y., Xu, C.-l.: Mixed Laguerre-Legendre pseudospectral method for incompressible fluid flow in an infinite strip. Math. Comput. 72, 95–125 (2003)
480
J Sci Comput (2012) 53:451–480
16. Guo, B.-y., Zhang, X.-y.: A new generalized Laguerre spectral approximation and its applications. J. Comput. Appl. Math. 181, 342–363 (2005) 17. Guo, B.-y., Zhang, X.-y.: Spectral method for differential equations of degenerate type by using generalized Laguerre functions. Appl. Numer. Math. 57, 455–471 (2007) 18. Guo, B.-y., Wang, L.-l., Wang, Z.-q.: Generalized Laguerre interpolation and pseudospectral method for unbounded domains. SIAM J. Numer. Anal. 43, 2567–2589 (2006) 19. Guo, B.-y., Sun, T., Zhang, C.: Jacobi and Laguerre quasi-orthogonal approximations and related interpolations. Math. Comput. (in press) 20. Karniadakis, G.E., Sherwin, S.J.: Spectral/hp Element Methods for CFD, 2nd edn. Oxford Univ. Press, Oxford (2005) 21. Shen, J.: Stable and efficient spectral methods in unbounded domains using Laguerre functions. SIAM J. Numer. Anal. 38, 1113–1133 (2000) 22. Wang, L.-l., Guo, B.-y.: Interpolation approximations based on Gauss-Lobatto-Legendre-Birkhoff quadrature. J. Approx. Theory 161, 142–173 (2009) 23. Xu, C.-l., Guo, B.-y.: Laguerre pseudospectral method for nonlinear partial differential equation. J. Comput. Math. 20, 413–428 (2002)
Domain Decomposition Spectral Method for Mixed Inhomogeneous Boundary Value Problems of High Order Differential Equations on Unbounded Domains Chao Zhang · Ben-yu Guo
Received: 30 December 2010 / Revised: 6 February 2012 / Accepted: 12 February 2012 / Published online: 28 February 2012 © Springer Science+Business Media, LLC 2012
Abstract In this paper, we develop domain decomposition spectral method for mixed inhomogeneous boundary value problems of high order differential equations defined on unbounded domains. We introduce an orthogonal family of new generalized Laguerre functions, with the weight function x α , α being any real number. The corresponding quasiorthogonal approximation and Gauss-Radau type interpolation are investigated, which play important roles in the related spectral and collocation methods. As examples of applications, we propose the domain decomposition spectral methods for two fourth order problems, and the spectral method with essential imposition of boundary conditions. The spectral accuracy is proved. Numerical results demonstrate the effectiveness of suggested algorithms. Keywords Laguerre quasi-orthogonal approximation · Laguerre-Gauss-Radau type interpolation · Domain decomposition spectral methods · High order problems with mixed inhomogeneous boundary conditions
The work of C. Zhang is supported in part by NSF of China No. 11171227 and Research Fund for young teachers of Jiangsu Normal University No. 11XLR27. The work of B.-y. Guo is supported in part by NSF of China No. 11171227, Fund for Doctor Authority of Chinese Educational Ministry No. 20080270001, Shanghai Leading Academic Discipline Project No. S30405, and Fund for E-institute of Shanghai Universities No. E03004. C. Zhang Jiangsu Normal University, Xuzhou 221116, China e-mail: [email protected] B.-y. Guo () Shanghai Normal University, Shanghai 200234, China e-mail: [email protected] B.-y. Guo Scientific Computing Key Laboratory, Division of Computational Science of E-institute, Shanghai Universities, Shanghai, China
452
J Sci Comput (2012) 53:451–480
1 Introduction Spectral and pseudospectral methods have been used successfully for scientific computation, see [3, 5–7, 9, 10, 20] and the references therein. During the past decade, more and more attentions were paid to spectral and pseudospectral methods for unbounded domains. Among those approaches, four direct and commonly used methods are the Hermite spectral and pseudospectral methods, and the Laguerre spectral and pseudospectral methods, see [8, 11, 13, 15, 21, 23] and the references therein. The existing results on the Laguerre approximation are mostly suitable for second order problems. However, it is also interesting to consider high order problems, cf. [4, 15]. Recently, Guo, Sun and Zhang [19] proposed the generalized Laguerre quasi-orthogonal approximation with the weight function x α e−βx , where α is any real number and β > 0. By adjusting the parameters α and β, it matches considered functions at the origin properly, and simulates the asymptotic behaviors more reasonably. This approach provides a tool for treating high order problems with solutions growing up at infinity. Whereas, in many practical cases, such as some problems arising in fluid dynamics and astrophysics, the solutions do not grow up at infinity. In those cases, the above framework is no longer the most appreciate. On the other hand, the solutions may vary rapidly near the fixed boundaries. For raising numerical accuracy and saving work, we might use the Legendre approximation near the boundaries and the Laguerre approximation on the remaining domains. It leads to domain decomposition spectral method, cf. [2, 12]. Furthermore, for solving certain problems in statistical physics, such as the Fokker-Planck equation in an infinite channel, we also need domain decomposition spectral method, cf. [14]. But the presence of weight x α e−βx brings difficulties in domain decomposition spectral method for high order problems. In this paper, we consider another spectral method on unbounded domains. We shall introduce a new family of generalized Laguerre functions, which are mutually orthogonal with the weight function x α , α being any real number. Then we develop the related orthogonal approximation, which is more suitable for approximated functions decaying or growing algebraically at infinity. We next propose the corresponding quasi-orthogonal approximation, with which we could deal with various mixed inhomogeneous boundary value problems easily. In particular, the Laguerre quasi-orthogonal approximation with negative integer α, serves as a basic tool in domain decomposition spectral method of high order problems, since it keeps certain continuity on all common boundaries of adjacent subdomains. Moreover, such new approximation is applicable to numerical solutions of some exterior problems. We also study the related Laguerre-Gauss-Radau type interpolation, which is the mathematical foundation of the new Laguerre pseudospectral and collocation methods. As examples of applications, we provide the domain decomposition spectral methods for two high order problems, and prove their spectral accuracy. We also propose spectral method with essential imposition of boundary conditions for a model problem, with the error analysis. Numerical results demonstrate the high accuracy of suggested algorithms. It is noted that Shen [21], and Guo and Ma [12] considered the case with α = 0, β = 1, and Guo and Zhang [17] studied the case with α > −1, β > 0. The related results are available for second order problems. Our new results improve and generalize the existing results, and so enlarge the applications of Laguerre approximation essentially. In particular, they are very useful for domain decomposition spectral method of high order problems with inhomogeneous boundary conditions, as well as exterior problems. This paper is organized as follows. The next section is for preliminaries. In Sect. 3, we develop the new Laguerre quasi-orthogonal approximation. In Sect. 4, we study the new Laguerre-Gauss-Radau type interpolation. Section 5 is devoted to domain decomposition
J Sci Comput (2012) 53:451–480
453
spectral methods and spectral method with essential imposition of boundary conditions. In Sect. 6, we present some numerical results. The final section is for concluding remarks.
2 Preliminaries In this section, we recall the recent results on the generalized Laguerre orthogonal approximation and the Laguerre-Gauss-Radau type interpolation. Let = {x | 0 < x < ∞} and χ (x) be certain a weight function. For integer r ≥ 0, we define the weighted Sobolev space Hχr () in the usual way, with the inner product (· , ·)r,χ , , the semi-norm | · |r,χ , and the norm · r,χ , . In particular, the inner product and the norm dk v of L2χ () are denoted by (· , ·)χ , and · χ , , respectively. For simplicity, we denote dx k by ∂xk v. For integer r ≥ 1, r 0 Hχ ()
= v ∈ Hχr () | ∂xk v(0) = 0, 0 ≤ k ≤ r − 1 .
We omit the subscript χ in notations whenever χ (x) ≡ 1. The scaled generalized Laguerre polynomials of degree l ≥ 0 were given in [16], as (α,β)
Ll
(x) =
1 −α βx l l+α −βx x e ∂x x e , l!
α > −1, β > 0, l ≥ 0.
Now, let α be any real number. Denote by [α] the largest integer ≤ α. Let l¯α = [−α] for α ≤ −1, and l¯α = 0 for α > −1. Meanwhile, lα = l − [−α] for α ≤ −1, and lα = l for α > −1. The new generalized Laguerre functions with real parameter α of degree l ≥ 0 were introduced in [19], namely, (α,β) Ll (x)
=
x −α Llα (x), (α,β) Ll (x), (−α,β)
α ≤ −1, l ≥ l¯α = [−α], α > −1, l ≥ l¯α = 0.
(2.1)
They satisfy the following Sturm-Liouville problem (see (4.6) of [19]), (α,β) (α,β) (α,β) ∂x x α+1 e−βx ∂x Ll (x) + λl x α e−βx Ll (x) = 0,
l ≥ l¯α ,
(2.2)
with the eigenvalues (α,β) λl
=
β(lα − α) = β(l − [−α] − α), βlα = βl,
for α ≤ −1, for α > −1.
(2.3)
Next, let ωα,β (x) = x α e−βx , and (α,β)
γl
=
(l + α + 1) , β α+1 l!
α > −1, β > 0, l ≥ 0.
According to (4.8) of [19], we have (α,β)
Ll
(α,β)
(x)Ll
(α,β)
(x)ωα,β (x)dx = ηl
δl,l
(2.4)
454
J Sci Comput (2012) 53:451–480
where δl,l is the Kronecker symbol, and (−α,β) (−α,β) = γl−[−α] , for α ≤ −1, l ≥ l α , γα (α,β) ηl = l(α,β) γl , for α > −1, l ≥ 0. (α,β)
The set of all Ll
(2.5)
(x) is complete in L2ωα,β (). Hence, for any v ∈ L2ωα,β (), we have v(x) =
∞
(α,β)
vˆl
(α,β)
Ll
(2.6)
(x)
l=l¯α
with (α,β) vˆl
=
1
(α,β)
v(x)Ll
(α,β) ηl
(2.7)
(x)ωα,β (x)dx.
Now, let N ≥ l¯α be non-negative integer, and (α,β) (α,β) (α,β) (α,β) QN () = span Ll¯ (x), Ll¯ +1 (x), . . . , LN (x) . α
α
(α,β)
The orthogonal projection PN,α,β, : L2ωα,β () → QN (PN,α,β, v − v, φ)ωα,β , = 0,
() is defined by (α,β)
∀φ ∈ QN
().
Due to Theorem 4.2 of [19], we have the following basic result. Lemma 2.1 If ∂xν v ∈ L2ω−m+ν,β () for ν = 0, k, r, and integers 1 ≤ m ≤ min(r, N ), then for 0 ≤ k ≤ r ≤ N + 1, k k−r ∂ (PN,−m,β, v − v) ≤ c(βN ) 2 ∂xr v ω . (2.8) x ω , , −m+r,β
−m+k,β
Hereafter, we denote by c a generic positive constant independent of any function and N . Next, we introduce the following space with integer r ≥ 0, Hωr −m,β ,A () = v | v is measurable on and vHωr
−1 and β > 0, we (α,β) (α,β) denote by ξG,N,j, (0 ≤ j ≤ N ) the zeros of the polynomial LN+1 (x), which are arranged in (α,β) ascending order. Meanwhile, ωG,N,j, (0 ≤ j ≤ N ) stand for the corresponding Christoffel numbers such that φ(x)ωα,β (x)dx =
N (α,β) (α,β) φ ξG,N,j, ωG,N,j, ,
∀φ ∈ P2N+1 ().
(2.11)
j =0
We know from (2.10) of [18] that (α,β)
ωG,N,j, =
1
(N + α + 2) , (α,β) (α,β) (α,β) α β (N + 2) ξG,N,j, [∂x LN+1 (ξG,N,j, )]2
0 ≤ j ≤ N.
Next, let m ≥ 0, and (m,β)
(m,β)
ξN,j, = ξG,N−m,j, , (m,β) −2m (m,β) (m,β) ωN,j, = ωG,N−m,j, ξG,N−m,j, ,
(2.12)
0 ≤ j ≤ N − m.
Thanks to (5.9) of [19], we have N−m
(m,β) (m,β) (m,β) φ ξN,j, ψ ξN,j, ωN,j, = (φ, ψ)ω−m,β , ,
(−m,β)
∀φ ∈ QN
(−m,β)
(), ψ ∈ QN+1 ().
j =0
(2.13) For any integer r ≥ 0, we denote by C r () the space consisting of all r-times differential functions. Further, for integer r ≥ m − 1, r () = v ∈ C r () | ∂xk v(0) = 0, 0 ≤ k ≤ m − 1 . C0,m,β m−1 () and m ≥ 1, the new Laguerre-Gauss-Radau type interpolation For any v ∈ C0,m,β (−m,β) IN,−m,β, v ∈ QN () is determined uniquely by
(m,β)
(m,β)
IN,−m,β, v ξN,j, = v ξN,j, ,
0 ≤ j ≤ N − m.
(2.14)
Thanks to Lemma 5.2 of [19], we have the following result. Lemma 2.3 If v ∈ 0 Hωm−m,β ,A () ∩ Hωk−m+k ,β (), ∂xr v ∈ L2ω−m+r,β () and integers 1 ≤ m ≤ min(r, N ), then for 0 ≤ k ≤ r ≤ N + 1, k ∂ (IN,−m,β, v − v) x
ω−m+k,β ,
1 1 k+1−r ≤ c β − 2 + 1 (ln N ) 2 (βN ) 2 ∂xr v ω
−m+r,β ,
In the end of this section, we derive an equality, which will be used in Sect. 5.
.
(2.15)
456
J Sci Comput (2012) 53:451–480
Proposition 2.1 If u ∈ 0 Hω20,β (), then 2 2 ∂ u x
2 4 2 − β x 2 ∂ e 2 u + β ∂x e− β2 x u 2 + β e− β2 x u 2 . = x ω0,β , 2 16
(2.16)
β
Proof Let v(x) = e− 2 x u(x). Then
β β2 ∂x2 u(x) = ∂x2 v(x) + β∂x v(x) + v(x) e 2 x . 4 Since u ∈ 0 Hω20,β (), we have v ∈ 0 H 2 (). Therefore,
v(x)∂x v(x)dx =
∂x v(x)∂x2 v(x)dx = 0,
v(x)∂x2 v(x)dx = −∂x v2 .
Consequently, 2 2 ∂ u x ω
2 β2 v(x) dx 4 4 2 β β2 = ∂x2 v + β 2 ∂x v2 + v2 + v(x)∂x2 v(x)dx 16 2
0,β ,
=
∂x2 v(x) + β∂x v(x) +
2 β2 β4 = ∂x2 v + ∂x v2 + v2 . 2 16
The proof is completed.
3 New Laguerre Quasi-Orthogonal Approximation In this section, we develop the new Laguerre quasi-orthogonal approximation. 3.1 New Laguerre Orthogonal Approximation For any real number α and β > 0, the new generalized Laguerre functions are given by L˜ l
(α,β)
β
(x) = e− 2 x Ll
(α,β)
(x),
l ≥ l¯α .
(α,β) Due to (2.2), L˜ l (x) is the lth eigenfunction of the following Sturm-Liouville problem,
β β (α,β) ∂x x α+1 e−βx ∂x e 2 x v(x) + λl x α e− 2 x v(x) = 0,
l ≥ l¯α .
(3.1)
(α,β) Let ω˜ α (x) = x α . By virtue of (2.4), the set of all L˜ l (x) conforms a complete L2ω˜ α ()orthogonal system, namely, (α,β) (α,β) (α,β) L˜ l (x)L˜ l (x)ω˜ α (x)dx = ηl δl,l . (3.2)
J Sci Comput (2012) 53:451–480
457
Thus, for any v ∈ L2ω˜ α (), we have v(x) =
∞
(α,β)
v˜l
L˜ l
(α,β)
(3.3)
(x)
l=l¯α
with (α,β)
v˜l
=
1
v(x)L˜ l
(α,β)
(α,β)
ηl
(x)ω˜ α (x)dx.
(3.4)
Next, let β (α,β) (α,β) Q˜ N () = e− 2 x ψ | ψ ∈ QN () . ˜ (α,β) The L2ω˜ α ()-orthogonal projection P˜N,α,β, : L2ω˜ α () → Q () is defined by N (P˜N,α,β, v − v, φ)ω˜ α , = 0,
(α,β) ∀φ ∈ Q˜ N (),
or equivalently, P˜N,α,β, v(x) =
N
(α,β)
v˜l
L˜ l
(α,β)
(3.5)
(x).
l=l¯α
We now estimate the approximation error. To do this, we introduce the Sturm-Liouville operator: β β A˜ α,β v(x) = −x −α e 2 x ∂x x α+1 e−βx ∂x e 2 x v(x) . Evidently, A˜ α,β L˜ l
(α,β)
(α,β)
(x) = λl
L˜ l
(α,β)
(x),
l ≥ l¯α .
(3.6)
By integration by parts, we obtain
(α,β)
v, L˜ l
ω˜ α ,
(α,β) −1 (α,β) −1 (α,β) (α,β) = λl = λl A˜ α,β v, L˜ l . v, A˜ α,β L˜ l ω˜ α , ω˜ α ,
(3.7)
Therefore, if u, v are in the domain of the operator A˜ α,β , then (A˜ α,β u, v)ω˜ α , = (u, A˜ α,β v)ω˜ α , . Accordingly, A˜ α,β is a positive definite and self-conjugate operator. Thus, we could define the following Sobolev-type spaces with integer r ≥ 0, D A˜ rα,β = v | A˜ kα,β v ∈ L2ω˜ α (), 0 ≤ k ≤ r , equipped with the following semi-norm and norm, |v|D(A˜ r
α,β
r ˜ ) = Aα,β v
ω˜ α ,
vD(A˜ r
,
α,β )
=
r k=0
12 |v|2D(A˜ k ) α,β
.
Theorem 3.1 For any v ∈ D(A˜ rα,β ) and integers 0 ≤ μ ≤ r ≤ N + 1, |P˜N,α,β, v − v|D(A˜ μ
α,β )
≤ c(βN )μ−r |v|D(A˜ rα,β ) .
(3.8)
458
J Sci Comput (2012) 53:451–480
Proof We use (3.5) and (3.6) to deduce that |P˜N,α,β, v − v|2D(A˜ μ
α,β )
=
∞
(α,β) 2μ (α,β) (α,β) 2 ηl v˜l
λl
l=N+1
=
∞
(α,β) 2μ (α,β) −1 (α,β) 2 . ηl v, L˜ l ω˜ α ,
λl
l=N+1
Further, using (3.7) repeatedly, yields
(α,β) v, L˜ l
ω˜ α ,
(α,β) −r r (α,β) = λl A˜ α,β v, L˜ l . ω˜ α ,
The above two equalities, together with (2.3), lead to that |P˜N,α,β, v − v|2D(A˜ μ
α,β )
∞
=
(α,β) 2μ−2r (α,β) −1 ˜ r (α,β) 2 Aα,β v, L˜ l ηl ω˜ α ,
λl
l=N+1 ∞ (α,β) 2μ−2r (α,β) −1 r (α,β) 2 ≤ c λN+1 A˜ α,β v, L˜ l ηl ω˜
α ,
l=N+1
≤ c(βN )2μ−2r
∞
(α,β) −1
ηl
(α,β) 2 . ω˜ α ,
A˜ rα,β v, L˜ l
l=l¯α
On the other hand, thanks to (3.2)–(3.4), we have ∞ r 2 (α,β) −1 r (α,β) 2 A˜ v = A˜ α,β v, L˜ l ηl α,β ω˜ α , ω˜ l=l¯α
α ,
.
Consequently, |P˜N,α,β, v − v|D(A˜ μ
α,β )
≤ c(βN )μ−r A˜ rα,β v ω˜ α , = c(βN )μ−r |v|D(A˜ rα,β ) .
The proof is completed. 3.2 Laguerre Orthogonal Approximation with Negative Integer α
We now focus on a special Laguerre orthogonal approximation, which is the mostly appreciated for numerical solutions of high order differential equations defined on unbounded domains. In the forthcoming discussions, we assume integer m ≥ 1. Clearly, ∂xk L˜ l
(−m,β)
(0) = 0,
0 ≤ k ≤ m − 1, l ≥ m.
Proposition 3.1 For any v ∈ 0 Hω1α,β () and α ≤ 0, vωα,β , ≤
2 ∂x vωα,β , . β
(3.9)
J Sci Comput (2012) 53:451–480
459
Proof For any x ∈ , we have v 2 (x)x α e−βx =
x
∂ξ v 2 (ξ )ξ α e−βξ dξ
0
x
=2
α −βξ
v(ξ )∂ξ v(ξ )ξ e 0
x
dξ + α
v 2 (ξ )ξ α−1 e−βξ dξ
0
x
−β
v 2 (ξ )ξ α e−βξ dξ.
(3.10)
0
Since v ∈ 0 Hω1α,β (), we have v 2 (x)x α e−βx → 0 as x → ∞. Letting x → ∞ in (3.10), and using the Cauchy-Schwarz inequality, we derive that v2ωα,β , ≤
4 2α 4 v2ωα−1,β , ≤ 2 ∂x v2ωα,β , . ∂x v2ωα,β , + 2 β β β
This is the desired result. We next introduce the following space with integer r ≥ 0, Hωr˜ −m ,A () = v | v is measurable on and vHωr˜
−m ,A
max(m, n), we introduce the spaces r r (I ) = v | v is measurable on I and vHm,n,A Hm,n,A (I ) < ∞ , r r H0,m,n,A (I ) = v ∈ Hm,n,A (I ) | ∂xk v(−2) = 0 for 0 ≤ k ≤ n − 1, and ∂xk v(0) = 0 for 0 ≤ k ≤ m − 1 ,
466
J Sci Comput (2012) 53:451–480
equipped with the following semi-norm and norm, r r |v|Hm,n,A (I ) = ∂x v χ (−m+r,−n+r) ,I ,
r vHm,n,A (I ) =
r k=0
12 |v|2H k (I ) m,n,A
.
Next, Jl(m,n) (x) stands for the standard Jacobi polynomial of degree l, and (m,n) Yl(m,n) (x) = χ (m,n) (x)Jl−m−n (x + 1), x ∈ I, l ≥ m + n, (I ) = span Yl(m,n) (x), m + n ≤ l ≤ N . Q(m,n) N μ,0 μ For integers max(m, n) ≤ μ ≤ m + n, the projection PN,m,n,I : H0,m,n,A (I ) → Q(m,n) (I ) is N defined by
μ,0 v , ∂xμ φ χ (−m+μ,−n+μ) = 0, ∂xμ v − PN,m,n,I
∀φ ∈ Q(m,n) (I ). N
− + We also introduce the polynomials qm,n,j,I (x), qm,n,j,I (x) ∈ Pm+n−1 (I ) as follows,
− (x) = qm,n,j,I
n−1−j (m + l − 1)! 1 m (−x) (2 + x)l+j , m l l!(m − 1)! 2 j! 2 l=0
m−1−j (n + l − 1)! (−1)j + (−x)l+j . (x) = n (2 + x)n qm,n,j,I l l!(n − 1)! 2 j! 2 l=0
(5.1)
It can be checked that − (−2) = δk,j , ∂xk qm,n,j,I + (−2) ∂xl qm,n,j,I
− ∂xl qm,n,j,I (0) = 0, + ∂xk qm,n,j,I (0)
= 0,
= δk,j ,
0 ≤ j, k ≤ n − 1, 0 ≤ l ≤ m − 1, 0 ≤ j, k ≤ m − 1, 0 ≤ l ≤ n − 1.
(5.2)
r For any v ∈ Hm,n,A (I ) and r ≥ max(m, n), we set
vm,n,b,I (x) =
n−1 j =0
− ∂xj v(−2)qm,n,j,I (x) +
m−1
+ ∂xj v(0)qm,n,j,I (x).
j =0
Let v(x) ˜ = v(x) − vm,n,b,I (x). The Jacobi quasi-orthogonal projection is defined by (cf. (2.18) of [19]) μ μ,0 PN,m,n,I v(x) = PN,m,n,I v(x) ˜ + vm,n,b,I (x).
(5.3)
By using the results (2.21) and (2.24) of [19], we reach the following result. μ Proposition 5.1 If v ∈ Hm,n,A (I ), ∂xr v ∈ L2χ (−m+r,−n+r) (I ), integers m, n, r ≥ 1, N ≥ m + n, 0 ≤ k ≤ r ≤ N + 1, max(m, n, k) ≤ μ ≤ m + n, and r ≥ m + n (or max(m, n) ≤ r ≤ m + n − 1, 1 ≤ m, n ≤ 4), then
k μ ∂ P x
N,m,n,I v
− v χ (−m+k,−n+k) ,I ≤ cN k−r ∂xr v χ (−m+r,−n+r) ,I .
(5.4)
J Sci Comput (2012) 53:451–480
467
Now, let = {x | −2 < x < ∞}. We denote the inner product and norm of L2χ () by (· , ·)χ and · χ , respectively. We omit χ in notations, whenever χ ≡ 1. Let W (), W () H 2 () and λ ≥ 0. We introduce the bilinear form Aλ (u, v) = ∂x2 u, ∂x2 v + λ(u, v), ∀u ∈ W (), v ∈ W (). Next, we set N = (N1 , N2 ), and ˜ (0,β) WN () = φ ∈ W ()|φ|I ∈ PN1 (I ), φ| ∈ Q N2 () , ˜ (0,β) W N () = φ ∈ W ()|φ|I ∈ PN1 (I ), φ| ∈ Q N2 () . 2 We define the operator ∗ PN,λ : W () → WN (), by
2 Aλ ∗ PN,λ v − v, φ = 0,
∀φ ∈ W N ().
(5.5)
2 Proposition 5.2 If v ∈ W (), w ∈ WN () and ∗ PN,λ v − w ∈ W N (), then
2 2 Aλ v − ∗ PN,λ v, v − ∗ PN,λ v ≤ Aλ (v − w, v − w).
(5.6)
Proof A direct calculation shows 2 2 2 2 Aλ (v − w, v − w) = Aλ v − ∗ PN,λ v, v − ∗ PN,λ v + Aλ ∗ PN,λ v − w, ∗ PN,λ v−w 2 2 + 2Aλ v − ∗ PN,λ v, ∗ PN,λ v−w . Thanks to (5.5), we have
2 2 Aλ v − ∗ PN,λ v, ∗ PN,λ v − w = 0.
Since λ ≥ 0, we have
2 2 Aλ ∗ PN,λ v − w, ∗ PN,λ v − w ≥ 0.
Then, the desired result (5.6) follows from the previous statements immediately.
5.2 Dirichlet Problem of High Order Equation We first consider the following simple model problem, ∂x4 U (x) + λU (x) = f (x), x ∈ , ∂x U (−2) = b, U (−2) = a,
(5.7)
where f ∈ L2 (), a, b and λ > 0 are given constants. In addition, we require that 1 1 x 2 ∂xk U (x) → 0, as x → ∞, 0 ≤ k ≤ 2, and x − 2 ∂x3 U (x) → 0, as x → ∞. Let 1 H˜ 2 () = H 2 () ∩ v | x − 2 ∂x3 v(x) → 0, as x → ∞ , V () = 0 H 2 (). V () = v | v ∈ H˜ 2 () and ∂x v(−2) = b, v(−2) = a , Multiplying (5.7) by v ∈ V () and integrating the resulting equation over , we obtain the weak formulation of (5.7). It is to look for U ∈ V () such that
468
J Sci Comput (2012) 53:451–480
Aλ (U, v) = (f, v),
∀v ∈ V ().
(5.8)
It can be shown that problem (5.8) admits a unique solution. We now construct the domain decomposition spectral scheme for (5.8). Let ˜ (0,β) VN () = φ ∈ V ()|φ|I ∈ PN1 (I ), φ| ∈ Q N2 () , ˜ (0,β) V N () = φ ∈ V ()|φ|I ∈ PN1 (I ), φ| ∈ Q N2 () . The domain decomposition spectral method for (5.8), is to seek uN ∈ VN () such that Aλ (uN , φ) = (f, φ),
∀φ ∈ V N ().
(5.9)
For the existence of solutions of (5.9), it suffices to check the uniqueness of solution. Assume (2) (2) that u(1) ˜ N (x) = u(1) N (x) and uN (x) are solutions of (5.9), and u N (x) − uN (x) ∈ V N (). Then Aλ (u˜ N , φ) = 0,
∀φ ∈ V N ().
(5.10)
Putting φ = u˜ N in (5.10), we obtain 2 2 ∂ u˜ N + λu˜ N 2 = Aλ (u˜ N , u˜ N ) = 0. x Thus, u˜ N (x) ≡ 0. This implies the existence and the uniqueness of solution of (5.9). For numerical analysis of the proposed method, we introduce the composite quasiorthogonal projection PN v(x), defined by PN21 ,2,2,I v(x), x ∈ I, PN v(x) = (5.11) P˜N22 ,−2,β, v(x), x ∈ . Obviously, PN U ∈ VN (). By using (5.4), we deduce that for integer r1 ≥ 2, k 2 2 ∂ U − P 2 ≤ c ∂ k U − P 2 (−2+k,−2+k) x N1 ,2,2,I U x N1 ,2,2,I U I χ ,I 2 ≤ cN1 2(k−r1 ) ∂xr1 U χ (−2+r1 ,−2+r1 ) ,I , k = 0, 1, 2.
(5.12)
On the other hand, according to (3.14) and (3.16), we set
β β U˜ (x) = U (x) − Ub,2, (x), Ub,2, (x) = e− 2 x U (0) + ∂x U (0) + U (0) x , 2 β β P˜N22 ,−2,β, U (x) = e− 2 x 0 PN22 ,−2,β, e 2 x U˜ (x) + Ub,2, (x). We use (2.16) and (2.10) to derive that for integer r2 ≥ 2, k 2 β x 2 ∂ U − P˜ 2 = ∂ k e− β2 x e β2 x U˜ − 0 P 2 ˜ 2 U x N2 ,−2,β, U x N2 ,−2,β, e β β 2 ≤ cβ 2k−4 ∂x2 e 2 x U˜ − 0 PN22 ,−2,β, e 2 x U˜ ω
0,β ,
β 2 ≤ cβ 2k−4 (βN2 )2−r2 ∂xr2 e 2 x U˜ ω β 2 ≤ cβ 2k−4 (βN2 )2−r2 ∂xr2 e 2 x U ω
−2+r2 ,β ,
−2+r2 ,β ,
,
k = 0, 1, 2. (5.13)
J Sci Comput (2012) 53:451–480
469
We now estimate the error of numerical solution. For this purpose, we introduce the 2 auxiliary projection P N,λ : V () → VN (), defined by
2
Aλ P N,λ v − v, φ = 0,
∀φ ∈ V N ().
(5.14)
∀φ ∈ V N ().
(5.15)
∀φ ∈ V N ().
(5.16)
We have from (5.8) and (5.14) that
2
Aλ P N,λ U, φ = (f, φ),
Subtracting (5.15) from (5.9), yields
2
Aλ uN − P N,λ U, φ = 0, 2
Taking φ = uN − P N,λ U ∈ V N () in (5.16), we obtain
2
2
Aλ uN − P N,λ U, uN − P N,λ U = 0. 2
This implies uN = P N,λ U . Let r = (r1 , r2 ). For simplicity, we shall use the following notations, 2 1 2 A(v, k) = ∂xk v χ (−2+k,−2+k) ,I + ∂xk v 2 , 2 B(v, N, r) = N1 4−2r1 ∂ r1 U (−2+r ,−2+r ) x
+ 1+β
−4
1
χ
1 ,I
β 2 (βN2 )2−r2 ∂xr2 e 2 x U ω
12 −2+r2 ,β ,
.
Furthermore, we put in Proposition 5.2 that W () = V (), v = U,
W () = V (), w = PN U,
2 ∗ PN,λ U
WN () = VN (),
W N () = V N (),
2
= P N,λ U.
Then, we use the Gagliardo-Nirenberg inequality, projection theorem with (5.14), (5.12) and (5.13) successively, to obtain that for integers r1 , r2 ≥ 2, U − uN H 2 () ≤ c ∂x2 (U − uN ) + λU − uN 2 2 = c ∂x2 U − P N,λ U + λ U − P N,λ U ≤ c ∂x2 (U − PN U ) + λU − PN U ≤ c A(U − PN U, 2) + λA(U − PN U, 0) ≤ cB(U, N, r).
(5.17)
5.3 Mixed Inhomogeneous Boundary Value Problem We now consider the following mixed inhomogeneous boundary value problem, ∂x4 U (x) + λU (x) = f (x), x ∈ , U (−2) = a, ∂x2 U (−2) = b,
(5.18)
470
J Sci Comput (2012) 53:451–480 1
where f ∈ L2 (), a, b and λ > 0 are given constants. We also require that x 2 ∂xk U (x) → 0, 1 as x → ∞, 0 ≤ k ≤ 2, and x − 2 ∂x3 U (x) → 0, as x → ∞. Let V () = v | v ∈ H˜ 2 () and v(−2) = a , V () = v | v ∈ H 2 () and v(−2) = 0 . The weak formulation of (5.18) is to look for U ∈ V () such that Aλ (U, v) + b∂x v(−2) = (f, v),
∀v ∈ V ().
(5.19)
It can be shown that problem (5.19) admits a unique solution. Next, let ˜ (0,β) VN () = φ ∈ V () | φ|I ∈ PN (I ), φ| ∈ Q () , N 1
2
˜ (0,β) V N () = φ ∈ V () | φ|I ∈ PN1 (I ), φ| ∈ Q N2 () .
The domain decomposition spectral method for (5.19), is to seek uN ∈ VN () such that Aλ (uN , φ) + b∂x φ(−2) = (f, φ),
∀φ ∈ V N ().
(5.20)
Like (5.9), the above problem has a unique solution. We now estimate the error of numerical solution. Let PN v(x) be the same as in (5.11). 2 We also introduce the auxiliary operator P N,λ : V () → VN (), defined by
2
Aλ P N,λ v − v, φ = 0,
∀φ ∈ V N ().
(5.21)
We have from (5.19) and (5.21) that
2
Aλ P N,λ U, φ + b∂x φ(−2) = (f, φ),
∀φ ∈ V N ().
(5.22)
Subtracting (5.22) from (5.20), yields
2
Aλ uN − P N,λ U, φ = 0,
∀φ ∈ V N ().
(5.23)
2
Taking φ = uN − P N,λ U ∈ V N () in (5.23), we reach
2
2
Aλ uN − P N,λ U, uN − P N,λ U = 0. 2
Thereby, uN = P N,λ U . Finally, by virtue of the Gagliardo-Nirenberg inequality, projection theorem with (5.21), (5.12), (5.13) and Proposition 5.2 with W () = V (), v = U,
W () = V (), w = PN U,
2 ∗ PN,λ U
we obtain that for integers r1 , r2 ≥ 2,
WN () = VN (), =
2 P N,λ U,
W N () = V N (),
J Sci Comput (2012) 53:451–480
471
U − uN H 2 () ≤ c ∂x2 (U − uN ) + λU − uN = c ∂x2 (U − P N,λ U ) + λU − P N,λ U ≤ c ∂x2 (U − PN U ) + λU − PN U ≤ c A(U − PN U, 2) + λA(U − PN U, 0) ≤ cB(U, N, r).
(5.24)
Remark 5.1 We may use spectral method without domain decomposition to solve (5.18). In this case, we let (α,β) () = span Ll (x + 2), l¯α ≤ l ≤ N , β (α,β) (α,β) Q˜ N () = span e− 2 (x+2) φ, φ ∈ QN () , (α,β)
QN
VN () = V () ∩ Q˜ N
(0,β)
(),
V N () = V () ∩ Q˜ N
(0,β)
().
The spectral method for (5.18) is to find uN ∈ VN () such that Aλ (uN , φ) + b∂x φ(−2) = (f, φ),
∀φ ∈ V N ().
(5.25)
∗ (x) = ωα,β (x + 2). By an argument similar to the derivation of (5.24), we obtain Let ωα,β
β 2 U − uN 2H 2 () ≤ c 1 + β −4 (βN )2−r ∂xr e 2 x U ω∗
−2+r,β
.
5.4 Spectral Method with Essential Imposition of Boundary Conditions The Laguerre quasi-orthogonal approximation also plays an important role in spectral method with essential imposition of mixed inhomogeneous boundary conditions. To do this, we set V () = v | v ∈ H˜ 2 () and v(−2) = a ,
V () = v | v ∈ H 2 () and v(−2) = 0 .
The weak formulation of (5.18) is to look for U ∈ V () such that Aλ (U, v) + b∂x v(−2) = (f, v),
∀v ∈ V ().
(5.26)
This problem admits a unique solution. Let () ∩ v | ∂x2 v(−2) = b , (0,β) V N () = V () ∩ Q˜ N () ∩ v | ∂x2 v(−2) = 0 .
˜N VN () = V () ∩ Q
(0,β)
The spectral method for (5.26), is to seek uN ∈ VN () such that Aλ (uN , φ) + b∂x φ(−2) = (f, φ),
∀φ ∈ V N ().
Like (5.9) and (5.20), problem (5.27) possess a unique solution.
(5.27)
472
J Sci Comput (2012) 53:451–480 2
We now deal with the convergence of (5.27). We introduce the projection P N,λ : V () → VN (), defined by 2 Aλ P N,λ v − v, φ = 0, ∀φ ∈ V N (). (5.28) We have from (5.26) and (5.28) that
2
Aλ P N,λ U, φ + b∂x φ(−2) = (f, φ),
∀φ ∈ V N ().
(5.29)
Subtracting (5.29) from (5.27), yields
2
Aλ uN − P N,λ U, φ = 0,
∀φ ∈ V N ().
(5.30)
2
Taking φ = uN − P N,λ U ∈ V N () in (5.30), we obtain
2
2
Aλ uN − P N,λ U, uN − P N,λ U = 0. 2
This implies uN = P N,λ U . For error estimation of numerical solution, we need the Laguerre quasi-orthogonal projection defined on . According to (3.14) and (3.16), we introduce the function U˜ (x) = U (x) − Ub,3, (x), with
β β Ub,3, (x) = e− 2 (x+2) U (−2) + ∂x U (−2) + U (−2) (x + 2) 2
1 2 β β2 2 + ∂ U (−2) + ∂x U (−2) + U (−2) (x + 2) . 2 x 2 8 ∗ (x) be the same as in Remark 5.1. We define the orthogonal projecLet QN () and ωα,β 3 tion 0 PN,−3,β, by (α,β)
3 3 ∂x 0 PN,−3,β, v − v , ∂x3 v ω∗ = 0, 0,β
(−3,β)
∀φ ∈ QN
().
Like (3.16), the Laguerre quasi-orthogonal projection defined on , is defined by β (x+2) β 3 3 U (x) = e− 2 (x+2) 0 PN,−3,β, U˜ (x) + Ub,3, (x). e2 P˜N,−3,β, Thanks to Lemma 2.2 of [16], for any v ∈ 0 Hω1∗ (), 0,β
v2ω∗ ≤ 0,β
4 ∂x v2ω∗ . 0,β β2
(5.31)
Accordingly, we use an equality like (2.16), (5.31) and an inequality like (2.10) successively to derive that for integer r ≥ 3, k 2 k − β (x+2) β (x+2) β (x+2) 2 3 ∂ U − P˜ 3 = ∂ e 2 U˜ − 0 PN,−3,β, U˜ e2 e2 x N,−3,β, U x β β (x+2) 2 3 ≤ cβ 2k−4 ∂x2 e 2 (x+2) U˜ − 0 PN,−3,β, U˜ ω∗ e2 0,β
≤
cβ 2k−6 ∂x3
e
β 2 (x+2)
U˜
3 − 0 PN,−3,β,
e
β 2 (x+2)
2 U˜ ∗
ω0,β
J Sci Comput (2012) 53:451–480
473
β 2 ≤ cβ 2k−6 (βN )3−r ∂xr e 2 (x+2) U˜ ω∗
−3+r,β
β 2 ≤ cβ 2k−6 (βN )3−r ∂xr e 2 (x+2) U ω∗
−3+r,β
,
k = 0, 1, 2. (5.32)
On the other hand, like Proposition 5.2, we can prove that
3 3 Aλ U − P N,λ U, U − P N,λ U ≤ Aλ U − P˜N,−3,β, U, U − P˜N,−3,β, U . 2
2
Finally, we use the Gagliardo-Nirenberg inequality, projection theorem with (5.28), and (5.32) successively, to obtain that for integer r ≥ 3, 2 U − uN 2H 2 () ≤ c ∂x2 (U − uN ) + λU − uN 2 2 = c ∂x2 (U − P N,λ U ) + λU − P N,λ U 2 2 2 3 3 ≤ c ∂x2 U − P˜N,−3,β, U + λ U − P˜N,−3,β, U
β 2 c λ . (5.33) ≤ 2 1 + 4 (βN )3−r ∂xr e 2 x U ω∗ −3+r,β β β For other spectral and collocation methods with essential imposition of boundary conditions, we refer to the work of [1, 22]. 6 Numerical Results In this section, we present some numerical results to illustrate the efficiency of proposed methods. 6.1 Spectral Scheme (5.9) We need six kinds of basis functions. The first kind of basis functions are local, and correspond to subdomains I and respectively, namely, 1 Y (2,2) (x), x ∈ I, (6.1) ψ1,l (x) = 4(l−2)(l−3) l 0, otherwise, 2 (−2,β) β L˜ (x), x ∈ , ψ2,l (x) = l(l−1) l (6.2) 0, otherwise. Next, we define the basis function for matching numerical solution at x = 0, as (x), x ∈ I, q+ σ (x) = 2,2,0,Iβ β (1 + 2 x)e− 2 x , x ∈ .
(6.3)
Similarly, we define the basis function for matching the first derivative of numerical solution at x = 0, as (x), x ∈ I , q+ (6.4) ρ(x) = 2,2,1,I β x ∈ . xe− 2 x ,
474
J Sci Comput (2012) 53:451–480
Obviously, σ (0) = ∂x ρ(0) = 1 and ∂x σ (0) = ρ(0) = 0. Furthermore, we take the basis function for matching numerical solution at x = −2, as (x), x ∈ I , q− (6.5) ζ (x) = 2,2,0,I 0, otherwise. We also take the basis function for matching the first order derivative of numerical solution at x = −2, as (x), x ∈ I , q− (6.6) ξ(x) = 2,2,1,I 0, otherwise. Clearly, ζ (−2) = ∂x ξ(−2) = 1 and ∂x ζ (−2) = ξ(−2) = 0. In actual computation, we expand the numerical solution as uN (x) =
N1
uˆ 1,l ψ1,l (x) +
l=4
N2
uˆ 2,l ψ2,l (x) + uN (0)σ (x) + ∂x uN (0)ρ(x)
l=2
+ aζ (x) + bξ(x) ∈ VN ().
(6.7)
Let g(φ) = (f, φ). Inserting (6.7) into the domain decomposition spectral scheme (5.9), we obtain Aλ (uN , φ) = g(φ),
∀φ ∈ V N ().
(6.8)
where the test functions φ stand for the basis functions given by (6.1)–(6.5). We now use algorithm (6.8) to solve problem (5.7) with λ = 1. Let xLe,N1 ,j and ωLe,N1 ,j be the nodes and weights of the Legendre-Gauss-Lobatto numerical quadrature on the interval (−1, 1), respectively. Meanwhile, xLa,N2 ,j and ωLa,N2 ,j stand for the nodes and weights of the Leguerre-Gauss-Radau numerical quadrature on the interval (0, ∞), respectively. We measure the numerical errors by the following quantities:
N 1 2 U (xLe,N1 ,j − 1) − uN (xLe,N1 ,j − 1) ωLe,N1 ,j U − uN N,ave = j =0
+
N2
2 U (xLa,N2 ,j ) − uN (xLa,N2 ,j ) ωLa,N2 ,j
12 ,
j =1
U − uN N,max = max max U (xLe,N1 ,j − 1) − uN (xLe,N1 ,j − 1), 0≤j ≤N1
max U (xLa,N2 ,j ) − uN (xLa,N2 ,j ) .
1≤j ≤N2
We first take the test function U (x) = e−θx sin γ x,
θ, γ > 0,
(6.9)
which decays exponentially as x goes to infinity. In Tables 1 and 2, we list the values of U − uN N,ave and U − uN N,max , with θ = γ = 1, N2 = 4N1 and β = 1, 2, v.s. the modes N1 . The numerical results demonstrate that the numerical errors decay fast as N1 increases. This confirms the theoretical analysis in Sect. 5.2.
J Sci Comput (2012) 53:451–480 Table 1 The values U − uN N,ave of algorithm (6.8)
Table 2 The values U − uN N,max of algorithm (6.8)
475 N1 = 5
N1 = 10
N1 = 15
β =1
3.31E–3
4.85E–9
1.81E–13
β =2
2.86E–3
4.84E–9
6.88E–14
N1 = 5
N1 = 10
N1 = 15
β =1
1.81E–3
2.62E–9
3.51E–14
β =2
1.95E–3
2.62E–9
2.35E–14
We also see that a suitable choice of parameter β could raise the numerical accuracy. However, so far, there is no theoretical result on the best choice of β. Generally speaking, it seems reasonable to take certain β so that the asymptotic behaviors of numerical solutions are more similar to those of exact solutions. Clearly, the test function (6.9) with θ = 1 decays like e−x , while the numerical solution with β = 2 does so also. Thus it is more accurate. For a given problem, the better choice of β depends on some parameters involved in the underlying differential equation, which affect the asymptotic behavior of solution. For instance, we consider the following simple equation, ∂x4 U (x) + λU (x) = f (x),
λ > 0, 0 < x < ∞.
In addition, ∂xm U (x) → 0 as x → ∞, m = 1, 2. If f (x) fulfills certain condition, then the √ √ 4 solution U (x) decays like e− λx as x → ∞. In this case, we may take β = 2 4 λ. On the other hand, if the solution varies faster near the fixed boundary, then it is better to take suitably bigger β, since the distances between two adjacent roots of the scaled Laguerre functions with bigger β are small than those with small β. We refer to [17, 18] for more discussions. We next use (6.8) to solve (5.7) with the test function U (x) =
sin x + 1 , (x + 3)k
k ≥ 1,
(6.10)
which decays algebraically as x goes to infinity. In Tables 3 and 4, we list the values of U − uN N,ave and U − uN N,max , with k = 3, N2 = 4N1 and β = 1, 2, v.s. the mode N1 . The numerical results show again the high accuracy of algorithm (6.8), and the affecting of parameter β. In Table 5, we present the values of U − uN N,ave and U − uN N,max , with k = 3, 4 and β = 2, v.s. the mode N1 . The numerical results indicate that the smoother the exact solution, the better the numerical results. This coincides with the theoretical analysis. In fact, for the test function U (x) given by (6.10), β the norm ∂xr (e 2 x U )2ω−2+r,β , is finite, as long as r < 2k + 1. Thereby, the estimate (5.17) implies U − uN H 2 () β ≤ c N1 2−2k ∂x2k U χ (−2+2k,−2+2k) ,I + 1 + β −2 (βN2 )1−k ∂x2k e 2 x U ω
−2+2k,β ,
.
476
J Sci Comput (2012) 53:451–480
Table 3 The values U − uN N,ave of algorithm (6.8)
Table 4 The values U − uN N,max of algorithm (6.8)
N1 = 5
N1 = 10
N1 = 15
N1 = 20
N1 = 25
β =1
1.74E–2
6.70E–5
5.89E–6
2.74E–6
1.48E–6
β =2
1.80E–2
6.55E–5
2.90E–6
1.50E–6
8.42E–7
N1 = 5
N1 = 10
N1 = 15
N1 = 20
N1 = 25
β =1
1.12E–2
3.63E–5
2.57E–6
1.20E–6
6.34E–7
β =2
1.12E–2
3.63E–5
1.53E–6
6.50E–7
3.48E–7
Table 5 The errors of algorithm (6.8) k=3
k=4
U − uN N,ave
U − uN N,max
U − uN N,ave
U − uN N,max
N1 = 15
2.90E–6
1.53E–6
7.41E–7
3.32E–7
N1 = 20
1.50E–6
6.50E–7
1.40E–8
6.35E–9
N1 = 25
8.42E–7
3.48E–7
6.36E–9
2.65E–9
6.2 Spectral scheme (5.20) We need the same six kinds of basis functions as in (6.1)–(6.6). In actual computation, we expand the numerical solution as uN (x) =
N1 l=4
uˆ 1,l ψ1,l (x) +
N2
uˆ 2,l ψ2,l (x) + uN (0)σ (x) + ∂x uN (0)ρ(x)
l=2
+ aζ (x) + ∂x uN (−2)ξ(x) ∈ VN ().
(6.11)
Let g(φ) = (f, φ) − b∂x φ(−2). Inserting (6.11) into (5.20), we obtain Aλ (uN , φ) = g(φ),
∀φ ∈ V N ().
(6.12)
where the test functions φ stand for the basis functions given by (6.1)–(6.6). We now use algorithm (6.12) to solve the mixed inhomogeneous boundary value problem (5.18) with λ = 1. We first take the test function (6.9). In Tables 6 and 7, we list the values of U − uN N,ave and U − uN N,max , with σ = γ = 1, N2 = 4N1 and β = 1, 2, v.s. the modes N1 . Evidently, the errors decay very fast as N1 increases. Also, a suitable choice of parameter β leads to more accurate numerical results. We next use (6.12) to solve (5.18) with the test function (6.10). In Tables 8 and 9, we list the values of U − uN N,ave and U − uN N,max , with k = 3, N2 = 4N1 and β = 1, 2, v.s. the mode N1 . In Table 10, we present the values of U − uN N,ave and U − uN N,max , with k = 3, 4 and β = 2, v.s. the mode N1 . This coincides with the theoretical analysis again. For comparison, we consider problem (5.18) with the test function (6.9) and θ = 4, γ = 2. Clearly, the solution oscillates seriously near the endpoint x = −2. We use (6.12) with N2 = 4N1 and (5.25) with N = 5N1 , to solve the same problem. The numerical errors of
J Sci Comput (2012) 53:451–480 Table 6 The values U − uN N,ave of algorithm (6.12)
Table 7 The values U − uN N,max of algorithm (6.12)
Table 8 The values U − uN N,ave of algorithm (6.8)
Table 9 The values U − uN N,max of algorithm (6.12)
477 N1 = 5
N1 = 10
N1 = 15
β =1
3.65E–3
4.85E–9
1.86E–13
β =2
2.86E–3
4.84E–9
7.55E–14
N1 = 5
N1 = 10
N1 = 15
β =1
1.84E–3
2.62E–9
3.51E–14
β =2
1.82E–3
2.62E–9
2.43E–14
N1 = 5
N1 = 10
N1 = 15
N1 = 20
N1 = 25
β =1
6.45E–2
6.74E–5
5.89E–6
2.74E–6
1.48E–6
β =2
6.76E–2
6.58E–5
2.90E–6
1.49E–6
8.43E–7
N1 = 5
N1 = 10
N1 = 15
N1 = 20
N1 = 25
β =1
4.00E–2
3.72E–5
2.57E–6
1.20E–6
6.34E–7
β =2
4.00E–2
3.72E–5
1.53E–6
6.50E–7
3.49E–7
Table 10 The errors of algorithm (6.12) k=3
k=4
U − uN N,ave
U − uN N,max
U − uN N,ave
U − uN N,max
N1 = 15
2.90E–6
1.53E–6
7.41E–7
3.31E–7
N1 = 20
1.49E–6
6.50E–7
1.40E–8
6.35E–9
N1 = 25
8.43E–7
3.49E–7
6.36E–9
2.65E–9
(5.25) are measured by
U − uN N,ave,∗ =
N
U (xLa,N,j − 2) − uN (xLa,N,j
2 − 2) ωLa,N,j
12 ,
j =0
U − uN N,max,∗ = max U (xLa,N,j − 2) − uN (xLa,N,j − 2). 0≤j ≤N
In Table 11, we list the values of U − uN N,ave and U − uN N,max of algorithm (6.12), and the values of U − uN N,ave,∗ and U − uN N,max,∗ of algorithm (5.25). The parameter β = 2. Obviously, the spectral method with domain decomposition provides better numerical results.
478
J Sci Comput (2012) 53:451–480
Table 11 The errors of (6.12) with N2 = 4N1 and (5.25) with N = 5N1 Algorithm (6.12)
Algorithm (5.25)
U − uN N,ave
U − uN N,max
U − uN N,ave,∗
U − uN N,max,∗
N1 = 20
3.88E–10
1.47E–10
6.33E–9
2.18E–9
N1 = 25
7.64E–11
1.65E–11
1.53E–8
5.30E–9
6.3 Spectral scheme (5.27) We define the basis functions as follows,
β β2 β ψ0 (x) = 1 + (x + 2) + (x + 2)2 e− 2 (x+2) , 2 8
β β ψ1 (x) = x + 2 + (x + 2)2 e− 2 (x+2) , 2
(6.13)
β 1 ψ2 (x) = (x + 2)2 e− 2 (x+2) , 2
ψl (x) =
β3 (−3,β) L˜ (x + 2), l(l − 1)(l − 2) l
3 ≤ l ≤ N.
Evidently, ∂xk ψl (−2) = δk,l for 0 ≤ k, l ≤ 2, and ∂xk ψl (−2) = 0 for 0 ≤ k ≤ 2, 3 ≤ l ≤ N . In actual computation, we expand the numerical solution as uN (x) =
N
uˆ N,l ψl (x) + bψ2 (x) + ∂x uN (−2)ψ1 (x) + aψ0 (x) ∈ VN ().
l=3
Let g(φ) be the same as in (6.12). Then, the scheme (5.27) with essential imposition of boundary condition, becomes Aλ (uN , φ) = g(φ),
∀φ ∈ V N ().
(6.14)
Here we take the test functions φ, as presented in (6.13). We now use (6.14) to solve (5.18) with λ = 1 and the test function (6.9) with θ = γ = 1. In Table 12, we present the values of U − uN N,ave,∗ and U − uN N,max,∗ of algorithms (5.25) and (6.14), with β = 2, v.s. the mode N . We find that like [1], the algorithm (6.14) with essential imposition of boundary condition, provides better numerical results, than scheme (5.25) without domain decomposition and essential imposition of boundary condition. But, the error estimate (5.33) is not as good as the estimate for (5.25), see Remark 5.1. How to improve the estimate (5.33) is still an open problem.
7 Concluding Remarks In this paper, we developed the generalized Laguerre quasi-orthogonal approximation by using the new basis functions, which are mutually orthogonal with the weight function x α , α being any real number. This new approximation matches the considered function at the origin properly, and fits the asymptotic behaviors (growing or decaying algebraically at the
J Sci Comput (2012) 53:451–480
479
Table 12 The errors of algorithms (5.25) and (6.14) Algorithm (5.25)
Algorithm (6.14)
U − uN N,ave,∗
U − uN N,max,∗
U − uN N,ave,∗
U − uN N,max,∗
N = 100
1.44E–12
5.80E–13
1.36E–12
5.07E–13
N = 125
2.65E–11
9.08E–12
1.15E–12
3.89E–13
infinity) reasonably. By using this approximation, it is more convenient to deal with various mixed inhomogeneous boundary value problems of high order differential equations. Moreover, it plays an important role in the domain decomposition spectral method of high order problems defined on unbounded domains. We also studied the corresponding LaguerreGauss-Radau type interpolation, which serves as an important tool for the new Laguerre pseudospectral and collocation methods. As examples of applications, we provided the domain decomposition spectral methods for two high order problems with mixed inhomogeneous boundary conditions. The numerical results demonstrated their high accuracy. We also developed the spectral method with essential imposition of boundary condition, which often leads to better numerical results. Although we only considered several model problems, we could use the techniques developed in this work, coupled with other tricks, to treat other problems, such as the stream function form of the Navier Stokes equations in an infinite strap or outside certain obstacles, the infinite beam equations and so on.
References 1. Auteri, F., Parolini, N., Quartapelle, L.: Essential imposition of Neumann condition in Galerkin-Legendre elliptic solver. J. Comput. Phys. 185, 427–444 (2003) 2. Azaiez, M., Shen, J., Xu, C., Zhuang, Q.: A Laguerre-Legendre spectral method for the Stokes problems in a semi-infinite channel. SIAM J. Numer. Anal. 47, 271–292 (2008) 3. Bernardi, C., Maday, Y.: Spectral methods. In: Ciarlet, P.G., Lions, J.L. (eds.) Handbook of Numerical Analysis. Techniques of Scientific Computing, vol. 5 pp. 209–486. Amsterdam, Elsevier (1997) 4. Bernardi, C., Coppoletta, G., Maday, Y.: Some spectral approximations of multi-dimensional fourthorder problems. Internal report 90021, Laboratoire d’Analyse Numérique, Université Pierre et Marie Curie, Paris (1990) 5. Boyd, J.P.: Chebyshev and Fourier Spectral Methods, 2nd edn. Dover, New York (2001) 6. Canuto, C., Hussaini, M.Y., Quarteroni, A., Zang, T.A.: Spectral Methods, Fundamentals in Single Domains. Springer, Berlin (2006) 7. Funaro, D.: Polynomial Approximations of Differential Equations. Springer, Berlin (1992) 8. Funaro, D., Kavian, O.: Approximation of some diffusion evolution equations in unbounded domains by Hermite function. Math. Comput. 57, 597–619 (1999) 9. Gottlieb, D., Orszag, S.A.: Numerical Analysis of Spectral Methods: Theory and Applications. SIAM/CBMS, Philadelphia (1977) 10. Guo, B.-y.: Spectral Methods and Their Applications. World Scientific, Singapore (1998) 11. Guo, B.-y.: Error estimation of Hermite spectral method for nonlinear partial differential equations. Math. Comput. 68, 1067–1078 (1999) 12. Guo, B.-y., Ma, H.-p.: Composite Legendre-Laguerre approximation in unbounded domains. J. Comput. Math. 19, 101–112 (2001) 13. Guo, B.-y., Shen, J.: Laguerre-Galerkin method for nonlinear partial differential equations on a semiinfinite interval. Numer. Math. 86, 635–654 (2000) 14. Guo, B.-y., Wang, T.-j.: Composite generalized Laguerre-Legendre spectral method with domain decomposition and its application to Fokker-Planck equation in an infinite channel. Math. Comput. 78, 129–151 (2009) 15. Guo, B.-y., Xu, C.-l.: Mixed Laguerre-Legendre pseudospectral method for incompressible fluid flow in an infinite strip. Math. Comput. 72, 95–125 (2003)
480
J Sci Comput (2012) 53:451–480
16. Guo, B.-y., Zhang, X.-y.: A new generalized Laguerre spectral approximation and its applications. J. Comput. Appl. Math. 181, 342–363 (2005) 17. Guo, B.-y., Zhang, X.-y.: Spectral method for differential equations of degenerate type by using generalized Laguerre functions. Appl. Numer. Math. 57, 455–471 (2007) 18. Guo, B.-y., Wang, L.-l., Wang, Z.-q.: Generalized Laguerre interpolation and pseudospectral method for unbounded domains. SIAM J. Numer. Anal. 43, 2567–2589 (2006) 19. Guo, B.-y., Sun, T., Zhang, C.: Jacobi and Laguerre quasi-orthogonal approximations and related interpolations. Math. Comput. (in press) 20. Karniadakis, G.E., Sherwin, S.J.: Spectral/hp Element Methods for CFD, 2nd edn. Oxford Univ. Press, Oxford (2005) 21. Shen, J.: Stable and efficient spectral methods in unbounded domains using Laguerre functions. SIAM J. Numer. Anal. 38, 1113–1133 (2000) 22. Wang, L.-l., Guo, B.-y.: Interpolation approximations based on Gauss-Lobatto-Legendre-Birkhoff quadrature. J. Approx. Theory 161, 142–173 (2009) 23. Xu, C.-l., Guo, B.-y.: Laguerre pseudospectral method for nonlinear partial differential equation. J. Comput. Math. 20, 413–428 (2002)
