# Solving the Linear Integral Equations Based on Radial Basis Function Interpolation

Authors: Huaiqing Zhang; Yu Chen; Xin Nie
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

The radial basis function (RBF) method, especially the multiquadric (MQ) function, was introduced in solving linear integral equations. The procedure of MQ method includes that the unknown function was firstly expressed in linear combination forms of RBFs, then the integral equation was transformed into collocation matrix of RBFs, and finally, solving the matrix equation and an approximation solution was obtained. Because of the superior interpolation performance of MQ, the method can acquire higher precision with fewer nodes and low computations which takes obvious advantages over thin plate splines (TPS) method. In implementation, two types of integration schemes as the Gauss quadrature formula and regional split technique were put forward. Numerical results showed that the MQ solution can achieve accuracy of1
E
-
5. So, the MQ method is suitable and promising for integral equations.


---

## Body

## 1. Introduction

In mechanics, electromagnetic theory, thermal, radiation science, and other fields, the mathematical models of science or engineering problems can be attributed to the integral equations. Comparing to the differential equation, the integral equation embodies the initial value or boundary information just in a brief form. Moreover, the numerical integral brings smaller relative error than the numerical differentiation. So it can obtain higher accuracy, and the research on numerical solution of integral equations is of great importance. The general form for linear integral equation is(1)
f
(
P
)
=
μ
∫
Ω
k
(
P
,
Q
)
f
(
Q
)
d
Q
+
g
(
P
)
P
,
Q
∈
Ω
.
In the above, g
(
P
) is free function, k
(
P
,
Q
) is kernel function, f
(
Q
) is the unknown function, Ω is integral region, and μ is the parameter where the g
(
P
), k
(
P
,
Q
), and μ are known. The numerical methods for solving integral equations mainly include the undetermined coefficient approximation method, the direct numerical integration method, the successive approximation method, and the approximation method for specific kernel function.Currently, a variety of new numerical methods were proposed, for instance, the combination method of Galerkin and collocation with Haar wavelet function for the first kind of linear Fredholm equation [1], the Galerkin method with Daubechies wavelet for the second kind of linear Volterra equations [2], however, the calculation error of wavelet methods is relatively large. The Sinc-collocation method [3] and the modified block pulse function [4] are for one dimension problems. The Chebyshev collocation method [5], block pulse method [6], and polynomial approximation with smooth kernel function [7] are for solving two-dimensional Fredholm or Volterra equations. In the above unknown function approximation methods, the polynomial, Chebyshev polynomials, or wavelet function, were adopted, respectively. However, the polynomial interpolation is unstable due to the strong rigidity of the polynomial function space, the interpolation error is relatively large for wavelets, the suitable interpolation basis functions for two or multidimensional problems are hard to construct, and so on. The interpolation basis function with high accuracy and convenience for multidimensions is particularly necessary.And, fortunately, the radial basis function (RBF) not only is insensitive to the space dimensions but also has high interpolation accuracy. The radial basis functions were first developed by Hardy [8] in 1971 as a multidimensional scattered interpolation method in modeling of the earth’s gravitational field. It was not recognized by most of the academic researchers until Franke [9] published a review paper in the evaluation of two-dimensional interpolation methods. The RBF method in solving integral equations was initially proposed in 2006 [10–13]. The related research attracted a lot of attention recently. Therefore, this paper attempted to adopt the RBF method for solving linear integral equations. The RBF and interpolation principle were described in Section 2. And, in Section 3, the RBF method for solving linear integral equation and the numerical integration scheme for coefficients matrix were proposed. And, finally, the RBF method was applied in one- and two-dimensional integral equations; its performance was analyzed and compared.


## 2. Radial Basis Function and Interpolation Principle

Radial basis functionФ
:
R
+
→
R (domain: R
d) is defined as the function of distance r
=
|
|
x
-
x
j
|
|. The commonly used RBF contains the thin plate splines function (
Ф
(
r
)
=
r
2
ln
⁡
r
), the Gaussian function (
Ф
(
r
)
=
exp
⁡
(
-
β
r
2
)
), Hardy’s multiquadric (MQ) functions, and so forth. Franke has proved that the MQ method has the superior comprehensive performance in 29 kinds of scattered data interpolation methods. So, in this paper, we choose the MQ as the basis function which is compared with thin plate splines (TPS) function. The widely used MQ function has the following expression:
(2)
ϕ
(
x
)
=
∥
x
-
c
∥
2
+
α
2
.
In the above, c means the center of basis function and α is the shape parameter which is generally associated with the distance between adjacent centers. So α
=
β
|
|
c
i
-
c
j
|
|, and β is also called the shape parameter.The principle of RBF interpolation regards the unknown function as linear combination of radial functions. So the approximation function can be obtained after calculating the coefficients. For a series of known data points(
x
i
,
f
(
x
i
)
), i
=
1,2
,
…
,
N, the approximation function can be constructed by RBF as
(3)
f
^
(
x
)
=
∑
j
=
1
N
λ
j
ϕ
j
(
x
)
=
∑
j
=
1
N
λ
j
(
x
-
c
j
)
2
+
α
2
.
Then, substituting the interpolation data points into the above equation, we obtain
(4)
f
(
x
i
)
=
∑
j
=
1
N
λ
j
ϕ
j
(
x
i
)
=
∑
j
=
1
N
λ
j
(
x
i
-
c
j
)
2
+
α
2
.
With its abstract matrix form,
(5)
[
Φ
d
]
[
λ
]
=
[
f
d
]
,
where [
f
d
]
=
[
f
(
x
1
)
,
f
(
x
1
)
,
…
,
f
(
x
N
)
]
T is the known data, [
λ
]
=
[
λ
1
,
λ
2
,
…
,
λ
N
]
T is the coefficients, and [
Φ
d
] is the coefficient matrix whose element calculation formula is Φ
m
n
=
(
x
m
-
c
n
)
2
+
α
2.The MQ function is globally defined which results in a full resultant coefficient matrix. And, consequently, the coefficient matrix in (5) is nonsingular and usually ill conditioned. The accuracy of the MQ solution depends heavily on the shape parameter. In general, for a fixed number of centers N, smaller shape parameters produce more accurate approximations, but they also are associated with a poorly conditioned interpolation matrix. Also, the condition number grows with N, for fixed values of the shape parameter α. The choice of this optimal value is still under intensive investigation. Therefore, solving (5) can acquire
(6)
[
λ
]
=
[
Φ
d
]
-
1
[
f
d
]
.
The approximation value f
^
(
x
) at any point x is
(7)
f
^
(
x
)
=
[
Φ
(
x
)
]
[
λ
]
=
[
Φ
(
x
)
]
[
Φ
d
]
-
1
[
f
d
]
.


## 3. Solving Linear Integral Equations Based on the RBF Interpolation

The basic idea of RBF method for solving integral equation is to employ a linear combination of RBF to approximate the unknown function; thus, the integral equation is transformed into the combination form of RBFs and their coefficients. Then, the weight coefficients of RBF can be calculated by using the collocation type weighted residual method. And, finally, an approximate representation of the unknown function can be obtained. Specifically, the RBF method for solving linear integral equations is proposed in the following.


### 3.1. RBF Interpolation to Approximate One-Dimensional Linear Integral Equation

The general form of one-dimensional linear integral equation is(8)
f
(
x
)
=
μ
∫
a
b
∣
x
k
(
x
,
t
)
f
(
t
)
d
t
+
g
(
x
)
,
x
∈
[
a
,
b
]
.
Here, b
∣
x represents the upper limit of integral. The constant b or variable x is, respectively, corresponding to Fredholm or Volterra equation. In the following discussion, we choose the upper limit b for demonstration. In RBF approximation, the unknown function f
(
x
) is expressed as a combination of RBF. Thus, the approximate formula (9) is obtained:
(9)
∑
i
=
1
N
λ
i
ϕ
i
(
x
)
=
μ
∫
a
b
k
(
x
,
t
)
∑
i
=
1
n
λ
i
ϕ
i
(
t
)
d
t
+
g
(
x
)
,
x
∈
[
a
,
b
]
.Furthermore, for any collocation pointx
d in [
a
,
b
], by using the linear property of integral, we obtain the specific expression of the above equation as
(10)
∑
i
=
1
N
λ
i
ϕ
i
(
x
d
)
-
μ
∑
i
=
1
N
λ
i
∫
a
b
k
(
x
d
,
t
)
ϕ
i
(
t
)
d
t
=
g
(
x
d
)
,
11111111111111111111111111111111
x
d
∈
[
a
,
b
]
,
when choosing M collocation points, hence forming the collocation equations in matrix form as
(11)
[
Φ
-
μ
Κ
]
[
λ
]
=
[
G
]
,
where [
G
]
=
[
g
(
x
1
)
,
g
(
x
1
)
,
…
,
g
(
x
M
)
]
T, and [
λ
]
=
[
λ
1
,
λ
2
,
…
,
λ
N
]
T. The elements in matrix [
Ф
] and [
K
] are defined:
(12)
Φ
m
n
=
(
x
m
-
c
n
)
2
+
α
2
K
m
n
=
∫
a
b
k
(
x
m
,
t
)
(
t
-
c
n
)
2
+
α
2
d
t
=
∫
a
b
h
(
t
)
d
t
.
The notations x
m and c
n represent correspondingly the collocation point and RBF center, while the element K
m
n is the definite integral of function h
(
t
) which can be calculated through the following Gauss quadrature formula:
(13)
∫
-
1
1
h
(
ξ
)
d
ξ
≈
∑
q
=
1
Q
W
q
h
(
ξ
q
)
.
The integral range can be transformed with formula t
=
(
b
+
a
)
/
2
+
(
b
-
a
)
ξ
/
2
=
p
(
ξ
), so the final calculation formula is
(14)
K
m
n
=
∫
a
b
k
(
x
m
,
t
)
ϕ
n
(
t
)
d
t
=
b
-
a
2
∫
-
1
1
k
[
x
m
,
p
(
ξ
)
]
ϕ
n
[
p
(
ξ
)
]
d
ξ
≈
b
-
a
2
∑
q
=
1
Q
W
q
k
[
x
m
,
p
(
ξ
q
)
]
ϕ
n
[
p
(
ξ
q
)
]
.


### 3.2. RBF Interpolation to Approximate Two-Dimensional Linear Integral Equation

Two-dimensional linear integral equation has the general form(15)
f
(
x
,
y
)
=
μ
∫
Ω
k
(
x
,
y
,
s
,
t
)
f
(
s
,
t
)
d
s
d
t
+
g
(
x
,
y
)
,
hhhhhhhhhhhhhhhhhhhhhhhhhhhh
(
x
,
y
)
∈
Ω
.
The corresponding RBF approximation is
(16)
∑
i
=
1
N
λ
i
ϕ
i
(
x
d
,
y
d
)
=
μ
∑
i
=
1
N
λ
i
∫
Ω
k
(
x
d
,
y
d
,
s
,
t
)
ϕ
i
(
s
,
t
)
d
s
d
t
+
g
(
x
d
,
y
d
)
.And, similarly, we have the interpolation matrix[
Φ
-
μ
Κ
]
[
λ
]
=
[
G
]. The differences lie in the distance formula and the definite integral of element K
m
n
=
∫
Ω
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t. According to the domain type, two integration schemes were put forward for regular and irregular situation.


#### 3.2.1. Definite Integral for Regular Domain

When the integral domainΩ is regular in [
a
,
b
]
×
[
c
,
d
], the two-dimensional Gauss quadrature formula is adopted with integral range transformation t
=
(
b
+
a
)
/
2
+
(
b
-
a
)
ξ
/
2
=
p
(
ξ
), s
=
(
d
+
c
)
/
2
+
(
d
-
c
)
η
/
2
=
q
(
η
). We have
(17)
K
m
n
=
∫
c
d
∫
a
b
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t
≈
(
b
-
a
)
(
d
-
c
)
4
∑
i
=
1
n
1
∑
j
=
1
n
2
W
i
W
j
k
[
x
m
,
y
m
,
p
(
ξ
i
)
,
q
(
η
j
)
]
ϕ
n
1111111111111111111
×
[
p
(
ξ
i
)
,
q
(
η
j
)
]
.


#### 3.2.2. Definite Integral for Irregular Domain

When the integral regionΩ is irregular, we can employ the similar numerical integration scheme as
(18)
K
m
n
=
∫
Ω
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t
=
∑
q
W
q
K
(
x
m
,
y
m
,
s
q
,
t
q
)
ϕ
n
(
s
q
,
t
q
)
.
However, the traditional regional split technique was adopted in this paper. Firstly, the irregular region was split geometrically. Then, the entire domain definite integral can be obtained by adding up all the integral of subregions which is indicated as the following:
(19)
K
m
n
=
∫
Ω
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t
=
∑
i
∫
Ω
i
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t
≈
∑
i
I
i
,
where Ω
i represents the ith subregion. Specifically, the triangular element is chosen because of its strong adaptability in conforming two-dimensional irregular domain. The triangulation can be implemented automatically by the PDE toolbox in Matlab. Then, let the vertex coordinates denote by (
x
i
j
,
y
i
j
), j
=
1,2
,
3, so the element area S
i and the three midpoints (
x
i
c
j
,
y
i
c
j
) of element boundary were easily obtained. And, finally, the quadratic interpolation formula was applied to calculate the definite integral as
(20)
I
i
=
S
i
3
∑
j
=
1
3
k
(
x
m
,
y
m
,
x
i
c
j
,
y
i
c
j
)
ϕ
n
(
x
i
c
j
,
y
i
c
j
)
.


### 3.3. The Implementation Procedure of RBF for Solving Linear Integral Equation

In summary, the implementation procedure of RBF for solving linear integral equation was listed as the following.(1)
Setting the centersc of RBFs and shape parameters α, then setting the collocation points x
d, the numbers of RBFs and collocation points are, respectively, N and M in general. Note that, to simplify and speed up the computations, we can identify the set of centers with the set of collocation points.
(2)
Substituting approximate solution based on RBF into the linear integral equation, then the collocation conditions are imposed to form the coefficient matrix.
(3)
Matrix elements calculation: the Gauss quadrature formula is adopted for one-dimensional and two-dimensional regular domain problems, while for irregular domain case, the numerical integration scheme based on regional split technique was employed.
(4)
Coefficient calculation: the equation[
Φ
-
μ
Κ
]
[
λ
]
=
[
G
] was solved by SVD method and the [
λ
] was obtained.
(5)
RBF solution for unknown function then can be achieved as(21)
f
(
x
)
≈
f
^
(
x
)
=
[
Φ
(
x
)
]
[
λ
]
=
[
Φ
(
x
)
]
[
Φ
-
μ
K
]
-
1
[
G
]
.


### 3.4. Error Analysis

Utilizing the same way in error analysis as [13], let the operator K
:
R
1
→
R
1 for one-dimensional case be defined as
(22)
K
[
f
(
x
)
]
=
∫
a
b
∣
x
k
(
x
,
t
)
f
(
t
)
d
t
.
Now, we can rewrite (8) in abstract form as
(23)
(
I
-
μ
K
)
f
=
g
.
Let K be a bounded operator with μ
|
|
K
|
|
<
1, and then the operator I
-
μ
K is a contraction operator; by the Banach contraction mapping principle, (8) has a unique solution. In such a case, the geometric series theorem implies that the operator I
-
μ
K has a bounded inverse with
(24)
∥
(
I
-
μ
K
)
-
1
∥
≤
1
1
-
∥
K
∥
.So, the error of the presented RBF method is mainly based upon the RBF interpolation error which emerges in (3) and the error caused by performing the numerical integration scheme over domain which appears in (18). For sufficiently large nodes N, the error of the RBF interpolation is dominated over the integration error and, thence, increasing the number of nodes in the numerical integration method has no significant effect on the error.


## 3.1. RBF Interpolation to Approximate One-Dimensional Linear Integral Equation

The general form of one-dimensional linear integral equation is(8)
f
(
x
)
=
μ
∫
a
b
∣
x
k
(
x
,
t
)
f
(
t
)
d
t
+
g
(
x
)
,
x
∈
[
a
,
b
]
.
Here, b
∣
x represents the upper limit of integral. The constant b or variable x is, respectively, corresponding to Fredholm or Volterra equation. In the following discussion, we choose the upper limit b for demonstration. In RBF approximation, the unknown function f
(
x
) is expressed as a combination of RBF. Thus, the approximate formula (9) is obtained:
(9)
∑
i
=
1
N
λ
i
ϕ
i
(
x
)
=
μ
∫
a
b
k
(
x
,
t
)
∑
i
=
1
n
λ
i
ϕ
i
(
t
)
d
t
+
g
(
x
)
,
x
∈
[
a
,
b
]
.Furthermore, for any collocation pointx
d in [
a
,
b
], by using the linear property of integral, we obtain the specific expression of the above equation as
(10)
∑
i
=
1
N
λ
i
ϕ
i
(
x
d
)
-
μ
∑
i
=
1
N
λ
i
∫
a
b
k
(
x
d
,
t
)
ϕ
i
(
t
)
d
t
=
g
(
x
d
)
,
11111111111111111111111111111111
x
d
∈
[
a
,
b
]
,
when choosing M collocation points, hence forming the collocation equations in matrix form as
(11)
[
Φ
-
μ
Κ
]
[
λ
]
=
[
G
]
,
where [
G
]
=
[
g
(
x
1
)
,
g
(
x
1
)
,
…
,
g
(
x
M
)
]
T, and [
λ
]
=
[
λ
1
,
λ
2
,
…
,
λ
N
]
T. The elements in matrix [
Ф
] and [
K
] are defined:
(12)
Φ
m
n
=
(
x
m
-
c
n
)
2
+
α
2
K
m
n
=
∫
a
b
k
(
x
m
,
t
)
(
t
-
c
n
)
2
+
α
2
d
t
=
∫
a
b
h
(
t
)
d
t
.
The notations x
m and c
n represent correspondingly the collocation point and RBF center, while the element K
m
n is the definite integral of function h
(
t
) which can be calculated through the following Gauss quadrature formula:
(13)
∫
-
1
1
h
(
ξ
)
d
ξ
≈
∑
q
=
1
Q
W
q
h
(
ξ
q
)
.
The integral range can be transformed with formula t
=
(
b
+
a
)
/
2
+
(
b
-
a
)
ξ
/
2
=
p
(
ξ
), so the final calculation formula is
(14)
K
m
n
=
∫
a
b
k
(
x
m
,
t
)
ϕ
n
(
t
)
d
t
=
b
-
a
2
∫
-
1
1
k
[
x
m
,
p
(
ξ
)
]
ϕ
n
[
p
(
ξ
)
]
d
ξ
≈
b
-
a
2
∑
q
=
1
Q
W
q
k
[
x
m
,
p
(
ξ
q
)
]
ϕ
n
[
p
(
ξ
q
)
]
.


## 3.2. RBF Interpolation to Approximate Two-Dimensional Linear Integral Equation

Two-dimensional linear integral equation has the general form(15)
f
(
x
,
y
)
=
μ
∫
Ω
k
(
x
,
y
,
s
,
t
)
f
(
s
,
t
)
d
s
d
t
+
g
(
x
,
y
)
,
hhhhhhhhhhhhhhhhhhhhhhhhhhhh
(
x
,
y
)
∈
Ω
.
The corresponding RBF approximation is
(16)
∑
i
=
1
N
λ
i
ϕ
i
(
x
d
,
y
d
)
=
μ
∑
i
=
1
N
λ
i
∫
Ω
k
(
x
d
,
y
d
,
s
,
t
)
ϕ
i
(
s
,
t
)
d
s
d
t
+
g
(
x
d
,
y
d
)
.And, similarly, we have the interpolation matrix[
Φ
-
μ
Κ
]
[
λ
]
=
[
G
]. The differences lie in the distance formula and the definite integral of element K
m
n
=
∫
Ω
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t. According to the domain type, two integration schemes were put forward for regular and irregular situation.


### 3.2.1. Definite Integral for Regular Domain

When the integral domainΩ is regular in [
a
,
b
]
×
[
c
,
d
], the two-dimensional Gauss quadrature formula is adopted with integral range transformation t
=
(
b
+
a
)
/
2
+
(
b
-
a
)
ξ
/
2
=
p
(
ξ
), s
=
(
d
+
c
)
/
2
+
(
d
-
c
)
η
/
2
=
q
(
η
). We have
(17)
K
m
n
=
∫
c
d
∫
a
b
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t
≈
(
b
-
a
)
(
d
-
c
)
4
∑
i
=
1
n
1
∑
j
=
1
n
2
W
i
W
j
k
[
x
m
,
y
m
,
p
(
ξ
i
)
,
q
(
η
j
)
]
ϕ
n
1111111111111111111
×
[
p
(
ξ
i
)
,
q
(
η
j
)
]
.


### 3.2.2. Definite Integral for Irregular Domain

When the integral regionΩ is irregular, we can employ the similar numerical integration scheme as
(18)
K
m
n
=
∫
Ω
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t
=
∑
q
W
q
K
(
x
m
,
y
m
,
s
q
,
t
q
)
ϕ
n
(
s
q
,
t
q
)
.
However, the traditional regional split technique was adopted in this paper. Firstly, the irregular region was split geometrically. Then, the entire domain definite integral can be obtained by adding up all the integral of subregions which is indicated as the following:
(19)
K
m
n
=
∫
Ω
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t
=
∑
i
∫
Ω
i
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t
≈
∑
i
I
i
,
where Ω
i represents the ith subregion. Specifically, the triangular element is chosen because of its strong adaptability in conforming two-dimensional irregular domain. The triangulation can be implemented automatically by the PDE toolbox in Matlab. Then, let the vertex coordinates denote by (
x
i
j
,
y
i
j
), j
=
1,2
,
3, so the element area S
i and the three midpoints (
x
i
c
j
,
y
i
c
j
) of element boundary were easily obtained. And, finally, the quadratic interpolation formula was applied to calculate the definite integral as
(20)
I
i
=
S
i
3
∑
j
=
1
3
k
(
x
m
,
y
m
,
x
i
c
j
,
y
i
c
j
)
ϕ
n
(
x
i
c
j
,
y
i
c
j
)
.


## 3.2.1. Definite Integral for Regular Domain

When the integral domainΩ is regular in [
a
,
b
]
×
[
c
,
d
], the two-dimensional Gauss quadrature formula is adopted with integral range transformation t
=
(
b
+
a
)
/
2
+
(
b
-
a
)
ξ
/
2
=
p
(
ξ
), s
=
(
d
+
c
)
/
2
+
(
d
-
c
)
η
/
2
=
q
(
η
). We have
(17)
K
m
n
=
∫
c
d
∫
a
b
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t
≈
(
b
-
a
)
(
d
-
c
)
4
∑
i
=
1
n
1
∑
j
=
1
n
2
W
i
W
j
k
[
x
m
,
y
m
,
p
(
ξ
i
)
,
q
(
η
j
)
]
ϕ
n
1111111111111111111
×
[
p
(
ξ
i
)
,
q
(
η
j
)
]
.


## 3.2.2. Definite Integral for Irregular Domain

When the integral regionΩ is irregular, we can employ the similar numerical integration scheme as
(18)
K
m
n
=
∫
Ω
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t
=
∑
q
W
q
K
(
x
m
,
y
m
,
s
q
,
t
q
)
ϕ
n
(
s
q
,
t
q
)
.
However, the traditional regional split technique was adopted in this paper. Firstly, the irregular region was split geometrically. Then, the entire domain definite integral can be obtained by adding up all the integral of subregions which is indicated as the following:
(19)
K
m
n
=
∫
Ω
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t
=
∑
i
∫
Ω
i
k
(
x
m
,
y
m
,
s
,
t
)
ϕ
n
(
s
,
t
)
d
s
d
t
≈
∑
i
I
i
,
where Ω
i represents the ith subregion. Specifically, the triangular element is chosen because of its strong adaptability in conforming two-dimensional irregular domain. The triangulation can be implemented automatically by the PDE toolbox in Matlab. Then, let the vertex coordinates denote by (
x
i
j
,
y
i
j
), j
=
1,2
,
3, so the element area S
i and the three midpoints (
x
i
c
j
,
y
i
c
j
) of element boundary were easily obtained. And, finally, the quadratic interpolation formula was applied to calculate the definite integral as
(20)
I
i
=
S
i
3
∑
j
=
1
3
k
(
x
m
,
y
m
,
x
i
c
j
,
y
i
c
j
)
ϕ
n
(
x
i
c
j
,
y
i
c
j
)
.


## 3.3. The Implementation Procedure of RBF for Solving Linear Integral Equation

In summary, the implementation procedure of RBF for solving linear integral equation was listed as the following.(1)
Setting the centersc of RBFs and shape parameters α, then setting the collocation points x
d, the numbers of RBFs and collocation points are, respectively, N and M in general. Note that, to simplify and speed up the computations, we can identify the set of centers with the set of collocation points.
(2)
Substituting approximate solution based on RBF into the linear integral equation, then the collocation conditions are imposed to form the coefficient matrix.
(3)
Matrix elements calculation: the Gauss quadrature formula is adopted for one-dimensional and two-dimensional regular domain problems, while for irregular domain case, the numerical integration scheme based on regional split technique was employed.
(4)
Coefficient calculation: the equation[
Φ
-
μ
Κ
]
[
λ
]
=
[
G
] was solved by SVD method and the [
λ
] was obtained.
(5)
RBF solution for unknown function then can be achieved as(21)
f
(
x
)
≈
f
^
(
x
)
=
[
Φ
(
x
)
]
[
λ
]
=
[
Φ
(
x
)
]
[
Φ
-
μ
K
]
-
1
[
G
]
.


## 3.4. Error Analysis

Utilizing the same way in error analysis as [13], let the operator K
:
R
1
→
R
1 for one-dimensional case be defined as
(22)
K
[
f
(
x
)
]
=
∫
a
b
∣
x
k
(
x
,
t
)
f
(
t
)
d
t
.
Now, we can rewrite (8) in abstract form as
(23)
(
I
-
μ
K
)
f
=
g
.
Let K be a bounded operator with μ
|
|
K
|
|
<
1, and then the operator I
-
μ
K is a contraction operator; by the Banach contraction mapping principle, (8) has a unique solution. In such a case, the geometric series theorem implies that the operator I
-
μ
K has a bounded inverse with
(24)
∥
(
I
-
μ
K
)
-
1
∥
≤
1
1
-
∥
K
∥
.So, the error of the presented RBF method is mainly based upon the RBF interpolation error which emerges in (3) and the error caused by performing the numerical integration scheme over domain which appears in (18). For sufficiently large nodes N, the error of the RBF interpolation is dominated over the integration error and, thence, increasing the number of nodes in the numerical integration method has no significant effect on the error.


## 4. Numerical Experiments and Analysis

In this section, the RBF method is applied to some numerical examples involving one-dimensional and two-dimensional integral equations on the regular and nonrectangular regions. Two types of RBF as multiquadric (MQ) function and thin plate splines (TPS) function were adopted simultaneously in the following simulations. Their parameters were set as the same except that the additional shape parameter was provided for MQ. In order to measure the accuracy of the method, the maximum absolute error (MAE), the maximum relative error (MRE), and the root mean squared error (RMSE) have been used, respectively. The RMSE has the following definition as(25)
RMSE
=
1
m
N
∑
k
=
1
m
N
[
f
e
(
x
k
)
-
f
(
x
k
)
]
2
,
where f
e
(
x
k
), f
(
x
k
) denote the exact and approximation value at point x
k and symbol m
N is the total testing points.Example 1.
Consider the following Fredholm integral equation in [1]:
(26)
∫
0
1
x
+
t
f
(
t
)
d
t
=
g
(
x
)
,
g
(
x
)
=
-
16
105
x
7
/
2
+
16
105
x
2
(
x
+
1
)
3
/
2
111111
-
8
35
x
(
x
+
1
)
3
/
2
+
2
7
(
x
+
1
)
3
/
2
,
111111111111111111111111
0
≤
x
≤
1
.
The exact solution for this equation is f
(
x
)
=
x
2.Setting RBF centers with intervalh
=
0.2, the collocation points are the same as the center points N
=
6. For MQ, the shape parameter β
=
12. Then, the 20-point Gauss quadrature formula is adopted in matrix element calculation. And, finally, the numerical result was shown in Table 1. Comparing to [1] where the Haar wavelet was adopted as base function, for the scaling j
=
4, the RMSE, MAE, and MRE were 4.1774
E
-
2, 1.007
E
-
1, and 380.46
%, and, for j
=
6, the corresponding errors were 2.0457
E
-
2, 5.1024
E
-
2, and 187.15%. And in TPS method, the errors were 1.3403
e
-
1, 2.2625
E
-
1, and 159.79%. But, for the MQ method, the accuracy was significantly improved as 5.7774
E
-
5, 1.0512
E
-
4, and 0.4112%, respectively, whose accuracy was improved nearly 3 orders of magnitude.Table 1
Comparison between RBF methods and Haar wavelet methods.
x
t
Exact
Haar wavelet withj
=
4
Haar wavelet withj
=
6
Thin plate splines
Multiquadric
0.1
0.01
0.04804645706772
0.02871498616505
0.01369148850547
0.01004112468149
0.2
0.04
0.02853148488475
0.04216759916577
−0.02391498932119
0.04003606824426
0.3
0.09
0.05136013092094
0.07818735894931
0.14584667811187
0.09004712350520
0.4
0.16
0.15024502143722
0.14464852362516
0.28759641578216
0.16005432363045
0.5
0.25
0.21756512590070
0.23821128159978
0.18330827120676
0.25004112179590
0.6
0.36
0.37282520080419
0.37175451578151
0.13743517284623
0.36001852111326
0.7
0.49
0.54342797321883
0.50029535575412
0.50098212437262
0.49001930600502
0.8
0.64
0.63155345618083
0.66060990603746
0.86624590306761
0.64006180352681
0.9
0.81
0.80994050257895
0.80342706047881
0.83621252785184
0.81008909789603
1.0
1.00
0.89929178048446
0.94897579805255
0.77598833938371
0.99989488295933
RMSE
4
.
1775
E
-
002
2
.
045
7
E
-
002
1
.
3
403
E
-
001
5.7774
E
-
005Then, the shape parameter was analyzed with discrete interval of 0.1; the RMSE and the condition number of matrix[
Φ
-
μ
Κ
] were calculated in case of h
=
0.2. It is obvious to see in Figure 1 that the optimal shape parameter lies in the interval of [
10
,
20
] showed in red dashed box. And, in the whole range of [0.1,100], the MQ method can maintain the RMSE at the level of 1
E
-
3. The curve of condition number was shown in Figure 2, and it fluctuates in the range of [1
E
+
16, 1
E
+
18] when β is greater than 20. So the MQ method is stable when the shape parameter varies in a relative wide range.Figure 1
The RMSE curve versus the shape parameterβ.Figure 2
The condition number curve versus the shape parameterβ.And, finally, the RBF’s center spacing was analyzed in case ofβ
=
12. The RMSE curve was shown in Figure 3. It can be seen clearly that the optimal center spacing was h
=
0.2 and 0.0833. But for mostly center spacing, the MQ method can maintain RMSE at the level of 1
E
-
3. As for condition number, whose curve also fluctuates in the range of [1
E
+
13, 1
E
+
20], it was shown in Figure 4.Figure 3
The RMSE curve versus the RBF’s center spacingh.Figure 4
The condition number curve versus the RBF’s center spacingh.Example 2.
Consider the one-dimensional Volterra integral equation in [4]:
(27)
∫
0
x
cos
⁡
(
x
-
t
)
f
(
t
)
d
t
=
x
sin
(
x
)
,
0
≤
x
≤
1
.
The exact solution is f
(
x
)
=
2
sin
⁡
(
x
).Firstly, setting the RBF center spacingh
=
0.1, the collocation points and centers are the same as N
=
11. The shape parameter β
=
10 for MQ and the 20-point Gauss quadrature formula is adopted. The total test points are 101 with equidistant h
t
=
0.01. We focus on the center spacing variation; the approximation solutions were shown in Figure 5 for TPS method. Generally, the TPS can acquire better approximation as the center spacing decreases which means that the solution of h
=
0.05 was better than that of h
=
0.1. However, numerical simulation also indicated that keeping on decreasing h may result in unstable solution.Figure 5
The TPS solutions for different RBF center spacingh.The comparison between MQ and TPS was shown in Table2.Table 2
Error comparison between MQ and TPS.
h
MQ
TPS
RMSE
MAE
MRE
Cond.
RMSE
MAE
MRE
Cond.
0.25
1.3784
e
-
4
1.3784
e
-
4
6.2050
e
-
4
3.6563
e
+
0
8
2.3237
8.0507
2
.
800
4
e
1
2.8490
e
+
03
0.20
1.6127
e
-
5
9.9905
e
-
5
5.9363
e
-
5
5.4846
e
+
0
9
1
.
1164
e
-
1
1
.
8176
e
-
1
2.0428
4
.
7439
e
+
02
0.10
2.1030
e
-
5
2.4501
e
-
5
1
.
225
1
e
-
3
1.8064
e
+
13
6
.
1529
e
-
2
2
.
2521
e
-
1
1.5646
2.9389
e
+
04
0.05
9.6970
e
-
5
1.8010
e
-
4
8
.
903
5
e
-
3
1.4434
e
+
16
1.9299
e
-
2
8
.
386
8
e
-
2
6.2500
e
-
1
8.2771
e
+
05
0.025
3.4048
e
-
4
4.9369
e
-
4
2.4357
e
-
2
3.1722
e
+
17
1
.
906
5
e
-
1
1.1247
6
.
7714
e
-
01
5.7935
e
+
09Therefore, we can conclude from Table2 that the MQ method takes obvious advantages over TPS whose accuracy is improved about 3 or 4 orders of magnitude. And the MQ method is also much stable for wide range of center spacing h. The MQ especially can acquire high precision even with larger center spacing as h
=
0.25. Compared to the method in [4], the MQ also showed its superiority; the highest precision of modified block pulse functions method is about the level of 1
E
-
2.Example 3.
Consider the two-dimensional integral equation in [14]:
(28)
f
(
x
,
y
)
=
u
(
x
,
y
)
+
∫
0
1
(
t
sin
(
s
)
+
1
)
f
(
s
,
t
)
d
s
d
t
,
u
(
x
,
y
)
=
x
cos
⁡
(
y
)
-
sin
⁡
(
1
)
(
3
+
sin
⁡
(
1
)
)
6
,
hhhhhhhhhhhhhhhhhhhh
0
≤
x
,
y
≤
1
.
The exact solution is f
(
x
,
y
)
=
x
cos
⁡
(
y
).Setting the RBF center spacingh
=
0.1, the collocation points and centers are N
=
121. The shape parameter is β
=
35 for MQ and the 20-point Gauss quadrature formula was also adopted in both x and y directions. The total testing points are 441 with equidistant h
t
=
0.05. The RMSE, MAE, MRE, and condition number for MQ are 2.6583
E
-
6, 5.5309
E
-
6, 1.6077
E
-
4, and 8.0241
E
+
18, respectively. However, the above values for TPS are 1.2273
E
-
3, 9.9312
E
-
3, 3.4146
E
-
1, and 1.9795
E
+
3. And comparing to the rational Haar function method in [14] whose accuracy was about level 1
E
-
2 and 1
E
-
3, the MQ acquires much higher accuracy than others. The error surfaces of MQ and TPS were shown in Figure 6.The error surfaces for MQ method and TPS method.
(a)
MQ method
(b)
TPS methodThen, analyzing the shape parameter of MQ with discrete step of 1.0, the RMSE and the condition number were calculated in case ofh
=
0.1. The local optimal shape parameters in this instance were β
=
7, 35, 58, and 64 revealed in Figure 7. The overall RMSE of MQ method was about the level of 1
E
-
4 and 1
E
-
5. And, correspondingly, the curve of condition number was similar to that in Figure 2 except that the curve fluctuates in the range of [
1
E
+
18
,
1
E
+
20
] when β is greater than 10.Figure 7
The RMSE curve versus the shape parameterβ.Example 4.
Consider the two-dimensional integral equation in an irregular domain in [13]:
(29)
f
(
x
,
y
)
-
∫
Ω
(
t
sin
(
s
)
+
1
)
f
(
s
,
t
)
d
s
d
t
=
x
cos
⁡
(
y
)
-
0.2804057157
111111
Ω
=
{
1
-
4
(
y
-
0.5
)
2
(
x
,
y
)
∈
R
2
:
0
≤
y
≤
1
,
1111111111111111111111
0.25
-
(
y
-
0.5
)
2
1111111111111111111111
≤
x
≤
1
-
4
(
y
-
0.5
)
2
}
.
The exact solution is f
(
x
,
y
)
=
x
cos
⁡
(
y
).Firstly, setting the RBF centers and collocation points as Figure8, the total number was N
=
35. The shape parameter was α
=
β
h
=
2.3 for MQ and the regional splitting with triangular.Figure 8
RBF centers and collocation points.Then, the PDE toolbox in Matlab was utilized in triangulation. The above irregular regionΩ was split into 2272 triangular elements with 1225 nodes. In simplicity, the testing points were chosen at the element nodes. And the quadratic form of interpolation was employed in calculating the definite integral of subregions in (20). The results showed that the RMSE, MAE, MRE, and condition number are 2.6260
E
-
5, 3.0345
E
-
5, 1.5783
E
-
3, and 3.3527
E
+
14. However, the above values for TPS method are 3.9563
E
-
3, 1.9925
E
-
2, 8.6492
E
-
1, and 4.3509
E
+
2. The error curves of MQ and TPS are shown in Figure 9. Comparing to [13] where 33 collocation points were selected, the RMSE is 1.96
E
-
4, and for GA weight function the RMSE is 1.65
E
-
4. So the RBF method proposed in this paper can obtain higher accuracy.The error curves for MQ method and TPS method.
(a)
MQ method
(b)
TPS method


## 5. Conclusion

In this paper, the RBF method, specifically, the MQ, was introduced in solving linear integral equations. The following conclusions can be drawn.(1)
The procedure of MQ method includes that the unknown function was firstly expressed in linear combination forms of RBFs, then the integral equation was transformed into collocation matrix of RBFs, and finally, solving the matrix equation and an approximation solution was obtained.
(2)
Due to the excellent interpolation performance of MQ, the method can acquire higher precision with fewer nodes and low computations. Numerical results showed that the MQ method can achieve solution accuracy of1
E
-
5. Numerical results showed that the MQ method takes obvious advantages over TPS. In addition, the method was very convenient for solving higher dimensional integral equations because the only modification in implementation is the distance formula.
(3)
One of the key issues in MQ method is the definite integral calculation. Two types of integration schemes as the Gauss quadrature formula and regional split technique were put forward. Specifically, the former is adopted in one-dimensional and two-dimensional regular domain problems and the latter is employed in irregular domain case.



---

*Source: 793582-2014-06-01.xml*