Computers and Mathematics with Applications (
)
–
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Numerical solutions of elliptic partial differential equations using Chebyshev polynomials B. Khatri Ghimire, H.Y. Tian ∗ , A.R. Lamichhane Department of Mathematics, The University of Southern Mississippi, Hattiesburg, MS 39406, USA
article
abstract
info
Article history: Received 2 September 2015 Received in revised form 20 April 2016 Accepted 11 June 2016 Available online xxxx Keywords: Chebyshev polynomials Particular solutions The method of fundamental solutions Collocation Trefftz method Poisson equation Modified Helmholtz equation
We present a simple and effective Chebyshev polynomial scheme (CPS) combined with the method of fundamental solutions (MFS) and the equilibrated collocation Trefftz method for the numerical solutions of inhomogeneous elliptic partial differential equations (PDEs). In this paper, CPS is applied in a two-step approach. First, Chebyshev polynomials are used to approximate a particular solution of a PDE. Chebyshev nodes which are the roots of Chebyshev polynomials are used in the polynomial interpolation due to its spectral convergence. Then the resulting homogeneous equation is solved by boundary type methods including the MFS and the equilibrated collocation Trefftz method. Numerical results for problems on various irregular domains show that our proposed scheme is highly accurate and efficient. © 2016 Elsevier Ltd. All rights reserved.
1. Introduction Over the past couple of decades, meshless methods have been gaining popularity in engineering and scientific computing. Many researchers in mathematics and engineering have successfully applied meshless methods to solve many challenging problems in their fields [1–4]. In recent years meshless methods have become a useful alternative to traditional methods such as finite element method (FEM) and finite difference method (FDM). There are two types of meshless methods: domain-type and boundary-type. The domain type methods are evolved from the FEM whereas the boundary-type methods are from the boundary element method (BEM). One of the common goals of developing meshless methods is to solve a given set of partial differential equations (PDEs) with minimum human and computational costs. Hence, other than the accuracy, the simplicity and efficiency of a meshless algorithm is also of great importance. Let us consider the following boundary value problem:
Lu(x, y) = f (x, y),
(x, y) ∈ Ω ,
(1)
with Dirichlet boundary condition u(x, y) = g (x, y),
(x, y) ∈ ∂ Ω1 ,
(2)
and Neumann boundary condition
∂ u(x, y) = h(x, y), ∂n ∗
(x, y) ∈ ∂ Ω2 ,
Corresponding author. E-mail address:
[email protected] (H.Y. Tian).
http://dx.doi.org/10.1016/j.camwa.2016.06.012 0898-1221/© 2016 Elsevier Ltd. All rights reserved.
(3)
2
B. Khatri Ghimire et al. / Computers and Mathematics with Applications (
)
–
where Ω ⊂ R2 is a simply connected domain bounded by a simple closed curve, ∂ Ω1 and ∂ Ω2 are parts of the boundary ∂ Ω with ∂ Ω1 ∪ ∂ Ω2 = ∂ Ω and ∂ Ω1 ∩ ∂ Ω2 = ∅, L is an elliptic differential operator, f , g and h are given functions, and the function f can be smoothly extended to a rectangular domain containing Ω . In this paper, we use an approach that combines a particular solution method and a boundary method to solve the problem (1)–(3). Let u = up + uh be the solution of (1)–(3) where up is a particular solution that satisfies (1), and its associated homogeneous solution uh (x, y) satisfies,
Luh (x, y) = 0, uh (x, y) = g (x, y) − up (x, y), ∂ up (x, y) ∂ uh (x, y) = h(x, y) − , ∂n ∂n
(x, y) ∈ Ω , (x, y) ∈ ∂ Ω1 ,
(5)
(x, y) ∈ ∂ Ω2 .
(6)
(4)
The basic idea is to transform the given inhomogeneous problem into a homogeneous one via a numerically obtained particular solution to the inhomogeneous equation. A particular solution satisfies the given differential equation (1) in the infinite domain, but it does not necessarily satisfy the given boundary conditions (2)–(3). In 1967, Fox et al. [5] proposed the method of particular solutions (MPS) for solving elliptic eigenvalue problems using the combination of Bessel functions and sine functions as basis functions. In the recent decades, many different numerical techniques of the method of particular solutions have been developed for solving partial differential equations [1,2,6–12]. In these approaches, a variety of basis functions have been used. Among them, radial basis functions (RBFs) have been very popular due to their simplicity and effectiveness in implementation. Typically, there are two approaches for constructing the particular solution Φ (r ) that satisfies the equation
LΦ (r ) = φ(r ).
(7)
One approach is to utilize the RBFs as the basis functions for constructing the particular solution Φ (r ) through collocation method. The other approach is to use the RBFs as the basis functions for approximating φ(r ) and for deriving a particular solution Φ (r ) from Eq. (7) by reverse differentiation process. We remark that the first approach using RBFs does not guarantee the matrix invertibility [13] while the matrix by the second approach possesses the positive definite property [14]. Kansa’s method [2] as a popular method in the meshless literature belongs in the first approach, while the method of approximate particular solutions (MAPS) in [1,9,12] belong in the second approach. In this paper, we employ the first approach to obtain a particular solution, but we use Chebyshev polynomials as basis functions. There exist different methods of using polynomials [8,11,15–17] for the evaluation of approximate particular solutions. Chen et al. in [8] obtained particular solutions in analytical form for 2-D Poisson equation when the forcing term f is a homogeneous polynomial. Golberg et al. in [11] implemented MPS using Chebyshev interpolants for 2-D and 3-D Helmholtz type equations when the forcing term is monomial. However, book keeping of the many monomial terms of the approximate polynomial f and the particular solutions corresponding to these terms becomes very tedious and inefficient. Chen et al. [15] used a finite geometric series expansion on a differential operator to directly obtain a particular solution corresponding to each Chebyshev polynomials. However, the procedure of the actual implementation is quite cumbersome since some tedious algebraic operations require symbolic manipulations that impede the computational efficiency. Karageorghis and Kyza [17] derived particular solutions for Poisson and bi-harmonic equations directly using Chebyshev polynomials as basis functions. In [17], they made use of special properties of the second derivatives of Chebyshev polynomials without the need for the expanded form of the polynomials. They used approach two for the MPS, and in the solution process they needed to solve a block matrix system, the dimension of which is big when the degree of the Chebyshev polynomials used is high. Ding et al. in [16] proposed a recursive formulation or a matrix free method to derive a particular solution. This method requires the expansion of each of the Chebyshev polynomials which comprises the approximate forcing term. To avoid the expansion of the Chebyshev polynomials, we propose a particular solution method that combines the direct collocation and the reduction of the second order derivative of Chebyshev polynomials. The method is easy in coding and implementation. We use the Gauss–Lobatto nodes to ensure the spectral convergence for the Chebyshev interpolation. Also there is no problem of ill-conditioning which occurs when RBF approximations are used. Once a particular solution is found, we continue to solve the associated homogeneous problem (4)–(6) by boundary methods. We use one of the two options: the method of fundamental solution and the equilibrated approach in collocation Trefftz method (CTM). The MFS was first proposed in 1964 by Kupradze and Aleksidze [18] and it has been used for solving various problems since then. The high accuracy of both MFS and the MPS of Chebyshev polynomials helps produce a highly accurate solution of (1)–(3). Due to the fact that MFS requires a fictitious boundary in the solution process, we use the CTM [19] as an alternative. The CTM is based on T-complete functions and does not require a fictitious boundary for solving the homogeneous problem. To treat the highly ill-conditioned linear system resulting from the direct implementation of CTM, we have adopted the equilibrated matrix concept [3,20] to determine the scales and to construct an equivalent linear algebraic problem with a leading matrix less ill-conditioned. The organization of this paper is as follows: In Section 2, we briefly discuss the properties of the Chebyshev polynomials. Then we use the Chebyshev collocation method to find an approximate particular solution with the change of variables and re-scaling of PDEs. In Section 3, we describe the standard boundary-type meshless methods MFS and CTM, and the incorporation of equilibrated matrix approach in CTM. The numerical results for verification of the accuracy of our proposed method is presented in Section 4. Conclusions drawn from the work and the numerical results are given in Section 5.
B. Khatri Ghimire et al. / Computers and Mathematics with Applications (
)
–
3
2. Chebyshev collocation method Chebyshev polynomials are a sequence of orthogonal polynomials which can be defined recursively. The roots of Chebyshev polynomials, also called Gauss–Lobatto nodes, are used in the polynomial interpolation for high order of accuracy. The resulting polynomial minimizes the problem of Runge’s phenomena. It is called mini–max polynomial of a function since the polynomial minimizes the maximum error. For the numerical work, we have made use of the following recursive relations [21–23]: For the Chebyshev polynomials of the first kind, T0 (x) = 1,
T1 (x) = x,
Tn+1 (x) = 2xTn (x) − Tn−1 (x),
n = 1, 2, . . . , N .
For the Chebyshev polynomials of the second kind, U0 (x) = 1,
U1 (x) = 2x,
Un+1 (x) = 2xUn (x) − Un−1 (x),
n = 1, 2, . . . , N .
The relations between Tn and Un are as follows: Tn (x) = Un (x) − xUn−1 (x), Un (x) = xUn−1 (x) + Tn (x). The first order derivatives of Tn are dTn
= nUn−1 , for n = 1, 2, . . . , N , dx and the second order derivatives of Tn are d2 Tn dx2
=
n x2 − 1
(n + 1)(Tn − Un ),
n4 − n2 = , dx2 (x=1) 3 n4 − n2 d2 Tn = (−1)n , 2 dx 3 d2 Tn
(x=−1)
for n = 1, 2, . . . , N. The important part of the method involves the approximation of a particular solution up as a linear combination of Chebyshev polynomials, i.e., up (x, y) ≈ uˆ p (x, y) =
m n
αij Ti (x)Tj (y),
(x, y) ∈ Ω ,
(8)
i=0 j=0
where {αij } are the coefficients to be determined. By collocating on the Gauss–Lobatto nodes {(xk , yk )}Nk=1 , αij can be obtained by solving the linear system, m n
αij LTi (xk )Tj (yk ) = f (xk , yk ),
1 ≤ k ≤ N.
(9)
i =0 j =0
Since particular solutions do not necessarily satisfy the boundary conditions, we have the flexibility in extending the domain Ω to a rectangle [a, b] × [c , d] that contains Ω . It is important to re-scale the domain Ω and the rectangle ˆ [a, b] × [c , d] by the change of variables ξ = (2x − a − b) / (b − a) , η = (2y − c − d) / (d − c ) so that the new domain Ω is embedded in a square [−1, 1] × [−1, 1]. The PDE will change accordingly in terms of the new variables ξ and η. For illustration, we consider the following Poisson problem
1u(x, y) = f (x, y), (x, y) ∈ Ω , u(x, y) = g (x, y), (x, y) ∈ ∂ Ω . By the change of variables, x = α1 ξ + β1 ,
y = α2 η + β2 ,
or equivalently
ξ=
x − β1
α1
,
η=
y − β2
α2
,
(10) (11)
4
B. Khatri Ghimire et al. / Computers and Mathematics with Applications (
)
–
with
α1 =
b−a
α2 =
d−c
2 2
,
β1 =
a+b
,
β2 =
c+d
2 2
,
(12)
,
(13)
the problem (10)–(11) is transformed into a problem in the (ξ , η) plane. Letting u(x, y) = U (ξ , η), we have 1 ∂ 2U ∂ 2u = , ∂ x2 α12 ∂ξ 2
∂ 2u 1 ∂ 2U = . ∂ y2 α22 ∂η2
Thus the original problem (10)–(11) is transformed into the problem, 1 ∂ 2U
α ∂ξ 2 1
2
+
1 ∂ 2U
α22 ∂η2
= f (α1 ξ + β1 , α2 η + β2 ),
U (ξ , η) = g (α1 ξ + β1 , α2 η + β2 ),
(ξ , η) ∈ Ω ,
(14)
(ξ , η) ∈ ∂ Ω ,
where Ω ∪ ∂ Ω = [−1, 1] × [−1, 1]. For Eq. (14) to be a Poisson equation, we would like to have α1 = α2 , which is true when b − a = d − c by (12)–(13). By embedding the original domain [a, b] × [c , d] into a larger domain [a1 , b1 ] × [c1 , d1 ] such that b1 − a1 = d1 − c1 , we would be able to achieve this purpose. Consequently the MFS or the equilibrated CTM are used for the differential operators on the new domain. 3. Boundary collocation methods 3.1. The method of fundamental solutions The method of fundamental solutions is a popular classical meshless boundary-type method. It is free of integration and can be extended to higher dimensions. The MFS requires only boundary discretization and it shares the same advantages as the BEM over domain discretization methods [24]. The singularities of MFS are placed outside the domain. Hence a fictitious boundary is necessary to show where the singularities are located. For its easy implementation and high accuracy, the MFS is applied here for the homogeneous problem. With also an accurate approximate particular solution up from CPS, we expect an accurate solution u for the boundary value problem (1)–(3). The MFS can be applied to differential operators whose fundamental solutions are known. Readers can refer to [24] and references cited therein for more details. The idea of MFS is to approximate the solution uh as a linear combination of fundamental solutions, uh (x) ≈
M
αj G(x, y; ξj , ηj ),
(15)
j =0
where G(x, y; ξj , ηj ) are the fundamental solutions of the differential operator L with singularity at ξj , ηj . To determine the optimal location of the source points, the leave-one-out cross validation (LOOCV) method [25,26] can be used for the MFS. Once the source points have been chosen, the coefficients {αj } in Eq. (15) can be obtained by collocation on the boundary ∂ Ω . The MFS approach has been adopted by many numerical methods. Different formulations of modified MFS have also been developed so that the fictitious boundary can coincide with the domain boundary. These include regularized meshless method [4], the modified method of fundamental solutions [27], and the boundary distributed source method [28]. In these methods, source points can be directly located on the domain boundary.
3.2. Collocation Trefftz method The Trefftz method was first proposed by Trefftz in 1926 [29]. It is a boundary-type method. In [19], Kita and Kamiya classified the Trefftz formulation into the indirect and the direct ones. In this paper we use the indirect formulation, which is the original one presented by Trefftz. The solution of the problem is approximated by the superposition of the functions satisfying the governing homogeneous equation. These functions are called T-complete functions. Then the unknown parameters are determined by the collocation on the boundary. Li et al. named this as collocation Trefftz method in their book [30]. The algorithm for the CTM requires the explicit form of the T-complete functions. Considering a disk S = {(r , θ )|0 ≤ r ≤ R, 0 ≤ θ ≤ 2π }, we list below the T-complete functions [30] of a few typical equations in 2D, which are often used in practice.
B. Khatri Ghimire et al. / Computers and Mathematics with Applications (
)
–
5
Table 1 List of T-complete functions. PDEs
T-complete function
1u = 0 1u − λ2 u = 0 1u + λ2 u = 0 ∆2 u = 0
1, r n einθ I0 (λr ), In (λr )einθ J0 (λr ), Jn (λr )einθ 1, r n einθ , r 2 , r n+2 einθ
Note that In (r ) and Jn (r ) are respectively modified Bessel functions and Bessel functions of the first kind, i is the imaginary unit, and r and θ are the polar coordinates. As an example, let us consider the modified Helmholtz boundary value problem,
1u(x, y) − λ2 u(x, y) = 0, u(x, y) = g (x, y), ∂ u(x, y) = h(x, y), ∂n
(x, y) ∈ Ω , (x, y) ∈ ∂ Ω1 ,
(17)
(x, y) ∈ ∂ Ω2 ,
(18)
(16)
where n denotes the unit outward normal vector on the boundary, and g and h are given functions. By CTM, the solution of the problem (16)–(18) is represented as uh ≈
n
ai uˆ i ,
i =0
where uˆ i are the T-complete functions given in Table 1 satisfying Eq. (16). For an interior modified Helmholtz problem, the solution can be expressed as u(r , θ ) = a0 I0 (λr ) +
m
ai Ii (λr ) cos(iθ ) + bi Ii (λr ) sin(iθ ),
(19)
i =1
where m is the order of the series solution, and a0 , a1 , b1 , . . . , am , bm are the unknown coefficients to be determined. Solving the problem using boundary collocation, we obtain the following linear system, Be = d
(20)
where B = Bij is the resulting collocation matrix, e is the vector of unknown coefficients and d is the right hand side. The dimension of B is M × (2m + 1) if the number of boundary points is M. 3.3. Multiple scale CTM and equilibrated matrix We would like to point out that the matrix in (20) from CTM is highly ill-conditioned. To alleviate the ill-conditioning of this linear system, we incorporate the multiple scale Trefftz method and the equilibrated matrix concept [3,20]. An equilibrated matrix is a resulting matrix of which all the column norms or row norms are the same, and under this treatment the matrix has a significantly reduced condition number. By this method, the series solution (19) is replaced with the following form: u(r , θ ) = a0
I0 (λr ) R0
+
m Ii (λr ) Ii (λr ) ai cos(iθ ) + bi sin(iθ ), i =1
R2i−1
R2i
where
a0 = a0 R0 ,
ai = ai R2i−1 ,
bi = bi R2i ,
(21)
and R0 , R1 , . . . , R2m are the multiple-scale characteristic lengths. The concept of characteristic length was first proposed by Liu [20,31] to reduce the ill-posedness of the Trefftz method and MFS for solving Laplace equations. The characteristic lengths are computed by
R j −1
M = α B2 , ij
j = 1, 2, . . . , 2m + 1.
i=1
Then the linear system (20) is reduced to the system Ac = d,
(22)
with A = Aij and Aij = Bij /Rj−1 . That is, all the column norms of matrix A are 1/α . So A is an equilibrated matrix. Since A is well conditioned, we proceed to compute the linear system (22).
6
B. Khatri Ghimire et al. / Computers and Mathematics with Applications (
)
–
Fig. 1. Original domain.
Fig. 2. Computational domain.
4. Numerical results The method of Chebyshev particular solutions is coupled with one of the two boundary methods MFS and the equilibrated CTM for solving the elliptic PDEs. To verify the effectiveness of the method, we carry out numerical experiments on five examples from second to fourth order PDEs. When a fourth order PDE can be reduced to second order PDEs, the Chebyshev polynomial technique for particular solutions is applicable. While implementing our numerical schemes, we assume that the particular solution to be approximated can be extended smoothly outside its domain Ω . To validate the numerical accuracy, we calculate the maximum absolute error (MAE) which is defined as follows: MAE = max (ˆuj − uj ) ,
j
j = 1 . . . q,
where q is the number of testing nodes chosen randomly in the domain, and uˆ j and uj denote respectively the approximate solution and the exact solution at the jth node. In the numerical results, Nx and Ny denote the number of Chebyshev nodes in x and y directions respectively. Each original domain is transformed into the standard square domain by the change of variables as discussed in Chapter 2, which is illustrated by Figs. 1 and 2. Example 1. We first consider a Poisson equation with a Dirichlet boundary condition:
1u(x, y) = f (x, y), (x, y) ∈ Ω , u(x, y) = g (x, y), (x, y) ∈ ∂ Ω . The forcing term f (x, y) and the boundary data g (x, y) are given such that the exact solution is: u(x, y) = cos(3x) − cos(3y).
B. Khatri Ghimire et al. / Computers and Mathematics with Applications (
)
–
7
Fig. 3. The profile of amoeba domain.
Fig. 4. The profile of Cassini domain.
The Amoeba-shaped domain Ω as shown in Fig. 3 is bounded by the curve defined by the following parametric equation:
∂ Ω = {(x, y)|x = ρ cos θ , y = ρ sin θ , 0 ≤ θ ≤ 2π }, where
ρ = esin(θ) sin2 (2θ ) + ecos(θ) cos2 (2θ ). The Cassini-shaped curve as shown in Fig. 4 is defined by the following parametric equation:
∂ Ω = {(x, y)|x = ρ cos θ , y = ρ sin θ , 0 ≤ θ ≤ 2π }, where
ρ = cos(3θ ) +
13 2 − sin (3θ ) .
2
We employ both MFS and equilibrated CTM for the homogeneous solution. When using MFS, we let the source points distribute on a circle of radius r. By LOOCV, the sub-optimal radius r for the Cassini and Amoeba are 1.398 and 1.529 respectively. For CTM, the number of basis functions is 62. For each domain, we use 250 boundary points and 120 test points. We list in Table 2 the numerical results for this problem with Ω being an Amoeba or a Cassini domain. Since Chebyshev polynomials provide spectral convergence on the domain and MFS and CTM are highly accurate boundary methods, the overall solution is highly accurate. From Table 2, we observe that the equilibrated CTM performs slightly better than the MFS on both irregular domains (see Fig. 5).
8
B. Khatri Ghimire et al. / Computers and Mathematics with Applications (
)
–
Fig. 5. The profile of the error plot in Example 1.
Table 2 Max Abs error in Example 1. Nx = Ny
10 12 14 16 18 20 22 24 26 28 30
CPS–MFS
CPS–CTM
Cassini
Amoeba
Cassini
Amoeba
1.640e−03 4.034e−05 9.067e−07 1.199e−08 1.601e−10 1.210e−12 8.099e−14 1.922e−13 1.842e−13 6.796e−13 1.523e−13
4.163e−02 2.832e−03 1.396e−04 6.705e−06 2.369e−07 3.761e−09 7.877e−11 1.213e−12 4.833e−13 9.527e−13 1.169e−12
5.671e−05 1.066e−06 1.639e−08 1.945e−10 2.399e−12 1.831e−14 1.110e−15 1.443e−15 1.720e−15 1.276e−15 8.579e−11
4.604e−02 3.227e−03 1.524e−04 4.482e−06 1.749e−07 3.916e−09 8.989e−11 9.743e−13 1.887e−14 6.022e−15 5.020e−13
Example 2. We consider a modified Helmholtz problem,
(∆ − λ2 )u(x, y) = f (x, y), (x, y) ∈ Ω , u(x, y) = g (x, y), (x, y) ∈ ∂ Ω , where f (x, y) and g (x, y) are given such that the analytical solution is u(x, y) = cos(3x) − cos(3y). We use Amoeba-shaped domain Ω as shown in Fig. 3 and the Star-shaped domain Ω as shown in Fig. 6 which is bounded by the curve defined by the parametric equation:
∂ Ω = {(x, y)| x = ρ cos θ , y = ρ sin θ , 0 ≤ θ ≤ 2π }, where
ρ = 1 + cos2 (4θ ). For MFS the fictitious boundary (xs , ys ) where the source points are located is chosen to have the same shape as that of the boundary. By LOOCV, the sub-optimal radius r for the Amoeba and the Star domains are 1.030 and 1.023 respectively. For CTM, the number of basis functions is 62. For each domain, we use 250 boundary points and 120 test points and wave number λ = 35. Table 3 shows the numerical results by the CPS with the MFS and the equilibrated CTM on both domains. These results show the fast convergence and accuracy of the proposed method. We notice that the CPS–CTM performs slightly better than the CPS–MFS. In addition, the CTM avoids the determination of a suitable location for source points, which is a great advantage of the CTM over the MFS (see Fig. 7).
B. Khatri Ghimire et al. / Computers and Mathematics with Applications (
)
Fig. 6. Star shape domain.
Fig. 7. Error plot.
Table 3 Max Abs error in Example 2. Nx = Ny
10 12 14 16 18 20 22 24 26 28 30
CPS–MFS
CPS–CTM
Amoeba
Star
Amoeba
Star
8.777e−02 1.044e−02 3.355e−04 2.189e−05 1.158e−06 1.739e−08 3.661e−10 1.089e−11 1.887e−12 2.032e−12 8.590e−12
3.559e−02 3.105e−03 1.788e−04 7.609e−06 2.553e−07 6.729e−09 1.478e−10 2.651e−12 4.041e−14 1.377e−14 1.659e−14
1.765e−02 1.396e−03 6.037e−05 2.286e−06 8.017e−08 1.888e−09 5.944e−11 2.201e−12 6.974e−14 7.327e−15 1.690e−14
3.338e−03 1.930e−04 7.078e−06 2.205e−07 4.462e−09 1.039e−10 2.002e−12 5.018e−14 1.998e−14 2.048e−14 1.328e−14
Example 3. We consider the following Poisson problem with a mixed boundary condition,
1u(x, y) = f (x, y), (x, y) ∈ Ω , u(x, y) = g1 (x, y), (x, y) ∈ ∂ Ω1 , ∂ u(x, y) = g2 (x, y), (x, y) ∈ ∂ Ω2 , ∂n
–
9
10
B. Khatri Ghimire et al. / Computers and Mathematics with Applications (
)
–
Table 4 MAE by CPS–MFS in Example 3. Nx = Ny
Cassini
Unit square
Amoeba
10 12 14 16 18 20 22 24 26 28 30
1.266e−02 1.958e−04 3.565e−06 9.362e−08 1.832e−10 3.204e−12 6.902e−12 9.390e−12 2.390e−11 3.613e−11 8.021e−10
7.975e−05 1.920e−06 1.946e−08 2.457e−10 3.198e−12 4.559e−12 6.130e−12 1.154e−11 5.784e−11 1.606e−10 1.922e−10
5.430e−01 7.188e−02 2.700e−03 2.683e−04 1.375e−05 8.764e−07 9.825e−09 2.340e−11 6.443e−12 2.490e−13 1.469e−13
Fig. 8. The profile of the error plot in Example 3.
πy
the forcing term and boundary data are given such that the analytical solution u(x, y) = sin(π x) cos 2 . This problem is solved on three different domains: the Cassini domain as shown in Fig. 4, unit square and the Amoeba-shaped curve as given by Fig. 3. For each domain, ∂ Ω1 and ∂ Ω2 each represents half of the boundary ∂ Ω with ∂ Ω1 ∪ ∂ Ω2 = ∂ Ω and ∂ Ω1 ∩ ∂ Ω2 = ∅. For the MFS we use a circle of radius r for the source location and the LOOCV method for a suitable value of the radius which is 1.48 for Amoeba, 1.96 for Cassini and unit square. The numbers of boundary and test points are 250 and 120 for Amoeba and Cassini whereas they are 276 and 64 for unit square. The computational domain and errors are plotted in Figs. 4 and 8 respectively. The results shown in Table 4 indicate the high accuracy of the proposed method. Example 4. We consider the cases when three Franke’s benchmark test functions [32] are the exact solutions for Poisson equation in Example 1 and the modified Helmholtz equation in Example 2 on the unit square. The three Franke’s test functions are given below:
1 3 − 1 (9x+1)2 − 10 1 − 1 (9x−7)2 − 14 (9y−3)2 1 3 − 1 [(9x−2)2 +(9y−2)2 ] 2 2 (9y+1)2 e 4 + e 49 + e 4 − e−[(9x−4) −(9y−7) ] 4 4 2 5 1 F 2(x, y) = (tanh(9y − 9x) + 1) 9 1.25 + cos(5.4y) F 3(x, y) = . 6[1 + (3x − 1)2 ]
F 1(x, y) =
The plots of the three Franke’s test functions are given in Figs. 9–11. The surfaces of these functions are more complex than those in previous examples. So more Chebyshev nodes are needed to achieve high accuracy. While using the MFS we choose a fictitious boundary which has the same shape as that of the boundary. Using LOOCV for Helmholtz equation, the sub-optimal radius r for F 1, F 2 and F 3 are 1.025, 1.024 and 1.026 respectively; and for Poisson equation, the sub-optimal radius r for F 1, F 2 and F 3 are 1.449, 1.326 and 1.382 respectively. We let λ = 150 in the modified Helmholtz equation. The numbers of boundary points and test points are 396 and 64 respectively. The numerical results for the Poisson and the modified Helmholtz problems are listed in Tables 5 and 6 respectively. The error plots of the three Franke’s test functions for the modified Helmholtz problem are given by Fig. 12.
B. Khatri Ghimire et al. / Computers and Mathematics with Applications (
Table 5 MAE by CPS–MFS for Poisson in Example 4. Nx = Ny
MAE(F1)
MAE(F2)
MAE(F3)
10 15 20 25 30 35 40 45 50
9.821e−02 4.641e−03 1.147e−03 4.095e−05 2.289e−06 2.346e−07 3.115e−08 1.679e−08 1.791e−09
1.695e−02 1.196e−03 3.131e−04 4.516e−05 1.198e−06 6.167e−07 2.328e−07 1.931e−08 3.318e−08
3.631e−03 1.832e−04 6.534e−05 1.152e−07 1.711e−08 4.611e−09 1.724e−09 4.693e−10 1.688e−11
Fig. 9. Surface plot F1.
Fig. 10. Surface plot F2.
Fig. 11. Surface plot F3.
)
–
11
12
B. Khatri Ghimire et al. / Computers and Mathematics with Applications (
)
–
Fig. 12. Error plot of Franke’s TFs. Table 6 MAE by CPS–MFS for Helmholtz in Example 4. Nx = Ny
MAE(F1)
MAE(F2)
MAE(F3)
10 15 20 25 30 35 40 45 50
2.7217e−02 1.006e−03 1.986e−04 4.210e−05 1.769e−06 1.581e−08 2.611e−10 2.986e−11 5.327e−13
5.343e−03 4.210e−04 3.615e−04 2.679e−05 2.005e−06 3.395e−06 7.221e−07 1.218e−08 3.710e−09
1.832e−03 3.766e−05 8.479e−06 1.807e−07 2.357e−09 3.433e−11 4.208e−12 3.302e−14 2.126e−15
Example 5. We now consider the following fourth order PDE with domain Ω being a unit square.
(∆2 − λ4 )u = f (x, y), u = g (x, y), 1u = h(x, y),
(x, y) ∈ Ω , ( x, y ) ∈ ∂ Ω , (x, y) ∈ ∂ Ω
(23)
where λ = 50, f , g and h are given functions such that the exact solution is u(x, y) = sin(π x) cos
πy 2
,
( x, y ) ∈ Ω .
Since ∆2 − λ4 = (∆ − λ2 )(∆ + λ2 ). we calculate the particular solution up of Eq. (23) in two steps,
(∆ + λ2 )up = vp
(24)
(∆ − λ )vp = f .
(25)
2
By CPS, we first find vp for Eq. (25) and then find up for Eq. (24). If we solve the fourth order partial differential equation without converting it into two of the second order equations, we would have to evaluate the fourth order derivatives of the Chebyshev polynomials. However with the conversion into Eqs. (24) and (25), we can directly apply the scheme of Section 2. We use 276 boundary points and 64 test points. The LOOCV is used to decide the sub-optimal source location for a fictitious circle of radius r. In this case a circle of radius 2.4 is chosen. From Table 7, the high accuracy of our proposed method is again observed. This shows that our method is equally applicable for solving fourth order PDEs. 5. Conclusion The main purpose of this paper is to propose a simple, direct, and efficient two-step method for solving elliptic PDEs. In step one we use a Chebyshev polynomial scheme, through which we approximate the particular solution of a PDE. The idea is to use Chebyshev polynomials as basis functions to directly represent the particular solution and its derivatives. The expansion of Chebyshev polynomials is not necessary. It is simple and yet highly accurate. In addition, this approach avoids the problems of ill-conditioning and the uncertainty of the shape parameter when using RBFs. In step two we use the MFS or the collocation Trefftz method for the corresponding homogeneous problem. For the MFS, the LOOCV technique can be
B. Khatri Ghimire et al. / Computers and Mathematics with Applications (
)
–
13
Table 7 MAE by CPS–MFS in Example 5. Nx = Ny
Maximum Abs error
10 12 14 16 18 20 22 25
1.249e−05 3.968e−07 5.573e−09 7.300e−11 1.308e−12 1.554e−12 1.537e−12 5.730e−11
used for identifying a suitable location for source points. For the CTM, an equilibrated matrix is formulated to reduce the condition number of the discretization matrix from the CTM. Numerical examples show that the results by CPS–MFS and CPS–CTM are comparable, they are extremely satisfactory in solving elliptic boundary value problems. Despite the simplicity and accuracy of the proposed method, it has its limitations like many other numerical techniques. To apply this method, we assume that the forcing term f (x, y) can be smoothly extended to a rectangular domain which contains the given domain. Otherwise, the basis functions that work for scattered data is more appropriate for finding particular solutions. In this paper, we have successfully employed our approach to solve second and fourth order elliptic PDEs on irregular domains. Besides its accuracy, another advantage of the method is that it can be easily extended for solving higher dimensional problems. For our future research we also plan to implement the scheme for axisymmetric problems and time-dependent problems. References [1] C.S. Chen, C.M. Fan, P.H. Wen, The method of approximate particular solutions for solving elliptic Problems with variable coefficients, Int. J. Comput. Methods 8 (3) (2011) 545–559. [2] E.J. Kansa, Multiquadrics—a scattered data approximation scheme with applications to computational fluid dynamics-I, Comput. Math. Appl. 19 (8–9) (1990) 127–145. [3] C.S. Liu, A multiple scale Trefftz method for an incomplete Cauchy problem of biharmonic equation, Eng. Anal. Bound. Elem. 37 (11) (2013) 17–42. [4] D.L. Young, K.H. Chen, C.W. Lee, Novel meshless method for solving the potential problems with arbitrary domain, J. Comput. Phys. 209 (2005) 290–321. [5] L. Fox, P. Henrici, C.B. Moler, Approximations and bounds for eigenvalues of elliptic operators, SIAM J. Numer. Anal. 4 (1967) 89–102. [6] K.E. Atkinson, The numerical evaluation of particular solutions for Poisson’s equation, IMA J. Numer. Anal. 5 (1985) 319–338. [7] C.S. Chen, Y.F. Rashed, Evaluation of thin plate spline based particular solutions of Helmholtz type operators for the dual reciprocity method, Mech. Res. Commun. 25 (1998) 195–201. [8] C.S. Chen, M.A. Golberg, A.S. Muleshkov, The numerical evaluation of particular solutions for Poisson’s equation-a revisit, in: C.A. Brebbia, H. Power (Eds.), BEM, Vol. XXI, WIT press, Southampton, 1999, pp. 313–322. [9] C.S. Chen, C.M. Fan, P.H. Wen, The method of approximated particular solutions for solving certain partial differential equations, Numer. Methods Partial Differential Equations 28 (2) (2012) 506–522. [10] A.H.D. Cheng, Particular solutions of Laplacian, Helmholtz-type, and polyharmonic operators involving higher order radial basis functions, Eng. Anal. Bound. Elem. 24 (2000) 531–538. [11] M.A. Golberg, A.S. Muleshkov, C.S. Chen, A.H.D. Cheng, Polynomial particular solutions for certain partial differential operators, Numer. Methods PDEs 19 (2003) 112–133. [12] G.M. Yao, J. Kolibal, C.S. Chen, A localized approach for the method of approximate particular solutions, Comput. Math. Appl. 61 (9) (2011) 2376–2387. [13] G.E. Fasshauer, Solving partial differential equations by collocation with radial basis functions, Proc. Chamonix (1996) 1–6. [14] W. Chen, Z.J. Fu, C.S. Chen, Recent Advances in Radial Basis Function Collocation Methods, Springer, 2013. [15] C.S. Chen, S.W. Lee, C.S. Huang, Derivation of particular solutions using Chebyshev polynomial based functions, Int. J. Comput. Methods (2007) 04–15. [16] J. Ding, H.Y. Tian, C.S. Chen, The recursive formulation of particular solutions for inhomogeneous elliptic PDEs with Chebyshev basis functions, Commun. Comput. Phys. 5 (5) (2009) 942–958. [17] A. Karageorghis, I. Kyza, Efficient algorithms for approximating particular solutions of elliptic equations using Chebyshev polynomials, Commun. Comput. Phys. 2 (2007) 501–521. [18] V.D. Kupradze, M.A. Aleksidze, The method of functional equations for the approximate solution of certain boundary value problems, USSR Comput. Math. Math. Phys. 4 (1964) 82–126. [19] E. Kita, N. Kamiya, Trefftz method: an overview, Adv. Eng. Softw. 24 (1995) 3–12. [20] C.S. Liu, An equilibrated method of fundamental solutions to choose the best source points for the Laplace equation, Eng. Anal. Bound. Elem. 36 (8) (2012) 1235–1245. [21] J.P. Boyd, Chebyshev and Fourier Spectral Methods, second ed., Dover publications, 2000. [22] L. Fox, I.B. Parker, Chebyshev Polynomials in Numerical Analysis, Oxford University press, London, 1968. [23] J.C. Mason, D.C. Handscomb, Chebyshev Polynomials, Chapman and Hall/CRC, Boca Raton, 2003. [24] G. Fairweather, A. Karageorghis, The method of fundamental solution for elliptic boundary value problems, Adv. Comput. Math. 9 (1998) 69–95. [25] C.S. Chen, A. Karageorghis, Y. Li, On choosing the location of the sources in the MFS, Numer. Algorithms 72 (2016) 107–130. [26] S. Rippa, An algorithm for selecting a good value for the parameter c in radial basis function interpolation, Adv. Comput. Math. Math. 11 (1999) 193–210. [27] B. Sarler, Solution of potential flow problems by the modified method of fundamental solutions: formulations with the single layer and the double layer fundamental solutions, Eng. Anal. Bound. Elem. 33 (2009) 1374–1382. [28] Y.J. Liu, A new boundary meshfree method with distributed sources, Eng. Anal. Bound. Elem. 34 (2010) 914–919. [29] E. Trefftz, Ein Gegenstuck zum Ritz’schen Verfahren, Proc. 2nd Int Cong Appl Mech, Zurch 1926, pp. 131–137. [30] Z.C. Li, T.T. Lu, H.Y. Hu, A.H.D. Cheng, Trefftz and Collocation Methods, WIT press, 2008. [31] C.S. Liu, A modified Trefftz method for two-dimensional Laplace equation considering the domain’s characteristic length, CMES-Comput. Model. Eng. Sci. 21 (1) (2007) 53–65. [32] R. Franke, Scattered data interpolation: tests of some methods, Math. Comp. 176 (1971) 1905–1915.