Computation of Volterra and Fredholm Integro-Differential Nonlinear
Boundary Value Problems by a Modified Regula Falsi-Bisection Shooting
Approach
OKEY OSELOKA ONYEJEKWE
Robnello Unit for Continuum Mechanics and Nonlinear Dynamics, Umuagu Oshimili South,
Asaba, Delta State, NIGERIA
Abstract: A modified Regula Falsi shooting method approach is deployed to solve a set of Volterra and
Fredholm integro- nonlinear boundary value problems. Efforts to solve this class of problems with
traditional shooting methods have generally failed and the outcomes from most domain based approaches
are often plagued by ill conditioning. Our modification is based on an exponential series technique
embedded in a shooting bracketing method. For the purposes of validation, we initially solved problems
with known closed form solutions before considering those that do not come with this property but are
singular, genuinely nonlinear, and are of practical interest. Although in most of these tests, convergence
was found to be super linear, the errors decreased monotonically after few iterations. This suggests that the
method is robust and can be trusted to yield faithful results and by far surpasses in simplicity various other
techniques that have been applied to solve similar problems. In order to buttress the effectiveness and utility
of this approach, we display both the graphical and error analysis outcomes. And for each test case, the
method can be seen to demonstrate the closeness of numerically generated results to the analytical
solutions.
Keywords: Integro differential equations; boundary value problems; Fredholm; Volterra; Regula Falsi;
exponential series; shooting method.
Received: October 26, 2022. Revised: September 14, 2023. Accepted: October 17, 2023. Published: November 2, 2023.
1. Introduction
An integro-differential equation (IDE)
contains both the derivative as well as the
integral of the unknown function. Equations
of this kind represent situations where there
is a strong link between the current behavior
of the system together with its history at any
starting point. They are applicable to several
fields including travelling waves, spread of
diseases, population spread, viscoelastic
flows, particle dynamics, elasticity,
mechanics of continuous media etc. Since
only a few of them are amenable to closed
form solutions, the majority is solved
numerically.
Research in this field has been ongoing.
Prominent among these are spectral
techniques. Spectral methods adopt a basis
function to construct an approximation to the
solution. An automatic Chebyshev spectral
approximation has been applied with great
success by Driscoll et al. [1] to solve linear
and nonlinear IDEs. Several of these
applications can be found in [2, 3 4].
Adomian decomposition method yields a
series expansion of the solution sought and
has been extended by Hadizadeh and
Maleknejad [5] to solve Volterra integral and
higher order integral equations. Despite the
fact that the method requires considerable
symbolic manipulations of its components,
sometimes only a few terms are needed to
arrive at a satisfactory solution [6].
The sinc-derivative collocation numerical
technique has been applied by Abdella and
Ross [7] to arrive at a solution of a general
International Journal on Applied Physics and Engineering
DOI: 10.37394/232030.2023.2.15
Okey Oseloka Onyejekwe
E-ISSN: 2945-0489
152
Volume 2, 2023
second order IDE They deployed suitable
transformations to reduce the non-
homogeneous conditions to their
homogeneous analogs. In addition they
developed special .interpolation techniques
to handle the derivatives of unknown
variables. Their approach was found to be
very competitive when the numerical results
were compared with literature results. The
method has been found to be very accurate
and sometimes displays exponentially
decaying errors[8] in addition to effectively
control and damp out errors oftentimes
associated with numerical differentiation
[9,10].
Jaiswel [11] introduced a Newton-Steffensen
method to effectively handle the derivative
term in the Newton’s method. He
successively applied his derivative free
technique to the solution of nonlinear
algebraic equations. Overall it can be
surmised that most of the research work done
in this field has been geared towards the
formulation of more elegant numerical
approaches, to arrive at a higher order
convergence , reduction in the level of
complexity, differentiability enhancement
and numerical reliability. Efficient ways of
handling variations of derivative based
Newton-Raphson approach are ubiquitously
explored [12,13,14]. In addition,
considerable has also been done on the
shooting method bracketing-based
algorithms. Noticeable among these are the
work done by Brent [15] Dekker[16] and
Wu[17].
So far, there is no single optimal shooting
method for handling nonlinear problems.
None is self sufficient ; each has its own
strength and weaknesses. Any one method
may outperform the other in a particular
problem and may produce poor results in a
different dataset. Hybrid techniques
developed to counter some of these problems
exist in literature. Sabharwal [18] developed
a new hybrid Newton-Raphson algorithm and
a blend of bisection, Regula Falsi and
Newton-Raphson techniques. His method
exploited the advantages of the three
algorithms in each iteration to arrive at a
better approximation of the sought result.
Moreover he showed that the complexity of
the method is far less than that of any of the
three component algorithms. Sabharwal [18]
considered the number of iterations and not
the running time as a valid criterion for
assessing. Overall performance. Badr et al.
[19] adopted a different approach. They
positioned their stand on the fact that there
are some methods that take a relatively small
number of iterations to solve a problem;
despite that, the execution time may be quite
large and vise versa. Hence both factors were
considered as realistic assessments of a
numerical algorithm. Their novel blended
hybrid technique incorporated the advantages
of the trisection and Regula Falsi methods
and was claimed to outperform that of
Sabharwal [18].
We hasten to comment that quite a good
number of these methods are aimed at
nonlinear algebraic equations and contain
attempts to attenuate handle numerical
challenges by extensive algebraic
manipulations. However considerable
computational efforts are needed to translate
and adapt some of these advantages to the
solution of nonlinear integral boundary value
problems (BVPs). Having noted the above,
the foremost objective of the work reported
herein is to adopt a blended bracketing
shooting method approach to provide
accurate results for a nonlinear integro-
differential equation. In line with this, two
types of shooting bracketing methods are
considered, namely: the Regula Falsi and the
bisection techniques. Each of them has its
own peculiar characteristic that provides
valuable insights into their strengths and
weaknesses. The bisection method, depends
solely on the continuity of the function being
investigated, as a consequence, its value at
the midpoint should lie between those at the
endpoints.However, if the function is
International Journal on Applied Physics and Engineering
DOI: 10.37394/232030.2023.2.15
Okey Oseloka Onyejekwe
E-ISSN: 2945-0489
153
Volume 2, 2023
differentiable, the Regula Falsi takes
advantage of this by using linear
interpolation (or extrapolation) to provide a
better estimate of the root within the range.
Both methods known as closed or bracketing
techniques require two guesses that bracket
or contain the sought root. With each
iteration, the range is further refined to
facilitate the location of the root.
The manner of approaching this quest varies
with the method. Whereas with the bisection
method, the endpoints of the function should
contain the root ( right in the middle of the
latest refinement); with the Regular Falsi, this
guaranteed method of locating the root is
compromised in exchange for a better guess
for the location of the root. It does this by
finding the secant line between two pints and
then locating the root at the point where the
line cuts the x-axis. Since the location of the
root is no longer situated at half the interval
space
, it will now lie anywhere between 0
and 1
01

. This is a very useful
information in a hybrid algorithm
formulation that includes a switching
mechanism. Keeping a tab on the value of
will yield a clear directive as to which of the
two methods to use. Several modifications of
the Regula Falsi have been reported in
literature consisting of several variants of the
classical approach [20-23]. We must
comment that relevant literature in this field
contains numerous applications to nonlinear
algebraic problems; and similar solutions of
nonlinear boundary value problems is often
lacking, it is this gap we aim to lessen.
Traditionally, problems with varying
degrees of nonlinearity and even singularity
have been handled with domain based
techniques such as finite difference, finite
element, boundary element and several
variants of collocation techniques. Their
limits and advantages have also been widely
discussed. We do not attempt to resort to such
methods. In this work, we introduce an
algorithm , that is a blend of two shooting
method bracketing techniques. It a does not
rely on the differentiability of a function or a
predictor-corrector approach. Instead the
iterative procedure is more heuristic and
proceeds in a manner that is in consonance
with the optimal ratio at which the iteration
interval is monotonically increasing or
decreasing.
.
2.1 Numerical Formulation
The hybrid technique reported here is a
considerable modification of Thota [24] where a
root-finding method was used to detect non-zero
real roots of transcendental nonlinear algebraic
equations. We deal with a family of techniques
which constrain the root within an interval. Let
us assume that [a,b] is such an interval, that
contains such a root, then an approximation of
such a root with the Regula Falsi method is given
as :
.1
a f b bf a
f b f a
Linear interpolation yields :
In fact if
.0f f a
, then the interval
,a
must contain a root. After each iteration,
more information is needed to refine the range of
the root identification. This can be said of both
the secant and the Regula Falsi methods. The
major difference between the two is that like the
bisection method , the Regula Falsi ensures that
the root is bracketed within range of the
function; whereas the secant method obtains a
better approximation by modifying the previous
estimate of the root.
Without any loss in generality, the linear
interpolation for both methods is written as :
International Journal on Applied Physics and Engineering
DOI: 10.37394/232030.2023.2.15
Okey Oseloka Onyejekwe
E-ISSN: 2945-0489
154
Volume 2, 2023
1
1
1
3
nn
n n n
nn
zz
z z f z
f z f z

Equation (3) is actually the first two terms of the
expansion:
1
1
1
2
1
1
3
1
2
1
1
2
1.... 4
6
n n n
nn
nn
n n n
n n n
n n n
n n n
z z f z
zz
f z f z
z z f z
z f z f z
z z f z
z f z f z








Adopting the standard form of exponential
expansion, equation (4) can be recast to read:
1
1
2
2
1
1
1
3
3
1
1
1
1
11 5
2
11 .....
6
n n n
n n n
n n n
nn
n n n
n n n
n n n
z z f z
z f z f z
z z f z
zz z f z f z
z z f z
z f z f z































hence:
1
1
1
exp 6
n n n
nn
n n n
z z f z
zz z f z f z




where the order of approximation order of
equation (6) is :
4
1
3
1
1
24
n n n
nnn
z z f z
zf z z









Essentially we are trying to achieve three things;
namely (i) develop a very simple numerical
algorithm to locate a root (ii) bracket the root and
(iii) speed up the whole process. The hybrid
method discussed herein should be able to
identify the root, but in reality, due to its heavy
bias towards the Regula Falsi method, it will
exchange the more reliable interval halving
bisection approach for a better guess of the root
location. To optimize performance, the interval
ratio
is monitored as iteration proceeds., and
whenever
0.5
, a switch is made to the
bisection method [25].
2.2 The Algorithm
Step 1: Supply the two guess values, error
tolerance and other parameters for computation.
Step 2: Set up an iteration loop and apply
equation 6
Step 3 : Compute the value of
for each
iteration
Step 4: Use an ‘if’ statement to determine the
appropriate shooting technique to deploy; based
on the value of
Step 5 : Terminate or continue with the
computation based on the magnitude of the error
criterion
International Journal on Applied Physics and Engineering
DOI: 10.37394/232030.2023.2.15
Okey Oseloka Onyejekwe
E-ISSN: 2945-0489
155
Volume 2, 2023
Fig. 1a : Algorithm Fowchart
Supply problem parameters:
12
,,xx
Solve Equation (6)
Is Interval greater/less than 0.5
Yes: Do
Bisection
No: Continue
with Equation
(6)
Is Error
satisfied
If No; Go to
step 1
If Yes; STOP
Is Error
Satisfied
If No; Go to
Step 1
If Yes; STOP
International Journal on Applied Physics and Engineering
DOI: 10.37394/232030.2023.2.15
Okey Oseloka Onyejekwe
E-ISSN: 2945-0489
156
Volume 2, 2023
3.Numerical Experiments :
For simulation purposes, we use a termination
error tolerance of
6
10
together with an
upper bound iteration of 50 . Examples of
Integro-differential equations with analytic
solutions where available are chosen from
literature. .
3.1.Example 1
The first example from Filipov et al.[26] is of a
Volterra integral type.
2
1
'' 3 ' 1 3 / , 1, 2
, 1 1, 2 4 7
t
u t u t G t u t
G t u s ds u u a
The analytic solution is given as :
27u t t b
Fig. 1b shows excellent comparison between
numerical and analytic results. A monotonically
decreasing error profile is displayed in Fig. 1c.
The convergence of the guessed slopes at the left
hand side (LHS) boundary as the shooting
method algorithm proceeds is displayed in Fig. 1d
Fig.1b Comparison of numerical and analytic
solutions
Fig. 1c : Guessed slopes per iteration at the
LHS boundary
Fig. 1d Absolute error profile
3.2.Example 2
Problem 2 is taken from Abdella and Ross [7] and
is of the Fredholm type
1
1
'' 2 cos 1,1
1 1, 1 1 8
t
u u te t udt x t
u u a
The solution is given as
International Journal on Applied Physics and Engineering
DOI: 10.37394/232030.2023.2.15
Okey Oseloka Onyejekwe
E-ISSN: 2945-0489
157
Volume 2, 2023
22
1 2 3 cos 0.5 8
tt
u t c e c e c t t b
where
12
3
0.4206057265; 0.3545613032;
0.2662525683
cc
c

Figures 1a and 1b show the profiles of both the
numerical solution and the absolute error
between the numerical and the analytic solution.
The magnitude of the errors displayed confirms
that the numerical results are accurate .
Fig. 2a : Numerical solution profile
Fig. 2b : Absolute error profile
3.3.Example 3
This is a Fredholm type Integro-differential
equation given as:
1
2
0
2
'' 1 (0) 1, (1) 0.5 9
1
u
u t t t udt u u a
t
The analytic solution is
1 1 9u t t b
International Journal on Applied Physics and Engineering
DOI: 10.37394/232030.2023.2.15
Okey Oseloka Onyejekwe
E-ISSN: 2945-0489
158
Volume 2, 2023
Fig. 3a Comparison of numerical and analytic
solution
Fig. 3b : Convergence of guessed slopes per
iteration
Fig. 3c. Absolute error profile
Fig. 3d : Scalar profile and gradient
3.4.Example 4
This example is a Volterra type Integro-
differential equation taken from Abdella and
Ross [7].
2
2
0
3
12
2
2
1
'' ( ) 10
1
0 1, 1 2 10
1
2 2 10
4
11
t
u
t
u te tsu ds f t a
u
u u b
t
f t te t c
t


The results displayed are quite reliable. The
boundary conditions are not only satisfied in
figure 4a, the reliability of the algorithm is also
confirmed by the closeness of analytic and
numerical results. Fig. 4b is a candid picture of
what obtains at the left end boundary as the
shooting iteration proceeds. Essentially,
convergence starts at the fourth iteration as the
difference between the ‘hits and target (specified
Dirichlet right end boundary condition)’ gets
smaller per iteration. Fig. 4c is the absolute error
profile.. A salient observation shows that the
magnitude of the errors increases monotonically
as we move away from the left boundary, but
starts decreasing as we move closer to the right
International Journal on Applied Physics and Engineering
DOI: 10.37394/232030.2023.2.15
Okey Oseloka Onyejekwe
E-ISSN: 2945-0489
159
Volume 2, 2023
boundary as the boundary conditions are satisfied
at both ends. This must have been caused by the
restrictive nonlinear forcing function that comes
with this problem. Figure 4d shows the
ascending solution profile from the left to the
right boundaries. as well as the accompanying
gradient.. This serves to validate the reliability of
the results as well as the physics of the problem.
Fig. 4a : Profiles of numerical and analytic results
Fig. 4b : Convergence of guessed slopes at left end
boundary
Fig. 4c: Absolute error profile
Fig. 4d: Scalar and gradient profiles
3.5.Example 5
Here we consider a Fredholm integro-differential
equation with a Neumann boundary condition at
the left-end of the problem interval definition as
well as a singularity. We have purposely included
this problem to show that the hybrid formulation
International Journal on Applied Physics and Engineering
DOI: 10.37394/232030.2023.2.15
Okey Oseloka Onyejekwe
E-ISSN: 2945-0489
160
Volume 2, 2023
is capable of dealing with different boundary
conditions as well as singularity.
1
0
2'
'' 2 11
1
uu
u sds a
tu
' 0 0, 1 1 11u u b
The analytic solution to this problem is unknown.
Fig. 5a shows the scalar solution profile. The
boundary conditions at both ends are satisfied. In
addition Fig. 5b illustrates the convergence
between ‘hits’ and ‘target’ as computation
proceeds. Based on these results, the method has
been able to cope with additional challenges
contributed by this type of problem.
Fig. 5a : Scalar profile
Fig. 5b: Error between strikes and target at
the right end boundary
4 Conclusion
This article describes a blended numerical
shooting algorithm which handles integro-
differential equations straightforwardly. The
effort to arrive at a simple and robust algorithm is
deliberate and a key motivation. Results
obtained herein confirm that the formulation is
robust and the overall convergence is satisfactory
as the iteration proceeds. The maximum iteration
steps for all the problems handled herein is less
than 10. Beyond being easy to implement, the
method requires only continuity and neither
derivatives nor a predictor-corrector component.
More importantly because of its bracketing
approach, it is better guaranteed to converge to a
root. We believe that further numerical
experimentation included theoretical analysis are
needed to arrive at firmer conclusions.
International Journal on Applied Physics and Engineering
DOI: 10.37394/232030.2023.2.15
Okey Oseloka Onyejekwe
E-ISSN: 2945-0489
161
Volume 2, 2023
References
[1] T.A. Driscoll, F. Bornemann, L.N.
Trfethen , The chebop system for
automatic solution of differential
equations, BIT Numerical Mathematics
Vol. 48 2008 pp. 701-723
[2] C-I Gheorghiu Chebfun solutions to a
class of 1D singular and nonlinear
boundary value problems , Computation,
Vol. 10, 2022 116.
https://doi.org/10.3390/computation100
70116
[3] J.P. Boyd and C.I. Georghiu All roots
spectral methods : Constraints, floating
point arithmetic and root exclusion :
Appld. Math. Letters Vol. 67 201 pp. 28-
32
[4] J.P. Boyd, Chebyshev polynomial
expansions for simultaneous
approximations of two branches of a
function with application to the one-
dimensional Bratu equation Appld. Math.
Comput. Vol. 143 2003 pp. 189-200
[5] M. Hadizadeh, K. Maleknejad, The
numerical analysis of Adomian
decomposition method for some
nonlinear turbulent diffusion problems,
Nonlinear Studies Vol. 6 1999 pp. 85-89
[6] G. Adomian, Solving frontier problems
of physics : The decomposition method
, Kluwer Acad. Publishers , Boston,
1994
[7] K. Abdella and G. Ross, Solving Integro-
differential boundary value problems
using sinc-derivative collocation,
Mathematics Vol. 8 2020 1637,
doi:10.3390/math8091637
[8] K. Abdella, Solving differential
equations using sinc-collocation
methods with derivative interpolations, J.
Comput. Mthds Sci. Engnr. Vol. 15 2015
pp. 305-315
[9] K. Parand, A. Pirkhedri, A sinc-
collocation metod for solving
astrophysics equation, New Astron. Vol.
15 2010 pp. 533-537
[10] S. Yganeh, Y. Ordokhani and A.
Saadatmandi, A sinc-collocation
method for second order boundary value
problems of nonlinear integral
differential equations, J. Inf. Comput.
Sci. Vo. 7 2012 pp. 151-160
[11] J.P. Jaiswal Third-order derivative free
method for solving nonlinear equations ,
Universal Jnl. Appld. Math. Vol. 1 2013
pp. 131-135
[12] R.P. Brent, Algorithms for minimization
without derivatives , Prentice Hall
Englewood Cliffs NJ 1973
[13] M. Frontini, E. Sormani, Modified
Newton’s method with third order
convergence and multiple roots . J.
Comput. Appld. Math. Vol. 156 2003
pp. 345-354
[14] X.Y. Wu, J.L. Xa, R. Shao, Quadratically
convergent multiple roots finding
method witout derivatives Comput.
Math. Appld. Vol. 142 2001 pp. 115-119
[15] R.P. Brent Algorithms for minimization
without derivatives , Prentice –Hall,
Englewood Cliffs Nj. 1973
[16] T.J. Dekker Finding a zero by means of
successive linear interpolation in : B.
Dejon, P. Henrici (Eds.) Constructive
aspects of the fundamental theorem of
Algebra , Wiley-Interscience , New York
1969
[17] X. Wu, Improved Muller method and
bisection method with global and
asymptotic superlinear convergence of
both point and interval for solving
nonlinear equations , Appld. Math.
Comput. Vol. 166 2015 pp. 299-311
International Journal on Applied Physics and Engineering
DOI: 10.37394/232030.2023.2.15
Okey Oseloka Onyejekwe
E-ISSN: 2945-0489
162
Volume 2, 2023
[18] C.L. Sabharwal An Iterative hybrid
algorithm for roots of linear equations,
Eng. Vol. 2 2021 pp. 80-98
[19] E. Badr, S. Almotairi and A. El. Ghamry,
A comparative study amomg new hybrid
root finding algorithms and traditional
methods , Mathematics Vol . 9 2021
1306.
https:’’doi.org/10.3390/math9111306
[20] R.G. Gottlieb and B.F. Thompson,
Bisected direct quadratic regula falsi,
Appld. Math. Sci. Vol. 4 2010 pp. 709-
718
[21] M. Dowel and P. Jarrat A modified
regular falsi method for computing the
root of an equation, BIT, Vol. 11 1971 pp.
168-174
[22] M. A. Hafiz, A new combined bracketing
method for solving nonlinear equations
Journal of Mathematical and
Computational Science, Vol. 1 2016 pp.
44-47
[23] A. Suhaldonik, Combined bracketing
methods for solving nonlinear equations,
Appld. Math. [24]Letters Vol. 25 2012
pp. 1755-1760
[24] S. Thota, A new root-finding algorithm
using exponential series URAE Math.
Jnl. Vol. 5 2019 pp. 83-90
[25] A. Suhaldonik, Combined bracketing
methods for solving nonlinear equations
, Appld. Math. Lettrs. Vo. 25 2012 pp.
1755-1760
[26] S. M. Filipov, I.D. Gospodinov J.
Angelova, Solving two-point boundary
value problems for Integro-differential
equations using the simple shooting-
projection method . Springer
International Publishing Switzerland I.
Dimov et al. (Eds.): NMA 2014 LNCS
8962, pp. 169-177 2015
International Journal on Applied Physics and Engineering
DOI: 10.37394/232030.2023.2.15
Okey Oseloka Onyejekwe
E-ISSN: 2945-0489
163
Volume 2, 2023
Contribution of Individual Authors to the
Creation of a Scientific Article (Ghostwriting
Policy)
The author contributed in the present research, at all
stages from the formulation of the problem to the
final findings and solution.
Sources of Funding for Research Presented in a
Scientific Article or Scientific Article Itself
No funding was received for conducting this study.
Conflict of Interest
The author has no conflict of interest to declare that
is relevant to the content of this article.
Creative Commons Attribution License 4.0
(Attribution 4.0 International, CC BY 4.0)
This article is published under the terms of the
Creative Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en
_US