# Analytical Particular Solutions of Multiquadrics Associated with Polyharmonic Operators

Authors: Chia-Cheng Tsai
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

We derive two- and three-dimensional analytical particular solutions of multiquadrics (MQ) associated with the polyharmonic operators, named as the polyharmonic multiquadrics (PMQs). The methods of undetermined coefficients are constructed by observing the first few orders of the PMQs which are obtained by the symbolic software,Mathematica. By expanding the PMQs into the Laurent series, the unknown coefficients of the PMQs can be determined. The homogeneous parts of the PMQs are suitably arranged so that the PMQs are hierarchically unique and infinitely differentiable. Mathematica codes are provided for obtaining the PMQs of arbitrary orders. The derived PMQs are validated by numerical solutions for  Poisson’s equation. Numerical results indicate that the solutions obtained by the PMQs are more accurate than those by the MQ.


---

## Body

## 1. Introduction

A radial basis function (RBF) is a real-valued function whose value depends only on the distance from a certain prescribed center. Traditionally, sums of RBFs are used for approximating given functions or reconstructing scattered data. Among the various RBFs, the Gaussian function, Hardy’s multiquadrics (MQ) [1], Duchon’s augmented polyharmonic spline (APS) [2], and Wendland’s compactly supported RBF (CSRBF) [3] are the most popular ones. In the early 1990s, Kansa [4, 5] made the first attempt to use the MQ for approximating the solutions of partial differential equations (PDEs). The Kansa method is meshless, simple and has been used for a wide range of PDEs.Alternatively, Nardini and Brebbia developed the dual reciprocity method (DRM) for approximating the particular of the considered PDE by the RBFs, and then solved the complementary problem by the boundary element method (BEM) [6]. In addition to the BEM, there are also other boundary-type numerical methods such as the method of fundamental solutions (MFS) [7–11] and the Trefftz methods [12, 13].In an early development of DRM, the ad hoc function,1+r, was exclusively used. In order to improve the accuracy of the computation, researchers applied the RBF theory to the DRM [14, 15]. For example, Golberg [9, 16], Chen [17], and Karur and Ramachandran [15] demonstrated the superiority of the APS over the ad hoc trial function. Then, Golberg et al. [18] further improved the accuracy of the approximated particular solution by utilizing the exponential convergence rate of the MQ. A comprehensive review on the DRM development in this period was given by Golberg and Chen [19].When the MQ is adopted in a DRM procedure, the applicability depends on the availability of the analytical particular solution associated with the partial differential operator of a given problem. Golberg et al. [18] derived the analytical particular solutions for the Laplace equation. Samaan and Rashed [20, 21] found the analytical particular solutions for the two-dimensional problems of elasticity and elastodynamics, which could be converted to a biharmonic equation using the Galerkin-Papkovich vector [22]. Basically, the applicability of MQ is limited to the harmonic and biharmonic operators as shown in Yao’s thesis [23]. However, Yao’s method cannot be easily generalized since the method she applied to cancel the singularity is not straightforward for polyharmonic operators.In the appendix of a recent study [24], the present author has found analytical particular solutions of the MQ associated with the fourth- and third-ordered polyharmonic operators of two and three dimensions, respectively. In this study, the theory for deriving these polyharmonic multiquadrics (PMQs) is presented. In other words, we derive the analytical particular solutions of multiquadrics (MQ) associated with the two- and three-dimensional polyharmonic operators. The methods of undetermined coefficients are constructed by observing the first few orders of the PMQs which are obtained by straight integrations via the symbolic software, Mathematica. Furthermore, the Laurent expansions are used to find the unknown coefficients and to make the PMQs infinitely differentiable. Numerical experiments are carried out to validate the derived PMQs and to demonstrate the superior accuracy of the PMQs over the MQ for solving Poisson’s equation.This paper is organized as follows: the problem is mathematically modeled in Section2. Then, the two- and three-dimensional PMQs are derived in Sections 3 and 4, respectively. Some numerical experiments are carried out to validate the PMQ in Section 5 and the conclusions are drawn in Section 6.


## 2. Definition of the Problem

In this study, we derive the two- and three-dimensional polyharmonic multiquadrics (PMQs)ψL(2D) and ψL(3D) of order L, which are governed by
(1)Δr(*)ψL(*)=ψL-1(*)
with
(2)ψ0(*)=c2+r2,(3)Δr(2D)=1rddr(rddr),Δr(3D)=1r2ddr(r2ddr).
In (1), (2), and the following, (*) stands for (2D) or (3D). If we make the following change of variables
(4)r=cR,(5)ψL(*)=c(2L+1)ΨL(*),
equations (1) and (2), respectively, become
(6)ΔR(*)ΨL(*)=ΨL-1(*),(7)Ψ0(*)=1+R2.
Here, ΔR(*) is defined as (3) while replacing r by R. Equations (7) and (6) can be deduced to
(8)(ΔR(*))LΨL(*)=1+R2.It is clear that the solutions of (6)-(7) are also the solutions of (8). However, the converse is not always true since the homogeneous solutions of (8) can be arbitrarily chosen. In the next two sections, a specific choice will be made to ensure the solutions of (8) are also the solutions of (6)-(7). Then, the PMQs governed by (1)-(2) can be obtained by using (5).


## 3. Two-Dimensional PMQ

It is clear that the solution of (8) can be found by straight integrations. For the two-dimensional case, it can be expressed as follows:
(9)ΨL(2D)=(IR-1IR)L1+R2,
where the integral operator I is defined as ∫dR. When observing the first few orders of the PMQs obtained by the symbolic software, Mathematica[25], we can find that
(10)ΨL(2D)=∑j=0LAL,j(2D)R2j1+R2+∑j=0L-1BL,j(2D)R2jln(1+1+R22)+∑j=0L-1CL,j(2D)R2j+∑j=0L-1DL,j(2D)R2jlnR,
where the first two series form a particular solution with 2L+1 undetermined coefficients AL,j(2D) and BL,j(2D) and the last two series give the homogeneous solutions with 2L integration constants CL,j(2D) and DL,j(2D). To make ΨL(2D) infinitely differentiable, we obviously need
(11)DL,j(2D)=0,
and then (10) becomes
(12)ΨL(2D)=∑j=0LAL,j(2D)R2j1+R2+∑j=0L-1BL,j(2D)R2jln(1+1+R22)+∑j=0L-1CL,j(2D)R2j.In (12), the coefficients AL,j(2D) and BL,j(2D) should be determined such that the PMQ ΨL(2D) in (12) satisfies (8) and the arbitrary coefficients CL,j(2D) should be chosen to ensure the solution ΨL(2D) satisfies the hierarchical relation (6) and (7). Therefore, we expand (12) in the Maclaurin series as
(13)ΨL(2D)=∑i=0∞∑j=0min{L,i}(-1)i-j+1(2i-2j-3)!!(i-j)!2i-jAL,j(2D)R2i+∑i=1∞∑j=0min{L-1,i-1}(-1)i-j+1(2i-2j-1)!!(i-j)!2i-j+1(i-j)BL,j(2D)R2i+∑i=0L-1CL,i(2D)R2i,
where (·)!! is the double factorial defined as
(14)n!!{n·(n-2)⋯5·3·1,nisanoddpositiveintegern·(n-2)⋯6·4·2,nisanevenpositiveinteger1,n=0,-1(-1)(n+1)/2(-n-2)!!,n≤-3isanoddnegativeinteger.
The detailed derivation is given in Appendix A. In order to ensure the solution ΨL(2D) satisfies the hierarchical relation (6) and (7), we can make the following selection:
(15)CL,i(2D)=∑j=0i(-1)i-j(2i-2j-3)!!AL,j(2D)(i-j)!2i-j+∑j=0i-1(-1)i-j(2i-2j-1)!!BL,j(2D)(i-j)!2i-j+1(i-j)fori≤L-1.
Then, (13) becomes
(16)ΨL(2D)=∑i=L∞∑j=0L(-1)i-j+1(2i-2j-3)!!(i-j)!2i-jAL,j(2D)R2i+∑i=L∞∑j=0L-1(-1)i-j+1(2i-2j-1)!!(i-j)!2i-j+1(i-j)BL,j(2D)R2i.
The PMQ ΨL(2D), given in (16) or equivalently in (12) with CL,j(2D) defined in (15), is the unique infinitely differentiable PMQ with its first L even coefficients of the Maclaurin series being zero if the undetermined coefficients AL,j(2D) and BL,j(2D) can be found by (8). The proposed condition can be alternatively stated as
(17)(ΔR(2D))jΨL(2D)(R⟶0)=0forj=0,1,…,L-1,
where the following relation has to be considered:
(18)(ΔR(2D))jR2i={4ji!2(i-L)!2R2i-2jfori≥j0fori<j.
Applying ΔR(2D) to the PMQ ΨL(2D) in (16) and using (18) indicate that ΨL-1(2D) also satisfies (17) and of the same form in (16). Therefore, it can be understood that the prescribed PMQ also satisfies the hierarchical relation (6) and (7).In order to determineAL,j(2D) and BL,j(2D), (16) is substituted into (8) which results in
(19)∑i=L∞∑j=0L(-1)i-j+1(2i-2j-3)!![i!2](i-j)!2i-j-2L[(i-L)!2]AL,j(2D)R2i-2L+∑i=L∞∑j=0L-1(-1)i-j+1(2i-2j-1)!![i!2](i-j)!2i-j-2L+1(i-j)[(i-L)!2]BL,j(2D)R2i-2L=∑i=0∞(-1)i-1(2i-3)!!i!2iR2i,
where (18) and the following Maclaurin series are used:
(20)1+R2=∑i=0∞(12i)R2i=∑i=0∞(-1)i-1(2i-3)!!i!2iR2i,
with (··) being the usual combination formula.Substitutingi by i+L in the left-hand side of (19) and comparing the coefficients of R2i, we have
(21)∑j=0L(-1)i+L-j+1(2i+2L-2j-3)!![(i+L)!2](i+L-j)!2i-L-j[i!2]AL,j(2D)+∑j=0L-1(-1)i+L-j+1(2i+2L-2j-1)!![(i+L)!2](i+L-j)!2i-L-j+1(i+L-j)[i!2]BL,j(2D)=(-1)i-1(2i-3)!!i!2i
or equivalently
(22)∑j=0L(-1)L-j2L+j(i+L)!i!︸Lthorders(i+L)!(i+L-j)!︸jthorders×(2i+2L-2j-3)!!(2i-3)!!︸(L-j)thordersAL,j(2D)+∑j=0L-1(-1)L-j2L+j-1(i+L)!i!(i+L-j)︸(L-1)thorders(i+L)!(i+L-j)!︸jthorders×(2i+2L-2j-1)!!(2i-3)!!︸(L+1-j)thordersBL,j(2D)=1.
The identity (22) should be satisfied for every nonnegative integer i to ensure the PMQ ΨL(2D) is the solution of (8). Since the left-hand side of (22) is a 2L-ordered polynomial of the variable i, it can be used to uniquely determine the 2L+1 unknown coefficients AL,j(2D) and BL,j(2D) if the system is not singular. AfterAL,j(2D) and BL,j(2D) are obtained, (15) can be used to have CL,j(2D). Therefore, the two-dimensional PMQ ΨL(2D) can be computed simply by (12). For the cases when R is very small, one can alternatively use (16).In practice, one can simply seti equal to 1,2,…,2L+1 in (22) which results in 2L+1 linear equations for solving the 2L+1 unknown coefficients AL,j(2D) and BL,j(2D) of the PMQ of order L. Mathematica codes for obtaining the 2D-PMQ are given in Algorithm 1. Alternatively, if a specific low order of 2D-PMQ is considered, one can use the explicit coefficients given in Tables 1, 2, and 3 and the results are plotted in Figure 1. Basically, the PMQs of higher orders are more flat in the near field and grow faster in the far field, as shown in the figure.Table 1
A
L
,
j
(
2
D
).
j
=
0
j
=
1
j
=
2
j
=
3
j
=
4
j
=
5
L
=
1
4
9
1
9
L
=
2
-
61
900
4
75
1
225
L
=
3
1019
352800
-
153
15680
8
3675
1
11025
L
=
4
-
3407
57153600
229
510300
-
33007
76204800
8
178605
1
893025
L
=
5
80939
110649369600
-
9535
983549952
735263
35407798272
-
3311969
354077982720
4
7203735
1
108056025Table 2
B
L
,
j
(
2
D
).
j
=
0
j
=
1
j
=
2
j
=
3
j
=
4
L
=
1
-
1
3
L
=
2
1
30
-
1
12
L
=
3
-
1
840
1
120
-
1
192
L
=
4
1
45360
-
1
3360
1
1920
-
1
6912
L
=
5
-
1
3991680
1
181440
-
1
53760
1
69120
-
1
442368Table 3
C
L
,
j
(
2
D
).
j
=
0
j
=
1
j
=
2
j
=
3
j
=
4
L
=
1
-
4
9
L
=
2
61
900
-
1
36
L
=
3
-
1019
352800
31
3600
1
1152
L
=
4
3407
57153600
-
599
1411200
1
3600
1
13824
L
=
5
-
80939
110649369600
2147
228614400
-
389
22579200
1
345600
1
589824Algorithm 1:
Mathematica codes for obtaining 2D-PMQs.
L
=
3;  (*set an arbitrary L*)
F
[
j
_
]
≔
(
-
1
)
L
-
j
2
L
+
j
(
2
i
+
2
L
-
2
j
-
3
)
!!
(
(
i
+
L
)
!
)
2
(
i
+
L
-
j
)
!
(
i
!
)
(
2
i
-
3
)
!!
A
[
L
,
j
];
G
[
j
_
]
≔
(
-
1
)
L
-
j
2
L
+
j
-
1
(
2
i
+
2
L
-
2
j
-
1
)
!!
(
(
i
+
L
)
!
)
2
(
i
+
L
-
j
)
!
(
i
+
L
-
j
)
(
i
!
)
(
2
i
-
3
)
!!
B
[
L
,
j
];
F
G
[
i
_
]
≔
∑
j
=
0
L
F
[
j
]
+
∑
j
=
0
L
-
1
G
[
j
];
sol = Solve[Table[FG[i]==1,{i,1,2L+1}],
Join[Table[A[L,j],{j,0,L}],  Table[B[L,j],{j,0,L-1}]]];
Do
[
A
A
[
L
,
j
]
=
(
A
[
L
,
j
]
/
.
s
o
l
)
[
[
1
]
]
,
{
j
,
0
,
L
}
];
Do
[
B
B
[
L
,
j
]
=
(
B
[
L
,
j
]
/
.
s
o
l
)
[
[
1
]
]
,
{
j
,
0
,
L
-
1
}
];
C
C
[
L
_
,
i
_
]
≔
∑j=0i(-1)i-j(2i-2j-3)!!AA[L,j](i-j)!2i-j+∑j=0i-1(-1)i-j(2i-2j-1)!!BB[L,j](i-j)!2i-j+1(i-j);
Phi
L
=
∑
j
=
0
L
A
A
[
L
,
j
]
R
2
j
1
+
R
2
+
∑
j
=
0
L
-
1
B
B
[
L
,
j
]
L
o
g
[
1
+
1
+
R
2
2
]
R
2
j
+
∑
j
=
0
L
-
1
C
C
[
L
,
j
]
R
2
j
Series
[
Phi
L
,
{
R
,
0,2
L
+
2
}
]
PhiTemp = PhiL;(*Check thePMQ*)
Do
[
P
h
i
0
=
∂
R
(
R
∂
R
P
h
i
T
e
m
p
)
R
;
P
h
i
T
e
m
p
=
P
h
i
0
,
{
j
,
1
,
L
}
];
P
rint["LaplaceL  PhiL  = ",  FullSimplify[Phi0]]Plots of the 2D-PMQs in (a) near field and (b) far field.
(a)
(b)


## 4. Three-Dimensional PMQ

Similarly, the three-dimensional PMQs can be obtained by the method of undetermined coefficients which begins with(23)ΨL(3D)=∑j=0LAL,j(3D)R2j1+R2+∑j=0L-1BL,j(3D)R2jsinh-1RR+∑j=0L-1CL,j(3D)R2j+∑j=0L-1DL,j(3D)R2j-1.
Here, the first two series form a particular solution with undetermined coefficients AL,j(3D) and BL,j(3D) and the last two series give the homogeneous solutions with 2L integration constants CL,j(3D) and DL,j(3D). Here, it should be noted that the three-dimensional PMQ involving the inverse hyperbolic sine function has been recently derived in the author’s previous study [24] and will be further extended in this study. To determine AL,j(3D), BL,j(3D), CL,j(3D), and DL,j(3D), we first expand (23) in terms of the Laurent series as
(24)ΨK,L(3D)=∑i=0∞∑j=0min{L,i}(-1)i-j-1(2i-2j-3)!!(i-j)!2i-jAL,j(3D)R2i+∑i=0∞∑j=0min{L-1,i}(-1)i-j(2i-2j-1)!!(i-j)!2i-j(2i-2j+1)BL,j(3D)R2i+∑i=0L-1CL,i(3D)R2i+∑i=0L-1DL,i(3D)R2i-1.
The detailed derivation is given in Appendix B. In order to make ΨL(3D) infinitely differentiable, it is required to cancel the singular term R-1 and all the odd coefficients in (24) [26]. In other words, we need
(25)DL,i(3D)=0.
Similar to (17), we additionally consider the following relation:
(26)(ΔR(3D))jΨL(3D)(R⟶0)=0forj=0,1,…,L-1
with
(27)(ΔR(3D))jR2i={(2i+1)!(2i-2j+1)!R2i-2jfori≥j0fori<j.
Then we use (24) and (26) to have the following relation:
(28)CL,i(3D)=∑j=0i(-1)i-j(2i-2j-3)!!AL,j(3D)(i-j)!2i-j+∑j=0i(-1)i-j-1(2i-2j-1)!!BL,j(3D)(i-j)!2i-j(2i-2j+1)fori≤L-1.
Substituting (28) into (24) gives
(29)ΨL(3D)=∑i=L∞∑j=0L(-1)i-j-1(2i-2j-3)!!(i-j)!2i-jAL,j(3D)R2i+∑i=L∞∑j=0L-1(-1)i-j(2i-2j-1)!!(i-j)!2i-j(2i-2j+1)BL,j(3D)R2i.
The PMQ ΨL(3D), given in (29) or equivalently in (23) with CL,i(3D) and DL,i(3D) defined in (28) and (25), is the only infinitely differentiable PMQ satisfying (26) if the undetermined coefficients AL,j(3D) and BL,j(3D) can be found by using (8). Then, we consider the following relation.Equation (27) with j=1 together with (29) is sufficient to infer that the prescribed PMQ also satisfies the hierarchical relation (6) and (7) as explained in the two-dimensional case. Substituting (29) into (8) and using (20) and (27) result in
(30)∑i=L∞∑j=0L(-1)i-j-1(2i-2j-3)!!(2i+1)!(i-j)!2i-j(2i-2L+1)!AL,j(3D)R2i-2L+∑i=L∞∑j=0L-1(-1)i-j(2i-2j-1)!!(2i+1)!(i-j)!2i-j(2i-2j+1)(2i-2L+1)!BL,j(3D)R2i-2L=∑i=0∞(-1)i-1(2i-3)!!i!2iR2i.
Substituting i by i+L for the left-hand side of (30) and comparing the coefficients of R2i, we have
(31)∑j=0L(-1)i+L-j-1(2i+2L-2j-3)!!(2i+2L+1)!AL,j(3D)(i+L-j)!2i+L-j(2i+1)!+∑j=0L-1(-1)i+L-j(2i+2L-2j-1)!!(2i+2L+1)!BL,j(3D)(i+L-j)!2i+L-j(2i+2L-2j+1)(2i+1)!=(-1)i-1(2i-3)!!i!2i
or equivalently
(32)∑j=0L(-1)L-j2j(2i+2L+1)!!(2i+1)!!︸Lthorders(i+L)!(i+L-j)!︸jthorders×(2i+2L-2j-3)!!(2i-3)!!︸(L-j)thordersAL,j(3D)+∑j=0L-1(-1)L-j+12j(2i+2L+1)!!(2i+1)!!(2i+2L-2j+1)︸(L-1)thorders(i+L)!(i+L-j)!︸jthorders×(2i+2L-2j-1)!!(2i-3)!!︸(L+1-j)thordersBL,j(3D)=1.
Similar to the two-dimensional cases, (32) can be used to uniquely determine the 2L+1 unknown coefficients AL,j(3D) and BL,j(3D) since its left-hand side is a 2L-ordered polynomial of the variable i. Similarly, Mathematica codes for obtaining the 3D-PMQ are given in Algorithm 2. Also, the first few orders of 3D-PMQ are given in Tables 4, 5, and 6 and plotted in Figure 2. In the figure, we can also observe that the PMQs of higher orders are more flat in the near field and grow faster in the far field.Table 4
A
L
,
j
(
3
D
).
j
=
0
j
=
1
j
=
2
j
=
3
j
=
4
j
=
5
L
=
1
5
24
1
12
L
=
2
-
9
160
7
360
1
360
L
=
3
919
322560
-
391
53760
3
4480
1
20160
L
=
4
-
3781
58060800
5777
14515200
-
1457
4838400
11
907200
1
1814400
L
=
5
122987
143061811200
-
406247
42918543360
926801
53648179200
-
55151
8941363200
13
95800320
1
239500800Table 5
B
L
,
j
(
3
D
).
j
=
0
j
=
1
j
=
2
j
=
3
j
=
4
L
=
1
1
8
L
=
2
-
1
96
1
16
L
=
3
1
3072
-
1
192
1
192
L
=
4
-
1
184320
1
6144
-
1
2304
1
5760
L
=
5
1
17694720
-
1
368640
1
73728
-
1
69120
1
322560Table 6
C
L
,
j
(
3
D
).
j
=
0
j
=
1
j
=
2
j
=
3
j
=
4
L
=
1
-
1
3
L
=
2
1
15
-
1
18
L
=
3
-
1
315
1
90
-
1
360
L
=
4
1
14175
-
1
1890
1
1800
-
1
15120
L
=
5
-
1
1091475
1
85050
-
1
37800
1
75600
-
1
1088640Algorithm 2:
Mathematica codes for obtaining 3D-PMQs.
L
=
3;  (*set an arbitrary L*)
F
[
j
_
]
:
=
(
-
1
)
L
-
j
2
-
j
(
2
i
+
2
L
+
1
)
!!
(
2
i
+
1
)
!!
(
i
+
L
)
!
(
i
+
L
-
j
)
!
(
2
i
+
2
L
-
2
j
-
3
)
!!
(
2
i
-
3
)
!!
A
[
L
,
j
];
G
[
j_]:=
(-1)L-j+12-j(2i+2L+1)!!(2i+1)!!(2i+2L-2j+1)(i+L)!(i+L-j)!(2i+2L-2j-1)!!(2i-3)!!B[L,j];
F
G
[
i
_
]
:
=
∑
j
=
0
L
F
[
j
]
+
∑
j
=
0
L
-
1
G
[
j
];
s
o
l
=
S
o
l
v
e
[
T
a
b
l
e
[
F
G
[
i
]
=
=
1
,
{
i
,
1,2
L
+
1
}
],
Join[Table[A[L,j],{j,0,L}],Table[B[L,j],{j,0,L-1}]]];
D
o
[
A
A
[
L
,
j
]
=
(
A
[
L
,
j
]
/
.
s
o
l
)
[
[
1
]
]
,
{
j
,
0
,
L
}
];
D
o
[
B
B
[
L
,
j
]
=
(
B
[
L
,
j
]
/
.
s
o
l
)
[
[
1
]
]
,
{
j
,
0
,
L
-
1
}
];
C
C
[
L_,i_]:=
∑j=0i(-1)i-j(2i-2j-3)!!AA[L,j](i-j)!2i-j+∑j=0i(-1)i-j-1(2i-2j-1)!!BB[L,j](i-j)!2i-j(2i-2j+1);
P
h
i
L
=
∑
j
=
0
L
A
A
[
L
,
j
]
R
2
j
1
+
R
2
+
∑
j
=
0
L
-
1
B
B
[
L
,
j
]
A
r
c
S
i
n
h
[
R
]
R
R
2
j
+
∑
j
=
0
L
-
1
C
C
[
L
,
j
]
R
2
j
Series
[
P
h
i
L
,
{
R
,
0,2
L
+
2
}
]
PhiTemp = PhiL;
(*Check the
P
MQ
*)
Do
[
P
h
i
0
=
∂
R
(
R
2
∂
R
P
h
i
T
e
m
p
)
R
2
;
P
h
i
T
e
m
p
=
P
h
i
0
,
{
j
,
1
,
L
}
];
P
rint["LaplaceL  PhiL  = ",  FullSimplify[Phi0]]Plots of the 3D-PMQs in (a) near field and (b) far field.
(a)
(b)


## 5. Validation

By the standard MFS-DRM, the PMQs derived in this paper have been validated for the cases of fourth and third orders, respectively, in two and three dimensions [24]. The results basically indicate that the solutions of the MQ are more accurate than those of the low-ordered APS. Numerical results for other orders of polyharmonic problems should be similar and thus will not be repeated here.Alternatively, we apply the PMQs to the standard MFS-DRM for solving Poisson’s equation(33)Δu=f(x);x∈Ω,
where Ω⊂R2orR3 is the considered domain, f is the given forcing term, and Δ denotes the Laplacian operator. For simplicity, we consider the boundary condition
(34)u=g(x);x∈Γ,
where Γ is the boundary of Ω and g is the given boundary data. In the formal MFS-DRM, the solution should be decomposed into two parts
(35)u=up+uh,
where the particular solution up satisfies
(36)Δup=f(x);x∈Ω,
but not necessarily the boundary condition (34). And the homogeneous solution uh is defined by
(37)Δuh=0;x∈Ω,uh=g(x)-up;x∈Γ.
The solution of (37) can be formally obtained by the MFS [7–11] and will not be additionally described.On the other hand, in order to solve the particular solution, we approximate the forcing term by(38)f(x)≈∑i=1MαiψL(*)(ri),
where ri=∥x-xi∥ is the Euclidean distance between the spatial coordinates x and M prescribed points xi in Ω and ψL(*) is the PMQ of order L. In order to solve the unknown coefficients αi, (38) should be collocated on the M prescribed points xi as
(39)f(xj)=∑i=1MαiψL(*)(rij)forj=1,2,…,M,
where rij=∥xj-xi∥. After the forcing term is approximated in (38), the particular solution can be correspondingly approximated as
(40)up(x)≈∑i=1MαiψL+1(*)(ri).ForL=0, the solution procedure is reduced to the formal MFS-DRM based on the MQ [18].


### 5.1. Two-Dimensional Case

We consider a two-dimensional problem governed by Poisson’s equation (33) with Ω=[-(1/2),1/2]2and(41)f(x)=f(x,y)=4cos(2x2+y)-(1+16x2)sin(2x2+y).
The boundary condition (34) is set according to the exact solution as
(42)u(x,y)=g(x,y)=sin(2x2+y).
Since the maximum value of the solution in Ω is equal to sin1≈0.84, we use the maximum error in this example. And the sources of the MFS are typically located on the boundary of a square with edge length equal to 5.We have hierarchically performed our tests forL=0,1,2,…,8. Among the similar results, the maximum errors for L=0,2,4,6,8 are addressed in Table 7. This should have validated the two-dimensional PMQs for L=0,1,2,…,9. For every individual test, a preliminary search is always performed for finding the optimal shape parameter c as demonstrated in Figure 3. It can be observed from Table 7 that the best maximum errors (marked by *) usually occur at L=6 and are about an order of magnitude better than those obtained by the MQ and the case of more collocation points (larger M) gives better accuracy. If we observe the PMQ coefficients, we will find that eight is the minimum order of the PMQ with some coefficients over the limit of 64-bit unsigned integer. To be more detailed, A8,5(2D) is equal to 1457447323/28298491779317760000 whose denominator is just larger than 264-1≈1.8×1019.Table 7
Maximum errors of the two-dimensional Poisson's problem.
M
=
9
×
9
M
=
11
×
11
M
=
13
×
13
L
=
0
1.07
E
-
06
1.64
E
-
07
6.38
E
-
08
L
=
2
8.23
E
-
07
1.22
E
-
07
4.99
E
-
08
L
=
4
3.50
E
-
07
2.85
E
-
08
1.47
E
-
08
L
=
6
2.08
E
-
07
*
1.12
E
-
08
*
5.06
E
-
09
L
=
8
2.53
E
-
07
1.43
E
-
08
2.25
E
-
09
*
*Means the best result.Search for the optimal shape parameterc for the 2D case with (a) M=9×9 and (b) M=13×13.
(a)
(b)


### 5.2. Three-Dimensional Case

Then, we consider a three-dimensional problem governed by Poisson’s equation (33) with Ω=[-(1/2),1/2]3and(43)f(x)=f(x,y,z)=6(4+x+y+z)3.The boundary condition (34) is set according to the exact solution as
(44)u(x,y,z)=g(x,y,z)=14+x+y+z.
Since the maximum value of the solution in the computational domain is equal to 0.25, we also use the maximum error in this example. And the sources of the MFS are typically located on the boundary of a cuboid with edge length equal to 5.Similarly, we also hierarchically perform numerical tests forL=0,1,2,…,8 and address the maximum errors for L=0,2,4,6,8 in Table 8. This should also validate the three-dimensional PMQs for L=0,1,2,…,9. Figure 4 gives the preliminary search for finding the optimal shape parameter c. It can also be observed from Table 8 that the best maximum errors basically occur at L=6 and are about an order of magnitude better than those obtained by the MQ. Similarly, we have A8,5(3D) equal to 977222021/25635104317734912000 whose denominator is larger than 264-1≈1.8×1019.Table 8
Maximum errors of the three-dimensional Poisson's problem.
M
=
5
×
5
×
5
M
=
7
×
7
×
7
M
=
9
×
9
×
9
L
=
0
1.23
E
-
06
9.22
E
-
08
1.73
E
-
08
L
=
2
2.03
E
-
06
6.67
E
-
08
1.66
E
-
08
L
=
4
4.61
E
-
07
1.07
E
-
08
1.02
E
-
09
L
=
6
2.63
E
-
07
*
3.68
E
-
08
6.07
E
-
10
*
L
=
8
1.09
E
-
02
2.03
E
-
08
*
7.55
E
-
10
*Means the best result.Search for the optimal shape parameterc for the 3D case with (a) M=5×5×5 and (b) M=9×9×9.
(a)
(b)It seems that the PMQ of higher order gives more accurate results if there is no numerical ill-conditioning. However for the case ofM=5×5×5 and L=8, the error goes worse as shown in Table 8. This result seems to indicate that the PMQ of higher order also requires more collocation points. In practice, one can safely choose L=4~6 and find some improvement over the MQ.


## 5.1. Two-Dimensional Case

We consider a two-dimensional problem governed by Poisson’s equation (33) with Ω=[-(1/2),1/2]2and(41)f(x)=f(x,y)=4cos(2x2+y)-(1+16x2)sin(2x2+y).
The boundary condition (34) is set according to the exact solution as
(42)u(x,y)=g(x,y)=sin(2x2+y).
Since the maximum value of the solution in Ω is equal to sin1≈0.84, we use the maximum error in this example. And the sources of the MFS are typically located on the boundary of a square with edge length equal to 5.We have hierarchically performed our tests forL=0,1,2,…,8. Among the similar results, the maximum errors for L=0,2,4,6,8 are addressed in Table 7. This should have validated the two-dimensional PMQs for L=0,1,2,…,9. For every individual test, a preliminary search is always performed for finding the optimal shape parameter c as demonstrated in Figure 3. It can be observed from Table 7 that the best maximum errors (marked by *) usually occur at L=6 and are about an order of magnitude better than those obtained by the MQ and the case of more collocation points (larger M) gives better accuracy. If we observe the PMQ coefficients, we will find that eight is the minimum order of the PMQ with some coefficients over the limit of 64-bit unsigned integer. To be more detailed, A8,5(2D) is equal to 1457447323/28298491779317760000 whose denominator is just larger than 264-1≈1.8×1019.Table 7
Maximum errors of the two-dimensional Poisson's problem.
M
=
9
×
9
M
=
11
×
11
M
=
13
×
13
L
=
0
1.07
E
-
06
1.64
E
-
07
6.38
E
-
08
L
=
2
8.23
E
-
07
1.22
E
-
07
4.99
E
-
08
L
=
4
3.50
E
-
07
2.85
E
-
08
1.47
E
-
08
L
=
6
2.08
E
-
07
*
1.12
E
-
08
*
5.06
E
-
09
L
=
8
2.53
E
-
07
1.43
E
-
08
2.25
E
-
09
*
*Means the best result.Search for the optimal shape parameterc for the 2D case with (a) M=9×9 and (b) M=13×13.
(a)
(b)


## 5.2. Three-Dimensional Case

Then, we consider a three-dimensional problem governed by Poisson’s equation (33) with Ω=[-(1/2),1/2]3and(43)f(x)=f(x,y,z)=6(4+x+y+z)3.The boundary condition (34) is set according to the exact solution as
(44)u(x,y,z)=g(x,y,z)=14+x+y+z.
Since the maximum value of the solution in the computational domain is equal to 0.25, we also use the maximum error in this example. And the sources of the MFS are typically located on the boundary of a cuboid with edge length equal to 5.Similarly, we also hierarchically perform numerical tests forL=0,1,2,…,8 and address the maximum errors for L=0,2,4,6,8 in Table 8. This should also validate the three-dimensional PMQs for L=0,1,2,…,9. Figure 4 gives the preliminary search for finding the optimal shape parameter c. It can also be observed from Table 8 that the best maximum errors basically occur at L=6 and are about an order of magnitude better than those obtained by the MQ. Similarly, we have A8,5(3D) equal to 977222021/25635104317734912000 whose denominator is larger than 264-1≈1.8×1019.Table 8
Maximum errors of the three-dimensional Poisson's problem.
M
=
5
×
5
×
5
M
=
7
×
7
×
7
M
=
9
×
9
×
9
L
=
0
1.23
E
-
06
9.22
E
-
08
1.73
E
-
08
L
=
2
2.03
E
-
06
6.67
E
-
08
1.66
E
-
08
L
=
4
4.61
E
-
07
1.07
E
-
08
1.02
E
-
09
L
=
6
2.63
E
-
07
*
3.68
E
-
08
6.07
E
-
10
*
L
=
8
1.09
E
-
02
2.03
E
-
08
*
7.55
E
-
10
*Means the best result.Search for the optimal shape parameterc for the 3D case with (a) M=5×5×5 and (b) M=9×9×9.
(a)
(b)It seems that the PMQ of higher order gives more accurate results if there is no numerical ill-conditioning. However for the case ofM=5×5×5 and L=8, the error goes worse as shown in Table 8. This result seems to indicate that the PMQ of higher order also requires more collocation points. In practice, one can safely choose L=4~6 and find some improvement over the MQ.


## 6. Conclusion

The polyharmonic multiquadrics (PMQs) are derived. The methods of undetermined coefficients are constructed and the PMQs are solved by the Laurent expansion. Some additional conditions are introduced to make the derived PMQs hierarchically unique and infinitely differentiable. Numerical experiments were carried out to validate the derived PMQs. Basically, the numerical results obtained by the PMQ showed some improvement over those by the multiquadrics.The method introduced in this study can be extended to derive the polyharmonic particular solutions of the generalized multiquadrics(r2+c2)K/2 with K=-1,1,3,… which is currently under investigation.



---

*Source: 613082-2013-05-22.xml*