††thanks:

# Normalized factorial moments of spatial distributions of particles in high
multiplicity events: A Toy model study

Sheetal Sharma, Salman K. Malik, Zarina Banoo  Ramni Gupta
<EMAIL_ADDRESS>Department of Physics, University of Jammu, Jammu,
India

###### Abstract

Scaling behavior of normalized factorial moments ($F_{q}$) of spatial
distributions of the particles comprising a system may be studied to probe and
to determine its characteristics. In heavy-ion collisions at ultra-
relativistic energies, a strongly interacting complex system of quarks and
gluons is created. The nature of the system created and multi particle
production mechanism in these collisions is predicted to be revealed by the
study of normalized factorial moments ($F_{{\rm{q}}}$) as function of various
parameters. In this work, observations from the Toy model study of the scaling
behaviour of $F_{{\rm{q}}}$ moments, resilience of these moments to detector
efficiencies and sensitivity towards fluctuations in the system will be
presented.

##  I Introduction

The Large Hadron Collider (LHC) at CERN, Switzerland [1], and the Relativistic
Heavy-Ion Collider (RHIC) at BNL, USA [2], serve the purpose of studying the
deconfined state of strongly interacting particles, the Quark-Gluon Plasma
(QGP) and unravelling its properties by colliding nuclei at ultra-relativistic
energies [3, 4]. To obtain the exotic state of matter QGP in the laboratory,
heavy nuclei are collided at high energies such that the energy density or
temperature is high enough so as to melt hadrons [5]. As the medium so created
cools down, a transition occurs back from the QGP state to the hadronic state.
One of the aims of these collider experiments is to probe phase diagram of
this strongly interacting matter and its properties [6, 7, 8].

In heavy-ion collisions, behavior of the system near critical point, which is
a point in the phase diagram, where the correlation length of the system
becomes infinite and the system becomes scale-invariant, is still not known
[9]. The study of critical behavior is vital to have a deeper understanding of
the properties of any system and knowledge about the nature of the phase
transition between different phases. At this point, the system exhibits large
fluctuations in various observables, such as particle multiplicity, transverse
momentum, pion to kaon ratio etc. [10, 11]. A study of the fluctuations in
heavy-ion collisions thus provides a powerful means to understand various
characteristics of the system.

The charged particle multiplicity in heavy-ion collisions is predicted to be a
function of the energy and temperature conditions of the collision.
Multiplicity fluctuations is an important observable of these experiments,
that promises to provide information on the characteristics of the state of
the system created in the collisions. In Ref.[12], it has been suggested that
the normalized factorial moments (NFMs) of the charged particle multiplicity
distributions be used to study multiplicity fluctuations from bin-to-bin and
event-to-event, thereby revealing some of the important features of the matter
at very high energy density[13].

A power-law behaviour of the normalized factorial moments of the particle
density fluctuations in spatial or momentum space, with increasing number of
bins, is termed as intermittency [14, 15, 16]. This analysis technique has
been performed for various systems at low energies in search for QGP and to
understand the quark-hadron phase-transition [13]. However, investigations of
heavy-ion collision data using this tool have led to no definitive conclusions
on the critical point, order of phase transition or nature of multiplicity
fluctuations [16] because of low bin multiplicities. The availability of data
with high charged particle density per bin has generated a renewed interest,
such as at STAR experiment at RHIC [17], ALICE at CERN [18], to use this
methodology to comprehend the multiparticle production processes.

In this work, Toy model events are studied to have baseline behavior of NFM in
the high multiplicity scenerio. In Section II, a brief introduction to the
event generation using the Toy model is given. This follows the details on
methodology of analysis in Section III. Observations and results of the
analysis are discussed in Section IV. Finally, Section V summarizes the whole
work.

##  II Event generation

In the field of heavy-ion collisions, event generators are extensively
employed by experimentalists to simulate experimental conditions and to get
the expected outcome from experimental conditions as per known physics. These
event generators are also widely used by phenomenologist to develop new models
and test predictions from existing models and theories. As multiparticle
production in heavy-ion collisions is still not completely understood, a
number of methodologies are used to investigate particle production mechanism
and to have insight into the quark-hadron phase transition. In this work, the
baseline behavior of the normalized factorial moments in the contours of
methodology similar to the intermittency analysis, proposed in [12], has been
investigated for its suitability to analyze heavy-ion collision data.

Intermittency analysis, first proposed in 1986 [14] to investigate
fluctuations in pseudorapidity distributions of some cosmic ray events, has
been suggested [19, 12, 20] to be performed for high multiplicity data
recorded by the detectors in the recent colliders, such as at the LHC, to
understand the multiparticle production and quark-hadron phase transition. The
ALICE experimental set up at LHC has central barrel detector that includes
TPC, ITS and TOF detectors. These together measure the charged particles
within a common angular phase space of pseudorapidity ($|\eta|\\!\leq$ 0.8)
and full azimuth ($0\\!\leq\\!\varphi\\!\leq\\!2\pi$) [21, 22].

Toy model events, with multiplicity distributions similar to the one recorded
for midrapidity region by ALICE detector for low $p_{\rm{T}}$ tracks
($<\\!2.0$ GeV/c), are generated. Four event samples with multiplicity
distributions, as shown in Figure 1, having mean 400, 900, 1300 and 1900 are
simulated using gRandom3() function of the TRandom class of ROOT [23] analysis
framework using system clock as seed. Particles are generated in the spatial
phase space $(\eta,\varphi)$ with uniform distribution, having $|\eta|\,\leq$
0.8 and $0\\!\leq\\!\varphi\\!\leq\\!2\pi$, in a way that there are no bin-to-
bin (spatial) fluctuations.

Figure 1: Multiplicity distributions for the four sets of Toy model events.

##  III Analysis methodology

Lattice QCD predicts large density fluctuations in a system at critical point
undergoing phase transition. It is argued [12] that behavior of NFM as
function of number of bins is a good measure of fluctuations in the spatial
configurations of produced particles in the heavy-ion collisions. This
methodology, similar to the intermittency analysis as proposed earlier [14,
24, 15, 25] has been tested for some of the models at LHC energies in the
recent work [26, 27] A two dimensional (2D) analysis is performed for the Toy
model events wherein for an event the kinematic phase space is partitioned
into a square lattice with $M_{{\rm\eta}}$ and $M_{\rm{\varphi}}$ numbers of
bins along $\eta$ and $\varphi$, respectively, as is depicted in Fig 2. With
$M_{{\rm{\eta}}}\\!=\\!M_{{\rm\varphi}}$ there are $M^{2}$ = $M$(say) number
of bins, such that $q^{\rm{th}}$ order vertically averaged horizontal moments
(NFM), for a total of $N$ events, is defined as

$F_{\rm{q}}(M)=\frac{\frac{1}{N}\displaystyle\sum_{e=1}^{N}\frac{1}{M}\sum_{i=1}^{M}f_{q}^{e}(n_{\rm{ie}})}{\left(\frac{1}{N}\displaystyle\sum_{e=1}^{N}\frac{1}{M}\sum_{i=1}^{M}f_{1}^{e}(n_{\rm{ie}})\right)^{q}}$
(1)

Figure 2: Graphic illustration of two dimensional ($\eta,\varphi$) phase space
partitioned in M$\,\times\,$M bins and mapping of tracks of an event onto this
phase space (grid).

$\rm{where}\,\,f_{\rm{q}}^{e}(n_{\rm{ie}})=\Pi_{j=0}^{q-1}(n_{\rm{ie}}-j),$
(2)

is the factorial moment and $n_{\rm{ie}}$ is the bin multiplicity in the
$i^{\rm{th}}$ bin of $e^{\rm{th}}$ event. $q$ takes positive integer values
$\geq 2$. By definition these normalized factorial moments ($F_{\rm{q}}$(M))
filter out statistical fluctuations. If fluctuations in the spatial
distributions of the particles is Poissonian, then $F_{q}$(M) = 1 [14]. A
power-law behaviour of $F_{q}$ for $q\geq 2$ as function of M

$F_{\rm{q}}(M)\propto M^{\phi_{\rm{q}}}$ (3)

then defines intermittency (also termed as M-scaling), where $\phi_{\rm{q}}>0$
is known as intermittency index. Another scaling where $F_{\rm{q}}(M)$ for q
$>$ 2 having power-law dependence on the second-order factorial moment
($F_{\rm{2}}(\rm{M})$) is proposed to exist [19, 28] such that

$F_{\rm{q}}\propto F_{\rm{2}}^{\beta_{\rm{q}}},$ (4)

Figure 3: (Left) log-log plot of $F_{\rm{q}}(M)$ vs $M^{2}$ (M-scaling) for q
= 2, 3, 4 and 5; (Right) log-log plot of $F_{\rm{q}}(M)$ vs $F_{\rm{2}}$
(F–scaling) for q = 3, 4 and 5.

where $\beta_{q}$ is the exponent which depends on bin resolution and number
of bins $M$. This scaling is termed as F-scaling. Relating $\beta_{q}$ to
order of the moment($q$), a scaling exponent ($\nu$) can be obtained from the
relation

$\beta_{\rm{q}}=(q-1)^{\nu}.$ (5)

$\nu$ is called the scaling exponent. The scaling behavior and sensitivity of
NFM of Toy model events is investigated. In the backdrop of various detector
effects which distort the original signal, the efficiency correction method
that is required to be applied to the experimental data for the extraction of
true behaviour has also been studied. In this work, $M$ is varied from 4 to
100 and $q$ is taken as 2, 3, 4, and 5.

##  IV Observations and Discussion

###  IV.1 Scaling of normalized factorial moments (NFM)

The analysis has been performed on event-by-event basis where each event is
mapped on to two dimensional phase space that is partitioned into $M^{2}$
bins. The number of particles that go into a bin is the bin multiplicity and
using Eq.(2) the event factorial moments are determined for q = 2, 3, 4, and
5. The process is repeated for all the events for a $M$ and thus to obtain
Normalized factorial moments for the whole event sample, as defined in Eq. 1
($F_{q=2,3,4,5}(M)$). To see the dependence of $F_{q}$ on $M^{2}$, a log-log
plot between the two is drawn as is shown in Fig. 3(a) for events with average
multiplicity 900. It is observed that $F_{\rm{q}}$(M)$>\\!1$ for the second
order NFM though the deviation is very small. As $q$ increase deviation of
$F_{q}$ from 1 also increases. However $F_{q}$ shows no dependence on M (no
power-law as in Eq. 3) that is as M increases there is negligible change in
the $F_{\rm{q}}$ for all $q{\rm{{}^{\prime}s}}$ and thus no intermittency in
the generated events. No significant rise in $\ln{F_{\rm{q}}}$ with
$\ln{M^{2}}$ seen for the Toy model indicates absence of self-similar
modelling of particle generation. The statistical uncertainties on the markers
are estimated using the sub-sampling method. For F-scaling as is defined in
Eq. 4, a $\ln\\!F_{q}(M)$ versus $\ln\\!F_{2}(M)$ plot is shown in Fig. 3 (b).
For $q>2$, $F_{q}$ shows linear dependence on $F_{\rm{2}}$ and line fit to
these graphs give slopes ($\beta_{q}$). A dimensionless scaling exponent is be
obtained (Eq. 5) which quantitatively characterizes the geometric spatial
configurations of the Toy model events (discussed in section IV.3).

###  IV.2 Sensitivity of NFM to gauge fluctuations

Independence of $F_{\rm{q}}$ on the number of bins ($M^{2}$) for the Toy model
events can be because of two possible reasons. One of the reasons is that the
observable $F_{\rm{q}}(M)$ are not sensitive to multiplicity fluctuations from
bin-to-bin and the second possible reason which can lead to this behaviour is
in case there are no fluctuations in the events that vary from bin-to-bin. The
events have been simulated (section II) for uniform spatial distributions and
thus it is important to see that in case there are fluctuations,
$F_{\rm{q}}(M)$ moments should be sensitive to the spatial fluctuations. To
investigate the sensitivity of the observable to gauge the density
fluctuations in the spatial configurations, fluctuations are added with hands
in each events and the scaling behaviour of normalized factorial moments is
studied.

The event sample that results in after the fluctuations are added with hand in
each Toy model event is termed as modified Toy model events. For adding
fluctuations with hand in an event, number of tracks equal to 5% of the
multiplicity of an event, are added in some region of the phase space of the
event. An equal number of tracks are removed from the rest of the region to
keep the multiplicity distribution of the modified Toy model events same as
that of the original Toy model events. This makes some of the bins more
populated than the rest of the bins thereby including multiplicity spikes in
the spatial configurations which replicates a system having critical
fluctuations. Fig.4, depicts ($\eta$, $\varphi$) lego plot of a Toy model
event, for the phase-space partitioned into 40 bins and the same event after
it is modified with fluctuations added with hand.

Figure 4: (Left plot) Spatial distribution of particles in the ($\eta,\phi$)
phase space of an event. Panel on right shows the spatial distribution of
particles of the same event with 5% tracks added in selected bins while
removing equal number of tracks from rest of the region (modified Toy model
events).

Figure 5:  log-log plot of (a) $F_{\rm{q}}$(M) vs $M^{2}$ (M-scaling) for q =
2, 3, 4 and 5 and (b) $F_{\rm{q}}$(M) vs $F_{\rm{2}}$(M) (F–scaling) for q =
3, 4 and 5

Analysis is performed on the modified Toy model events, and normalized
factorial moments are determined for q = 2, 3, 4 and 5 with M varying from 4
to 100. The M-scaling behaviour is studied where in Fig.5(a) log-log plot
between $F_{\rm{q}}$ and $M^{2}$ is given. A power-law dependence of
$F_{\rm{q}}$ on $M^{2}$ is observed with a linear dependence for high $M$
values. However, in this case, $F_{\rm{q}}$ is observed to be $>>$ 1 and
$F_{\rm{q}}(M)$ deviates from 1 as M increases. A linear dependence of
$F_{\rm{q}}$ on $M$ at higher values indicate the scale invariant
fluctuations. These results show that NFM are a suitable measure to gauge
fluctuations when there are bin-to-bin multiplicity fluctuations. This study,
in addition to what is observed in Fig. 3, makes the NFM calculations and the
study of scaling behaviour with increasing bin resolution an observable that
is suitable for gauging fluctuations in spatial distributions.

###  IV.3 Scaling exponent and its dependence on multiplicity

Scaling exponent ($\nu$) obtained from linear fit to the $\ln\beta_{q}$ vs
$\ln(q-1)$ plot (Eq. 5) provides valuable insight into the underlying physics
of the heavy-ion collision[13]. The scaling exponent for the Toy model events,
with no dynamical fluctuations, is extracted to be $1.603\\!\pm\\!0.016$. For
the modified Toy model events, where some fluctuations are added with hands,
the value of $\nu$ is found to be $0.998\\!\pm\\!0.004$ as is shown in Fig. 6.
Toy model events having no preferred bins for track generation, a little
addition of tracks in small region alters bin multiplicity and gives large
fluctuations in spatial patterns such that $\nu\approx 1.00$, a value obtained
for the two dimensional Ising model with large fluctuations [29, 30, 13].

To study the dependence of scaling exponent value on the mean multiplicity of
the events, Toy events with different mean values are generated. Analysis is
performed on these sets of events and scaling behavior is studied for each
case as in Fig. 3. Scaling exponent is calculated and observed to have no
dependence on average multiplicity as is shown in Fig. 7. Scaling exponent
values predicted from models and theory are also shown with lines.

Figure 6:  Scaling exponent ($\nu$) obtained from linear fit to
$ln\beta_{\rm{q}}$ vs $ln(q-1)$ plot for the Toy model events and modified Toy
model events. Figure 7: Scaling exponent as function of average
multiplicities. Scaling exponent values predicted from the other models are
also shown.

###  IV.4 Resilience of NFM to efficiency corrections

Particles produced in heavy-ion collisions are recorded by detectors that may
have some inefficiencies in detecting all the particles of interest. The
possible reasons may be, not so perfect detector resolution and inefficiencies
in tracking routines. The measurements and calculations of the observable are
thus affected by these inefficiencies and need to be corrected before drawing
any reasonable conclusions. Monte Carlo event generators are used to calculate
overall efficiencies by passing the generated events through the detector
geometry, recording and applying track reconstruction routines as in the real
setup. Detector efficiency is defined as the ratio of the number of tracks
detected by the detector geometry to the number of generated tracks in the
acceptance region. Constrained by different set of conditions, detector
efficiencies in the acceptance region can be binomial or non-binomial. The
inefficiencies may alter true values of the normalized factorial moments and
also their dependence on the resolution or number of bins $M^{2}$. If
$\epsilon_{i}$ denotes efficiency of the $i^{\rm{th}}$ bin in the phase space,
that has been divided in $M^{2}$ bins, then the normalized factorial moments
as defined in Eq.( 1) are corrected for the efficiencies as

$F_{\rm{q}}(M)=\frac{\frac{1}{N}\displaystyle\sum_{e=1}^{N}\frac{1}{M}\sum_{i=1}^{M}\frac{f_{q}^{e.rec}(n_{\rm{ie}})}{{\epsilon_{i}^{q}}}}{\left(\frac{1}{N}\displaystyle\sum_{e=1}^{N}\frac{1}{M}\sum_{i=1}^{M}\frac{f_{1}^{e.rec}(n_{\rm{ie}})}{\epsilon_{i}}\right)^{q}}$
(6)

For the uniform (binomial) efficiencies and the non-uniform (non-binomial)
efficiencies in the acceptance region, how efficiency corrections affect the
NFM has been studied. In general considering $80\%$ overall efficiency in
detecting the particles, to create events with uniform efficiencies, 20
percent of the tracks are uniformly removed from the Toy model events. These
events termed here as the reconstructed-uniform (recU) events, that have 80%
of the tracks that are there in the Toy events. Analysis is performed on these
events which gives NFM (Eq.1) $F_{q}^{recU}$(M). To get the efficiency
corrected NFM, two-dimensional efficiency maps are calculated for all $M$,
where for each $i^{{\rm th}}$ bin in the phase space divided in $M^{2}$ bins,
there is an efficiency value $\epsilon_{i}$. For recU events corrected
normalized factorial moments $F_{q}^{corrU}$(M) are calculated using Eq.6. In
Fig.8 it is observed that $F_{q}$(M) $\approx$ $F_{q}^{recU}$(M) $\approx$
$F_{\rm{q}}^{corrU}$(M). Ratio plots in the lower panel show that, for all
three data sets, the moments have the same value as the ratio is equal to 1.
Thus, for the binomial efficiencies, the true values of NFM are not affected
by these effects or in other words normalized factorial moments are robust
against binomial (uniform) efficiencies.

Events with 80 $\%$ non-uniform efficiencies are obtained from Toy model
events by removing 20 percent of the tracks from the selected phase space
region. This results in a set of events with non-binomial efficiencies, and
the resultant set of events is termed as reconstructed-non uniform (recNU)
events. Analysis is performed on recNU events and normalized factorial moments
($F_{q}^{recNU}$(M)) are calculated using Eq. 1. It is observed that
$F_{q}{\rm(M)}\\!\not\equiv\\!F_{q}^{recNU}$(M) that is if the efficiency
values are not binomial then there the value of the NFM is not the same as
true sample of events. The dependence of NFM on the number of bins will also
be there. Efficiency maps are obtained for these events for each $M$ value and
$i^{th}$ bin. Using efficiency values ($\epsilon_{i}$) from the efficiency
maps for the events with 80$\%$ non-binomial efficiencies, the corrected NFM,
$F_{q}^{corrNU}$(M) are calculated using Eq.6. As shown in Fig.9 it is
observed that $F_{q}{\rm(M)}\\!\equiv\\!F_{q}^{corrNU}$(M), that is applying
efficiency corrections the true value of the normalized factorial moments is
recovered. For non-binomial (non-uniform) detector efficiencies, Eq. 6 where
the efficiency values for the bins are used is resilient enough to obtain the
true values for NFM.

Figure 8: $F_{q=2}$ as function of $\ln M^{2}$ for Toy model events,
reconstructed-Uniform (recU) and efficiency corrected recU events. The lower
panel shows the ratio plots. Figure 9: $F_{q=2}$ as function of $\ln M^{2}$
for Toy model events, reconstructed-non Uniform (recNU) and efficiency
corrected recNU events. The lower panel shows the ratio plots.

##  V Summary

Intermittency analysis is performed on a set of high multiplicity Toy model
events. This methodology to search for the dynamical fluctuations is observed
to be sensitive to the density fluctuations in the multiplicity distributions.
A baseline value of the scaling exponent ($\nu$) from the Toy model events is
found to be $\,1.603\pm 0.016\\!$, a value greater than the value predicted
using second-order phase transition formalism in the framework of the
Ginzburg-Landau theory. The observable, normalized factorial moments
$F_{q}$(M), are observed to be robust against the uniform efficiencies and for
non-uniform efficiencies, the calculations using efficiency corrected NFM
recover the true value of the observable. In addition observations and results
obtained here give a baseline behavior of the normalized factorial moments for
the experimental studies.

##  VI Acknowledgement

This work is partially supported by RUSA 2.0 grant sanctioned in favour of one
of the authors by the Ministry of Education, Government of India.

## References

  * [1] K Aamodt et al.  Elliptic flow of charged particles in Pb-Pb collisions at 2.76 TeV.  Phys. Rev. Lett., 105:252302, 2010.
  * [2] Miklos Gyulassy and Larry McLerran.  New forms of QCD matter discovered at RHIC.  Nucl. Phys. A, 750:30–63, 2005.
  * [3] L. Van Hove.  Two Problems Concerning Hot Hadronic Matter and High-Energy Collisions (Equilibrium Formation, Plasma Deflagration).  Z. Phys. C, 21(1):93–98, 1983.
  * [4] John W. Harris and Berndt Muller.  The Search for the quark - gluon plasma.  Ann. Rev. Nucl. Part. Sci., 46:71–107, 1996.
  * [5] Wit Busza, Krishna Rajagopal, and Wilke van der Schee.  Heavy Ion Collisions: The Big Picture, and the Big Questions.  Ann. Rev. Nucl. Part. Sci., 68:339–376, 2018.
  * [6] Simon Hands.  The Phase diagram of QCD.  Contemp. Phys., 42:209–225, 2001.
  * [7] Yasuyuki Akiba et al.  The Hot QCD White Paper: Exploring the Phases of QCD at RHIC and the LHC.  2 2015.
  * [8] Edward V. Shuryak.  Quantum Chromodynamics and the Theory of Superdense Matter.  Phys. Rept., 61:71–158, 1980.
  * [9] Mikhail A. Stephanov.  QCD phase diagram and the critical point.  Prog. Theor. Phys. Suppl., 153:139–156, 2004.
  * [10] Marcus Bluhm et al.  Dynamics of critical fluctuations: Theory – phenomenology – heavy-ion collisions.  Nucl. Phys. A, 1003:122016, 2020.
  * [11] A. Aduszkiewicz et al.  Multiplicity and transverse momentum fluctuations in inelastic proton–proton interactions at the CERN Super Proton Synchrotron.  Eur. Phys. J. C, 76(11):635, 2016.
  * [12] Rudolph C. Hwa and C. B. Yang.  Local Multiplicity Fluctuations as a Signature of Critical Hadronization at LHC.  Phys. Rev. C, 85:044914, 2012.
  * [13] W. Kittel and E. A. De Wolf.  Soft multihadron dynamics.  2005\.
  * [14] A. Bialas and Robert B. Peschanski.  Moments of Rapidity Distributions as a Measure of Short Range Fluctuations in High-Energy Collisions.  Nucl. Phys. B, 273:703–718, 1986.
  * [15] A. Bialas.  Intermittency ’90.  Nucl. Phys. A, 525:345–360, 1991.
  * [16] E.A. De Wolf et al.  Phys. Rep., 270:1–141, 1996.
  * [17] The STAR Collaboration.  Energy dependence of intermittency for charged hadrons in au+au collisions at rhic, 2023.
  * [18] Sheetal Sharma and Ramni Gupta.  Local multiplicity fluctuations in Pb$-$Pb collisions at $\sqrt{s_{\rm{NN}}}$ = 2.76 TeV with ALICE at the LHC, 2023.
  * [19] Rudolph C. Hwa and M. T. Nazirov.  Intermittency in second order phase transition.  Phys. Rev. Lett., 69:741–744, 1992.
  * [20] Rudolph C. Hwa.  Recognizing Critical Behavior amidst Minijets at the Large Hadron Collider.  Adv. High Energy Phys., 2015:526908, 2015.
  * [21] K. Aamodt et al.  The ALICE experiment at the CERN LHC.  JINST, 3:S08002, 2008.
  * [22] Jaroslav Adam et al.  Pseudorapidity dependence of the anisotropic flow of charged particles in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV.  Phys. Lett. B, 762:376–388, 2016.
  * [23] R. Brun and F. Rademakers.  ROOT: An object oriented data analysis framework.  Nucl. Instrum. Meth. A, 389:81–86, 1997.
  * [24] A. Bialas and Robert B. Peschanski.  Intermittency in Multiparticle Production at High-Energy.  Nucl. Phys. B, 308:857–867, 1988.
  * [25] A. Bialas and K. Zalewski.  Phase Structure of Selfsimilar Multiparticle Systems and Experimental Determination of Intermittency Parameters.  Phys. Lett. B, 238:413–416, 1990.
  * [26] Ramni Gupta and Salman Khurshid Malik.  Intermittency study of charged particles generated in Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}\text{= 2.76 TeV}$ using EPOS3.  Adv. High Energy Phys., 2020:5073042, 2020.
  * [27] Sheetal Sharma and Ramni Gupta.  Intermittency Analysis of Toy Monte Carlo Events.  SciPost Phys. Proc., 10:024, 2022.
  * [28] Wolfgang Ochs and Jacek Wosiek.  Intermittency and Jets.  Phys. Lett. B, 214:617–620, 1988.
  * [29] Rudolph C. Hwa and Jicai Pan.  Intermittency in the Ginzburg-Landau theory.  Phys. Lett. B, 297:35–38, 1992.
  * [30] Rudolph C. Hwa.  Scaling exponent of multiplicity fluctuation in phase transition.  Phys. Rev. D, 47:2773–2781, 1993.


A New family of methods for solving delay differential equations
Yogita Mahatekar, Pallavi S. Scindia

Department of Mathematics, College of Engineering Pune, Pune - 411005, India,
<EMAIL_ADDRESS><EMAIL_ADDRESS>

Abstract

In the present paper, we introduce a new family of $\theta-$methods for
solving delay differential equations. New methods are developed using a
combination of decomposition technique viz. new iterative method proposed by
Daftardar Gejji and Jafari and existing implicit numerical methods. Using
Butcher tableau, we observed that new methods are non Runge-Kutta methods.
Further, convergence of new methods is investigated along with its stability
analysis. Applications to variety of problems indicates that the proposed
family of methods is more efficient than existing methods.
Keywords: New iterative method (NIM), Delay differential equation.

##  1 Introduction

A delay differential equation (DDE) is a differential equation in which state
function is given in terms of value of the function at some previous times.
Introduction of delay term in modelling allows better representation of real
life phenomenon and enriches its dynamics. Due to presence of delay terms in
the model, Delay differential equations (DDEs) are infinite dimentional and
hence are difficult to analyse. Hence now a days, solving delay differential
equations is an important area of research. Every DDE cannot be integrated
analytically and hence there is a need to be dependant on numerical methods to
solve DDEs. To develop efficient, stable and accurate numerical algorithms is
primarily important task of research.

DDEs are receiving increasing importance in many areas of science and
engineering like biological processes, population growth and decay models,
epidemiology, physiology, neural networks etc [20], [6], [3]. The classical
numerical methods such as Eulers method, trapezoidal method, Runge-Kutta
methods are discussed in [8]. Enright and Hu [14] have developed a tool to
solve DDEs with vanishing delay with iteration and interpolation technique.
Karoui and Vaillancourt [19] presneted a SYSDEL code to solve DDEs using Runge
kutta methods of desired convergence order. In [21], a new MATLAB program
dde23 has been developed to solve wide range of DDEs with constant delays. A
new Adomian decomposition method is given in [15] to solve DDEs. Further in
[17], two point predictor-corrector block method for solving DDEs is
described. Recently New iterative method (NIM) developed by Daftardar-Gejji
and Jafari is used to develop many efficient numerical methods to solve
fractional differential equations [12], partial differential equations [9],
boundary value problems [10]. In [2], a solution of DDE is acheived by using
Aboodh transformation method. Recently in [13] , a new numerical technique for
solving fractional generalised pantograph-delay differential equations by
using fractional order hybrid Bessel functions is developed.

In the present article, we use the powerful technique of NIM to generate new
efficient numerical tools to solve DDEs which are reducible to solve ODEs.

The paper is organized as follows. In section 2, important preliminaries like
Delay differential equations (DDEs), New iterative method (NIM) etc. are
reviewed. In next section 3, we developed a new family of numerical methods to
solve DDEs and system of DDEs. In section 4, error analysis of newly proposed
methods is done along with its stability analysis in section 5. In section 6,
some illustrative examples are solved to check the accuracy of new methods
practically. In last section 7, some important observations are made on the
basis of theoretical stability analysis and error analysis of newly proposed
methods.

##  2  Preliminaries

###  2.1 Delay differential equations

Consider a general form of an initial value problem (IVP) representing a time
delay differential equation:

$y^{\prime}=f(t,y(t),y(t-\tau));\,t\in[0,T],$ (1) $y(t)=\phi(t),\,\,-\tau\leq
t\leq 0,$ (2)

where $\phi(t):[-\tau,0]\rightarrow\mathbb{R}^{n};\,n\in\mathbb{N}$ is a real
valued function which represents histroy of the solution in the past.
Existence and uniqueness conditions for the solution of DDEs are given in
[23]. For existence and uniqueness of the solution of the IVP (1-2), we assume
that $f$ is non-linear, bounded and continuous real valued operator defined on
$[0,T]\times\mathbb{R}^{n}\times C^{1}(\mathbb{R},\mathbb{R}^{n})$ to
$\mathbb{R}^{n}$ and fulfils Lipschitz conditions with respect to the second
and third arguments:

$|f(t,x_{1},u)-f(t,x_{2},u)|\leq L_{1}|x_{1}-x_{2}|,$ (3)
$|f(t,x,u_{1})-f(t,x,u_{2})|\leq L_{2}|u_{1}-u_{2}|,$ (4)

where $L_{1}$ and $L_{2}$ are positive constants.

###  2.2 Approximation to the delay term $y(t-\tau)$

The approximation to the delay term $y(t_{n}-\tau)$ is denoted by $\nu_{n}$
[12].
When $\tau$ is constant $(t_{n}-\tau)$ may not be a grid point $t_{n}$ for any
n. Suppose $(m+\delta)h=\tau,m\in\mathbb{N}$ and $0\leq\delta<1.$
If $\delta=0,\,mh=\tau$ and $y(t_{n}-\tau)$ is approximated as

$y(t_{n}-\tau)\approx\nu_{n}=\left\\{\begin{array}[]{ccc}y_{n-m}&\mbox{
if}&n\geq m;\\\ \phi({t_{n}-\tau})&\mbox{if}&n<m.\\\ \end{array}\right.$

###  2.3 New iterative method (NIM)

Daftardar-Gejji and Jafari [11] have proposed a new iterative method (NIM) for
solving linear/non-linear functional equations of the form

$u=f+L(u)+N(u)$ (5)

where $f$ is a known function, $L$ is a linear operator and $N$ a non linear
operator. This decomposition technique is used by many researchers for solving
variety of problems such as fractional differential equations [12], boundary
value problems [10], system of non-linear functional equations [5] etc.
Recently NIM is used to develop efficient numerical algorithms to solve
ordinary differential equations and time delay fractional differential
equations too [22],[18].

In this method, We assume that eq.(5) has a series solution of the form
$u=\displaystyle\sum^{\infty}_{i=0}u_{i}.$ Since $L$ is a linear operator, we
have
$L\left(\displaystyle\sum^{\infty}_{i=0}u_{i}\right)=\displaystyle\sum^{\infty}_{i=0}L(u_{i})$
and non linear operator $N$ is decomposed as

$N(u)=N(u_{0})+[N(u_{0}+u_{1})-N(u_{0})]+[N(u_{0}+u_{1}+u_{2})-N(u_{0}+u_{1})]+\dots$

Let
$G_{0}=N(u_{0})\,\text{and}\,G_{i}=N\left(\displaystyle\sum_{n=0}^{i}u_{n}\right)-N\left(\displaystyle\sum_{n=0}^{i-1}u_{n}\right),i=1,2,3,\dots$.
Observe that $N(u)=\displaystyle\sum_{i=0}^{\infty}G_{i}$.
Putting series solution $u=\displaystyle\sum^{\infty}_{i=0}u_{i}$ in eq.(5),
we get

$\displaystyle\sum^{\infty}_{i=0}u_{i}=f+\displaystyle\sum^{\infty}_{i=0}L(u_{i})+\displaystyle\sum_{i=0}^{\infty}G_{i}.$

Taking $u_{0}=f$, and $u_{n}=L(u_{n-1})+G_{n-1},n=1,2,3,\dots$
Note that

$u=u_{0}+u_{1}+u_{2}+\dots=f+L(u_{0})+N(u_{0})+L(u_{1})+[N(u_{0}+u_{1})-N(u_{0})]+\dots=f+L(u)+N(u).$

Hence $u$ satisfies the functional eq.(5). $k$-term NIM solution is given by
$u=\displaystyle\sum^{k-1}_{i=0}u_{i}.$

##  3 New methods to solve DDEs

In this section, we represent a family of new numerical methods based on NIM
and implicit numerical methods to solve delay differential equations. In [22],
we developed a new algorithm to solve DDEs and ODEs by improving existing
trapezoidal rule of integration. In present paper, we generalize this task by
improving general $\theta-$methods to solve delay differential equations as
follows.

For solving eq.(1-2) on $[0,T]$, consider the uniform grid
$t_{n}=nh,\,n=-m,-m+1,....,-1,0,1,....N,$ where $m$ and $N$ are integers such
that $N=T/h$ and $m=\tau/h.$ We integrate eq.(1) from the node $t_{n}$ to
$t_{n+1}$ on both sides, which gives us

$y(t_{n+1})=y(t_{n})+\displaystyle\int^{t_{n+1}}_{t_{n}}f(t,y(t),y(t-\tau))\,\,dt.$
(6)

Now we approximate integration on R.H.S. in above equation as;

$\displaystyle\int^{t_{n+1}}_{t_{n}}f(t,y(t),y(t-\tau))dt\approx
h\left[(1-\theta)f(t_{n},y(t_{n}),\nu_{n})+\theta
f(t_{n+1},y(t_{n+1}),\nu_{n+1})\right],$ (7)

where $\theta$ is a parameter which lies in the closed interval $[0,1]$ and
$\nu_{n}$ is approximation to delay term at the node $t=t_{n}.$ Applying the
approximation from eq.(7) in eq.(6), we get

$y_{n+1}=y_{n}+h\left[(1-\theta)f(t_{n},y_{n},\nu_{n})+\theta
f(t_{n+1},y_{n+1},\nu_{n+1})\right],\,n=0,1,\ldots,N-1.$ (8)

This is called as a family of $\theta-$ methods to solve DDEs.
In particular, when $\theta=1,$ eq.(8) reduces to

$y_{n+1}=y_{n}+hf(t_{n+1},y_{n+1},\nu_{n+1}),\,n=0,1,\ldots,N-1.$ (9)

This eq.(9) is referred as Implicit Euler’s method. Note that when $\theta=0,$
eq.(8) yields Explicit Euler’s method and when
$\theta=\displaystyle\frac{1}{2}$ we get implicit Trapezoidal rule of
integration to solve DDEs.

###  3.1 New Reults

To generate a new family of $\theta-$methids for solving DDEs, we note that
eq.(8) is of the form $u=f+N(u),$ and hence NIM can be employed as follows:

We write,

$u=y_{n+1},$ $f=y_{n}+h(1-\theta)f(t_{n},y_{n},\nu_{n}),$ $N(u)=h\theta
f(t_{n+1},y_{n+1},\nu_{n+1}).$

For simplicity we denote, $f_{n}=f(t_{n},y_{n},\nu_{n}).$ Three term NIM
solution of eq.(8) gives us

$\displaystyle u$ $\displaystyle=$ $\displaystyle u_{0}+u_{1}+u_{2}$
$\displaystyle=$ $\displaystyle u_{0}+N(u_{0})+N(u_{0}+u_{1})-N(u_{0})$
$\displaystyle=$ $\displaystyle u_{0}+N(u_{0}+u_{1})$ $\displaystyle=$
$\displaystyle u_{0}+N(u_{0}+N(u_{0})).$

That is
$y_{n+1}=y_{n}+h(1-\theta)f_{n}+N(u_{0}+u_{1})$
$=y_{n}+h(1-\theta)f_{n}+h\theta f(t_{n+1},u_{0}+u_{1},\nu_{n+1})$
Therefore,

$y_{n+1}=y_{n}+h(1-\theta)f_{n}+h\theta
f(t_{n+1},y_{n}+h(1-\theta)f_{n}+h\theta
f(t_{n+1},y_{n}+h(1-\theta)f_{n},\nu_{n+1}),\nu_{n+1})$ (10)

Eq.(10) represents new family of $\theta-$methods to solve DDEs which can be
expressed in the following more simple form as follows:

$\begin{rcases}\begin{array}[]{r@{}l}k_{1}=f(t_{n},y_{n},\nu_{n})\\\
k_{2}=f(t_{n+1},y_{n}+h(1-\theta)k_{1},\nu_{n+1})\end{array}\\\
\begin{array}[]{r@{}l}k_{3}=f(t_{n+1},y_{n}+h(1-\theta)k_{1}+h\theta
k_{2},\nu_{n+1})\\\ \text{where,\,\,}y_{n+1}=y_{n}+h(1-\theta)k_{1}+h\theta
k_{3}.\end{array}\\\ \end{rcases}\begin{array}[]{l}\end{array}$ (11)

Case 1: When $\theta=1,$

$\begin{rcases}\begin{array}[]{r@{}l}k_{1}=f(t_{n},y_{n},\nu_{n})\\\
k_{2}=f(t_{n+1},y_{n},\nu_{n+1})\end{array}\\\
\begin{array}[]{r@{}l}k_{3}=f(t_{n+1},y_{n}+hk_{2},\nu_{n+1})\\\
\text{where,\,\,}y_{n+1}=y_{n}+hk_{3}.\end{array}\\\
\end{rcases}\begin{array}[]{l}\end{array}$ (12)

Eqs.(12) represents new improved implicit Eulers method to solve DDEs.
Case 2: When $\theta=\frac{1}{2},$

$\begin{rcases}\begin{array}[]{r@{}l}k_{1}=f_{n}=f(t_{n},y_{n},\nu_{n})\\\
k_{2}=f(t_{n+1},y_{n}+\frac{hk_{1}}{2},\nu_{n+1})\end{array}\\\
\begin{array}[]{r@{}l}k_{3}=f(t_{n+1},y_{n}+\frac{hk_{1}}{2}+\frac{hk_{2}}{2},\nu_{n+1})\\\
\text{where,\,\,}y_{n+1}=y_{n}+\frac{hk_{1}}{2}+\frac{hk_{3}}{2}.,\end{array}\\\
\end{rcases}\begin{array}[]{l}\end{array}$ (13)

Eqs.(13) represents improved trapezoidal rule to solve DDEs which is developed
in [22].
Case 3: When $\theta=0,$

$\begin{rcases}\begin{array}[]{r@{}l}k_{1}=f_{n}=f(t_{n},y_{n},\nu_{n})\\\
k_{2}=f(t_{n+1},y_{n}+hk_{1},\nu_{n+1})\end{array}\\\
\begin{array}[]{r@{}l}k_{3}=f(t_{n+1},y_{n}+hk_{1},\nu_{n+1})\\\
\text{where,\,\,}y_{n+1}=y_{n}+hk_{1}.,\end{array}\\\
\end{rcases}\begin{array}[]{l}\end{array}$ (14)

Eqs.(14) represents a method which is usual explicit Eulers method to solve
DDEs.
Case 4: When $\theta=\frac{3}{4},$

$\begin{rcases}\begin{array}[]{r@{}l}k_{1}=f_{n}=f(t_{n},y_{n},\nu_{n})\\\
k_{2}=f(t_{n+1},y_{n}+\frac{hk_{1}}{4},\nu_{n+1})\end{array}\\\
\begin{array}[]{r@{}l}k_{3}=f(t_{n+1},y_{n}+\frac{hk_{1}}{4}+\frac{3hk_{2}}{4},\nu_{n+1})\\\
\text{where,\,\,}y_{n+1}=y_{n}+\frac{hk_{1}}{4}+\frac{3hk_{3}}{4}.,\end{array}\\\
\end{rcases}\begin{array}[]{l}\end{array}$ (15)

Eqs.(15) is a new method to solve DDEs.

In above family of methods to solve DDEs, when delay term $y(t_{n}-\tau)$ and
its approximation $\nu_{n}$ are zero then DDE (1) reduces to ODE without
delay. In accordance with this, above family of methods gets reduced to the
family of numerical methods for solving ODEs (without delay) which are
developed in [1].

###  3.2 Non-Runge Kutta methods

General form of Runge-Kutta method is given by

$y_{n+1}=y_{n}+h\displaystyle\sum^{3}_{i=1}b_{i}k_{i}$ and

$\begin{rcases}\begin{array}[]{r@{}l}k_{1}=f_{n}=f(t_{n},y_{n},\nu_{n})\\\
k_{2}=f(t_{n}+c_{2}h,y_{n}+ha_{21}k_{1},\nu_{n})\end{array}\\\
\begin{array}[]{r@{}l}k_{3}=f(t_{n}+c_{3}h,y_{n}+h(a_{31}k_{1}+a_{32}k_{2}),\nu_{n})\\\
\end{array}\\\ \end{rcases}\begin{array}[]{l}\end{array}$ (16)

In newly proposed $\theta-$methods represented by eqs.(11), we have
$b_{1}(\theta)=1-\theta,\,b_{2}(\theta)=0,\,b_{3}(\theta)=\theta,\,a_{21}=1-\theta,\,a_{31}=1-\theta,\,a_{32}=\theta,\,c_{2}=1,\,c_{3}=1.$
Therefore new Family of $\theta-$methods for solving DDEs can be stated in the
form of Butcher tableu as follows.

0 |  |  |  |
---|---|---|---|---
1 | 1-$\theta$ |  |  |
1 | 1-$\theta$ | $\theta$ |  |
| 1-$\theta$ | 0 | $\theta$ |

For a Runge Kutta method it is necessary to satisfy that
$\displaystyle\sum^{i-1}_{j=1}a_{ij}=c_{i}\,\forall\,i=2,3\,\,(cf.$[4]). From
the above tableau, for $i=2,\,$ $\displaystyle\sum^{1}_{j=1}a_{2j}=c_{2}$ if
and only if $\theta=0.$ This shows that, newly proposed family of numerical
methods for solving DDEs is different from Runge Kutta methods except for
$\theta=0.$

###  3.3 New $\theta-$methods for solving a system of delay differential
equations

The numerical algorithm presented in above section can be generalized for
solving the following system of DDEs:

$y^{\prime}_{1}(t)=f_{1}(t,\overline{y}(t),\overline{y}(t-\tau)),$
$y^{\prime}_{2}(t)=f_{2}(t,\overline{y}(t),\overline{y}(t-\tau)),$ $\vdots$
$y^{\prime}_{m}(t)=f_{m}(t,\overline{y}(t),\overline{y}(t-\tau)),$

with initial condition

$\overline{y}(t)=(y_{1}(t),y_{2}(t),\ldots,y_{m}(t))=(\phi_{1}(t),\phi_{2}(t),...,\phi_{m}(t));-\tau\leq
t\leq 0.$ (17)

We let $\overline{y_{n}}$ be a vector of independent variables representing
values of $(y_{1},y_{2},\ldots,y_{m})$ at node $t=t_{n}$ and
$\overline{\nu_{n}}$ is a vector approximation to
$(y_{1}(t-\tau),y_{2}(t-\tau),\ldots,y_{m}(t-\tau))$ at $t=t_{n}.$ We obtain a
new family of $\theta-$methods to solve a system of DDEs as follows:
$k_{1,y_{i}}=f_{i}(t_{n},\overline{y_{n}},\overline{\nu_{n}})$
$k_{2,y_{i}}=f_{i}(t_{n+1},\overline{y_{n}}+h(1-\theta)k_{1,y_{i}},\overline{\nu_{n+1}})$
$k_{3,y_{i}}=f_{i}(t_{n+1},\overline{y_{n}}+h(1-\theta)k_{1,y_{i}}+h\theta
k_{2,y_{i}},\overline{\nu_{n+1}})$
Where,

$\overline{y_{n+1}}=\overline{y_{n}}+h(1-\theta)k_{1,y_{i}}+h\theta
k_{3,y_{i}},\,i=1,2,\ldots,m.$

##  4 Error analysis

###### Theorem 4.1

The new family of $\theta-$methods given by eqs.(11) forms a second order
numerical method for $\theta=\frac{1}{2}$ and has a first order convergence
for any other value of $\theta\in[0,1].$

Proof: Using Taylor’s series expansion,

$y(t_{n+1})=y(t_{n}+h)=y_{n}+hf_{n}+\frac{h^{2}}{2}(f_{t}+ff_{y}+f_{\nu})+O(h^{3}).$

Now consider $\theta-$method and Taylors expansions applied in
$k_{2},\,k_{3}.$

$\displaystyle k_{1}$ $\displaystyle=$ $\displaystyle f(t_{n},y_{n},\nu_{n})$
(18) $\displaystyle k_{2}$ $\displaystyle=$ $\displaystyle
f(t_{n+1},y_{n}+h(1-\theta)k_{1},\nu_{n+1})$ $\displaystyle=$ $\displaystyle
f(t_{n}+h,y_{n}+h(1-\theta)k_{1},\nu_{n}+h)$
$k_{2}=f_{n}+\left(hf_{t}+h(1-\theta)k_{1}f_{y}+hf_{\nu}\right)+\\\
\frac{1}{2}\left(hf_{tt}+h^{2}(1-\theta)^{2}k^{2}_{1}f_{yy}+h^{2}f_{\nu\nu}+2h^{2}(1-\theta)k_{1}f_{ty}+2h^{2}f_{t\nu}+2h^{2}(1-\theta)k_{1}f_{y\nu}\right)+O(h^{3}).$
(19) $\displaystyle k_{3}$ $\displaystyle=$ $\displaystyle
f(t_{n+1},y_{n}+h(1-\theta)k_{1}+h\theta k_{2},\nu_{n+1})$ $\displaystyle=$
$\displaystyle f(t_{n}+h,y_{n}+h(1-\theta)k_{1}+h\theta k_{2},\nu_{n}+h)$
$k_{3}=f_{n}+\left(hf_{t}+h(1-\theta)k_{1}f_{y}+h\theta
k_{2}f_{y}+hf_{\nu}\right)+\frac{1}{2}\left(h^{2}f_{tt}+(h^{2}(1-\theta)^{2}k^{2}_{1}+h^{2}\theta^{2}k^{2}_{2}+2h^{2}(1-\theta)\theta
k_{1}k_{2})f_{yy}\\\ +h^{2}f_{\nu\nu}+(2h^{2}(1-\theta)k_{1}+2h^{2}\theta
k_{2})f_{yt}+2h^{2}f_{t\nu}+2(h^{2}(1-\theta)k_{1}+h^{2}\theta
k_{2})f_{y\nu}\right)+O(h^{3}).$ (20)

Putting (18), (19), (20) in $\theta-$method given by eqs.(11), error $e_{n+1}$
of the method is given by

$e_{n+1}=y(t_{n+1})-y_{n+1}\\\
=y_{n}+hf_{n}+\frac{h^{2}}{2}f_{t}+\frac{h^{2}ff_{y}}{2}+\frac{h^{2}f_{\nu}}{2}+O(h^{3})-\left(y_{n}+h(1-\theta)f_{n}+h\theta
f_{n}\\\ h^{2}\theta
f_{t}+h^{2}\theta(1-\theta)k_{1}f_{y}+h^{2}\theta^{2}k_{2}f_{y}+h^{2}\theta
f_{\nu}\right)+O(h^{3}).$ (21)

Clearly, for $\theta=0,1,3/4$ method has a linear convergence and for
$\theta=1/2$ its order of convegence is quadratic.

###  4.1 Error analysis of new improved implicit Euler’s method for solving
DDEs

Refer to Eqs.(12), For $\theta=1,$ we get new improved implicit Eulers method
to solve DDEs, which is given by the formula:

$y_{n+1}=y_{n}+hf[t_{n+1},y_{n}+hf(t_{n+1},y_{n},\nu_{n+1}),\nu_{n+1}].$ (22)

By inserting analytical solution in the above equation, we get the truncation
errorv $T_{n}$ as

$\displaystyle\frac{y(t_{n+1})-y(t_{n})}{h}-f[t_{n+1},y(t_{n})+hf(t_{n+1},y(t_{n}),\nu_{n+1}),\nu_{n+1}]=T_{n}.$
(23)

and by eq.(23), we get

$\displaystyle\frac{y_{n+1}-y_{n}}{h}-f[t_{n+1},y_{n}+hf(t_{n+1},y_{n},\nu_{n+1}),\nu_{n+1}]=0.$
(24)

Subtracting eq.$(\ref{em23})$ and eq.$(\ref{em24})$, we get

$T_{n}=\displaystyle\frac{e_{n+1}-e_{n}}{h}-f[t_{n+1},y(t_{n})+hf(t_{n+1},y(t_{n}),\nu_{n+1}),\nu_{n+1}]+f[t_{n+1},y_{n}+hf(t_{n+1},y_{n},\nu_{n+1}),\nu_{n+1}].$
(25)

This implies that,

$hT_{n}=e_{n+1}-e_{n}-h\left(f[t_{n+1},y(t_{n})+hf(t_{n+1},y(t_{n}),\nu_{n+1}),\nu_{n+1}]\\\
-f[t_{n+1},y_{n}+hf(t_{n+1},y_{n},\nu_{n+1}),\nu_{n+1}]\right).$

Therefore,

$|e_{n+1}|\leq|e_{n}|+hL_{1}|(y(t_{n})+hf(t_{n+1},y(t_{n}),\nu_{n+1})-(y_{n}+hf(t_{n+1},y_{n},\nu_{n+1})|+hT_{n}\\\
\leq|e_{n}|+hL_{1}|e_{n}|+h^{2}L_{1}^{2}|e_{n}|+|T_{n}|h\\\
\leq(1+hL_{1}+h^{2}L_{1}^{2})|e_{n}|+|T_{n}|h;\,\,n=0,1,\ldots N.$

Let $T=\displaystyle\max_{0\leq n\leq(N-1)}|T_{n}|$ Therefore,

$|e_{n+1}|\leq(1+hL+h^{2}L^{2})|e_{n}|+Th$

Now by induction,

$|e_{n}|\leq(1+hL+h^{2}L^{2})^{n}|e_{0}|+\left(\frac{(1+hL+h^{2}L^{2})^{n-1}-1}{(1+hL+h^{2}L^{2})-1}\right)Th\\\
\leq(1+hL+h^{2}L^{2})^{n}|e_{0}|+\left(\frac{e^{nhL}-1}{(Lh+1)Lh}\right)Th\\\
\leq e^{nhL}|e_{0}|+\left(\frac{e^{nhL}-1}{L}\right)T\\\ \leq
e^{(t_{n}-t_{0})L}|e_{0}|+\left(\frac{e^{(t_{n}-t_{0})L}-1}{L}\right)T$

Since $nh=t_{n}-t_{0}.$ Noting that, Truncation error
$T\leq\frac{1}{2}hy^{\prime\prime}(\zeta)$ and if $y^{\prime\prime}\leq M_{2}$
then $T\leq\frac{hM_{2}}{2}.$ Therefore,

$|e_{n}|\leq
e^{(t_{n}-t_{0})L_{1}}|e_{0}|+\left(\displaystyle\frac{(e^{(t_{n}-t_{0})L_{1}}-1)hM_{2}}{2L_{1}}\right)\\\
\leq\left(\displaystyle\frac{(e^{(t_{n}-t_{0})L_{1}}-1)hM_{2}}{2L_{1}}\right),\text{noting
that $|e_{0}|$ is zero.}$

hence as $h\rightarrow 0,$ $|e_{n}|\rightarrow 0.$ Here $L_{1}$ is Lipschitz
constant as given in eq.(3). Hence method is convergent. Similary, for
$\theta=\frac{1}{2},$ we proved the convergence of the new improved
trapezoidal rule for solving DDEs in [22] and for other values of $\theta$
convergenrce can be proved on similar lines.

##  5 Stability Analysis

###### Definition 5.1

The numerical method for solving IVP eq.(1-2) is said to be zero-stable if
small perturbation in the initial condition of IVP do not cause the numerical
approximation to diverge from the exact solution, provided the exact solution
of the IVP is bounded.

Consider the IVP eq.(1-2) and let $\epsilon{y}(0)=\epsilon{y_{0}}$ be the new
initial value (perturbed initial condition) obtained by making a small change
in $y(0)=y_{0}.$

###### Theorem 5.1

Let $y_{n}$ be the solution obtained by new improved implicit Euler’s method
(Case 1 in section3) for solving DDE at the node $t_{n}$ with the initial
condition $y(0)=y_{0}$ and let $\epsilon y_{n}$ be the solution obtained by
the same numerical method with perturbed initial condition
$\epsilon{y_{0}}=y_{0}+\epsilon_{0};\,\epsilon_{0}>0.$ We assume that $f(t,y)$
satisfies Lipschitz condition with respect to second variable and third
variable with Lipschitz constant $L_{1},\,L_{2}$ then $\exists$ positive
constants $k$ and $\epsilon_{1}$ such that $|y_{n}-\epsilon y_{n}|\leq
k\epsilon,\,\,\forall nh\leq T,\,\,h\in(0,\epsilon_{1})\,$ whenever
$|\epsilon_{0}|\leq\epsilon.$

Proof: We prove the result by induction.

$|y_{n}-\epsilon
y_{n}|=|y_{n-1}+hf\left(t_{n},y_{n-1}+hf\left(t_{n},y_{n-1},\nu_{n}\right),\nu_{n}\right)-\epsilon
y_{n-1}-hf\left(t_{n},\epsilon y_{n-1}+hf\left(t_{n},\epsilon
y_{n-1},\nu_{n}\right),\nu_{n}\right)|\\\ \leq|y_{n-1}-\epsilon
y_{n-1}|+hL_{1}|y_{n-1}+hf\left(t_{n},y_{n-1},\nu_{n}\right)-\epsilon
y_{n-1}-hf\left(t_{n},\epsilon y_{n-1},\nu_{n}\right)|\\\
\leq|y_{n-1}-\epsilon y_{n-1}|+hL_{1}|y_{n-1}-\epsilon
y_{n-1}|+h^{2}L_{1}^{2}|y_{n-1}-\epsilon y_{n-1}|\\\
=(1+hL_{1}+h^{2}L_{1}^{2})|y_{n-1}-\epsilon y_{n-1}|$

Therefore by induction,

$|y_{n}-\epsilon y_{n}|\leq(1+hL_{1}+h^{2}L_{1}^{2})^{n}|y_{0}-\epsilon
y_{0}|\\\ \leq e^{nhL_{1}}|\epsilon_{0}|\\\
=e^{T_{M}L_{1}}\epsilon=k\epsilon,\,\,\text{where}\,k=e^{T_{M}L_{1}}>0.$

This proves that the new method for ($\theta=1$) is stable.

##  6 Illustrative examples

To demonstrate applicability of newly proposed methods, We present here some
illustrative examples which are solved using Mathematica 12.

###### Example 6.1

Consider the delay logistic differential equation

$\displaystyle y^{\prime}(t)$ $\displaystyle=$ $\displaystyle
0.3y(t)(1-y(t-1)),\,y(t\leq 0)=0.1.$ (26)

(a) Solution of (26) by new method ($\theta=1$)

(b) Error in solution of (26) when solved using new method ($\theta=1$)

Figure 1: Dashed graph: Solution by new method ($\theta=1$), Dotted graph:
Exact solution. It is noted that in Fig.1 (a) exact solution and approximate
solution by the new method ($\theta=1$) overlaps on each other.

Step length is taken as $h=0.01.$

###### Example 6.2

Consider the differential equation without delay [7]

$\displaystyle y^{\prime}(t)$ $\displaystyle=$ $\displaystyle
2-e^{-4t}-2y,\,y(0)=1,\,0\leq t\leq 10.$ (27)

Exact solution of the differential equation (27) is
$y(t)=1+\displaystyle\frac{e^{-4t}-e^{-2t}}{2}.$

(a) Solution of (27) by new method ($\theta=1$)

(b) Error in solution of (27) when solved using $\theta$ method ($\theta=1$)

Figure 2: Dashed graph: Solution by new method, Dotted graph: Exact solution.
It is noted that in Fig.2 (a) exact solution and approximate solution by the
new method overlaps on each other.

Step length is taken as $h=0.01.$

The error in the new method is shown in Fig.2 (b). Thus, the new method is
accurate. In Table 1, error in implicit backward Euler’s method $E_{1}$ and
error in new method ($\theta=1$) $E_{2}$ (obtained by taking 3-term NIM
solution) while solving (27) are compared for various values of $t$. It is
noteworthy that $E_{2}$ is always smaller than $E_{1}$ in all cases. So the
new method is more accurate than implicit backward Euler’s method.

$t_{n}$ | $S_{1}$ | $S_{2}$ | $S_{3}$ | S | $e_{1}$ | $e_{2}$ | $e_{3}$
---|---|---|---|---|---|---|---
$0$ | 1 | 1 | $1$ | $1$ | $0$ | $0$ | $0$
$0.01$ | $0.99016$ | $0.9902$ | $0.990196$ | $0.990295$ | $0.000099$ | $0.0000954$ | $0.0000994$
$0.02$ | $0.980969$ | $0.980976$ | $0.980969$ | $0.981163$ | $0.000194$ | $0.000187$ | $0.0001948$
$0.03$ | $0.972292$ | $0.97230$ | $0.972292$ | $0.972578$ | $0.000286$ | $0.000275$ | $0.0002863$
$0.04$ | $0.96414$ | $0.9641$ | $0.96414$ | $0.964514$ | $0.0003739$ | $0.00036$ | $0.0003741$
$0.05$ | $0.95648$ | $0.9565$ | $0.956488$ | $0.956947$ | $0.000458$ | $0.0004413$ | $0.0004584$
$0.06$ | $0.949315$ | $0.949334$ | $0.949315$ | $0.949854$ | $0.000538$ | $0.0005194$ | $0.0005391$
$0.1$ | $0.924966$ | $0.924993$ | $0.924966$ | $0.925795$ | $0.0008283$ | $0.00080$ | $0.0008289$
$0.2$ | $0.980969$ | $0.980976$ | $0.980969$ | $0.981163$ | $0.000194$ | $0.000187$ | $0.0001948$
$0.3$ | $0.8745$ | $0.874563$ | $0.874528$ | $0.87619$ | $0.0016622$ | $0.001628$ | $0.001663$
$0.4$ | $0.87447$ | $0.874498$ | $0.87447$ | $0.87628$ | $0.001813$ | $0.001785$ | $0.001814$
$0.5$ | $0.881873$ | $0.881893$ | $0.881872$ | $0.883728$ | $0.001855$ | $0.001834$ | $0.0018554$
Table 1:  Ex.(8.1)
$S_{1}:$ Solution by backward Euler’s method
$S_{2}:$ Solution by new method ($\theta=1$) obtained by taking 3-term NIM
solution
$S_{3}:$ Solution by new method ($\theta=1$) obtained by taking 4-term NIM
solution
$S:$ Exact solution
$e_{1}:$ Error in solution by backward Euler’s method
$e_{2}:$ Error in solution by new method ($\theta=1$) obtained by taking
3-term NIM solution
$e_{3}:$ Error in solution by new method ($\theta=1$) obtained by taking
4-term NIM solution

Observation: In this example, it is observed that (4-term NIM solution) new
method ($\theta=1$) method and backward Euler’s method gives same error and
error in these two methods is greater than (3-term NIM solution) new method
($\theta=1$). Hence, new method with three term NIM solution gives better
accuracy than implicit backward euler method and new method with 4-term NIM
solution.

###  6.1 R$\ddot{o}$ssler System with delay

Consider the R$\ddot{o}ssler$ system [16] with delay given by the following
system of differential equations.

$\displaystyle\dot{x}$ $\displaystyle=-y(t)-z(t),$ $\displaystyle\dot{y}$
$\displaystyle=x(t)+ay(t-1),$ $\displaystyle\dot{z}$
$\displaystyle=b+z(t)(x(t)-c).$ (28)

Let $a=b=0.2,\,\,x(0)=y(0)=z(0)=0.0001.$ Note that $c$ is a control parameter.
Solving above system by new method for ODE we obtain the x-waveforms and x-y
phase potraits which are depicted in Figs.3(a)-(b), Figs.4(a)-(b) and in
Figs.5(a)-(b) for $c=2.3,\,c=2.9,\,c=7.9$ respectively.

(a) x-waveform of the system (28)

(b) xy-phase diagram of the system (28)

Figure 3: c=2.3 in system(28)

(a) x-waveform of the system (28)

(b) xy-phase diagram of the system (28)

Figure 4: c=2.9 in system(28)

(a) x-waveform of the system (28)

(b) xy-phase diagram of the system (28)

Figure 5: c=7.9 in system(28)

##  7 Conclusions

In the present work, a new family of methods for solving delay differential
equations (DDEs) has been proposed which are reducible to solve ordinary
differential equations too (without delay). Newly proposed methods are then
compared with existing methods with respect to stability and accuracy. New
methods formed found to be stable and accurate. Further error analysis and
stability analysis of new methods is carried out. Numerous illustrative
examples are solved using Mathematica 12 to demonstrate the efficiency of the
method. It is observed that, new methods are non-Runge Kutta methods and are
more accurate than existing numerical methods for solving DDEs.

## References

  * [1] O. Y. Ababneh.  New numerical methods for solving differential equations, 2019.
  * [2] K. Aboodh, R. Farah, I. Almardy, and A. Osman.  Solving delay differential equations by aboodh transformation method.  International Journal of Applied Mathematics & Statistical Sciences, 7(2):55–64, 2018.
  * [3] C. Baker, G. Bocharov, C. Paul, and F. Rihan.  Modelling and analysis of time-lags in some basic patterns of cell proliferation.  Journal of mathematical biology, 37(4):341–371, 1998.
  * [4] A. Bellen and M. Zennaro.  Numerical methods for delay differential equations.  Oxford University Press, 2013.
  * [5] S. Bhalekar and V. Daftardar-Gejji.  Solving a system of nonlinear functional equations using revised new iterative method.  International Journal of Mathematical and Computational Sciences, 6(8):968–972, 2012.
  * [6] G. A. Bocharov and F. A. Rihan.  Numerical modelling in biosciences using delay differential equations.  Journal of Computational and Applied Mathematics, 125(1-2):183–199, 2000.
  * [7] T. Bui.  Explicit and implicit methods in solving differential equations.  2010\.
  * [8] J. C. Butcher.  Numerical methods for ordinary differential equations in the 20th century.  Journal of Computational and Applied Mathematics, 125(1-2):1–29, 2000.
  * [9] V. Daftardar-Gejji and S. Bhalekar.  Solving fractional diffusion-wave equations using a new iterative method.  Fractional Calculus and Applied Analysis, 11(2):193p–202p, 2008\.
  * [10] V. Daftardar-Gejji and S. Bhalekar.  Solving fractional boundary value problems with dirichlet boundary conditions using a new iterative method.  Computers & Mathematics with Applications, 59(5):1801–1809, 2010\.
  * [11] V. Daftardar-Gejji and H. Jafari.  An iterative method for solving nonlinear functional equations.  Journal of Mathematical Analysis and Applications, 316(2):753–763, 2006.
  * [12] V. Daftardar-Gejji, Y. Sukale, and S. Bhalekar.  Solving fractional delay differential equations: A new approach.  Fractional Calculus and Applied Analysis, 18(2):400–418, 2015.
  * [13] H. Dehestani, Y. Ordokhani, and M. Razzaghi.  Numerical technique for solving fractional generalized pantograph-delay differential equations by using fractional-order hybrid bessel functions.  International Journal of Applied and Computational Mathematics, 6(1):1–27, 2020.
  * [14] W. H. Enright and M. Hu.  Interpolating runge-kutta methods for vanishing delay differential equations.  Computing, 55(3):223–236, 1995.
  * [15] D. J. Evans and K. Raslan.  The adomian decomposition method for solving delay differential equation.  International Journal of Computer Mathematics, 82(1):49–54, 2005\.
  * [16] K. Ibrahim, R. Jamal, and F. Ali.  Chaotic behaviour of the rossler model and its analysis by using bifurcations of limit cycles and chaotic attractors.  In J. Phys. Conf. Ser, volume 1003, page 012099, 2018.
  * [17] F. Ishak, M. Suleiman, and Z. Omar.  Two-point predictor-corrector block method for solving delay differential equations.  MATEMATIKA: Malaysian Journal of Industrial and Applied Mathematics, 24:131–140, 2008.
  * [18] A. Jhinga and V. Daftardar-Gejji.  A new numerical method for solving fractional delay differential equations.  Computational and Applied Mathematics, 38(4):1–18, 2019.
  * [19] A. Karoui and R. Vaillancourt.  A numerical method for vanishing-lag delay differential equations.  Applied Numerical Mathematics, 17(4):383–395, 1995.
  * [20] F. A. Rihan, D. Abdelrahman, F. Al-Maskari, F. Ibrahim, and M. A. Abdeen.  Delay differential model for tumour-immune response with chemoimmunotherapy and optimal control.  Computational and mathematical methods in medicine, 2014, 2014.
  * [21] L. F. Shampine, S. Thompson, and J. Kierzenka.  Solving delay differential equations with dde23.  URL http://www. runet. edu/~ thompson/webddes/tutorial. pdf, 2000\.
  * [22] Y. Sukale and V. Daftardar-Gejji.  New numerical methods for solving differential equations.  International Journal of Applied and Computational Mathematics, 3(3):1639–1660, 2017.
  * [23] E. Süli and D. F. Mayers.  An introduction to numerical analysis.  Cambridge university press, 2003.