- Email: [email protected]
A spectral collocation method based on Haar wavelets for Poisson equations and biharmonic equations
Mathematical and Computer Modelling 54 (2011) 2858–2868
Contents lists available at SciVerse ScienceDirect
Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm
A spectral collocation method based on Haar wavelets for Poisson equations and biharmonic equations Shi Zhi ∗ , Cao Yong-yan School of Science, Xi’an University of Architecture and Technology, Xi’an 710055, China
article
info
Article history: Received 27 February 2011 Received in revised form 4 July 2011 Accepted 7 July 2011 Keywords: Haar wavelet Poisson equation Biharmonic equation
abstract In this work, we present a computational method for solving Poisson equations and biharmonic equations which are based on the use of Haar wavelets. The first transform the spectral coefficients into the nodal variable values. The second use Kronecker products to construct the approximations for derivatives over a tensor product grid of the horizontal and vertical blocks. Finally, solve the obtained system of algebraic equations. The efficiency of the method is demonstrated by four numerical examples. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction The numerical solution of Poisson equations is an important problem in numerical analysis. In the past two decades, a great deal of research work has been published on the development of numerical solution of Poisson equations [1–6]. A number of alternatives to direct integration of the domain integral have been proposed, notably among them the important papers of Atkinson [7] and Golberg and Chen [8]. In mathematics, a biharmonic equation is a fourth-order partial differential equation which arises in areas of continuum mechanics, including linear elasticity theory and the solution of Stokes flows. The concept of using a radial basis function to solve Poisson equations and biharmonic equations has been researched by Mai-Duy et al. [9–18]. Haar wavelets are made up of pairs of piecewise constant functions and mathematically the simplest orthonormal wavelets with a compact support. Due to mathematical simplicity the Haar wavelet method has turned out to be an effective tool for solving differential and integral equations. Lepik [19–21] uses Haar wavelets to solve differential and integral equations. In this paper, we replace the radial basis function with a Haar wavelet to solve Poisson equations and biharmonic equations. The efficiency of the method is demonstrated by four numerical examples. The paper is organized as follows. In Section 2, Haar wavelets and their integrals are introduced. In Section 3, Haar wavelet solutions of Poisson equations are presented. In Section 4, the Haar wavelet solution of biharmonic equations on a rectangle is discussed. A conclusion is made in Section 5. 2. Haar wavelets and their integrals [19] Let us define M = 2J , where J is the maximal of resolution. Next two parameters are introduced: the dilatation parameter j = 0, 1, . . . , J and the translation parameter k = 0, 1, . . . , 2j − 1. The Haar wavelets family is defined for [A, B] as follows
hi ( x ) =
∗
1, −1, 0,
x ∈ [ξ1 (i), ξ2 (i)), x ∈ [ξ2 (i), ξ3 (i)), elsewhere,
Corresponding author. Tel.: +86 029 82205353. E-mail address: [email protected] (Z. Shi).
0895-7177/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2011.07.006
(2.1)
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
2859
where
ξ1 (i) = A + 2kµ1x, ξ2 (i) = A + (2k + 1)µ1x, ξ3 (i) = A + 2(2k + 1)µ1x, µ = M /2j , i = 2j + k + 1, 1x = (B − A)/(2M ). The case i = 1 corresponds to the scaling function h1 (x) =
1, 0,
A ⩽ x < B, elsewhere.
(2.2)
Any function y(x) ∈ L2 [0, 1) can be decomposed as y(x) =
+∞ −
ai hi (x),
i=1
where the coefficients ai are determined by ai = 2j
B
∫
y(x)hi (x)dx. A
B
Specially, a0 = A y(x)dx. The series expansion of y(x) contains infinite terms. If y(x) is piecewise constant by itself, or may be approximated as piecewise constant during each subinterval, then y(x) will be terminated with finite terms, that is y(x) =
2M −
ci hi (x).
i=1
In Section 3–4, the integrals pn,i (x) =
∫ x∫ A
x
... A
∫
x
hi (t )dt n =
A
n-times
1
(n − 1)!
x
∫
(x − t )n−1 hi (t )dt A
are needed, where i = 1, 2, . . . , 2M. The case n = 0 corresponds to function hi (t ). Taking account of (2.1) and (2.2) these integrals can be calculated analytically, by doing it we obtain
pn,i (x) =
0, 1 n n! [x − ξ1 (i)] ,
x < ξ1 (i), x ∈ [ξ1 (i), ξ2 (i)],
1
(2.3)
(x − A)n .
(2.4)
{[x − ξ1 (i)]n − 2[x − ξ2 (i)]n }, x ∈ [ξ2 (i), ξ3 (i)], n ! 1 {[x − ξ (i)]n − 2[x − ξ (i)]n + [x − ξ (i)]n }, x > ξ (i). 1 2 3 3 n! These formulas hold for i > 1. In the case i = 1 we have ξ1 = A, ξ2 = ξ3 = B and pn,1 (x) =
1 n!
In the present paper, the collocation method for solving partial differential equations is applied. Consider a two-dimensional domain: A ⩽ x, y ⩽ B. The domain of interest is represented by the collocation points in each coordinate direction xk = A + (k − 0.5)1x,
yk = A + (k − 0.5)1y,
k = 1, 2, . . . , 2M .
(2.5)
where 1x = 1y = (B − A)/(2M ). 3. Haar wavelets for solving Poisson equations on a rectangle In mathematics, a Poisson equation is a partial differential equation of elliptic type with broad utility in electrostatics. In two-dimensional Cartesian coordinates, a Poisson equation takes the form
∂ 2u ∂ 2u + 2 = f (x, y), ∂ x2 ∂y where f and u are real or complex-valued functions on a manifold.
2860
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
Problem 1. Consider Poisson equation
∂ 2u ∂ 2u + 2 = sin(π x) sin(π y), 0 ⩽ x, y ⩽ 1, ∂ x2 ∂y u(x, y) = 0 along the boundaries,
(3.1)
which has the exact solution ue (x, y) = −1/(2π 2 ) sin(π x) sin(π y). Let n = 2M. The variable v(x) = u(x, y) (u as a function of one variable at a time, y fixed) and its derivative along a grid line that runs parallel to the x-axis can be approximated by d2 v(x)
n −
=
dx2
ai hi (x).
(3.2)
i=1
By integrating (3.2), we obtain dv(x) dx
=
n −
ai p1,i (x) +
dv(0) dx
i =1 n −
v(x) =
ai p2,i (x) +
dv(0)
i =1
dx
,
(3.3)
x + v(0),
(3.4)
where p1,i and p2,i see (2.3)–(2.4). The presence of two integration constants allow the addition of two extra equations. The evaluation of (3.4) at collocation points leads to
[ ]
a , b
v = P1
(3.5)
the two additional equations are written as
] [ ] v(0) a = = P2 , v(1) b [
vextra
(3.6)
where v = (v(x1 ), v(x2 ), . . . , v(xn ))T , p2,1 (x1 ) p2,1 (x2 )
P1 =
.. . p2,1 (xn )
p2,1 (0) p2,1 (1)
p2,n (x1 ) p2,n (x2 )
··· ··· .. .
.. . p2,2 (xn )
p2,2 (0) p2,2 (1)
[ P2 =
p2,2 (x1 ) p2,2 (x2 )
a = (a1 , a2 , . . . , an )T ,
.. . p2,n (xn )
···
p2,n (0) p2,n (1)
··· ···
x1 x2
b=
dv(0) dx
, v(0)
T
,
1 1
.. , .
.. .
xn
1
]
0 1
1 . 1
According to (3.5) and (3.6), we obtain
[
]
v
vextra
[ ][ ] =
P1 P2
[ ]
a a =C , b b
where C = (P1 , P2 )T . Therefore
[ ]
[
]
a v = C −1 . b vextra
(3.7)
By substituting (3.7) into (3.2), the second derivative of the variable v can be expressed in terms of nodal values d2 v dx2
= HC −1
[
v
fextra
]
= B1 v + B2 vextra ,
(3.8)
where B1 and B2 are matrices that are formed by the first n columns and the last two columns of the HC −1 ,
(
d2 v(x1 ) dx2
,...,
d2 u(xn ) T dx2
),
h1 ( x 1 ) h1 (x2 )
H =
.. . h1 (xn )
h2 (x1 ) h2 (x2 )
.. . h2 ( x n )
··· ··· .. . ···
hn (x1 ) hn (x2 )
.. . hn (xn )
0 0
0 0
0
0
.. .
.. . .
d2 v dx2
=
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
2861
Fig. 1. Haar wavelet solution for Problem 1, J = 6.
Table 1 Error estimates for Problem 1. J
eJ
3 4 5 6
3.0827e−004 8.0287e−005 2.0276e−005 5.0820e−006
2
Now we can extend the approximation for ddx2v to a 2D domain. Assuming that the grids are numbered from bottom to top and from left to right, we can write the values of the derivatives of u over the whole domain by using the Kronecker tensor products as follows
∂2 u = (B1 ⊗ Iy )u + (B2 ⊗ Iy )f = Dx u + kx , ∂x2 where Dx = (B1 ⊗ Iy ), kx = (B2 ⊗ Iy )f , f is a zero vector of length 2n, Iy is the identity matrix.
(3.9)
Similarly, we have
∂2 u = (Ix ⊗ B1 )u + (Ix ⊗ B2 )f = Dy u + ky ∂y 2
(3.10)
where Dy = (Ix ⊗ B1 ), ky = (Ix ⊗ B2 )f , f is a zero vector of length 2n, Ix is the identity matrix. Applying the Eq. (3.1) at the inter points(ip) yields
(Dx + Dy )ip u = sin(π xk ) sin(π yk ) − kx − ky . The solution is plotted in Fig. 1. The accuracy of the results can be estimated by error function eJ =
max |u(xk , xl ) − ue (xk , xl )|.
1⩽k,l⩽2M
The eJ is shown in Table 1. Problem 2. Consider Poisson equation
∂ 2u ∂ 2u + 2 = −8 sin(2π x) sin(2π y), 0 ⩽ x, y ⩽ 1, ∂ x2 ∂y f1 (x = 0, y) = ux (0, y) = 2π sin(2π y), f2 (x = 1, y) = ux (1, y) = 2π sin(2π y), f3 (y = 0, x) = 0, f4 (y = 1, x) = 0. The exact solution is given by ue (x, y) = sin(2π x) sin(2π y).
(3.11)
2862
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868 Table 2 Error estimates for Problem 2. J
eJ
3 4 5 6
0.0175 0.0051 0.0013 3.6099e−004
Fig. 2. Haar wavelet solution for Problem 2.
Similar to method of the Problem 1, along a horizontal line, the additional equation can be written as ux (0, y) a = P1 , ux (1, y) b1
[ fextra1 =
]
[ ]
along a vertical line, the additional equation can be written as u(x, 0) a = P2 , u(x, 1) b2
[ fextra2 =
]
[ ]
where
[
p1,1 (0) P1 = p1,1 (1)
p1,2 (0) p1,2 (1)
··· ···
p1,n (0) p1,n (1)
1 1
0 , 0
p2,1 (0) P2 = p2,1 (1)
p2,2 (0) p2,2 (1)
··· ···
p2,n (0) p2,n (1)
0 1
1 , 1
[
a = (a1 , a2 , . . . , an )T ,
] ]
b1 = (ux (0, y), u(0, y))T ,
b2 = (uy (x, 0), u(x, 0))T .
The final form can be written as
∂2 u = (B1 ⊗ Iy )u + (B2 ⊗ Iy )f1 = Dx u + kx , ∂ x2 ∂2 u = (Ix ⊗ B1 )u + (Ix ⊗ B2 )f2 = Dy u + ky , ∂y 2 where Dx = B1 ⊗ Iy , Dy = Ix ⊗ B1 , kx = (B2 ⊗ Iy )f1 , ky = (Ix ⊗ B2 )f2 , f1 = (2π sin(2π x1 ), . . . , 2π sin(2π xn ), 2π sin(2π x1 ), . . . , 2π sin(2π xn ))T , f2 is a zero vector of length 2n, Applying the Eq. (3.11) at the inter points(ip) yields
(Dx + Dy )ip u = −8 sin(2π xk ) sin(2π yk ) − kx − ky . The solution is plotted in Fig. 2. The accuracy of the results can be estimated by error function eJ =
max |u(xk , xl ) − ue (xk , xl )|.
1⩽k,l⩽2M
The eJ is shown in Table 2.
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
2863
4. Haar wavelets solution of biharmonic equations on a rectangle Consider the 2D biharmonic equation with Dirichlet boundary condition of the first kind, i.e.
∂ 4u ∂ 4u ∂ 4u + 2 2 2 + 4 = f (x, y), 4 ∂x ∂x ∂y ∂y ∂u u = f1 (x, y), = f2 (x, y), 0 ⩽ x, y ⩽ 1. ∂n ∆4 u =
Problem 3. We consider the 2D biharmonic equation
∂ 4u ∂ 4u ∂ 4u + 2 2 2 + 4 = 64π 4 sin(2π x) sin(2π y), 4 ∂x ∂x ∂y ∂y u = 0 along the boundaries, ux = 2π sin(2π y), x = 0, 1, uy = 2π sin(2π x), y = 0, 1. ∆4 u =
(4.1)
The exact solution is given by ue (x, y) = sin(2π x) sin(2π y). Let n = 2M. The variable v(x) = u(x, y) (u as a function of one variable at a time, y fixed) and its derivative along a grid line that runs parallel to the x- axis can be approximated by using Haar wavelets d4 v(x)
=
dx4
n −
ai hi (x),
(4.2)
i =1
d3 v(x)
=
dx3
n −
ai p1,i (x) +
d3 v(0) dx3
i =1
d2 v(x)
=
dx2
n −
ai p2,i (x) +
d3 v(0) dx3
i =1
dv(x) dx
=
n −
ai p3,i (x) +
i=1 n −
v(x) =
ai p4,i (x) +
i=1
, x+
1 d3 v(0) dx3
2
1 d3 v(0) dx3
6
(4.3) d2 v(0)
x2 +
x3 +
,
dx2
(4.4)
d2 v(0) dx2
1 d2 v(0) 2
dx2
x+
dv(0) dx
x2 +
,
dv(0) dx
(4.5)
x + v(0).
(4.6)
The evaluation of (4.6) at collocation points leads to
[ ]
a , b
v = P1
(4.7)
the four additional equations can be written as
[ ] vextra = P2
a , b
(4.8)
where v = (v(x1 ), v(x2 ), . . . , v(xn ))T ,
b= vextra
a = (a1 , a2 , . . . , an )T , T d3 v(0) d2 v(0) dv(0) , , , v( 0 ) . 3 2 dx
dx
dx
T dv(0) dv(1) = v(0), v(1), , , dx
p4,1 (x1 )
p (x ) 4,1 2 P1 = .. . p4,1 (xn )
p4,2 (x1 )
dx
···
p4,n (x1 )
p4,2 (x2 )
···
p4,n (x2 )
.. .
..
.. .
p4,2 (xn )
···
.
p4,n (xn )
1 6 1 6 1 6
x31 x32
.. . x3n
1 2 1 2 1 2
x21 x22
.. . x2n
x1 x2
.. . xn
1
1 , .. . 1
2864
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
p4,1 (0)
p4,2 (0)
···
p4,n (0)
p (1) 4,1 P2 = p3,1 (0) p3,1 (1)
p4,2 (1)
p4,n (1)
p3,2 (0)
··· ···
p3,n (0)
p3,2 (1)
···
p3,n (1)
0 1
0 1
6 0 1
2 0
1
1
1
2
0
1
1
1
. 0 0
According to (4.7) and (4.8), we have
[
]
v
vextra
[ ][ ]
[ ]
a a =C . b b
P = 1 P2
Therefore
[ ]
[
]
a v = C −1 . b vextra
(4.9)
Substituting (4.9) into (4.2), the four derivatives of the variable v can be expressed in terms of nodal values d4 v dx4
= HC
[
−1
]
v
vextra
= B1 v + B2 vextra ,
where B1 and B2 are matrices that are formed by the first n columns and the last four columns of the HC −1 ,
(
d4 v(x1 ) dx4
d4 v(xn ) T dx4
,...,
d4 v dx4
=
) and
h1 ( x 1 ) h1 (x2 )
h2 (x1 ) h2 (x2 )
H =
.. . h1 (xn )
hn (x1 ) hn (x2 )
··· ··· .. .
.. . h2 ( x n )
.. . hn (xn )
···
0 0
.. .
0
0 0
.. . .
(4.10)
0
Similar to method of the Problem 1, we have
∂4 u = (B1 ⊗ Iy )u + (B2 ⊗ Iy )f1 = Dx u + kx , ∂ x4
(4.11)
where Dx = (B1 ⊗ Iy ), kx = (B2 ⊗ Iy )f1 , f1 = (0, . . . , 0, 0, . . . , 0, 2π sin(2π x1 ), . . . , 2π sin(2π xn ), 2π sin(2π x1 ), . . . , 2π sin(2π xn ))T . Iy is the identity matrix. Similarly
∂4 u = Dy u + ky , ∂y 4
(4.12)
where Dy = (Ix ⊗ B1 ), ky = (Ix ⊗ B2 )f2 , Ix is the identity matrix, f2 = (0, 0, 2π sin(2π y1 ), 2π sin(2π y1 ), 0, 0, 2π sin(2π y2 ), 2π sin(2π y2 ), . . . , 0, 0, 2π sin(2π yn ), 2π sin(2π yn ))T . Making use of the following expression 2
∂ 4u ∂2 = 2 2 2 ∂x ∂y ∂x
∂ 2u ∂ y2
+
∂2 ∂ y2
∂ 2u , ∂ x2
(4.13)
the mixed fourth-order partial derivative can be computed. The variable k(x) = ∂∂ y2u ( ∂∂ y2u as a function of one variable at a time, y fixed) and its derivative along a grid line that runs parallel to the x-axis can be approximated by using Haar wavelets, that is 2
d2 k(x) dx2 dk(x) dx
=
n −
2
ci hi (x),
i=1
=
k(x) =
n −
ci p1,i (x) +
i =1 n − i =1
ci p2,i (x) +
dk(0) dx
dk(0) dx
,
x + k(0).
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
2865
The evaluation of k(x) at collocation points leads to k(x1 ) k(x2 )
p2,1 (x1 ) p2,1 (x2 )
p2,2 (x1 ) p2,2 (x2 )
.. = .. . . k(xn ) p2,1 (xn ) [ ] c = Q1 ,
k =
.. . p2,2 (xn )
p2,n (x1 ) p2,n (x2 )
··· ··· .. .
.. . p2,n (xn )
···
x1 x2
1 1 [ c ]
xn
1
.. .
.. d .
(4.14)
d
the two additional equations can be written as k(0)
]
[ kextra =
p2,1 (0) p2,1 (1)
p2,2 (0) p2,2 (1)
[ =
k(1)
p2,n (0) p2,n (1)
··· ···
0 1
1 1
][ ] c d
[ ]
c , d
= Q2
(4.15) dk(0)
where c = (c1 , c2 , . . . , cn )T , d = ( dx , k(0))T . According to (4.14) and (4.15), we obtain
[
]
k
[
Q1 Q2
=
kextra
][ ] c d
[ ]
c . d
= C2
(4.16)
Therefore
[ ] c d
−1
[
= C2
]
k kextra
,
(4.17)
and
d2 k(x ) 1 2dx2 d k(x2 ) h1 (x1 ) d2 k dx2 .. = = . .. dx2 . h1 (xn ) 2 d k(xn )
h2 (x1 )
.. . h2 (xn )
··· .. . ···
hn ( x 1 )
.. . hn ( x n )
0
0
0
0
[ ] .. c . d
.. .
dx2
= H C2 −1
[
k kextra
]
= A1 k + A2 kextra ,
(4.18)
where A1 and A2 are matrices that are formed by the first n columns and the last two columns of the HC2−1 . Let l(y) = u(x, y) (u as a function of one variable at a time, x fixed). Similar to (4.14), the variable l(y) and its derivative along a grid line that runs parallel to the y-axis can be approximated by using a Haar wavelet, that is l(y1 ) l(y2 )
p2,1 (y1 ) p2,1 (y2 )
l = .. = .
l(yn )
.. . p2,1 (yn )
p2,2 (y1 ) p2,2 (y2 )
.. . p2,2 (yn )
··· ··· .. . ···
p2,n (y1 ) p2,n (y2 )
.. . p2,n (yn )
y1 y2
1 1 [a]
yn
1
.. .
.. b .
[ ] = P1
a , b
(4.19) dl(0)
where a = (a1 , . . . , an )T , b = ( dy , l(0))T . The two additional equations can be written as
[ lextra =
l(0)
p (0) = 2,1 p2,1 (1) l(1)
]
[
p2,2 (0) p2,2 (1)
··· ···
p2,n (0) p2,n (1)
0 1
][ ]
1 1
a b
[ ] = P2
a . b
(4.20)
2866
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
Fig. 3. Haar wavelet solution for Problem 3.
According to (4.19)–(4.20), we have
[
]
l
lextra
[ ][ ] P1 P2
=
[ ]
a a = C1 b b
(4.21)
where C1 = [P1 , P2 ]T . Similar to (4.18), we have
d2 l(y1 )
dy2 d2 l(y2 ) h1 (y1 ) 2 d l .. 2 = dy = . . dy 2 .. h1 (yn ) d2 l(y ) n
h2 (y1 )
.. . h2 (yn )
··· .. . ···
hn (y1 )
.. . hn (yn )
0
0
0
0
.. .
[ ] .. a . b
dy2
= H C1
−1
[
l
lextra
]
= B1 l + B2 lextra ,
(4.22)
where B1 and B2 are matrices that are formed by the first n columns and the last two columns of the H C1 −1 . According (4.18) and (4.22), the final form can be written as
∂2 ∂x2
∂2 u ∂y 2
= (A1 ⊗ Iy )
∂2 u + (A2 ⊗ Iy )g1 ∂y 2
= (A1 ⊗ Iy ) ((Ix ⊗ B1 )u + (Ix ⊗ B2 )g2 ) + (A2 ⊗ Iy )g1 = (A1 ⊗ Iy )(Ix ⊗ B1 )u + (A1 ⊗ Iy )(Ix ⊗ B2 )g2 + (A2 ⊗ Iy )g1 = Dxy u + mxy , where Dxy = (A1 ⊗ Iy )(Ix ⊗ B1 ), mxy = (A1 ⊗ Iy )(Ix ⊗ B2 )g2 + (A2 ⊗ Iy )g1 , g1 and g2 are zero vectors of length 2n. For
∂2 ∂2 u ( ), ∂y 2 ∂x2
we can obtain a similar expression. Applying the Eq. (4.1) at the inter points(ip) yields
(Dx + 2Dxy + Dy )ip u = 64π 4 sin(2π xk ) sin(2π yk ) − kx − ky − mxy − myx . The solution is plotted in Fig. 3. The accuracy of the results can be estimated by error function eJ =
max |u(xk , xl ) − ue (xk , xl )|.
1⩽k,l⩽2M
The eJ is shown in Table 3.
(4.23)
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
2867
Table 3 Error estimates for Problem 3. J
eJ
3 4 5 6
0.0261 0.0051 0.0010 2.2875e−004
Fig. 4. Haar wavelet solution for Problem 4.
Problem 4. Consider a biharmonic equation in the form
∂ 4u ∂ 4u ∂ 4u + 2 + = 0, ∂ x4 ∂ x2 ∂ y2 ∂ y4
(x, y) = [−1, 1] × [−1, 1],
(4.24)
with the boundary conditions u(x, −1) = u(x, 1) =
sin(π x)
π2
,
(4.25)
∂ 2 u(x, 1) sin(π x)(2π 2 − π sinh(2π )) ∂ 2 u(x, −1) = = . 2 2 ∂y ∂y 2π 2 + π sinh(2π )
(4.26)
The exact condition of Eqs. (4.24)–(4.26) is ue (x, y) = sin(π x)(h1 cosh(π y) + h2 y sinh(π y)), where h1 =
2 sinh(π)+2π cosh(π) 2π 3 +pi2 sinh(2π)
, h2 =
−2 sinh(pi) 2π 2 +π sinh(2π)
.
Let x → x − 1, y → y − 1, Eqs. (4.24)–(4.26) becomes
∂ 4u ∂ 4u ∂ 4u + 2 2 2 + 4 = 0, 4 ∂x ∂x ∂y ∂y
(x, y) = [0, 2] × [0, 2],
with the boundary conditions u(x, 0) = u(x, 2) =
sin(π (x − 1))
π2
,
∂ u(x, 0) ∂ u(x, 2) sin(π (x − 1))(2π 2 − π sinh(2π )) = = . ∂ y2 ∂ y2 2π 2 + π sinh(2π ) 2
2
Consider a wavelet with support over 0 ⩽ x ⩽ 2, using the method of Problem 3, we obtain an approximate solution. The result is plotted in Fig. 4. The accuracy of the results can be estimated by error function eJ =
max |u(xk , xl ) − ue (xk , xl )|.
1⩽k,l⩽2M
The eJ is shown in Table 4.
2868
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868 Table 4 Error estimates for Problem 4. J
eJ
3 4 5 6
0.0022 0.0010 3.2241e−004 9.0296e−005
5. Conclusion In the paper, we replace radial basis functions with Haar wavelets to solve Poisson equations and biharmonic equations. The reason for using Haar wavelets are sparse representation, fast transformation, and the possibility of implementing fast algorithms. An extension of the proposed method to the case of non-rectangular domains and using other wavelets (Legendre wavelet, Shannon wavelet, Harmonic wavelet) to solve biharmonic equations will be reported in future work. Acknowledgment This work is supported by the Specialized Research Fund of Shaanxi Provincial Department of Education of China under grant No. 09JK539 and Natural Science Foundation of Shanxi Province of China under Grant No. 2009JM1002. References [1] Y kwon, J.W. Stephenson, Single cell finite difference approximation for Poisson equation in three variables, Appl. Math. Notes 2 (1982) 13–17. [2] M.M. Gupta, A fourth-order Poisson solver, J. Comput. Phys. 55 (1985) 166–172. [3] J. Wang, W. Zhong, J. Zhang, A general meshsize fourth-order compact difference discretization scheme for 3D Poisson equation, Appl. Math. Comput. 183 (2006) 804–812. [4] Xuan Wang, Zhifeng Yang, Gordon Huang, Bin Chen, A high-order compact difference scheme for 2D Laplace and Poisson equations in non-uniform grid systems, Commun. Nonlinear Sci. Numer. Simul. 14 (2) (2009) 379–398. [5] Yongbin Ge, Multigrid method and fourth-order compact difference discretization scheme with unequal mesh sizes for 3D poisson equation, J. Comput. Phys. 229 (18) (2010) 6381–6391. [6] E. Perrey-Debain, H.G. ter Morsche, B-Spline approximation and fast wavelet transform for an efficient evaluation of particular solutions for Poisson’s equation, Eng. Anal. Bound. Elem. 26 (1) (2002) 1–13. [7] K.E. Atkinson, The numerical evaluation of particular solution for Poisson’s equation, IMA J. Numer. Anal. 5 (1985) 319–338. [8] M.A. Golberg, C.S. Chen, The theory of radial basis functions applied to The BEM for inhomogeneous partial differential equations, Bound Elem. Commun. 5 (1994) 57–61. [9] Mehdi Dehghan, Akbar Mohebbi, Multigrid solution of high order discretisation for three-dimensional biharmonic equation with Dirichlet boundary conditions of second kind, Appl. Math. Comput. 180 (2006) 575–593. [10] Youngmok Jeon, New indirect scalar boundary integral equation formulas for the biharmonic equation, J. Comput. Appl. Math. 135 (2) (2001) 313–324. [11] Quang A. Dang, Iterative method for solving the Neumann boundary value problem for biharmonic type equation, J. Comput. Appl. Math. 196 (2) (2006) 634–643. [12] N. Mai-Duy, R.I. Tanner, A spectral collocation method based on integrated Chebyshev polynomials for two-dimensional biharmonic boundary-value problems, J. Comput. Appl. Math. 201 (2007) 30–47. [13] N. Mai-Duy, R.I. Tanner, Solving high-order partial differential equations with indirect radial basis function networks, Int. J. Numer. Meth. Engng. 63 (2005) 1636–1654. [14] N. Mai-Duy, T. Tran-Cong, Solving biharmonic problems with scattered-point discretization using indirect radial-basis-function networks, Eng. Anal. Bound. Elem. 30 (2) (2006) 77–87. [15] N. Mai-Duy, R.I. Tanner, An effective high order interpolation scheme in BIEM for biharmonic boundary value problems, Eng. Anal. Bound. Elem. 29 (3) (2005) 210–223. [16] Xiaolin Li, Jialin Zhu, A Galerkin boundary node method for biharmonic problems, Eng. Anal. Bound. Elem. 33 (6) (2009) 858–865. [17] N. Mai-Duy, H. See, T. Tran-Cong, A spectral collocation technique based on integrated Chebyshev polynomials for biharmonic problems in irregular domains, Appl. Math. Model. 33 (1) (2009) 284–299. [18] Bernard Bialecki, A fast solver for the orthogonal spline collocation solution of the biharmonic Dirichlet problem on rectangles, J. Comput. Phys. 191 (2) (2003) 601–621. [19] ü. Lepik, Haar wavelet method for solving higher order differential equation, Internat. J. Math. Comput. 1 (8) (2008) 84–94. [20] ü. Lepik, Solving integral and differential equations by the aid of non-uniform Haar wavelets, Appl. Math. Comput. 198 (1) (2008) 326–332. [21] ü. Lepik, Solving fractional integral equations by the Haar wavelet method, Appl. Math. Comput. 214 (2) (2009) 468–478.
Contents lists available at SciVerse ScienceDirect
Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm
A spectral collocation method based on Haar wavelets for Poisson equations and biharmonic equations Shi Zhi ∗ , Cao Yong-yan School of Science, Xi’an University of Architecture and Technology, Xi’an 710055, China
article
info
Article history: Received 27 February 2011 Received in revised form 4 July 2011 Accepted 7 July 2011 Keywords: Haar wavelet Poisson equation Biharmonic equation
abstract In this work, we present a computational method for solving Poisson equations and biharmonic equations which are based on the use of Haar wavelets. The first transform the spectral coefficients into the nodal variable values. The second use Kronecker products to construct the approximations for derivatives over a tensor product grid of the horizontal and vertical blocks. Finally, solve the obtained system of algebraic equations. The efficiency of the method is demonstrated by four numerical examples. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction The numerical solution of Poisson equations is an important problem in numerical analysis. In the past two decades, a great deal of research work has been published on the development of numerical solution of Poisson equations [1–6]. A number of alternatives to direct integration of the domain integral have been proposed, notably among them the important papers of Atkinson [7] and Golberg and Chen [8]. In mathematics, a biharmonic equation is a fourth-order partial differential equation which arises in areas of continuum mechanics, including linear elasticity theory and the solution of Stokes flows. The concept of using a radial basis function to solve Poisson equations and biharmonic equations has been researched by Mai-Duy et al. [9–18]. Haar wavelets are made up of pairs of piecewise constant functions and mathematically the simplest orthonormal wavelets with a compact support. Due to mathematical simplicity the Haar wavelet method has turned out to be an effective tool for solving differential and integral equations. Lepik [19–21] uses Haar wavelets to solve differential and integral equations. In this paper, we replace the radial basis function with a Haar wavelet to solve Poisson equations and biharmonic equations. The efficiency of the method is demonstrated by four numerical examples. The paper is organized as follows. In Section 2, Haar wavelets and their integrals are introduced. In Section 3, Haar wavelet solutions of Poisson equations are presented. In Section 4, the Haar wavelet solution of biharmonic equations on a rectangle is discussed. A conclusion is made in Section 5. 2. Haar wavelets and their integrals [19] Let us define M = 2J , where J is the maximal of resolution. Next two parameters are introduced: the dilatation parameter j = 0, 1, . . . , J and the translation parameter k = 0, 1, . . . , 2j − 1. The Haar wavelets family is defined for [A, B] as follows
hi ( x ) =
∗
1, −1, 0,
x ∈ [ξ1 (i), ξ2 (i)), x ∈ [ξ2 (i), ξ3 (i)), elsewhere,
Corresponding author. Tel.: +86 029 82205353. E-mail address: [email protected] (Z. Shi).
0895-7177/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2011.07.006
(2.1)
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
2859
where
ξ1 (i) = A + 2kµ1x, ξ2 (i) = A + (2k + 1)µ1x, ξ3 (i) = A + 2(2k + 1)µ1x, µ = M /2j , i = 2j + k + 1, 1x = (B − A)/(2M ). The case i = 1 corresponds to the scaling function h1 (x) =
1, 0,
A ⩽ x < B, elsewhere.
(2.2)
Any function y(x) ∈ L2 [0, 1) can be decomposed as y(x) =
+∞ −
ai hi (x),
i=1
where the coefficients ai are determined by ai = 2j
B
∫
y(x)hi (x)dx. A
B
Specially, a0 = A y(x)dx. The series expansion of y(x) contains infinite terms. If y(x) is piecewise constant by itself, or may be approximated as piecewise constant during each subinterval, then y(x) will be terminated with finite terms, that is y(x) =
2M −
ci hi (x).
i=1
In Section 3–4, the integrals pn,i (x) =
∫ x∫ A
x
... A
∫
x
hi (t )dt n =
A
n-times
1
(n − 1)!
x
∫
(x − t )n−1 hi (t )dt A
are needed, where i = 1, 2, . . . , 2M. The case n = 0 corresponds to function hi (t ). Taking account of (2.1) and (2.2) these integrals can be calculated analytically, by doing it we obtain
pn,i (x) =
0, 1 n n! [x − ξ1 (i)] ,
x < ξ1 (i), x ∈ [ξ1 (i), ξ2 (i)],
1
(2.3)
(x − A)n .
(2.4)
{[x − ξ1 (i)]n − 2[x − ξ2 (i)]n }, x ∈ [ξ2 (i), ξ3 (i)], n ! 1 {[x − ξ (i)]n − 2[x − ξ (i)]n + [x − ξ (i)]n }, x > ξ (i). 1 2 3 3 n! These formulas hold for i > 1. In the case i = 1 we have ξ1 = A, ξ2 = ξ3 = B and pn,1 (x) =
1 n!
In the present paper, the collocation method for solving partial differential equations is applied. Consider a two-dimensional domain: A ⩽ x, y ⩽ B. The domain of interest is represented by the collocation points in each coordinate direction xk = A + (k − 0.5)1x,
yk = A + (k − 0.5)1y,
k = 1, 2, . . . , 2M .
(2.5)
where 1x = 1y = (B − A)/(2M ). 3. Haar wavelets for solving Poisson equations on a rectangle In mathematics, a Poisson equation is a partial differential equation of elliptic type with broad utility in electrostatics. In two-dimensional Cartesian coordinates, a Poisson equation takes the form
∂ 2u ∂ 2u + 2 = f (x, y), ∂ x2 ∂y where f and u are real or complex-valued functions on a manifold.
2860
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
Problem 1. Consider Poisson equation
∂ 2u ∂ 2u + 2 = sin(π x) sin(π y), 0 ⩽ x, y ⩽ 1, ∂ x2 ∂y u(x, y) = 0 along the boundaries,
(3.1)
which has the exact solution ue (x, y) = −1/(2π 2 ) sin(π x) sin(π y). Let n = 2M. The variable v(x) = u(x, y) (u as a function of one variable at a time, y fixed) and its derivative along a grid line that runs parallel to the x-axis can be approximated by d2 v(x)
n −
=
dx2
ai hi (x).
(3.2)
i=1
By integrating (3.2), we obtain dv(x) dx
=
n −
ai p1,i (x) +
dv(0) dx
i =1 n −
v(x) =
ai p2,i (x) +
dv(0)
i =1
dx
,
(3.3)
x + v(0),
(3.4)
where p1,i and p2,i see (2.3)–(2.4). The presence of two integration constants allow the addition of two extra equations. The evaluation of (3.4) at collocation points leads to
[ ]
a , b
v = P1
(3.5)
the two additional equations are written as
] [ ] v(0) a = = P2 , v(1) b [
vextra
(3.6)
where v = (v(x1 ), v(x2 ), . . . , v(xn ))T , p2,1 (x1 ) p2,1 (x2 )
P1 =
.. . p2,1 (xn )
p2,1 (0) p2,1 (1)
p2,n (x1 ) p2,n (x2 )
··· ··· .. .
.. . p2,2 (xn )
p2,2 (0) p2,2 (1)
[ P2 =
p2,2 (x1 ) p2,2 (x2 )
a = (a1 , a2 , . . . , an )T ,
.. . p2,n (xn )
···
p2,n (0) p2,n (1)
··· ···
x1 x2
b=
dv(0) dx
, v(0)
T
,
1 1
.. , .
.. .
xn
1
]
0 1
1 . 1
According to (3.5) and (3.6), we obtain
[
]
v
vextra
[ ][ ] =
P1 P2
[ ]
a a =C , b b
where C = (P1 , P2 )T . Therefore
[ ]
[
]
a v = C −1 . b vextra
(3.7)
By substituting (3.7) into (3.2), the second derivative of the variable v can be expressed in terms of nodal values d2 v dx2
= HC −1
[
v
fextra
]
= B1 v + B2 vextra ,
(3.8)
where B1 and B2 are matrices that are formed by the first n columns and the last two columns of the HC −1 ,
(
d2 v(x1 ) dx2
,...,
d2 u(xn ) T dx2
),
h1 ( x 1 ) h1 (x2 )
H =
.. . h1 (xn )
h2 (x1 ) h2 (x2 )
.. . h2 ( x n )
··· ··· .. . ···
hn (x1 ) hn (x2 )
.. . hn (xn )
0 0
0 0
0
0
.. .
.. . .
d2 v dx2
=
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
2861
Fig. 1. Haar wavelet solution for Problem 1, J = 6.
Table 1 Error estimates for Problem 1. J
eJ
3 4 5 6
3.0827e−004 8.0287e−005 2.0276e−005 5.0820e−006
2
Now we can extend the approximation for ddx2v to a 2D domain. Assuming that the grids are numbered from bottom to top and from left to right, we can write the values of the derivatives of u over the whole domain by using the Kronecker tensor products as follows
∂2 u = (B1 ⊗ Iy )u + (B2 ⊗ Iy )f = Dx u + kx , ∂x2 where Dx = (B1 ⊗ Iy ), kx = (B2 ⊗ Iy )f , f is a zero vector of length 2n, Iy is the identity matrix.
(3.9)
Similarly, we have
∂2 u = (Ix ⊗ B1 )u + (Ix ⊗ B2 )f = Dy u + ky ∂y 2
(3.10)
where Dy = (Ix ⊗ B1 ), ky = (Ix ⊗ B2 )f , f is a zero vector of length 2n, Ix is the identity matrix. Applying the Eq. (3.1) at the inter points(ip) yields
(Dx + Dy )ip u = sin(π xk ) sin(π yk ) − kx − ky . The solution is plotted in Fig. 1. The accuracy of the results can be estimated by error function eJ =
max |u(xk , xl ) − ue (xk , xl )|.
1⩽k,l⩽2M
The eJ is shown in Table 1. Problem 2. Consider Poisson equation
∂ 2u ∂ 2u + 2 = −8 sin(2π x) sin(2π y), 0 ⩽ x, y ⩽ 1, ∂ x2 ∂y f1 (x = 0, y) = ux (0, y) = 2π sin(2π y), f2 (x = 1, y) = ux (1, y) = 2π sin(2π y), f3 (y = 0, x) = 0, f4 (y = 1, x) = 0. The exact solution is given by ue (x, y) = sin(2π x) sin(2π y).
(3.11)
2862
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868 Table 2 Error estimates for Problem 2. J
eJ
3 4 5 6
0.0175 0.0051 0.0013 3.6099e−004
Fig. 2. Haar wavelet solution for Problem 2.
Similar to method of the Problem 1, along a horizontal line, the additional equation can be written as ux (0, y) a = P1 , ux (1, y) b1
[ fextra1 =
]
[ ]
along a vertical line, the additional equation can be written as u(x, 0) a = P2 , u(x, 1) b2
[ fextra2 =
]
[ ]
where
[
p1,1 (0) P1 = p1,1 (1)
p1,2 (0) p1,2 (1)
··· ···
p1,n (0) p1,n (1)
1 1
0 , 0
p2,1 (0) P2 = p2,1 (1)
p2,2 (0) p2,2 (1)
··· ···
p2,n (0) p2,n (1)
0 1
1 , 1
[
a = (a1 , a2 , . . . , an )T ,
] ]
b1 = (ux (0, y), u(0, y))T ,
b2 = (uy (x, 0), u(x, 0))T .
The final form can be written as
∂2 u = (B1 ⊗ Iy )u + (B2 ⊗ Iy )f1 = Dx u + kx , ∂ x2 ∂2 u = (Ix ⊗ B1 )u + (Ix ⊗ B2 )f2 = Dy u + ky , ∂y 2 where Dx = B1 ⊗ Iy , Dy = Ix ⊗ B1 , kx = (B2 ⊗ Iy )f1 , ky = (Ix ⊗ B2 )f2 , f1 = (2π sin(2π x1 ), . . . , 2π sin(2π xn ), 2π sin(2π x1 ), . . . , 2π sin(2π xn ))T , f2 is a zero vector of length 2n, Applying the Eq. (3.11) at the inter points(ip) yields
(Dx + Dy )ip u = −8 sin(2π xk ) sin(2π yk ) − kx − ky . The solution is plotted in Fig. 2. The accuracy of the results can be estimated by error function eJ =
max |u(xk , xl ) − ue (xk , xl )|.
1⩽k,l⩽2M
The eJ is shown in Table 2.
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
2863
4. Haar wavelets solution of biharmonic equations on a rectangle Consider the 2D biharmonic equation with Dirichlet boundary condition of the first kind, i.e.
∂ 4u ∂ 4u ∂ 4u + 2 2 2 + 4 = f (x, y), 4 ∂x ∂x ∂y ∂y ∂u u = f1 (x, y), = f2 (x, y), 0 ⩽ x, y ⩽ 1. ∂n ∆4 u =
Problem 3. We consider the 2D biharmonic equation
∂ 4u ∂ 4u ∂ 4u + 2 2 2 + 4 = 64π 4 sin(2π x) sin(2π y), 4 ∂x ∂x ∂y ∂y u = 0 along the boundaries, ux = 2π sin(2π y), x = 0, 1, uy = 2π sin(2π x), y = 0, 1. ∆4 u =
(4.1)
The exact solution is given by ue (x, y) = sin(2π x) sin(2π y). Let n = 2M. The variable v(x) = u(x, y) (u as a function of one variable at a time, y fixed) and its derivative along a grid line that runs parallel to the x- axis can be approximated by using Haar wavelets d4 v(x)
=
dx4
n −
ai hi (x),
(4.2)
i =1
d3 v(x)
=
dx3
n −
ai p1,i (x) +
d3 v(0) dx3
i =1
d2 v(x)
=
dx2
n −
ai p2,i (x) +
d3 v(0) dx3
i =1
dv(x) dx
=
n −
ai p3,i (x) +
i=1 n −
v(x) =
ai p4,i (x) +
i=1
, x+
1 d3 v(0) dx3
2
1 d3 v(0) dx3
6
(4.3) d2 v(0)
x2 +
x3 +
,
dx2
(4.4)
d2 v(0) dx2
1 d2 v(0) 2
dx2
x+
dv(0) dx
x2 +
,
dv(0) dx
(4.5)
x + v(0).
(4.6)
The evaluation of (4.6) at collocation points leads to
[ ]
a , b
v = P1
(4.7)
the four additional equations can be written as
[ ] vextra = P2
a , b
(4.8)
where v = (v(x1 ), v(x2 ), . . . , v(xn ))T ,
b= vextra
a = (a1 , a2 , . . . , an )T , T d3 v(0) d2 v(0) dv(0) , , , v( 0 ) . 3 2 dx
dx
dx
T dv(0) dv(1) = v(0), v(1), , , dx
p4,1 (x1 )
p (x ) 4,1 2 P1 = .. . p4,1 (xn )
p4,2 (x1 )
dx
···
p4,n (x1 )
p4,2 (x2 )
···
p4,n (x2 )
.. .
..
.. .
p4,2 (xn )
···
.
p4,n (xn )
1 6 1 6 1 6
x31 x32
.. . x3n
1 2 1 2 1 2
x21 x22
.. . x2n
x1 x2
.. . xn
1
1 , .. . 1
2864
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
p4,1 (0)
p4,2 (0)
···
p4,n (0)
p (1) 4,1 P2 = p3,1 (0) p3,1 (1)
p4,2 (1)
p4,n (1)
p3,2 (0)
··· ···
p3,n (0)
p3,2 (1)
···
p3,n (1)
0 1
0 1
6 0 1
2 0
1
1
1
2
0
1
1
1
. 0 0
According to (4.7) and (4.8), we have
[
]
v
vextra
[ ][ ]
[ ]
a a =C . b b
P = 1 P2
Therefore
[ ]
[
]
a v = C −1 . b vextra
(4.9)
Substituting (4.9) into (4.2), the four derivatives of the variable v can be expressed in terms of nodal values d4 v dx4
= HC
[
−1
]
v
vextra
= B1 v + B2 vextra ,
where B1 and B2 are matrices that are formed by the first n columns and the last four columns of the HC −1 ,
(
d4 v(x1 ) dx4
d4 v(xn ) T dx4
,...,
d4 v dx4
=
) and
h1 ( x 1 ) h1 (x2 )
h2 (x1 ) h2 (x2 )
H =
.. . h1 (xn )
hn (x1 ) hn (x2 )
··· ··· .. .
.. . h2 ( x n )
.. . hn (xn )
···
0 0
.. .
0
0 0
.. . .
(4.10)
0
Similar to method of the Problem 1, we have
∂4 u = (B1 ⊗ Iy )u + (B2 ⊗ Iy )f1 = Dx u + kx , ∂ x4
(4.11)
where Dx = (B1 ⊗ Iy ), kx = (B2 ⊗ Iy )f1 , f1 = (0, . . . , 0, 0, . . . , 0, 2π sin(2π x1 ), . . . , 2π sin(2π xn ), 2π sin(2π x1 ), . . . , 2π sin(2π xn ))T . Iy is the identity matrix. Similarly
∂4 u = Dy u + ky , ∂y 4
(4.12)
where Dy = (Ix ⊗ B1 ), ky = (Ix ⊗ B2 )f2 , Ix is the identity matrix, f2 = (0, 0, 2π sin(2π y1 ), 2π sin(2π y1 ), 0, 0, 2π sin(2π y2 ), 2π sin(2π y2 ), . . . , 0, 0, 2π sin(2π yn ), 2π sin(2π yn ))T . Making use of the following expression 2
∂ 4u ∂2 = 2 2 2 ∂x ∂y ∂x
∂ 2u ∂ y2
+
∂2 ∂ y2
∂ 2u , ∂ x2
(4.13)
the mixed fourth-order partial derivative can be computed. The variable k(x) = ∂∂ y2u ( ∂∂ y2u as a function of one variable at a time, y fixed) and its derivative along a grid line that runs parallel to the x-axis can be approximated by using Haar wavelets, that is 2
d2 k(x) dx2 dk(x) dx
=
n −
2
ci hi (x),
i=1
=
k(x) =
n −
ci p1,i (x) +
i =1 n − i =1
ci p2,i (x) +
dk(0) dx
dk(0) dx
,
x + k(0).
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
2865
The evaluation of k(x) at collocation points leads to k(x1 ) k(x2 )
p2,1 (x1 ) p2,1 (x2 )
p2,2 (x1 ) p2,2 (x2 )
.. = .. . . k(xn ) p2,1 (xn ) [ ] c = Q1 ,
k =
.. . p2,2 (xn )
p2,n (x1 ) p2,n (x2 )
··· ··· .. .
.. . p2,n (xn )
···
x1 x2
1 1 [ c ]
xn
1
.. .
.. d .
(4.14)
d
the two additional equations can be written as k(0)
]
[ kextra =
p2,1 (0) p2,1 (1)
p2,2 (0) p2,2 (1)
[ =
k(1)
p2,n (0) p2,n (1)
··· ···
0 1
1 1
][ ] c d
[ ]
c , d
= Q2
(4.15) dk(0)
where c = (c1 , c2 , . . . , cn )T , d = ( dx , k(0))T . According to (4.14) and (4.15), we obtain
[
]
k
[
Q1 Q2
=
kextra
][ ] c d
[ ]
c . d
= C2
(4.16)
Therefore
[ ] c d
−1
[
= C2
]
k kextra
,
(4.17)
and
d2 k(x ) 1 2dx2 d k(x2 ) h1 (x1 ) d2 k dx2 .. = = . .. dx2 . h1 (xn ) 2 d k(xn )
h2 (x1 )
.. . h2 (xn )
··· .. . ···
hn ( x 1 )
.. . hn ( x n )
0
0
0
0
[ ] .. c . d
.. .
dx2
= H C2 −1
[
k kextra
]
= A1 k + A2 kextra ,
(4.18)
where A1 and A2 are matrices that are formed by the first n columns and the last two columns of the HC2−1 . Let l(y) = u(x, y) (u as a function of one variable at a time, x fixed). Similar to (4.14), the variable l(y) and its derivative along a grid line that runs parallel to the y-axis can be approximated by using a Haar wavelet, that is l(y1 ) l(y2 )
p2,1 (y1 ) p2,1 (y2 )
l = .. = .
l(yn )
.. . p2,1 (yn )
p2,2 (y1 ) p2,2 (y2 )
.. . p2,2 (yn )
··· ··· .. . ···
p2,n (y1 ) p2,n (y2 )
.. . p2,n (yn )
y1 y2
1 1 [a]
yn
1
.. .
.. b .
[ ] = P1
a , b
(4.19) dl(0)
where a = (a1 , . . . , an )T , b = ( dy , l(0))T . The two additional equations can be written as
[ lextra =
l(0)
p (0) = 2,1 p2,1 (1) l(1)
]
[
p2,2 (0) p2,2 (1)
··· ···
p2,n (0) p2,n (1)
0 1
][ ]
1 1
a b
[ ] = P2
a . b
(4.20)
2866
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
Fig. 3. Haar wavelet solution for Problem 3.
According to (4.19)–(4.20), we have
[
]
l
lextra
[ ][ ] P1 P2
=
[ ]
a a = C1 b b
(4.21)
where C1 = [P1 , P2 ]T . Similar to (4.18), we have
d2 l(y1 )
dy2 d2 l(y2 ) h1 (y1 ) 2 d l .. 2 = dy = . . dy 2 .. h1 (yn ) d2 l(y ) n
h2 (y1 )
.. . h2 (yn )
··· .. . ···
hn (y1 )
.. . hn (yn )
0
0
0
0
.. .
[ ] .. a . b
dy2
= H C1
−1
[
l
lextra
]
= B1 l + B2 lextra ,
(4.22)
where B1 and B2 are matrices that are formed by the first n columns and the last two columns of the H C1 −1 . According (4.18) and (4.22), the final form can be written as
∂2 ∂x2
∂2 u ∂y 2
= (A1 ⊗ Iy )
∂2 u + (A2 ⊗ Iy )g1 ∂y 2
= (A1 ⊗ Iy ) ((Ix ⊗ B1 )u + (Ix ⊗ B2 )g2 ) + (A2 ⊗ Iy )g1 = (A1 ⊗ Iy )(Ix ⊗ B1 )u + (A1 ⊗ Iy )(Ix ⊗ B2 )g2 + (A2 ⊗ Iy )g1 = Dxy u + mxy , where Dxy = (A1 ⊗ Iy )(Ix ⊗ B1 ), mxy = (A1 ⊗ Iy )(Ix ⊗ B2 )g2 + (A2 ⊗ Iy )g1 , g1 and g2 are zero vectors of length 2n. For
∂2 ∂2 u ( ), ∂y 2 ∂x2
we can obtain a similar expression. Applying the Eq. (4.1) at the inter points(ip) yields
(Dx + 2Dxy + Dy )ip u = 64π 4 sin(2π xk ) sin(2π yk ) − kx − ky − mxy − myx . The solution is plotted in Fig. 3. The accuracy of the results can be estimated by error function eJ =
max |u(xk , xl ) − ue (xk , xl )|.
1⩽k,l⩽2M
The eJ is shown in Table 3.
(4.23)
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868
2867
Table 3 Error estimates for Problem 3. J
eJ
3 4 5 6
0.0261 0.0051 0.0010 2.2875e−004
Fig. 4. Haar wavelet solution for Problem 4.
Problem 4. Consider a biharmonic equation in the form
∂ 4u ∂ 4u ∂ 4u + 2 + = 0, ∂ x4 ∂ x2 ∂ y2 ∂ y4
(x, y) = [−1, 1] × [−1, 1],
(4.24)
with the boundary conditions u(x, −1) = u(x, 1) =
sin(π x)
π2
,
(4.25)
∂ 2 u(x, 1) sin(π x)(2π 2 − π sinh(2π )) ∂ 2 u(x, −1) = = . 2 2 ∂y ∂y 2π 2 + π sinh(2π )
(4.26)
The exact condition of Eqs. (4.24)–(4.26) is ue (x, y) = sin(π x)(h1 cosh(π y) + h2 y sinh(π y)), where h1 =
2 sinh(π)+2π cosh(π) 2π 3 +pi2 sinh(2π)
, h2 =
−2 sinh(pi) 2π 2 +π sinh(2π)
.
Let x → x − 1, y → y − 1, Eqs. (4.24)–(4.26) becomes
∂ 4u ∂ 4u ∂ 4u + 2 2 2 + 4 = 0, 4 ∂x ∂x ∂y ∂y
(x, y) = [0, 2] × [0, 2],
with the boundary conditions u(x, 0) = u(x, 2) =
sin(π (x − 1))
π2
,
∂ u(x, 0) ∂ u(x, 2) sin(π (x − 1))(2π 2 − π sinh(2π )) = = . ∂ y2 ∂ y2 2π 2 + π sinh(2π ) 2
2
Consider a wavelet with support over 0 ⩽ x ⩽ 2, using the method of Problem 3, we obtain an approximate solution. The result is plotted in Fig. 4. The accuracy of the results can be estimated by error function eJ =
max |u(xk , xl ) − ue (xk , xl )|.
1⩽k,l⩽2M
The eJ is shown in Table 4.
2868
Z. Shi, Y.-y. Cao / Mathematical and Computer Modelling 54 (2011) 2858–2868 Table 4 Error estimates for Problem 4. J
eJ
3 4 5 6
0.0022 0.0010 3.2241e−004 9.0296e−005
5. Conclusion In the paper, we replace radial basis functions with Haar wavelets to solve Poisson equations and biharmonic equations. The reason for using Haar wavelets are sparse representation, fast transformation, and the possibility of implementing fast algorithms. An extension of the proposed method to the case of non-rectangular domains and using other wavelets (Legendre wavelet, Shannon wavelet, Harmonic wavelet) to solve biharmonic equations will be reported in future work. Acknowledgment This work is supported by the Specialized Research Fund of Shaanxi Provincial Department of Education of China under grant No. 09JK539 and Natural Science Foundation of Shanxi Province of China under Grant No. 2009JM1002. References [1] Y kwon, J.W. Stephenson, Single cell finite difference approximation for Poisson equation in three variables, Appl. Math. Notes 2 (1982) 13–17. [2] M.M. Gupta, A fourth-order Poisson solver, J. Comput. Phys. 55 (1985) 166–172. [3] J. Wang, W. Zhong, J. Zhang, A general meshsize fourth-order compact difference discretization scheme for 3D Poisson equation, Appl. Math. Comput. 183 (2006) 804–812. [4] Xuan Wang, Zhifeng Yang, Gordon Huang, Bin Chen, A high-order compact difference scheme for 2D Laplace and Poisson equations in non-uniform grid systems, Commun. Nonlinear Sci. Numer. Simul. 14 (2) (2009) 379–398. [5] Yongbin Ge, Multigrid method and fourth-order compact difference discretization scheme with unequal mesh sizes for 3D poisson equation, J. Comput. Phys. 229 (18) (2010) 6381–6391. [6] E. Perrey-Debain, H.G. ter Morsche, B-Spline approximation and fast wavelet transform for an efficient evaluation of particular solutions for Poisson’s equation, Eng. Anal. Bound. Elem. 26 (1) (2002) 1–13. [7] K.E. Atkinson, The numerical evaluation of particular solution for Poisson’s equation, IMA J. Numer. Anal. 5 (1985) 319–338. [8] M.A. Golberg, C.S. Chen, The theory of radial basis functions applied to The BEM for inhomogeneous partial differential equations, Bound Elem. Commun. 5 (1994) 57–61. [9] Mehdi Dehghan, Akbar Mohebbi, Multigrid solution of high order discretisation for three-dimensional biharmonic equation with Dirichlet boundary conditions of second kind, Appl. Math. Comput. 180 (2006) 575–593. [10] Youngmok Jeon, New indirect scalar boundary integral equation formulas for the biharmonic equation, J. Comput. Appl. Math. 135 (2) (2001) 313–324. [11] Quang A. Dang, Iterative method for solving the Neumann boundary value problem for biharmonic type equation, J. Comput. Appl. Math. 196 (2) (2006) 634–643. [12] N. Mai-Duy, R.I. Tanner, A spectral collocation method based on integrated Chebyshev polynomials for two-dimensional biharmonic boundary-value problems, J. Comput. Appl. Math. 201 (2007) 30–47. [13] N. Mai-Duy, R.I. Tanner, Solving high-order partial differential equations with indirect radial basis function networks, Int. J. Numer. Meth. Engng. 63 (2005) 1636–1654. [14] N. Mai-Duy, T. Tran-Cong, Solving biharmonic problems with scattered-point discretization using indirect radial-basis-function networks, Eng. Anal. Bound. Elem. 30 (2) (2006) 77–87. [15] N. Mai-Duy, R.I. Tanner, An effective high order interpolation scheme in BIEM for biharmonic boundary value problems, Eng. Anal. Bound. Elem. 29 (3) (2005) 210–223. [16] Xiaolin Li, Jialin Zhu, A Galerkin boundary node method for biharmonic problems, Eng. Anal. Bound. Elem. 33 (6) (2009) 858–865. [17] N. Mai-Duy, H. See, T. Tran-Cong, A spectral collocation technique based on integrated Chebyshev polynomials for biharmonic problems in irregular domains, Appl. Math. Model. 33 (1) (2009) 284–299. [18] Bernard Bialecki, A fast solver for the orthogonal spline collocation solution of the biharmonic Dirichlet problem on rectangles, J. Comput. Phys. 191 (2) (2003) 601–621. [19] ü. Lepik, Haar wavelet method for solving higher order differential equation, Internat. J. Math. Comput. 1 (8) (2008) 84–94. [20] ü. Lepik, Solving integral and differential equations by the aid of non-uniform Haar wavelets, Appl. Math. Comput. 198 (1) (2008) 326–332. [21] ü. Lepik, Solving fractional integral equations by the Haar wavelet method, Appl. Math. Comput. 214 (2) (2009) 468–478.