# The Approximate Solution of Fredholm Integral Equations with Oscillatory Trigonometric Kernels

Authors: Qinghua Wu
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

A method for approximating the solution of weakly singular Fredholm integral equation of the second kind with highly oscillatory trigonometric kernel is presented. The unknown function is approximated by expansion of Chebychev polynomial and the coefficients are determinated by classical collocation method. Due to the highly oscillatory kernels of integral equation, the discretised collocation equation will give rise to the computation of oscillatory integrals. These integrals are calculated by using recursion formula derived from the fundamental recurrence relation of Chebyshev polynomial. The effectiveness and accuracy of the proposed method are tested by numerical examples.


---

## Body

## 1. Introduction

The Fredholm integral equations of second kind,(1)
y
(
x
)
+
λ
∫
a
b
‍
K
(
x
,
t
)
|
x
-
t
|
α
y
(
t
)
d
t
=
f
(
x
)
,
y
(
x
)
+
λ
y
y
x
∈
[
a
,
b
]
,
0
<
α
<
1
,
where K
(
x
,
t
) is a continuous function and f
(
x
) is a given function, have many applications in mathematical physics and engineering, such as heat conduction problem, potential problems, quantum mechanics, and seismology image processing [1–4]. Particularly, when K
(
x
,
x
)
≠
0 and 0
<
α
<
1, (1) is called weakly singular.In most of the cases, the integral equation cannot be done analytically and one has to resort to numerical methods. Many numerical methods, such as collocation method and Galerkin method, have been developed to solve (1); for details see [2, 5, 6]. These methods are well-established numerical algorithms; however, standard version of these classical methods may suffer from difficulty for computation of (1), containing highly oscillatory kernels since the computation of the highly oscillatory integrals by standard quadrature methods is exceedingly difficult and the cost steeply increases with the frequency. Furthermore, Galerkin method requires many double integrals when approximating the solution of integral equation. Specially, when the kernel is highly oscillatory, it requires the evaluation of many highly oscillatory double integrals, which can become computationally expensive.Recently, for weakly singular Volterra integral equations of the second kind with highly oscillatory Bessel kernels, it was found that the collocation methods are much more easily implemented and can get higher accuracy than discontinuous Galerkin methods under the same piecewise polynomials space; for details see [7–10]. In addition, collocation method only involves single integrals which are a little easier to evaluate. Motivated by this fact, here we investigate the application of collocation method for the solution of Fredholm integral equation
(2)
y
(
x
)
+
λ
∫
a
b
‍
e
i
ω
(
x
-
t
)
|
x
-
t
|
α
y
(
t
)
d
t
=
f
(
x
)
,
(
x
)
+
i
λ
x
∈
[
a
,
b
]
,
α
<
1
,
ω
≫
1
.
Here y
(
x
) is the unknown function and f
(
x
) is a given function and assume that λ is not the eigenvalue of integral equation, since for any finite interval [
a
,
b
] it can be transformed to interval [
-
1,1
] by linear transformation. In this paper, we only consider the case of [
a
,
b
]
=
[
-
1,1
].On the other hand, when the unknown function is approximated by Chebychev polynomial, the discretization of the integral equation (2) by collocation method will give rise to highly oscillatory integral which can be formulated as
(3)
∫
a
b
‍
e
(
i
ω
(
x
-
t
)
)
|
x
-
t
|
α
T
j
(
t
)
d
x
,
x
∈
[
a
,
b
]
,
0
<
α
<
1
,
where T
j
(
x
) denotes the Chebychev polynomials of the first kind.In last few years many efficient methods have been devised for the evaluation of oscillatory integral, such as asymptotic method [11], Filon-type method [12, 13], Levin’s collocation method [14], modified Clenshaw-Curtis method, Clenshaw-Curtis-Filon-type method [15], and generalized quadrature rule [16], although some of these methods do not involve Chebyshev polynomial. Piessens and Poleunis [17] consider a simpler case of α
=
0 and use Chebyshev polynomials to evaluate the integral by somewhat indirect method involving a truncated infinite series of Bessel functions. The evaluation of this infinite summation suggests a defect of this method. An alternative procedure which avoids this infinite series is the original Bakhvalo and Vasil'eva-Legendre work involved in the evaluation of the integral
(4)
M
i
(
ω
)
=
∫
-
1
1
‍
P
i
(
x
)
e
i
ω
x
d
x
,
by recurrence relation, where P
i
(
x
) denotes Legendre polynomial. Sadly, this relation proved to be unstable in the forwards direction for small ω. In [18], Alaylioglu et al proposed a simple alternative approach analogous to Newton-Cotes based formulae. This method avoids the instability to a certain degree compared with these two methods mentioned above. We used the same idea and generalized it to the case of weakly singular and calculated oscillatory integral (3).This paper is organized as follows: in Section2 we derive some basic formulae and introduce some mathematical preliminaries of the proposed method. In Section 3 we discuss the evaluation of the integrals occurring in collocation equation. In Section 4 numerical experiments are conducted to illustrate the performance of the proposed method.


## 2. Fundamental Relations

Firstly, by separation of real and imaginary part ofy
(
x
) and e
i
w
(
x
-
t
), we transform the integral equation (2) into equivalent systems of two linear integral equations of Fredholm in the forms
(5)
y
1
(
x
)
+
λ
∫
-
1
1
‍
cos
⁡
⁡
(
ω
(
x
-
t
)
)
|
x
-
t
|
α
y
1
(
t
)
d
t
-
λ
∫
-
1
1
‍
sin
⁡
(
ω
(
x
-
t
)
)
|
x
-
t
|
α
y
2
(
t
)
d
t
=
f
1
(
x
)
,
y
2
(
x
)
+
λ
∫
-
1
1
‍
sin
⁡
(
ω
(
x
-
t
)
)
|
x
-
t
|
α
y
1
(
t
)
d
t
+
λ
∫
-
1
1
‍
cos
⁡
⁡
(
ω
(
x
-
t
)
)
|
x
-
t
|
α
y
2
(
t
)
d
t
=
f
2
(
x
)
,
where y
(
x
)
=
y
1
(
x
)
+
i
y
2
(
x
), f
(
x
)
=
f
1
(
x
)
+
i
f
2
(
x
), and i
=
-
1.Assume thaty
(
x
)
=
(
y
1
(
x
)
,
y
2
(
x
)
)
T and f
(
x
)
=
(
f
1
(
x
)
,
f
2
(
x
)
)
T; then, systems of linear integral equations (5) can be written in the matrix form
(6)
I
y
(
x
)
+
λ
∫
-
1
1
‍
K
(
x
,
t
)
y
(
t
)
d
t
=
f
(
x
)
,
where
(7)
I
=
(
1
0
0
1
)
,
K
(
x
,
t
)
=
(
cos
⁡
⁡
(
ω
(
x
-
t
)
)
|
x
-
t
|
α
-
sin
⁡
(
ω
(
x
-
t
)
)
|
x
-
t
|
α
sin
⁡
(
ω
(
x
-
t
)
)
|
x
-
t
|
α
cos
⁡
⁡
(
ω
(
x
-
t
)
)
|
x
-
t
|
α
)
.
Set y
i
(
x
)
=
T
(
x
)
b
i, i
=
1,2, where T
(
x
)
=
[
T
0
(
x
)
,
T
1
(
x
)
,
…
,
T
N
(
x
)
], b
i
=
[
b
i
0
,
b
i
1
,
…
,
b
i
N
]
T.Hence, unknown functions can be expressed by(8)
y
(
x
)
=
𝒯
B
,
where
(9)
y
(
x
)
=
(
y
1
(
x
)
y
2
(
x
)
)
,
𝒯
(
x
)
=
(
T
(
x
)
0
0
T
(
x
)
)
,
B
=
(
b
1
b
2
)
.Then, the aim is to find Chebyshev coefficients, that is, the matrixB. We first substitute the Chebyshev collocation points, which are defined by x
i
=
cos
⁡
⁡
(
i
π
/
N
),  i
=
0
,
…
,
N, into (6) and then rearrange a new matrix form to determine B:
(10)
ℐ
𝒴
+
λ
𝒦
=
ℱ
in which 𝒦 is the integral part of (6) and
(11)
ℐ
=
(
I
0
⋯
0
0
I
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
I
)
,
𝒴
=
(
y
(
x
0
)
y
(
x
1
)
⋮
y
(
x
N
)
)
,
ℱ
=
(
f
(
x
0
)
f
(
x
1
)
⋮
f
(
x
N
)
)
,
𝒦
=
(
I
K
(
x
0
)
I
K
(
x
1
)
⋮
I
K
(
x
N
)
)
.By substituting (8) into (10), the unknown coefficients can be easily computed from this linear algebraic equations and therefore we find the solution of integral equation (2).


## 3. Evaluation of the IntegralI
[
ω
,
α
,
j
,
x
]
=
∫
-
1
1
‍
(
e
(
i
ω
(
x
-
t
)
)
/
|
x
-
t
|
α
)
T
j
(
t
)
d
t

The discretization of integral equation will lead to the calculation of integralI
[
ω
,
α
,
j
,
x
]. In this section, we will derive the recursion formula to compute it efficiently from the fundamental recurrence relation of Chebyshev polynomial.The Chebyshev polynomials are of the form(12)
T
0
(
x
)
=
1
T
1
(
x
)
=
x
⋮
T
n
(
x
)
=
2
x
T
n
-
1
(
x
)
-
T
n
-
2
(
x
)
,
n
=
2,3
,
…
,
and the coefficients C
i
,
j of x
j in T
i
(
x
) can be easily calculated by the means of the recurrence relation
(13)
C
i
,
j
=
2
C
i
-
1
,
j
-
1
-
C
i
-
2
,
j
,
i
≥
2
,
j
≤
i
.
By expanding the Chebyshev polynomial T
i
(
x
) in terms of powers of x, the integral I
[
ω
,
α
,
j
,
x
] would be transformed into the form of
(14)
I
[
ω
,
α
,
j
,
x
]
=
∑
k
=
0
j
‍
C
j
,
k
I
1
[
ω
,
α
,
k
,
x
]
,
where I
1
[
ω
,
α
,
k
,
x
]
=
∫
-
1
1
‍
(
e
(
i
ω
(
x
-
t
)
)
/
|
x
-
t
|
α
)
t
k
d
t.By separation of real and imaginary part ofe
i
w
(
x
-
t
), (14) is transformed into
(15)
I
1
c
[
ω
,
α
,
j
,
x
]
=
∫
-
1
1
‍
t
j
|
x
-
t
|
α
cos
⁡
⁡
(
ω
(
x
-
t
)
)
d
t
,
I
1
s
[
ω
,
α
,
j
,
x
]
=
∫
-
1
1
‍
t
j
|
x
-
t
|
α
sin
⁡
(
ω
(
x
-
t
)
)
d
t
,
whose main difficulty now is turning to the evaluation of the following basic integrals:
(16)
M
1
c
[
ω
,
α
,
j
,
x
]
=
∫
-
1
1
‍
t
j
|
x
-
t
|
α
cos
⁡
⁡
(
ω
t
)
d
t
,
M
1
s
[
ω
,
α
,
j
,
x
]
=
∫
-
1
1
‍
t
j
|
x
-
t
|
α
sin
⁡
(
ω
t
)
d
t
.With substitution for definite integral and binomial theorem, (16) is equivalent to
(17)
M
1
c
[
ω
,
α
,
j
,
x
]
=
∑
k
=
0
j
‍
C
j
k
x
j
-
k
[
(
-
1
)
k
∫
0
x
+
1
‍
t
k
-
α
cos
⁡
⁡
(
ω
(
x
-
t
)
)
d
t
i
+
∫
0
1
-
x
‍
t
k
-
α
cos
⁡
⁡
(
ω
(
x
-
t
)
)
d
t
]
,
where C
j
k is number of combination and M
1
s
[
ω
,
α
,
j
,
x
] can be formulated analogy.Efficient evaluation of (18) is based on accurate calculation of integrals in the formula which can be computed explicitly by the incomplete Gamma function Γ
(
z
,
α
) [19]; that is,
(18)
∫
t
μ
-
1
cos
⁡
⁡
(
ω
t
)
d
t
=
1
2
[
(
i
ω
)
-
μ
Γ
(
μ
,
i
ω
t
)
+
(
-
i
ω
)
-
μ
Γ
(
μ
,
-
i
ω
t
)
]
,
(19)
∫
t
μ
-
1
sin
⁡
(
ω
t
)
d
t
=
i
2
[
(
i
ω
)
-
μ
Γ
(
μ
,
i
ω
t
)
-
(
-
i
ω
)
-
μ
Γ
(
μ
,
-
i
ω
t
)
]
,
in which μ
>
0
,
x
>
0; for details see ([20], pp 215).Once the integralI
[
ω
,
α
,
j
,
x
] is obtained, by substituting into (10), the coefficient matrix is derived; then, we can compute the solution of integral by solving these linear algebraic equations.


## 4. Numerical Examples

In this section, we give some numerical examples to illustrate the performance of proposed method. In all the following examples,N
+
1 is the number of mesh points, y
^
i
N
(
x
) denotes the approximate solution, where N is the number of terms of the Chebyshev series, and y
(
x
) denotes the exact solution, respectively. All the computations have been performed by using Matlab R2012a on a 2.5 GHz PC with 2 GB of RAM.Example 1.
We considery
(
x
)
+
∫
-
1
1
(
e
i
ω
(
x
-
t
)
/
|
x
-
t
|
α
)
‍
y
(
t
)
d
t
=
f
(
x
) and set f
(
x
)
=
1
+
∫
-
1
1
(
e
i
ω
(
x
-
t
)
/
|
x
-
t
|
α
)
‍
d
t; then, the solution of which is
(20)
y
(
x
)
=
1
.
We plot the absolute error |
y
(
x
)
-
y
^
(
x
)
| by Matlab 2012a internal functionezplot. Specially, for N
=
4, with ω
=
1
0
3 and α
=
0.5, solve linear algebraic equations (10); it can be found that the approximate solutions are
(21)
y
^
1
(
x
)
=
0.9999999996
-
0.110
×
1
0
-
15
x
+
0.281
×
1
0
-
9
x
2
+
0.323
×
1
0
-
15
x
3
+
0.563
×
1
0
-
10
x
4
,
y
^
2
(
x
)
=
0.460
×
1
0
-
14
-
0.329
×
1
0
-
9
x
+
0.554
×
1
0
-
11
x
2
+
0.322
×
1
0
-
9
x
3
-
0.555
×
1
0
-
11
x
4
.It is easy to see from Figure1 that the proposed method is converging and we only need N
=
4 so we could achieve high accuracy.Example1: the absolute error versus the x-coordinate for N
=
4 (a), N
=
8 (b), and N
=
16 (c) with negative data ignored.
(a)
(b)
(c)Example 2.
We considery
(
x
)
+
∫
-
1
1
(
e
i
ω
(
x
-
t
)
/
|
x
-
t
|
α
)
‍
y
(
t
)
d
t
=
f
(
x
) and set f
(
x
)
=
e
x
+
∫
-
1
1
(
e
i
ω
(
x
-
t
)
/
|
x
-
t
|
α
)
‍
e
t
d
t; then, the solution of which is
(22)
y
(
x
)
=
e
x
.In this case, forN
=
4, with ω
=
1000 and α
=
1
/
3, the approximate solutions are
(23)
y
^
1
(
x
)
=
1.000000001
+
0.9956819990
x
+
0.4992866904
x
2
+
0.1795192648
x
3
+
0.04379395375
x
4
,
y
^
2
(
x
)
=
-
0.530
×
1
0
-
7
+
0.189
×
1
0
-
7
x
+
0.426
×
1
0
-
6
x
2
-
0.246
×
1
0
-
7
x
3
-
0.420
×
1
0
-
6
x
4
.It is obvious to see from Figure2 that the proposed method is efficient and could achieve high accuracy with little number of collocation points.Example2: the absolute error versus the x-coordinate for N
=
4 (a), N
=
8 (b), and N
=
16 (c) with negative data ignored.
(a)
(b)
(c)Example 3.
For a general case, we considery
(
x
)
+
∫
-
1
1
‍
(
e
i
ω
(
x
-
t
)
/
|
x
-
t
|
α
)
y
(
t
)
d
t
=
e
x and set α
=
0.5. Although the solution of this equation is unknown, we can verify the calculation precision of the method to a certain degree by comparing the approximate evaluation at mesh points.In this case, forN
=
4, with ω
=
1000 and α
=
0.5, the approximate solutions are
(24)
y
^
1
(
x
)
=
0.9242239091
+
0.8836218495
x
+
0.4292470859
x
2
+
0.2458272504
x
3
+
0.1276741404
x
4
,
y
^
2
(
x
)
=
-
0.001077278807
+
0.06120358877
x
+
0.05269959347
x
2
-
0.1180433250
x
3
-
0.09433340997
x
4
.It is also easy to see from Tables1, 2, 3, and 4 that the presented method is efficient and accurate, although the exact solution is unknown.Table 1
Approximations fory
(
x
)
+
∫
-
1
1
(
e
i
ω
(
x
-
t
)
/
|
x
-
t
|
0.5
)
y
(
t
)
d
t
=
e
x with ω
=
1
0
3.
x
cos
⁡
(
π
)
cos
⁡
(
3
π
/
4
)
cos
⁡
(
0
)
y
^
1
4
0.351696035495271
0.459037927481918
0.924223909091303
y
^
1
8
0.351696947581043
0.459033811940511
0.924224547405338
y
^
1
16
0.351741610883315
0.459036527929146
0.924226969620880
y
^
2
4
0.0141286409694892
0.000146310591227839
-
0.00107727880701281
y
^
2
8
0.0141193392298137
0.000145019893750479
-
0.00107593754108643
y
^
2
16
0.0140583682955762
0.000144178555797097
-
0.00106881635211410Table 2
Approximations fory
(
x
)
+
∫
-
1
1
(
e
i
ω
(
x
-
t
)
/
|
x
-
t
|
0.5
)
y
(
t
)
d
t
=
e
x with ω
=
1
0
3.
x
cos
⁡
(
π
/
4
)
1
y
^
1
4
1.88249404678421
2.61059423521207
y
^
1
8
1.88248275737555
2.61059120470589
y
^
1
16
1.88248361661370
2.61073031676337
y
^
2
4
0.00323202027752064
-
0.0995508315886191
y
^
2
8
0.00321634589900894
-
0.0994917961267987
y
^
2
16
0.00319500282444326
-
0.0993143176430846Table 3
Approximations fory
(
x
)
+
∫
-
1
1
(
e
i
ω
(
x
-
t
)
/
|
x
-
t
|
0.5
)
y
(
t
)
d
t
=
e
x with ω
=
1
0
6.
x
cos
⁡
(
π
)
cos
⁡
(
3
π
/
4
)
cos
⁡
(
0
)
y
^
1
4
0.367419631536920
0.491835233010361
0.997500711679274
y
^
1
8
0.367419631541538
0.491835233006225
0.997500711679218
y
^
1
16
0.367419632527633
0.491835233064998
0.997500711678921
y
^
2
4
0.000458466865252363
0.150404190
×
1
0
-
5
-
0.219327280
×
1
0
-
5
y
^
2
8
0.000458466860662437
0.150403636
×
1
0
-
5
-
0.219326739
×
1
0
-
5
y
^
2
16
0.000458463698666909
0.150390826
×
1
0
-
5
-
0.219326736
×
1
0
-
5Table 4
Approximations fory
(
x
)
+
∫
-
1
1
(
e
i
ω
(
x
-
t
)
/
|
x
-
t
|
0.5
)
y
(
t
)
d
t
=
e
x with ω
=
1
0
6.
x
cos
⁡
(
π
/
4
)
1
y
^
1
4
2.02304108762491
2.71487514654707
y
^
1
8
2.02304108760947
2.71487514654071
y
^
1
16
2.02304108780835
2.71487515460854
y
^
2
4
0.384042582
×
1
0
-
5
-
0.00339813438891857
y
^
2
8
0.384040621
×
1
0
-
5
-
0.00339813431478428
y
^
2
16
0.3840356665
×
1
0
-
5
-
0.00339813322237100


## 5. Conclusion

In this paper, we explore quadrature methods for weakly singular Fredholm integral equation of the second kind with oscillatory trigonometric kernels and present collocation methods with Chebyshev series for calculation of the solution. For integral equation with highly oscillatory kernels, the standard collocation methods with classical quadrature methods are not suitable for the numerical approximation of the solution of integral equation, since the computation of the highly oscillatory integrals by standard quadrature methods is exceedingly difficult and the cost steeply increases with the frequency. Based on the recursion formula derived in Section3, we compute the highly oscillatory integrals occurring in collocation equation, directly and efficiently. Numerical examples demonstrate the performance of algorithm.



---

*Source: 172327-2014-04-17.xml*