EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
Vol. 14, No. 3, 2021, 969-979
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Generalized iteration method for the solution of fourth
order BVP via Green′s function

Fatma Aydin Akgun1,∗, Zaur Rasulov2

1 Yildiz Technical University, Mathematical Engineering Department, Davutpasa Campus,
34210, Istanbul, Turkey

Abstract. The aim of this paper is to extend and generalize Picard-Green’s fixed point iteration
method for the solution of fourth-order Boundary Value Problems. Several numerical applications
to linear and nonlinear fourth-order Boundary Value Problems are discussed to illustrate the main
results.

2020 Mathematics Subject Classifications: 47J26, 41A25

Key Words and Phrases: Boundary value problem, fixed point iteration method, Green’s
function, integral operator, rate of convergence

1. Introduction

An iterative method is an important tool for solving linear and nonlinear Boundary
Value Problems (BVPs). It has been used in the research areas of mathematics and
several branches of science and other fields. In recent years, the study of generalization
for the solution of third and fourth-order BVPs has attracted many scientists’ interest.
For instance, in mathematical modeling applications of a two-dimensional channel with
porous walls, deformation of elastic beams, fluid dynamics. These problems have been
studied by several mathematicians(see [11], [1], [4], [7], [2], [3], [8], [5], [15], [14] and the
references therein. [11] converted the BVP to Initial Value Problem (IVP) and solved the
linear equation while [13] solved the problem by transforming the fourth-order BVP to
third-order BVP and applied by contracting mapping principle. Meanwhile, [6] applied
the contracting mapping principle to the fourth-order BVP. Both of them showed the
existence, uniqueness, and convergence of the proposed methods.

[12] introduced his method for IVP by showing the existence and uniqueness for ordi-
nary differential equations, whereas [10] introduced his method.

[9] considered a cantilever beam problem in equilibrium position as follows:

u
′′′′

(t) = f(t, u(t), u
′
(t), u

′′
(t), u

′′′
(t)) (1)

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v14i3.4020

Email addresses: fakgun@yildiz.edu.tr (F.A. Akgun), rasulovzaurr@gmail.com (Z. Rasulov)

http://www.ejpam.com 969 © 2021 EJPAM All rights reserved.



F.A.Akgun, Z. Rasulov / Eur. J. Pure Appl. Math, 14 (3) (2021), 969-979 970

subject to the boundary conditions:

u(0) = 0, u
′
(0) = 0, u

′′
(1) = 0, u

′′′
(1) = 0 (2)

where t ∈ (0, 1) and f : [0, 1] × R3
+ × R− → R+ is continuous. By assuming that

f(t, u, x, y, z) is superlinear on u, x, y, z and satisfying a Nagumo-type condition, the
existence and uniqueness of the problem Eq. (1) − (2) are provided in the same arti-
cle. However, Quanq in [6] showed that the theorems in [9] are not sufficient for solving
resembling problems and provided some examples.

In this paper, we extend and generalize Picard-Green’s known iteration method to
fourth-order BVP. We proved the existence and uniqueness of theorems satisfying neces-
sary conditions for convergence. Furthermore, we give some numerical examples, including
linear and nonlinear BVPs, to illustrate the main results. Finally, we confirm better ap-
proximation with minimum errors than existing solutions by MATLAB.

2. Green’s function and methodology

Consider the following fourth-order BVP,

L[u] = p(t)u
′′′′

(t) + q(t)u
′′′

(t) + r(t)u
′′
(t) + h(t)u

′
(t) + g(t)u(t) = f(t) (3)

with the boundary conditions

Ba[u] = α1u(a) + α2u
′
(a) + α3u

′′
(a) + α4u

′′′
(a) = α

Bb[u] = β1u(b) + β2u
′
u(b) + β3u

′′
(b) + β4u

′′′
(b) = β

Bc[u] = γ1u(c) + γ2u
′
(c) + γ3u

′′
(c) + γ4u

′′′
(c) = γ

Bd[u] = ω1u(c) + ω2u
′
(c) + ω3u

′′
(c) + ω4u

′′′
(c) = ω (4)

where t ∈ (a, b) , α, β, γ and ω are constants and either c = a or c = b and either d = a or
d = b. The existence and uniqueness results for the solution of the problems Eq. (3)-(4)
are given in [13],[6] and [3].

The Green’s function G(t, s) corresponding to linear term in Eq. (3)

G(t, s) =

{
a1u1 + b1u2 + c1u3 + d1u4, a < t < s

a2u1 + b2u2 + c2u3 + d2u4, s < t < b
(5)

where t 6= s , u1, u2, u3 and u4 are linearly independent solutions of L[u] and ai, bi, ci
and di, (i = 1, 2) are constants. We follow five conditions to find the constants and final
version of Green’s function of the fourth-order boundary value problem.

1. G satisfies the homogeneous boundary conditions;

Ba[G(t, s)] = Bb[G(t, s)] = Bc[G(t, s)] = Bd[G(t, s)] = 0, (6)



F.A.Akgun, Z. Rasulov / Eur. J. Pure Appl. Math, 14 (3) (2021), 969-979 971

2. G is continuous at t = s;

a1u1(s) + b1u2(s) + c1u3(s) + d1u4(s) = a2u1(s) + b2u2(s) + c2u3(s) + d2u4(s), (7)

3. G
′

is continuous at t = s;

a1u
′
1(s) + b1u

′
2(s) + c1u

′
3(s) + d1u

′
4(s) = a2u

′
1(s) + b2u

′
2(s) + c2u

′
3(s) + d2u

′
4(s), (8)

4. G
′′

is continuous at t = s;

a1u
′′
1(s) + b1u

′′
2(s) + c1u

′′
3(s) + d1u

′′
4(s) = a2u

′′
1(s) + b2u

′′
2(s) + c2u

′′
3(s) + d2u

′′
4(s). (9)

5. G
′′′

has jumping discontinuous at t = s;

a1u
′′′
1 (s)+b1u

′′′
2 (s)+c1u

′′′
3 (s)+d1u

′′′
4 (s)+1/p(s) = a2u

′′′
1 (s)+b2u

′′′
2 (s)+c2u

′′′
3 (s)+d2u

′′′
4 (s).
(10)

As a consequence of these calculations, the Green’s function can be written in the following
form

up =

∫ b

a
G(t, s)f(s)ds, (11)

where un is the particular solution of Eq. (3).

3. Picard-Green’s fixed-point iterative method (PGEM)

Let us define the following linear integral operator to implement the proposed Picard-
Green’s fixed point iteration method.

T [u] = uh +

∫ b

a
G(t, s)(p(s)u

′′′′
(s)+ q(s)u

′′′
(s)+ r(s)u

′′
(s)+h(s)u

′
(s)+g(s)u(s))ds, (12)

By using Eq. (12), we obtain

T [u] = uh +

∫ b

a
G(t, s)(p(s)u

′′′′
(s) + q(s)u

′′′
(s) + r(s)u

′′
(s) + h(s)u

′
(s) + g(s)u(s)− f(s))ds

+

∫ b

a
G(t, s)f(s)ds.(13)

Using Eq. (11), we get

T [u] = u+

∫ b

a
G(t, s)(p(s)u

′′′′
(s) + q(s)u

′′′
(s) + r(s)u

′′
(s) +h(s)u

′
(s) + g(s)u(s)− f(s))ds.

(14)
Let the starting function u0 be the homogeneous solution of L[u] = 0 and un+1 = T [un],
for all n ≥ 0, then Picard-Green’s fixed point iteration method for Eq. (3) is defined as:

un+1 = un+

∫ b

a
G(t, s)(p(s)u

′′′′
n (s)+q(s)u

′′′
n (s)+r(s)u

′′
n(s)+h(s)u

′
n(s)+g(s)un(s)−f(s))ds.

(15)



F.A.Akgun, Z. Rasulov / Eur. J. Pure Appl. Math, 14 (3) (2021), 969-979 972

4. Mann-Green’s fixed point iterative method (MGEM)

The operator for Mann-Green’s fixed point iteration method is un+1 = (1 − αn)un +
αnT [un], for n ≥ 0. We generalize the Mann-Green’s method by applying the solutions of
Eq. (12)-(14) as follows:

un+1 = (1− αn)un + αn

∫ b

a

G(t, s)(p(s)u
′′′′

n (s) + q(s)u
′′′

n (s) + r(s)u
′′

n(s) + h(s)u
′

n(s) + g(s)un(s)− f(s))ds

= un + αn

∫ b

a

G(t, s)(p(s)u
′′′′

n (s) + q(s)u
′′′

n (s) + r(s)u
′′

n(s) + h(s)u
′

n(s) + g(s)un(s)− f(s))ds. (16)

Here, the choice of starting point is similar to the PGEM.

5. Convergence analysis

In this section, we analyze the convergence and determine the convergence rate. The
convergence analysis will be provided by using nonlinear differential equations and by the
contraction principle.

Consider the fourth-order BVP

u
′′′′

(t) = f(t, u(t), u
′
(t), u

′′
(t), u

′′′
(t)), (17)

with the boundary conditions

u(0) = α, u
′
(0) = β, u

′′
(1) = γ, u

′′′
(1) = ω. (18)

In particular, the solution of the problem Eq. (17) − (18) is as follows:

up =

∫ 1

0
G(t, s)f(s, up, u

′
p, u

′′
p , u

′′′
p )ds. (19)

Note that Eq. (19) cannot be replaced in Eq. (13). Due to this, Eq. (12) should be
modified as

T [up] =

∫ b

a
G(t, s)(p(s)u

′′′′
p (s) + q(s)u

′′′
p (s) + r(s)u

′′
p(s) + h(s)u

′
p(s) + g(s)up(s))ds. (20)

From Eq. (19) and Eq. (20), we get

T [up] = up +

∫ b

a
G(t, s)(p(s)u

′′′′
p (s) + q(s)u

′′′
p (s) + r(s)u

′′
p(s) + h(s)u

′
p(s) (21)

+g(s)up(s)− f(s, up, u
′
p, u

′′
p , u

′′′
p ))ds.

Let up = u. The Green’s function of Eq. (1) − (2) is as follows:

G(t, s) =

{
−t3
6 + st2

2 , 0 < t < s,
−s3
6 + ts2

2 , s < t < 1,
(22)



F.A.Akgun, Z. Rasulov / Eur. J. Pure Appl. Math, 14 (3) (2021), 969-979 973

and the adjoint of Eq. (22) is as follows:

G∗(t, s) = −

{
−s3
6 + ts2

2 , 0 < s < t,
−t3
6 + st2

2 , t < s < 1.
(23)

By using Eq. (23) and Eq. (15), we have

un+1 = un +

∫ b

a
G∗(t, s)(u

′′′′
n (s)− f(s, un(s), u

′
n(s), u

′′
n(s), u

′′′
n (s))ds. (24)

In particular, we have

un+1 = un −
∫ t

0
(
−s3

6
+
ts2

2
)(u

′′′′
n (s)− f(s, un(s), u

′
n(s), u

′′
n(s), u

′′′
n (s))ds (25)

−
∫ 1

t
(
−t3

6
+
st2

2
)(u

′′′′
n (s)− f(s, un(s), u

′
n(s), u

′′
n(s), u

′′′
n (s))ds.

Theorem 1. Let f(t, u, u′, u′′, u′′′) be a continuous function with bounded derivative u and

K := 1/8Lc < 1,

where

Lc = max
[0,1]×R3

|∂(f)

∂(u)
|.

Then, the iterative sequence un(t)∞(n=1) given in Eq. (26) converges uniformly to the solu-

tion of Eq. (1) − (2) for any bounded function f(t, u, u′, u′′, u′′′) on [0,1].

Proof. Let ‖u‖ = max0≤ t≤1 |u(t)| in C[0, 1] . By integrating (26) by parts, we get

un+1 (t) = un (t) +G∗(t, 1)u
′′′
n (1)−G∗(t, 0)u

′′′
n (0)−G∗s(t, 1)u

′′
n(1) +G∗s(t, 0)u

′′
n(0)

+G∗ss(t, 1)u
′
n(1)−G∗ss(t, 0)u

′
n(0)−G∗sss(t, 1)un(1) +G∗sss(t, 0)un(0)

+

∫ 1

0
G∗ssss(t, s)un(s)ds−

∫ 1

0
G∗(t, s)f(s, un(s), u

′
n(s), u

′′
n(s), u

′′′
n (s))ds. (26)

By using (18) and the conditions of Green’s function in Eq. (26), we obtain

un+1 (t) =
t3ω

6
+
γ + ω

2
t2 + βt− α−

∫ 1

0
G∗(t, s)(f(s, un(s), u

′
n(s), u

′′
n(s), u

′′′
n (s))ds. (27)

Let un+1 = TG(u) where TG : C[0, 1]→ C[0, 1], then

TG(u) ≡ ω

6
t3 +

γ + ω

2
t2 + βt− α−

∫ 1

0
G∗(t, s)(f(s, un(s), u

′
n(s), u

′′
n(s), u

′′′
n (s))ds. (28)

Theorem 2. Let TG(u) be a contractive mapping, then the fixed-point iteration method is
convergent.



F.A.Akgun, Z. Rasulov / Eur. J. Pure Appl. Math, 14 (3) (2021), 969-979 974

From Eq. (28),

|TG(u)− TG(z)| = |
∫ 1

0
G∗(t, s)(f(s, u, u

′
, u
′′
, u
′′′

)− f(s, z, z
′
, z
′′
, z
′′′

))ds|, (29)

and by using the fact that

max
0≤t≤1

|
∫ 1

0
G∗(t, s)ds| = 1/8, (30)

we get

|TG(u)− TG(z)| ≤ 1

8

∫ 1

0
|f(s, u, u

′
, u
′′
, u
′′′

)− f(s, z, z
′
, z
′′
, z
′′′

)|ds. (31)

By applying Mean Value Theorem, we obtain the following inequality:

|TG(u)− TG(z)| ≤ 1

8
max
[0,1]
|f(t, u(t), u

′
(t), u

′′
(t), u

′′′
(t)) (32)

−f(t, z(t), z
′
(t), z

′′
(t), z

′′′
(t))| ≤ 1

8
Lc‖u− z‖,

where ‖u− z‖ = max0≤t≤1 |u(t)− z(t)| and Lc = max[0,1]×R3 |∂f(t,u,u
′,u′′,u′′′)

∂u |. From The-
orem 1, we get

‖TG(u)− TG(z)‖ ≤ K‖u− z‖, (33)

where K ∈ (0, 1). This proves that TG is a contraction mapping. On the other hand, the
rate of convergence can be estimated as follows. Consider

‖un+1 − un‖ = ‖TG(un)− TG(un−1)‖ ≤ K‖un − un−1‖ ≤ Kn‖u1 − u0‖. (34)

If m > n > 0, from Eq. (34), we obtain

‖um − un‖ = ‖um − um−1‖+ ...+ ‖un+1 − un‖ ≤ (Km−1 + ...+Kn)‖u1 − u0‖

≤ Kn(1 +K +K2 +K3 + ...)‖u1 − u0‖ =
Kn

1−K
‖u1 − u0‖. (35)

The estimated error is

‖u∗ − un‖ ≤
Kn

1−K
‖u1 − u0‖ (36)

while n→∞.

We establish the convergence of MGEM analogously.



F.A.Akgun, Z. Rasulov / Eur. J. Pure Appl. Math, 14 (3) (2021), 969-979 975

6. Numerical Examples

In this section, we give numerical examples to confirm the applicability of the main
results.

Example 1. Consider the fourth order BVP

u
′′′′

(t) = −3
u
′
u
′′′

1152
+

(u
′′
)2

576
+
t

4
+

95

4
, (37)

with the boundary conditions

u(0) = u
′
(0) = u

′′
(1) = u

′′′
(1) = 0. (38)

The exact solution of Eq. (37) − (38) is

u(t) = t4 − 4t3 + 6t2, (39)

and the Green’s function is

G(t, s) =

{
−t3
6 + st2

2 , 0 < t < s
−s3
6 + ts2

2 , s < t < 1.
(40)

By applying PGEM, we get

un+1 = un −
∫ t

0
(
−s3

6
+
ts2

2
)(u

′′′′
n (s) + 3

u
′
nu
′′′
n

1152
− (u

′′
n)2

576
− s

4
− 95

4
)ds (41)

−
∫ 1

t
(
−t3

6
+
st2

2
)(u

′′′′
n (s) + 3

u
′
nu
′′′
n

1152
− (u

′′
n)2

576
− s

4
− 95

4
)ds,

where the starting function is u0 = 0. The absolute error of the problem is estimated by

Err = |u(t)− un(t)|. (42)

Table 1 shows the maximum errors obtained by PGEM for several iterations, which
indicate high accuracy. Meanwhile, Table 2 lists the errors of the results obtained by
PGEM, MGEM, and contraction mapping for comparison. Here, it is assumed that α =
0.99. From Table 2 it is clear that PGEM has a better rate of convergence than other
methods.

Table 1: Maximum Errors of Exercise 4.1.

Number of Iterations Maximum Error

3 1.75E − 07
5 3.42E − 12
7 6.68E − 17



F.A.Akgun, Z. Rasulov / Eur. J. Pure Appl. Math, 14 (3) (2021), 969-979 976

Table 2: Error comparisons of Exercise 4.1.

Picard Contraction Mapping Mann

t Error (7) Error (7) Error(7)

0.1 1.64E − 18 4.44E − 16 9.76E − 15
0.2 5.85E − 18 1.33E − 15 3.51E − 14
0.3 1.17E − 17 1.89E − 15 7.10E − 14
0.4 1.87E − 17 3.55E − 15 1.14E − 13

Example 2. Consider the following fourth-order BVP

u
′′′′

(t) =
u(t)

6
(u(t) + u

′
(t) + u

′′
(t)− u′′′(t)) + 1, (43)

with the boundary conditions

u(0) = u
′
(0) = u

′′
(1) = u

′′′
(1) = 0. (44)

Note that problem Eq. (43) − (44) does not have an exact solution. On the other hand,
the Green’s function of the problem is

G(t, s) =

{
−t3
6 + st2

2 , 0 < t < s
−s3
6 + ts2

2 , s < t < 1.
(45)

Using Eq. (45)and PGEM,

un+1 = un −
∫ t

0
(
−s3

6
+
ts2

2
)(u

′′′′
n (s)− un(s)

6
(un(s) + u

′
n(s) + u

′′
n(s)− u′′′n (s))− 1)ds

−
∫ 1

t
(
−t3

6
+
st2

2
)(u

′′′′
n (s)− un(s)

6
(un(s) + u

′
n(s) + u

′′
n(s)− u′′′n (s))− 1)ds (46)

is obtained. Here, the starting function is u0 = 0. In Table 3, the errors obtained by
PGEM and MGEM are illustrated. For MGEM, α = 0.99 and α = 0.80 are chosen,
respectively.From Table 3, it is clear that the closer α to 1, the higher accuracy for MGEM
is obtained. Moreover, Picard-Green’s method is the special case of Mann-Green’s method,
when α = 1.

Example 3. Consider the fourth-order BVP

u
′′′′

(t) + 4u(t) = 1, (47)

with the boundary conditions

u(−1) = u(1) = u
′
(−1) = u

′
(1) = 0. (48)

The Green’s function of the problem Eq.(47) − (48)is

G(t, s) =

{
s3−3s+2

24 (−1− t)3 + s3−s2−s+1
8 (−1− t)2, −1 < s < t

s3−3s−2
24 (1− t)3 + −s3−s2+s+1

8 (1− t)2, t < s < 1.
(49)



F.A.Akgun, Z. Rasulov / Eur. J. Pure Appl. Math, 14 (3) (2021), 969-979 977

Table 3: Error comparisons of Exercise 4.2.

Picard Mann

t Error (7) Error (7) Error (7)

0.1 5.89E − 15 2.63E − 13 1.55E − 07
0.3 4.72E − 14 1.89E − 12 1.22E − 06
0.5 1.16E − 13 4.63E − 12 2.96E − 06
0.7 1.98E − 13 7.94E − 12 5.04E − 06
0.9 2.86E − 13 1.15E − 11 7.26E − 06

When Eq.(49) is embedded into PGEM, the iterative method is obtained as

un+1 = un −
∫ t

−1
(
s3 − 3s− 2

24
(1− t)3 +

−s3 − s2 + s+ 1

8
(1− t)2)(u4n(s) + 4un(s)− 1)ds

−
∫ 1

t
(
s3 − 3s+ 2

24
(−1− t)3 +

s3 − s2 − s+ 1

8
(−1− t)2)(u4n(s) + 4un(s)− 1)ds.(50)

In Eq.(50) the starting function u0 = 0. The numerical results obtained by PGEM at 20th
iteration are presented with the errors in Table 4. In this case, to apply MGEM, α = 0.95
and α = 0.70 were chosen. The Figure 1 below represents the 20th iteration by PGEM
which is very close to the exact solution of the problem Eq. (47) − (48) .

Table 4: Error comparisons of Exercise 4.3.

Picard Mann

t Numerical Value Error (20) Error(20) Error (20)

−1.0 0.0 3.49E − 18 3.49E − 18 3.49E − 18
−0.8 0.004829301786589 9.19E − 19 9.13E − 19 1.12E − 14
−0.6 0.015199277038320 5.81E − 19 6.01E − 19 1.80E − 14
−0.4 0.026098892616535 1.14E − 18 1.17E − 18 8.75E − 15
−0.2 0.034019288419593 8.72E − 19 9.19E − 19 6.17E − 15
0.0 0.036887761992940 1.01E − 19 5.05E − 20 1.31E − 14

7. Conclusion

In this paper, we present the Picard-Green’s method, one of the most popular methods,
generalized and extended for the fourth-order nonlinear and linear BVP by embedding
Green’s function into the Picard-Green’s Mann-Green’s fixed point iteration methods.
Different examples were solved to demonstrate its accuracy. The obtained results were
compared to other numerical and analytical solutions earlier found in the literature during
the solution. Moreover, we proved the convergence and found the rate of convergence.



REFERENCES 978

Figure 1: 20th iteration of Example 3

Acknowledgements

The authors thank the referees for their positive recommendations.

References

[1] Waleed Al-Hayani and Luis Casasus. Approximate analytical solution of fourth order
boundary value problems. Numerical Algorithms, 40(1):67–78, 2005.

[2] Arshed Ali. Numerical solution of fourth order boundary-value problems using Haar
wavelets. Applied Mathematical Sciences, 5(63):3131–3146, 2011.

[3] Shih-Hsiang Chang. Existence-uniqueness and fixed-point iterative method for general
nonlinear fourth-order boundary value problems. Journal of Applied Mathematics and
Computing, pages 1–11, 2021.

[4] Gilbert Choudury and Philip Korman. Computation of solutions of nonlinear bound-
ary value problems. Computers & Mathematics with Applications, 22(8):49–55, 1991.

[5] Alicia Cordero, Eva G Villalba, Juan R Torregrosa, and Paula Triguero-Navarro.
Convergence and Stability of a Parametric Class of Iterative Schemes for Solving
Nonlinear Systems. Mathematics, 9(1):86, 2021.

[6] Quang A Dang and Thi Kim Quy Ngo. Existence results and iterative method for
solving the cantilever beam equation with fully nonlinear term. Nonlinear Analysis:
Real World Applications, 36:56–68, 2017.



REFERENCES 979

[7] A Golbabai and M Javidi. Application of homotopy perturbation method for solv-
ing eighth-order boundary value problems. Applied Mathematics and Computation,
191(2):334–346, 2007.

[8] Sana Keita, Abdelaziz Beljadid, and Yves Bourgault. Efficient second-order semi-
implicit finite element method for fourth-order nonlinear diffusion equations. Com-
puter Physics Communications, 258:107588, 2021.

[9] Yongxiang Li. Existence of positive solutions for the cantilever beam equations with
fully nonlinear terms. Nonlinear Analysis: Real World Applications, 27:221–237,
2016.

[10] W Robert Mann. Mean value methods in iteration. Proceedings of the American
Mathematical Society, 4(3):506–510, 1953.

[11] Tsung Yen Na. Computational methods in engineering boundary value problems.
Academic press, 1980.

[12] Emile Picard. Memoir on the theory of partial derivative equations and the method
of successive approximations. Journal of Pure and Applied Maths, 6:145–210, 1890.

[13] Yongfang Wei, Qilin Song, and Zhanbing Bai. Existence and iterative method for
some fourth order nonlinear boundary value problems. Applied Mathematical Letters,
87:101–107, 2019.

[14] Ting-Ting Wu, Yu-Fei Yang, and Zhi-Feng Pang. A modified fixed-point iterative
algorithm for image restoration using fourth-order pde model. Applied numerical
mathematics, 62(2):79–90, 2012.

[15] Zhiyong Xing, Liping Wen, and Hanyu Xiao. A fourth-order conservative difference
scheme for the Riesz space-fractional Sine-Gordon equations and its fast implemen-
tation. Applied Numerical Mathematics, 159:221–238, 2021.