EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS
2025, Vol. 18, Issue 2, Article Number 5704
ISSN 1307-5543 – ejpam.com
Published by New York Business Global

Development of the Nyström Method for Weakly
Singular Functional Integral Equations

B. H. Alrikabi1, P. Darania1,∗, S. Pishbin1

1 Department of Mathematics, Faculty of Science, Urmia University, P.O.Box 165,
Urmia-Iran

Abstract. In this research, we apply the standard product integration method (Nyström method)
for solving the delay nonlinear weakly singular Volterra integral equations. Typically, in weakly
singular integral equations, the singularity of the kernel leads to the derivatives of the solution be-
coming singular at the boundary of the domain. The Chelyshkov polynomials serving as orthogonal
polynomials, find application in numerical integration. Here, we use roots of these polynomials to
make Lagrange interpolating polynomial for approximating the kernel functions in weakly singu-
lar functional integral equation. The proposed method’s convergence analysis is developed, and
numerical examples demonstrate the method’s reliability and efficiency.

2020 Mathematics Subject Classifications: 65R20, 65L20, 65L80, 34K28

Key Words and Phrases: Functional equations, Nyström method, Non-vanishing delays, Van-
ishing delays, Convergence

1. Introduction

The focal issue involves delay nonlinear weakly singular Volterra integral equations{
y(t) = f(t) + (Vαy)(t) + (Vα,θy)(t), t ∈ (t0, T ],
y(t) = ϕ(t), t ∈ [θ(t0), t0].

(1)

where

(Vαy)(t) =

∫ t

t0

P1,h(t, w)k1(t, w, y(w))dw, (2)

with k1 ∈ C(D × R), D = {(t, w) : t0 ≤ w ≤ t ≤ T}, k2 ∈ C(Dθ × R), Dθ =
{(t, w) : θ(t0) ≤ w ≤ θ(t)} and

(Vα,θy)(t) =

∫ θ(t)

t0

P2,h(t, w)k2(t, w, y(w))dw, (3)

∗Corresponding author.
DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.5704

Email addresses: Bahaahussainmath@gmail.com (B. H. Alrikabi),
p.darania@urmia.ac.ir (P. Darania), s.pishbin@urmia.ac.ir (S. Pishbin)

https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0)



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 2 of 19

ϕ, f as given functions are at least continuous on their respective domains. Normal forms
of weakly singular kernels P.,h(t, w), are

P.,h(t, w) =


1

|t−w|µ , 0 < µ < 1, h = 1,

log |t− w|, h = 2.

(4)

Here, the delay function θ will be constrained by the following conditions:

(O1) θ(t) = t− η(t), τ ∈ Cν(I), I = [t0, T ], for some ν ≥ 0,

(O2)


For vanishing delay: η(t0) = 0, and η(t) > t0 > 0, for t0 < t ≤ T,

For non-vanishing delay: η(t) ≥ η0 > 0, for all t ∈ I,

(O3) θ is strictly increasing on I.

Due to the assumption (O2) that the delay η(t) does not become zero on I, having
smooth data in (1) typically does not result in a solution y that is globally smooth. It
is widely recognized that these equations often exhibit discontinuity in the solution or its
derivatives at the initial point of the integration domain. This discontinuity propagates
along the integration interval, giving rise to subsequent points referred to as singular
points. Determining these singular points in advance is challenging, and the solution
derivatives at these points tend to smooth out along the interval. Many existing numeri-
cal methods for such equations are highly sensitive to these singular points. Therefore, it is
crucial for these methods to incorporate a process that detects and includes these points in
the mesh to ensure the desired accuracy (For further details see [1].) By virtue of Volterra
integral equations and also their functional counterparts along with delay terms whether
vanishing or non-vanishing, a wide spectrum of science subjects, namely in biology, ecol-
ogy, physics, and chemistry have been mathematically well-formulated in order to analyse
and study underlying phenomena. Particularly, these classes of mathematical modellings
can be found in fluid dynamics, viscoelasticity of materials, population growth dynam-
ics, heat conduction, epidemiology, controlled liquidation in obsolete production units,
and renovation in economic systems, see [2–6] and references therein. Some numerical
methods have been proposed to obtain approximate solutions of weakly singular Volterra
integral equation, such as spectral method [7–9], collocation method [10, 11], radial basis
function method [12], Bernstein and Genocchi polynomial method [13, 14] and successive
approximation method [15].

Delay IEs have been solved approximately by many authors, such as, spectral methods
[16], two-point multistep block method with constant step-size [17], collocation and con-
tinuous implicit Runge-Kutta methods [18], spline collocation method [19, 20], piecewise
collocation method [21], and linear multistep methods [22]. Also, one- step polynomial
collocation method [1] has been applied to find numerical solution of delay IEs and delay
weakly singular IEs. In [23], a category of Volterra delay integral equations (VDIEs) in-
volving noncompact operators is estimated using collocation methods. The study explores



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 3 of 19

the characteristics of the associated operators and delves into the discussions regarding
the existence, uniqueness, and regularity of the exact solution. The paper establishes the
existence and uniqueness of collocation solutions, specifically under two distinct graded
meshes. Additionally, the convergence conditions and order of convergence are presented.
To validate the theoretical orders of convergence, numerical examples are provided.

The paper is organized as follows: Section 2 is dedicated to proposing some preliminar-
ies about the product rules using the roots of certain types of generalized Jacobi orthogonal
polynomials and specific properties of orthogonal polynomials. In Section 3, we construct
and analyze the product integration method to solve the equation (1) numerically and
in section 4 convergence of numerical solutions is investigated. This is succeeded by the
discussion of three test problems in Section 5 to validate the theoretical results. Finally, in
Section 6, we conclude the paper and suggest potential future avenues for research, which
are currently less explored.

2. Preliminaries

Consider the product rules∫ 1

0
P.,h(t, w)G(w)dw ≈

m∑
i=1

wm,i(t)G(tm,i) = IN (G, t), (5)

derived from the roots of particular classes of generalized Jacobi orthogonal polynomials.
We can consider the error term of (5) as

RN (G; t) =

∫ 1

0
P.,h(t, w)G(w)dw − IN (G, t). (6)

From[24, 25], we have RN (G, t) = O(N−m) when G ∈ Cm[0, 1]. If G is a polynomial of
degree N , then RN (G, t) = 0, so from [24], for all polynomial PN of degree N , we have

RN (G, t) =

∫ 1

0
P.,h(t, w)(G(w)− PN (w))dw −

∫ 1

0
P.,h(t, w)ΩN (G− PN , w)dw,

where ΩN (G, t) represents the Lagrange interpolation polynomial that interpolates G at
the points {ti}Ni=0, the expression is defined as follows:

ΩN (G, t) =
N∑
i=0

G(ti)lN,j(t),

where

lN,j(t) =
N∏
i=0
i ̸=j

t− ti
tj − ti

, j = 0, 1, ..., N. (7)



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 4 of 19

With a careful selection of the polynomial sequence {ΩN}, it becomes possible to establish
upper bounds as:

R1,N (G, t) =

∫ 1

0
| P.,h(t, w) || G(w)− PN (w) | dw,

R2,N (G, t) =

∫ 1

0
| P.,h(t, w) || ΩN (G− PN , w) | dw,

where P.,h(t, w) is defined in (4).

Theorem 2.1. (From [8]) Assume that G(t) = (1− t)µ where µ > −1 is not an integer
and τ > −1, with µ+ τ > −1, then∫ 1

0
| G(t)− ΩN (G, t) || t− w |τ dt ≤ C

{
N−2−2µ−2τ logN, | w |≤ 1, τ < 0,
N−2−2µ logN, | w |≤ 1, τ ≥ 0,

(8)

such that C represents a constant that remains unaffected by both w and N .

Corollary 2.2. (From [8]) Assume that G(t) = (1− t)µ where µ > 0, is not an integer,
then we have∫ 1

0
| G(t)− ΩN (G, t) || log | t− w || dt ≤ C


N−2−2µ log2N, | w |≤ 1,

N−2−2µ logN, 0 ≤ w < 1,
(9)

where C represents a constant that remains unaffected by both w and N .

2.1. Orthogonal polynomials

Chelyshkov has introduced polynomials in [26], specifically designed to be orthogonal
over the interval [0, 1] which the Rodrigues type representation as

Pn,l(t) =
1

(n− l)!

1

tl+1

dn−l

dtn−l
(tn+l+1(1− t)n−l), l = 0, 1, . . . , n,

and are explicitly characterized by

Pn,l(t) =
n−l∑
j=0

(−1)j
(
n− l

j

)(
n+ l + 1 + j

n− l

)
tl+j , l = 0, 1, ..., n. (10)

Given a fixed value for n, the orthogonality property within the interval [0, 1] estab-
lishes an immediate link between the polynomials Pn,l(t) and a set of Jacobi polynomials

P
(α,β)
m (t) [26] as:

Pn,l(t) = (−1)n−ltlP
(0,2l+1)
n−l (2t− 1), l = 0, 1, . . . , n. (11)



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 5 of 19

From [27], Jacobi polynomials are the polynomial eigenfunctions of the singular Sturm-
Liouville problem. An explicit formula is given by

P (α,β)
m (t) =

1

2m

m∑
j=0

(
m+ α

j

)(
m+ β

m− j

)
(t− 1)m−j(t+ 1)j .

An important consequence of the symmetry of the weight function w(t) = (1− t)α(1+ t)β,
and the orthogonality of the Jacobi polynomials, is the symmetry relation

P (α,β)
m (t) = (−1)mP (α,β)

m (−t).

The ultraspherical polynomials are simply Jacobi polynomials with α = β, and normalized
differently:

P (α)
m (t) =

Γ(α+ 1)Γ(m+ 2α+ 1)

Γ(2α+ 1)Γ(m+ α+ 1)
P (α,α)
m (t),

where Γ(·) is the gamma function. The relation between Legendre and Chebyshev poly-
nomials with the ultraspherical polynomials is

Pm(t) = P (0)
m (t), Tm(t) =

P
(− 1

2
,− 1

2
)

m (t)

P
(− 1

2
,− 1

2
)

m (1)
.

Also for considering the other orthogonal polynomials, we can refer to [28, 29]. Accord-
ing to [26], for each value of n, the polynomial Pn,0(t) has precisely n distinct roots within
the interval (0, 1) and this set of polynomials exhibits all the typical characteristics found
in other widely recognized orthogonal polynomial families such as Legendre or Chebyshev
polynomials. Yon can see that distribution of roots of these polynomials in Fig. 1 and
Fig. 2.

3. Nyström method

Here, we outline the Nyström method employed for the numerical solution of equation
(1). In equation (1), for the sake of convenience and without any loss of generality, we
assume certain conditions

y(t) = f(t) + (Vαy)(t) + (Vα,θy)(t), t ∈ (0, 1],

y(t) = ϕ(t), t ∈ [θ(0), 0].
(12)

By chosen (N + 1) distinct points {ti}Ni=1

⋃
{t0 = 0}, in the interval Ih ⊆ [0, 1] and

collocate the equation (12) at the underlying local mesh {ti}Ni=0, we have

y(ti) = f(ti) + (Vαy)(ti) + (Vα,θy)(ti), i = 0, 1, ..., N, (13)

where

(Vαy)(ti) =

∫ ti

0
P1,h(ti, w)k1(ti, w, y(w))dw, (14)



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 6 of 19

0.0 0.2 0.4 0.6 0.8

0.00

0.05

0.10

0.15

0.20

Roots of Chebyshev polynomial

Roots of Legendre polynomial

Roots of Chelyshkov polynomial

Figure 1: Plot of distribution of roots of Chelyshkov, Legendre and Chebyshev polynomials with m = 4.

0.0 0.2 0.4 0.6 0.8

0.00

0.05

0.10

0.15

0.20

Roots of Chebyshev polynomial

Roots of Legendre polynomial

Roots of Chelyshkov polynomial

Figure 2: Plot of distribution of roots of Chelyshkov, Legendre and Chebyshev polynomials with m = 6.

(Vα,θy)(ti) =

∫ θ(ti)

0
P2,h(ti, w)k2(ti, w, y(w))dw. (15)



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 7 of 19

Subsequently, we employ the Lagrange interpolation polynomial

ΩN (kh, w) =

N∑
j=0

lN,j(w)kh(ti, tj , y(tj)), h = 1, 2, (16)

to approximate kh(ti, w, y(w)) and get

yN (ti) = f(ti) +

N∑
j=0

[W1,i,j k1(ti, tj , yN (tj))

+ W2,i,j k2(ti, tj , yN (tj))] ,

(17)

with

W1,i,j =

∫ ti

0
P1,h(ti, w)lN,j(w)dw,

W2,i,j =

∫ θ(ti)

0
P2,h(ti, w)lN,j(w)dw, i, j = 0, 1, 2, · · · , N.

(18)

Note that the relationship given in equation (17) constitutes a nonlinear system of equa-
tions with dimensions (N + 1) × (N + 1). This system possesses a unique solution, as
demonstrated in sources such as [30] and [31]. The resolution of this nonlinear system
yields the values of yN (ti) for i = 0, 1, 2, ..., N , representing the solutions of (12) at the
points {ti}Ni=0.

Remark 3.1. Using the relation (16), approximate the integrals in (12), obtaining a new
equation:

yN (t) = f(t) +
N∑
j=0

(
W1,t,j k1(t, tj , yN (tj)) +W2,t,j k2(t, tj , yN (tj))

)
, (19)

where for j = 0, 1, 2, . . . , N ,

W1,t,j =

∫ t

0
P1,h(t, w)lN,j(w)dw, W2,t,j =

∫ θ(t)

0
P2,h(t, w)lN,j(w)dw,

We write this as an exact equation with a new unknown function yN (t). To find the
solution at the node points, let y(t) run through the quadrature node points ti. This yields
the nonlinear system (17) of order N+1. Each solution yN (t) of (19) furnishes a solution
to (17): merely evaluate yN (t) at the node points. The converse is also true. To each
solution (yN (t0), . . . , yN (tN ))T of (17), there is a unique solution of (19). Therefore, by
inserting (yN (t0), . . . , yN (tN ))T in (19), we can obtain the approximate solution in each
desired point t. In fact, the relation (19) is the Nyström interpolation formula.



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 8 of 19

With these symbols in place, we can encapsulate the steps in the following algorithm:

Algorithm 1.

Input:
N ;
begin
For l=0,1,2. . . ,N:

Compute tl as simple roots of N + 1 st-degree orthogonal polynomial in [0, 1];
For i=0,1. . . ,N:
For j=0,1,. . . ,N:

For r=1,2:
Compute Wr,i,j from (18);
Solve nonlinear system (17) and Compute yN (ti);

end.

4. Convergence Theorem

Here, we investigate the convergence analysis of the equation

y(t) = f(t) + (Vαy)(t) + (Vα,θy)(t), t ∈ [0, 1]. (20)

where

(Vαy)(t) =

∫ t

0
P1(t, s)k1(t, s)y(s)ds,

(Vα,θy)(t) =

∫ θ(t)

0
P2(t, s)k2(t, s)y(s)ds.

(21)

Assuming that the function f(t) belongs to the space C([0, 1]) and the kernels Ph, where
h = 1, 2, exhibit weak singularity in the forms (4), the equation (20) possesses a distinct
solution y in the interval C[0, 1]. It is anticipated that this solution may have unbounded
derivatives at the endpoints.

If we utilize the method (17) on the test problem (20) for a given mesh {ti}Ni=1

⋃
{0},

the resulting approximate solution yN (t) takes the following form

yN (t) = f(t) +

N∑
j=0

(ω1,j(P1,h, t) + ω2,j(P2,h, θ(t)))yN (tj), (22)

where

ω1,j(P1,h, t) =

∫ t

0
P1,h(ti, w)lN,j(w)dw,

ω2,j(P2,h, θ(t)) =

∫ θ(t)

0
P2,h(ti, w)lN,j(w)dw.

(23)



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 9 of 19

To assess the uniform convergence of yN (t) to y(t) as solution of (20), we rewrite

y(t)− yN (t) =

N∑
j=0

ω1,j(P1,h, t)(y(tj)− yN (tj))

+
N∑
j=0

ω2,j(P2,h, θ(t))(y(tj)− yN (tj))

+T1,N (P1,h, y; t) + T2,N (P2,h, y; θ(t)),

(24)

where

T1,N (P1,h, y; t) =

∫ t

0
P1,h(t, w)y(w)dw −

N∑
j=0

ω1,j(P1,h, t)y(tj),

T2,N (P2,h, y; θ(t)) =

∫ θ(t)

0
P2,h(t, w)y(w)dw −

N∑
j=0

ω2,j(P2,h, θ(t))y(tj).

(25)

If we set

ωj(P1,h,P2,h; t, θ(t)) = ω1,j(P1,h, t) + ω2,j(P2,h, θ(t)),

TN (P1,h,P2,h, y; t, θ(t)) = T1,N (P1,h, y; t) + T2,N (P2,h, y; θ(t)),
(26)

then equation (24) reduces to the following equation

y(t)− yN (t) =
N∑
j=0

ωj(P1,h,P2,h; t, θ(t))(y(tj)− yN (tj))

+TN (P1,h,P2,h, y; t, θ(t)).

(27)

Now, define linear operator AN as:
AN : C[0, 1] −→ C[0, 1]

AN (f(t)) =
N∑
j=0

ωj(P1,h,P2,h; t, θ(t))f(tj), f ∈ C[0, 1],
(28)

then

∥ y(t)− yN (t) ∥∞= ∥ AN(y(t)− yN (t))+ TN (P1,h,P2,h, y; t, θ(t)) ∥∞

≤ ∥ AN ∥∞∥ y(t)− yN (t) ∥∞ + ∥ TN (P1,h,P2,h, y; t, θ(t)) ∥∞,
(29)



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 10 of 19

and by considering (27), we have

∥ y(t)− yN (t) ∥∞≤∥ (I −AN )−1 ∥∞∥ TN ∥∞ . (30)

Our ultimate objective is to establish an upper bound for (30). We begin by preparing a
set of auxiliary theorems and lemmas concerning kernels of types (4).

Theorem 1. Consider {ti}Ni=0 as the zeros of the (N + 1)-th-degree member of a set of
polynomials that are orthogonal on Ih ⊆ [0, 1] with respect to the weight function

W (t) = g(t)(1− t)µ̄(1 + t)ν̄ , −1 < µ̄ ≤ 3

2
, ν̄ ≥ −1

2
, (31)

where g(t) is a positive and continuous function within the interval Ih, and the modulus of

continuity u for g fulfills the condition

∫ 1

0
u(g, s)

ds

s
< ∞. Also, let ΩN (y, w) represent the

interpolating polynomial of degree at most N that matches the function y at {ti}Ni=1∪{t0 =
0}. Additionally, assume that P1,h(t, w) and P2,h(t, w) are kernels of type (4). Then, for
any function y ∈ C(Ih), we have

lim
N→∞

∥
∫ t

0
P1,h(t, w)(y(w)− ΩN (y, w))dw +

∫ θ(t)

0
P2,h(t, w)(y(w)− ΩN (y, w))dw∥ = 0.

(32)
Especially, for 0 < µ < 1, the limits are as follows

∥TN (|t− w|−µ, |t− w|−µ, y; t, θ(t))∥∞ = O{h1}, (33)

∥TN (log |t− w|, log |t− w|, y; t, θ(t))∥∞ = O{h2}, (34)

∥TN (|t− w|−µ, log |t− w|, y; t, θ(t))∥∞ = O{h3}, (35)

where h1 = (N+1)2µ−2κ−2 log(N+1), h2 = (N+1)−2−2κ log2(N+1) and h3 = min{h1, h2}.



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 11 of 19

Proof. By using the equation (26), we have

| TN (P1,h,P2,h, y; t, θ(t)) |= |
∫ t

0
P1,h(t, w)y(w)dw −

N∑
j=0

ω1j(P1,h, t)y(tj)

+

∫ θ(t)

0
P2,h(t, w)y(w)dw −

N∑
j=0

ω2j(P2,h, θ(t))y(tj) |

≤
∫ t

0
| P1,h(t, w) || y(w)− ΩN (y, w) | dw

+

∫ θ(t)

0
| P2,h(t, w) || y(w)− ΩN (y, w) | dw

≤
∫ 1

0
| P1,h(t, w) || y(w)− ΩN (y, w) | dw

+

∫ 1

0
| P2,h(t, w) || y(w)− ΩN (y, w) | dw.

(36)

The proof is readily derived from an outcome of Theorem 2.2. The inequalities (33) and
(35) provide a measure of the convergence rate.

Now, let’s examine the characteristics of the initial term ∥ (I −AN )−1 ∥∞

Theorem 2. For a given set of nodes {ti}Ni=0 defined as Theorem 1 with the restriction
−1

2 < µ̄, ν̄ < 3
2 , let lN,j(w) show the corresponding j−th Lagrange polynomial. Fix a

subinterval [a, b] ⊆ [0, 1]. Then, there is a positive number C and a value of q1 greater
than 1, such that

sup

N∑
j=0

|
∫ b

a
(P1,h(b, w) + P2,h(b, w))lN,j(w)dw |≤ C(∥ P1,h(b, w) ∥q1 + ∥ P2,h(b, w) ∥q1),

(37)

for all P.,h ∈ Lq1 with ∥ P.,h ∥q1= (

∫ 1

0
| P.,h(b, w) |q1 dw)

1
q1 .

Proof. In the same manner in [32], we set B = {f ∈ C[0, 1] : ∥f∥∞ = 1}. Then we



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 12 of 19

have

N∑
j=0

|
∫ b

a
[P1,h(b, w) + P2,h(b, w)]lN,j(w)dw |

= sup
f∈B

|
∫ b

a

N∑
j=0

(P1,h(b, w) + P2,h(b, w))lN,j(w)f(tj)dw |

≤ sup
f∈B

∫ b

a
|

N∑
j=0

lNj(w)f(tj) | {| P1,h(b, w) | + | P2,h(b, w) |}dw

≤ sup
f∈B

{[
∫ b

a
|

N∑
j=0

lNj(w)f(tj) |q2 dw]
1
q2 {[

∫ b

a
| P1,h(b, w) |q1 dw]

1
q1 + [

∫ b

a
| P2,h(b, w) |q1 dw]

1
q1 }}

≤ sup
f∈B

{[
∫ 1

0
|

N∑
j=0

lNj(w)f(tj) |q2 dw]
1
q2 {[

∫ b

a
| P1,h(b, w) |q1 dw]

1
q1 + [

∫ b

a
| P2,h(b, w) |q1 dw]

1
q1 }}

≤ sup
f∈B

{∥
N∑
j=0

lNj(w)f(tj) ∥q2 [∥ P1,h(b, w) ∥q1 + ∥ P2,h(b, w) ∥q1 ]},

(38)
for all P.,h ∈ Lq1 , h = 1, 2 with q1, q2 > 1 such that 1

q1
+ 1

q2
= 1. Now, from [25], given

the conditions specified in this theorem, we obtain:

sup
f∈B

∥
N∑
j=0

lN,j(w)f(tj) ∥q2≤ C ∥ f ∥∞, (39)

for any bounded function f and 0 < q2 < ∞. Hence, the bound (37) follows.

Lemma 1. If the kernels Ph satisfy{
P.,h ∈ Lq1 , q1 > 1, h = 1, 2,
lim
t1→t

∥ P.,h(t1, w)− P.,h(t, w) ∥q1= 0, ∀t ∈ Ih.
(40)

Then

lim
(t1,θ(t1))→(t,θ(t))

sup
N

N∑
j=0

| ωj(P1,h,P2,h; t1, θ(t1))− ωj(P1,h,P2,h; t, θ(t)) |= 0, (41)

for all t ∈ Ih.



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 13 of 19

Proof. By using ωj = ω1j + ω2j we get

sup
N

{
N∑
j=0

| ωj(P1,h,P2,h; t1, θ(t1))− ωj(P1,h,P2,h; t, θ(t)) |}

= sup
N

{
N∑
j=0

| [ω1j(P1,h, t1)− ω1j(P1,h, t)] + [ω2j(P2,h, t1)− ω2j(P2,h, t)] |}

= sup
N

{
N∑
j=0

|
∫ t1

0
P1,h(t1, w)lN,j(w)dw −

∫ t

0
P1,h(t, w)lN,j(w)dw

+

∫ θ(t1)

0
P2,h(t1, w)lN,j(w)dw −

∫ θ(t)

0
P2,h(t, w)lN,j(w)dw |}

= sup
N

{
N∑
j=0

|
∫ t

0
[P1,h(t1, w)− P1,h(t, w)]lN,j(w)dw +

∫ t1

t
P1,h(t1, w)lN,j(w)dw

+

∫ θ(t)

0
[P2,h(t1, w)− P2,h(t, w)]lN,j(w)dw +

∫ θ(t1)

θ(t)
P2,h(t1, w)lN,j(w)dw |}

≤ sup
N

{
N∑
j=0

|
∫ t1

0
P1,h(t1, w)lN,j(w)dw | −

N∑
j=0

|
∫ t

0
P1,h(t, w)lN,j(w)dw |

+
N∑
j=0

|
∫ θ(t1)

0
P2,h(t1, w)lN,j(w)dw | −

N∑
j=0

|
∫ θ(t)

0
P2,h(t, w)lN,j(w)dw |}

≤ C{[
∫ t1

0
| P1,h(t1, w)− P1,h(t, w) |q1 dw]

1
q1 + [

∫ t1

t
| P1,h(t1, w) |q1 dw]

1
q1

+[

∫ t1

0
| P2,h(t1, w)− P2,h(t, w) |q1 dw]

1
q1 + [

∫ θ(t1)

θ(t)
| P2,h(t1, w) |q1 dw]

1
q1 }

≤ C{[
∫ 1

0
| P1,h(t1, w)− P1,h(t, w) |q1 dw]

1
q1 + [

∫ t1

t
| P1,h(t1, w) |q1 dw]

1
q1

+[

∫ 1

0
| P2,h(t1, w)− P2,h(t, w) |q1 dw]

1
q1 + [

∫ θ(t1)

θ(t)
| P2,h(t1, w) |q1 dw]

1
q1 }

≤ C{∥ P1,h(t1, w)− P1,h(t, w) ∥q1 +[

∫ θ(t1)

θ(t)
| P2,h(t1, w) |q1 dw]

1
q1

+ ∥ P2,h(t1, w)− P2,h(t, w) ∥q1 +[

∫ θ(t1)

θ(t)
| P2,h(t1, w) |q1 dw]

1
q1 }.

(42)



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 14 of 19

The Lemma now follow from these statements.

Theorem 3. Consider the operator AN defined as in (28), and let the nodes {ti}Ni=0 be
selected as outlined in Theorem 1. If conditions (32), (37), and (41) are satisfied, then
for all sufficiently large values of N , there exists a constant C > 0 that is independent of
N , so that

∥ (I −AN )−1 ∥≤ C. (43)

Proof. The proof immediately stems from an implication of Theorem 2 in [32] and
Lemma 1 in [33].

The outcomes of our efforts in this section lead to the following principal theorem:

Theorem 4. Consider y(t) and yN (t) as the exact and approximated solutions, respec-
tively, to the equation (20). These solutions are constructed based on a collection of distinct
nodes {ti}Ni=1∪{t0 = 0}. If the nodes {ti} are the zeroes of the orthogonal polynomial in Ih
and ph(t, w) is the kernel function of the form (4), then yN (t) converges uniformly to y(t).
Furthermore, the convergence rate aligns with the product integration quadrature that we
select to approximate the integral term (20).

5. Numerical results

In this section, we present the numerical outcomes of various test problems solved using
the method proposed in this article. Different forms of kernels are taken into account for
computational purposes in the following test problems. The discretization algorithm relies
on nodes that coincide with the zeros of the specified orthogonal polynomials of the (N)-th
degree, along with t = 0. Additionally, a product integration method described in Section
3 is employed. All calculations were conducted using Mathematica software.

Example 5.1. The non-linear, weakly singular Volterra functional integral equation

y(t) = f(t) +

∫ t

0
ln |t− w| (tw2 − y2(w))dw +

∫ t
10

0
|t− w|−

1
2 y2(w)dw, t ∈ [0, 1]

with f(t) such that possesses the exact solution y(t) = t
13
2 .

Example 5.2. The nonlinear weakly singular Volterra functional integral equation in [0, 1]

y (t) = f (t) +

∫ t

0
|t− w|−

1
2 y2 (w) dw +

∫ t− 1
10

0
|t− w|−

1
2 y2 (w) dw,

with f(t) such that possesses the analytical solution y(t) = e2t.

We consider the other example which the exact solution has the low regularity. The
exact solution is y(t) =

√
t, then derivative of this solution is unbounded at t = 0 and

reflects the general qualitative regularity behaviour of the solution near t = 0+.



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 15 of 19

N Chelyshkov Polynomials Legendre Polynomials Chebyshev Polynomials

4 2.81× 10−2 3.87× 10−2 7.31× 10−2

5 9.50× 10−3 1.48× 10−2 1.77× 10−2

6 2.48× 10−3 3.98× 10−3 5.10× 10−3

7 6.34× 10−4 1.00× 10−3 1.38× 10−3

8 1.25× 10−4 2.10× 10−4 2.57× 10−4

9 2.02× 10−5 3.44× 10−5 4.30× 10−5

10 2.40× 10−6 3.90× 10−6 5.06× 10−6

Table 1: Maximum errors by using different types of orthogonal polynomials in Example 5.1.

N Chelyshkov Polynomials Legendre Polynomials Chebyshev Polynomials

4 1.70× 10−2 1.72× 10−2 2.13× 10−2

5 4.28× 10−3 4.42× 10−3 5.20× 10−3

6 8.35× 10−4 8.90× 10−4 9.93× 10−4

7 1.30× 10−4 1.45× 10−4 1.53× 10−4

8 1.61× 10−5 1.97× 10−5 1.98× 10−5

9 1.47× 10−6 2.31× 10−6 2.51× 10−6

10 1.64× 10−7 2.31× 10−7 6.71× 10−7

Table 2: Maximum errors by using different types of orthogonal polynomials in Example 5.2.

Example 5.3. The non-linear, weakly singular Volterra functional integral equation

y(t) = f(t)−
∫ t

0
|t− w|−

1
2 y3(w)dw +

∫ 0.9t

0
|t− w|−

1
2 y2(w)dw, t ∈ [0, 1]

with f(t) such that possesses the exact solution y(t) = t
1
2 .

We give the numerical solution of the Examples, at the root of N -st-degree orthogonal
polynomials. The maximum errors obtained using the presented method are compared
with the exact solution in Tables 1, 2 and 3. In the case of examples, the obtained nonlinear
systems are solved using Newton’s method.

N Chelyshkov Polynomials Legendre Polynomials Chebyshev Polynomials

4 6.16× 10−4 6.69× 10−4 1.16× 10−3

5 2.91× 10−4 3.44× 10−4 5.50× 10−4

6 1.54× 10−4 2.23× 10−4 3.19× 10−4

7 8.87× 10−5 1.44× 10−4 2.06× 10−4

8 5.45× 10−5 1.00× 10−4 1.40× 10−4

9 3.53× 10−5 6.93× 10−5 9.53× 10−5

10 2.38× 10−5 4.91× 10−5 4.67× 10−5

Table 3: Maximum errors by using different types of orthogonal polynomials in Example 5.3.



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 16 of 19

N Example 5.1 Example 5.2 Example 5.3

4 6.03× 10−4 8.60× 10−4 2.72× 10−5

6 3.60× 10−4 9.54× 10−6 4.07× 10−6

8 4.79× 10−6 3.41× 10−7 8.67× 10−7

10 9.95× 10−8 1.81× 10−7 2.86× 10−7

Table 4: Absolute errors at the point tN+1 = 1 for different values of N .

In Tables 1, 2 and 3, we report maximum errors in Examples 5.1, 5.2 and 5.3. We con-
sider Chelyshkov, Legendre and Chebyshev polynomials as orthogonal polynomials and
observe that our proposed numerical method works well for theses different types of or-
thogonal polynomials. Also, according to the dispersion of the roots of these polynomials,
it seems that the errors reported by Chelyshkov polynomials is slightly better than the
other polynomials.

Remark 5.4. To get the approximate solution at any given point η within the interval
Ih ⊆ [0, 1], our focus lies on a rule that relies on both β and y(η)∫ 1

0
y(w)dw ≈

N∑
j=0

βjy(tj) + βy(η).

From [34] and [35], we know that this method is exact for all polynomials whose degree does
not exceed 2N . Nevertheless, when considering quadrature rules that incorporate the point
η as one of their nodes, serving as a collocation point, we will be faced with a nonlinear
system of equations of size (N + 1) × (N + 1). The solutions to this system provide the
values at our grid points, particularly at tN+1 = η([8]).

In Table 4, we set tN+1 = 1 and report absolute value errors at this points for the
examples 5.1, 5.2 and 5.3.

6. Conclusion

Using incorporating the Nyström method, the non linear weakly singular Volterra
functional can be converted into a nonlinear system. This system can be solved by some
classical techniques. The method’s effectiveness and accuracy in solving nonlinear equa-
tions have been assessed through various problems. Another significant advantage of the
proposed method is that the unknown coefficients can be determined quite easily using
computer programs. In our future work, we will study the Nyström method to solve delay



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 17 of 19

weakly singular integral-algebraic equations in the following form:

B(t)X(t) +

∫ t

0
P1,h(t, s)K1(t, s)X(s)ds+

∫ θ(t)

0
P2,h(t, w)K2(t, s)X(s)ds = F (t), t ∈ I,

subject to det(B(t)) = 0, ∀t ∈ I.

References

[1] H. Brunner. Collocation Methods for Volterra Integral and Related Functional Equa-
tions. Cambridge University Press, Cambridge, 2004.

[2] H. Brunner and Y. Yatsenko. Spline collocation methods for non- linear volterra
integral equations with unknown delay. Journal of computational and applied math-
ematics, 71(1):67–81, 1996.

[3] A. Cardone and D. Conte. Multistep collocation methods for volterra integro-
differential equations. Applied Mathematics and Computation, 221:770–7856, 2013.

[4] W. Ming and C. Huang. Collocation methods for volterra functional integral equations
with non-vanishing delays. Applied mathematics and computation, 269:198–214, 2017.

[5] C. Huang W. Ming and L. Zhao. Optimal super- convergence results for volterra
functional integral equations with proportional vanishing delays. Applied Mathematics
and Computation, 320:292–301, 2018.

[6] M. Wang. Multistep collocation method for fredholm integral equations of the second
kind. Applied Mathematics and Computation, 420(126870):1–16, 2022.

[7] T. Tang C. Huang and Z. Zhang. Supergeometric convergence of spectral collocation
methods for weakly singular Volterra and Fredholm integral equations with smooth
solutions. Journal of Computational Mathematics, 29(6):698–719, 2011.

[8] M. Rasty and M. Hadizadeh. A product integration approach on new orthogonal
polynomials for nonlinear weakly singular integral equations. Acta Applicandae Math-
ematicae, 109:861–873, 2010.

[9] M.A. Zaky and I. G. Ameen. A novel jacobi spectral method for multi-dimensional
weakly singular nonlinear volterra integral equations with nonsmooth solutions. En-
gineering with Computers-Germany, 37:2623–2631, 2021.

[10] H. Cai and Y. Chen. A fractional order collocation method for second kind volterra
integral equations with weakly singular kernels. Journal of Scientific Computing,
75:970–992, 2018.

[11] T. Herdman Y. Cao and Y. Xu. A hybrid collocation method for volterra integral
equations with weakly singular kernels. SIAM Journal on Numerical Analysis, 41:364–
381, 2003.

[12] P. Assari and M. Dehghan. A meshless method for the numerical solution of nonlinear
weakly singular integral equations using radial basis functions. The European Physical
Journal Plus, 132:1–23, 2017.

[13] M. A. Ebadi E. Hashemizadeh and S. Noeiaghdam. Matrix method by genocchi
polynomials for solving nonlinear volterra integral equations with weakly singular
kernels. Symmetry, 12:1–16, 2020.



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 18 of 19

[14] J. Huang Y. B. Pan and Y.Y Ma. Bernstein series solutions of multidimensional
linear and nonlinear volterra integral equations with fractional order weakly singular
kernels. Applied Mathematics and Computation, 374:149–161, 2019.

[15] A. M. Wazwaz. Linear and Nonlinear Integral Equations. Springer, Heidelberg, 2011.
[16] H. Brunner I. Ali and T. Tang. Spectral methods for pantograph-type differential and

integral equations with multiple delays. Frontiers of Mathematics in China, 4:49–61,
2009.

[17] N. Senu N. A. Baharum, Z.A. Majid and H. Rosali. Numerical Approach for Delay
Volterra Integro-Differential Equation. Sains Malaysiana, 51(12):4125–4144, 2002.

[18] H. Brunner. Collocation and continuous implicit runge-kutta methods for a class of
delay volterra integral equations. Journal of Computational and Applied Mathematics,
53:61–72, 1994.

[19] E. Marchetti F. Calio and R. Pavani. About the deficient spline collocation method
for particular differential and integral equations with delay. Rendiconti del Seminario
Matematico - Politecnico di Torino, 61:287–300, 2003.

[20] R. Pavani F. Calio, E. Marchetti and G. Micula. About some volterra problems solved
by a particular spline collocation. Studia Universitatis Babes-Bolyai, 48:45–52, 2003.

[21] F. Ghoreishi M. Khasi and M. Hadizadeh. Numerical analysis of a high order method
for state-dependent delay integral equations. Numerical Algorithms, 66:177–201, 2013.

[22] S. Gan S. Wu. Errors of linear multistep methods for singularly perturbed volterra
delay integro-differential equations. Mathematics and Computers in Simulation,
79:3148–3159, 2009.

[23] Y. Xiao H. Song and M. Chen. Collocation methods for third-kind volterra integral
equations with proportional delays. Applied Mathematics and Computation, 388:1–11,
2020.

[24] G. Mastroianni G. Criscuolo and G. Monegato. Convergence properties of a class of
product formulas for weakly singular integral equations. Mathematics of Computation,
55:213–230, 1990.

[25] P. Nevai. Mean convergence of lagrange interpolation. Trans. Am. Math. Soc.,
282:669–689, 1984.

[26] V.S. Chelyshkov. Alternative orthogonal polynomials and quadratures. Electronic
Transactions on Numerical Analysis, 25(7):17–26, 2006.

[27] S. Gottlieb J. Hesthaven and D. Gottlieb. Spectral Methods for Time-dependent
Problems. Cambridge University Press, Cambridge, 2007.

[28] C. Cesarano. Generalized Chebyshev polynomials. Hacettepe Journal of Mathematics
and Statistics, 43(5):731–740, 2014.

[29] B. Germano C. Cesarano and P.E. Ricci. Laguerre-type Bessel functions. Integral
Transforms and Special Functions, 16(4):315–322, 2005.

[30] H. Kaneko and Y. Xu. Gauss-type quadratures for weakly singular integrals and their
application to fredholm integral equations of second kind. Mathematics of Computa-
tion, 62:739–753, 1994.

[31] L. Tao and H. Yong. Extrapolation method for solving weakly singular nonlinear
volterra integral equations of second kind. Journal of Mathematical Analysis and



B. H. Alrikabi, P. Darania, S.Pishbin / Eur. J. Pure Appl. Math, 18 (2) (2025), 5704 19 of 19

Applications, 324:225–237, 2006.
[32] A.P. Orsi. Product integration for volterra integral equations of the second kined

with weakly singular kernels. Mathematics of Computation, 212:1201–1212, 1996.
[33] I.H. Sloan. Analysis of genral quadrature methods for integral equations with con-

tinuous or discontinuous terms. Journal of the Institute of Mathematics and its Ap-
plications, 26:175–186, 1980.

[34] V.I. Krylov. Approximate Calculation of Integrals. Macmillan Company, New York,
1962.

[35] P.K. Kythe and P. Puri. Computational Methods for Linear Integral Equations.
Birkhuser, Boston, 2002.