Walsh function-based numerical approach for nonlinear stochastic integral equations: Application to stochastic logistic models

Abstract

The current research study proposes an efficient numerical method for obtaining an approximate solution to nonlinear stochastic integral equations implementing the collocation method and the Walsh operational matrices. Given the complexity of solving each integral equation exactly, we use a numerical approach that converts the equation into a system of algebraic equations, yielding an approximate solution. Through error analysis, the method is found to exhibit linear convergence, underscoring its efficiency. We offer numerical tests to show the accuracy and usefulness of the proposed approach. This work also demonstrates the applicability of the method to dealing with the stochastic logistic model.

Introduction

Over the past decades, stochastic calculus and stochastic differential equations (SDEs) [1, 2] have emerged as fundamental tools across diverse scientific disciplines, driven by their ability to model systems subject to random fluctuations. The field traces its origins to 20th century developments in incorporating white noise processes into differential equations and analyzing their solution structures. While deterministic calculus provides exact solutions for many problems, the stochastic counterpart presents unique challenges - exact solutions remain elusive for most SDEs, particularly nonlinear formulations [3].

This mathematical framework has proven particularly transformative in biological sciences. Researchers have employed SDEs to model biochemical networks [4], simulate cellular processes [5], and develop hybrid stochastic systems for biological modeling [6]. Recent work by [7] demonstrates their value in epidemiological modeling through stochastic SIR frameworks for zoonotic diseases. Physical sciences similarly benefit, with applications ranging from molecular fluid dynamics [8] to constrained parameter estimation in physical systems [9] and neuro-physical systems with delayed interactions [10].

Chemical sciences have adopted stochastic approaches for diverse applications including: parameter estimation in chemical reaction networks [11], quantum chemical calculations via stochastically-optimized basis sets [12], and neuronal signaling pathway modeling [13]. In economics and finance, SDEs underpin modern theories, from multivariate volatility modeling [14] to jump-diffusion processes in financial mathematics [15]. The inherent complexity of nonlinear SDEs has spurred development of specialized numerical methods. While early work focused on theoretical foundations [1], contemporary research emphasizes computational techniques for approximate solutions [3]. This methodological evolution continues to drive interdisciplinary applications, making stochastic calculus an increasingly vital tool across the scientific spectrum.

Among the diverse numerical approaches for solving stochastic differential equations, function approximation methods using orthogonal bases have emerged as a powerful tool. Due to the orthogonality property of the function, the calculation of the coefficient becomes simpler which results in simplifying the process of approximation. Compared to other approximation methods, the orthogonal function approximation becomes more accurate especially when the number of elements in the basis is larger. These techniques approximate solutions either directly for differential equations or via their integral formulations, leveraging the inherent properties of orthogonal systems to transform complex problems into tractable algebraic systems. The Walsh function, introduced by [16] as a complete set of normal orthogonal functions, has proven particularly effective in this context. Recent work by [17] and [18] demonstrates its success in solving multi-dimensional stochastic Volterra and Volterra-Fredholm integral equations, showcasing its computational efficiency and spectral accuracy.

Beyond Walsh functions, other orthogonal systems offer unique advantages such as block-pulse functions provide simplicity in discretization, enabling efficient numerical solutions for nonlinear SDEs [19] and hybrid methods for Itô-Volterra equations [20], wavelet-based methods, such as Haar wavelets [21] and Taylor wavelet [22], combine localization properties with multi-resolution analysis, making them ideal for problems with sharp gradients or singularities. Moreover, polynomial bases like Legendre polynomials [23] and modified forms (e.g., QR factorized bases [24]) excel in spectral convergence for smooth solutions, as seen in epidemiological modeling [25]. Pell polynomials and their extensions [26, 27] have shown promise for fractional stochastic systems and high-order nonlinear problems due to their recursive structure and numerical stability. Also, Vieta-Lucas polynomials have been used to solve nonlinear stochastic Itô-Volterra integral equation and provide solutions without solving the algebraic equation [28].

Some other methods based on interpolation including meshless barycentric rational interpolation [29], Floater-Hormann interpolation [30], least-square method [31], and cubic B-spline [32] used to obtain an approximation to the solution of nonlinear SDEs.

In this work, we focus on Walsh functions for solving nonlinear stochastic integral equations (NLSIEs). The Walsh functions are simple to generate and manipulate, which become an efficient algorithm for solving complex problems. Their binary-valued orthogonality reduces computational overhead, while their ability to handle discontinuities aligns well with the irregular solutions typical of stochastic systems. This choice builds on their proven success in recent literature while addressing gaps in handling nonlinearities efficiently. As a key application, we employ Walsh functions to solve the stochastic logistic equation, a fundamental model in population dynamics that incorporates randomness to account for environmental fluctuations.

The stochastic logistic equation extends the classical deterministic model of population growth, originally developed by [33] and later refined in ecological studies by [34], to incorporate random environmental fluctuations. Although early work by [35] established foundational principles of population dynamics, recent advances have focused on fractional and stochastic formulations that better capture complex biological realities. The authors in [36] introduced a novel discretization approach for the nonlinear fractional logistic equation, demonstrating improved accuracy in modeling population systems with memory effects. Comparative studies by [37] systematically evaluated various fractional basis functions for solving these models, identifying optimal approximation schemes. Further refinements were achieved through Legendre-collocation methods [38], while [39] developed sophisticated numerical approximations using fractional-order Bessel and Legendre bases, significantly advancing our ability to simulate chaotic population behavior. These mathematical innovations provide powerful tools for analyzing stochastic population dynamics across ecological and biological systems, bridging the gap between classical theory and modern computational techniques.

This article is divided into seven major sections, in which Sect. 2 is dedicated to the definitions and properties of the Walsh function, which is used to provide an approximation for any suitable function, which becomes the prerequisite to the future section. In Sect. 3, a numerical solution to the general NLSIE is proposed in which the integral equation is approximated by the Walsh function and reduced to an algebraic system of equations that can further be computed to obtain an approximate solution of the problem. In Sect. 4, we analyze the convergence of the method and establish that its numerical solution converges to the true solution at a linear rate, specifically with an error of order $\mathcal{O}(h)$. In Sect. 5, a couple of examples are provided in support of the suggested approximated approach.

In Sect. 6, we illustrate the practical utility of our Walsh function approximation method by applying it to a biological system. Specifically, we analyze a logistic growth model under both deterministic conditions and stochastic perturbations (white noise), thereby showcasing the method’s effectiveness in handling real-world biological scenarios with varying noise intensities. This comparative analysis serves to validate the robustness and applicability of our approach in biological modeling. The final Sect. 7 synthesizes our key findings and presents concluding remarks. Additionally, we outline several productive directions for future investigation, including potential extensions of the method to more complex biological systems and higher-dimensional stochastic problems.

Walsh function approximation

In 1923, Walsh developed a complete orthonormal set of functions later known as Walsh function [16]. We begin this section by defining the Walsh functions as

Definition 2.1

(Walsh Function)

The $n^{\text{th}}$ Walsh function, denoted $w_n(t)$, is defined as the product $w_n(t)=\prod _{k=1}^s{(r_k(t))}^{b_k},$ where:

• $s=\lfloor {log}_{2}n\rfloor +1$ with $\lfloor \cdot \rfloor$ denoting the floor function,

• The coefficients $b_s,b_{s-1},\dots ,b_1\in \{0,1\}$ represent the binary digits of n, satisfying $n=\sum _{k=1}^s{b}_k{2}^{k-1}.$

Here, $r_i(t)$ is called the ith Rademacher’s function which is defined in [16].

The Walsh functions satisfy the following properties:

Orthonormality

The set of Walsh functions is orthonormal. i.e., ${\int }_0^1{w}_i(t)w_j(t)dt=\{\begin{array}{rl}1& \text{if}j=j\\ 0& \text{if}i\ne j.\end{array}$

Completeness

For every $f\in L^2[0,1)$${\int }_0^1{f}^2(t)dt=\sum _{i=0}^{\mathrm{\infty }}{f}_i^2{\parallel w_i(t)\parallel }^2,$ where $f_i={\int }_0^{1}f(t)w_i(t)dt$.

Function approximation

We start this subsection by referring to the most important theorem for orthogonal function approximation, which is as follows:

Theorem 2.2

[40] If$\{u_{\alpha }\}$be an orthonormal basis in Hilbert spaceH, then for every$f\in H$, we have$f=\sum _n⟨f,u_n⟩u_n,$where$\{u_1,u_2,\dots \}=\{u_{\alpha }:⟨f,u_{\alpha }⟩\ne 0\}$is an orthonormal basis forH.

Since $\{w_n(t)\}$ is a set of complete orthonormal functions [41], hence we can obtain a Walsh approximation of any $f(t)\in L^2[0,1)$ by $f(t)\approx f_m(t)=\sum _{i=0}^{m-1}{c}_i{w}_i(t),$ where $c_i={\int }_0^{1}f(s)w_i(s)ds$. Fixing $m=2^k$ for $k=0,1,\dots$, and $h=\frac{1}{m}$ and further simplification leads to

2.1

$f(t)\approx f_m(t)=F^{\mathrm{\top }}{T}_{W}W(t),$ where $\begin{array}{rl}W(t)& ={\left[\begin{array}{ccc}{w}_0(t)& w_1(t)& w_2(t)\dots w_{m-1}(t)\end{array}\right]}^{\mathrm{\top }},\\ F& ={\left[\begin{array}{ccc}{f}_0& f_1& f_2\dots f_{m-1}\end{array}\right]}^{\mathrm{\top }},\\ f_i& ={\int }_{ih}^{(i+1)h}f(s)ds.\end{array}$ Also, $T_W$ is a $m\times m$ non-singular matrix such that $T_W{T}_W^{\mathrm{\top }}=mI$ (I is identity matrix) and $T_W=[w_i({\eta }_j)]$, ${\eta }_j\in ((j-1)h,jh)$, $i,j=1,2,\dots ,m$. The following results are direct: Here’s a refined presentation of your statement with improved mathematical clarity and notation:

Theorem 2.3

[17] Let$\mathbf{W}(t)={[w_1(t),w_2(t),\dots ,w_m(t)]}^{\mathrm{\top }}$be a Walsh function vector. The integral of$\mathbf{W}(t)$with respect totsatisfies${\int }_0^t\mathbf{W}(s)ds=\mathrm{\Lambda }\mathbf{W}(t),$where Λ is the operational matrix of integration given by$\mathrm{\Lambda }=\frac{1}{m}{T}_{W}P{T}_W.$Here, $T_W$is the Walsh transform matrix and the integration operational matrixPfor the Walsh functions is given by$P=\frac{h}{2}{\left[\begin{array}{ccccc}1& 2& 2& \cdots & 2\\ 0& 1& 2& \cdots & 2\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 1\end{array}\right]}_{m\times m}.$Also, mis the dimension of the Walsh function system.

The above Theorem leads directly to the following conclusion: $\begin{array}{rcl}{\int }_0^{t}f(s)ds& =& {\int }_0^t{F}^{\mathrm{\top }}{T}_{W}W(s)ds=F^{\mathrm{\top }}{T}_W{\int }_0^{t}W(s)ds\\ & =& F^{\mathrm{\top }}P{T}_{W}W(t).\end{array}$ Next, we have the following result:

Theorem 2.4

[17] Let$\mathbf{W}(t)={[w_1(t),\dots ,w_m(t)]}^{\mathrm{\top }}$be a Walsh function vector and$B(t)$a Brownian motion. The Itô integral of$\mathbf{W}(t)$satisfies: ${\int }_0^t\mathbf{W}(s)dB(s)={\mathrm{\Lambda }}_S\mathbf{W}(t),$where the operational matrix${\mathrm{\Lambda }}_S$is given by: ${\mathrm{\Lambda }}_S=\frac{1}{m}{T}_W{P}_S{T}_W,$and$P_S$is the Stochastic integration matrix with entries: $P_S=\left[\begin{array}{cccc}B\left(\frac{h}{2}\right)& B(h)& \cdots & B(h)\\ 0& B\left(\frac{3h}{2}\right)-B(h)& \cdots & B(2h)-B(h)\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & B\left(\frac{(2m-1)h}{2}\right)-B((m-1)h)\end{array}\right].$

Hence, from the above Theorems we obtain that

2.2

${\int }_0^{t}W(s)W^{\mathrm{\top }}(s)T_{W}Xds=m{T}_W\tilde{X}P\frac{1}{m}{T}_{W}W(t)=T_W\tilde{X}P{T}_{W}W(t),$ and

2.3

${\int }_0^{t}W(s)W^{\mathrm{\top }}(s)T_{W}XdB(s)=m{T}_W\tilde{X}{P}_S\frac{1}{m}{T}_{W}W(t)=T_W\tilde{X}{P}_S{T}_{W}W(t),$ where $\tilde{X}=diag(X)$ denotes the diagonal matrix with entries ${\tilde{X}}_{ii}=x_i$ for $i=0,2,3,\dots ,m-1$ with $x_i={\int }_{ih}^{(i+1)h}x(p)dp$. The detailed descriptions can be found in [17].

Approximate solution of nonlinear stochastic integral equations

We consider following NLSIE

3.1

$x(t)=f(t)+{\int }_0^t{k}_1(s)\rho (x(s))ds+{\int }_0^t{k}_2(s)\sigma (x(s))dB(s),$ where the stochastic processes $x(t)$, $f(t)$, $k_1(t)$, $k_2(t)$ form a collection of stochastic processes indexed by $t\in [0,T)$, all adapted to the filtered probability space $(\mathrm{\Omega },F_t,P)$ and $\rho (\cdot )$, $\sigma (\cdot )$ being measurable functions of the state process $x(t)$. Here, $B(t)$ is a standard Brownian motion and ${\int }_0^t{k}_2(s)x(s)dB(s)$ is the Itô Integral. Here, $x(t)$ is to be evaluated.

Letting $z_1(t)=\rho (x(t))$ and $z_2=\sigma (x(t))$ in (3.1) we get

3.2

$x(t)=f(t)+{\int }_0^t{k}_1(s)z_1(s)ds+{\int }_0^t{k}_2(s)z_2(s)dB(s).$ Using equation (2.1) in (3.2), we arrive at $\begin{array}{rcl}{X}^{\mathrm{\top }}{T}_{W}W(t)& =& F^{\mathrm{\top }}{T}_{W}W(t)+{\int }_0^t{K}_1^{\mathrm{\top }}{T}_{W}W(s)W^{\mathrm{\top }}(s)T_W{Z}_{1}ds\\ & & +{\int }_0^t{K}_2^{\mathrm{\top }}{T}_{W}W(s)W^{\mathrm{\top }}(s)T_W{Z}_{2}dB(s)\\ & =& F^{\mathrm{\top }}{T}_{W}W(t)+K_1^{\mathrm{\top }}{T}_W{\int }_0^{t}W(s)W^{\mathrm{\top }}(s)T_W{Z}_{1}ds\\ & & +K_2^{\mathrm{\top }}{T}_W{\int }_0^{t}W(s)W^{\mathrm{\top }}(s)T_W{Z}_{2}dB(s).\end{array}$ Substituting (2.2) and (2.3) in the above equation, we get $X^{\mathrm{\top }}{T}_{W}W(t)=F^{\mathrm{\top }}{T}_{W}W(t)+K_1^{\mathrm{\top }}{T}_W{T}_W\widetilde{Z_1}P{T}_{W}W(t)+K_2^{\mathrm{\top }}{T}_W{T}_W\widetilde{Z_2}{P}_S{T}_{W}W(t),$ which gives us $X^{\mathrm{\top }}{T}_{W}W(t)=F^{\mathrm{\top }}{T}_{W}W(t)+m{K}_1^{\mathrm{\top }}\widetilde{Z_1}P{T}_{W}W(t)+m{K}_2^{\mathrm{\top }}\widetilde{Z_2}{P}_S{T}_{W}W(t).$ Further simplifications reveal that Hence can be solved to provide the unique solution to the integral equation specified in (3.1). Here $Z_1$ and $Z_2$ is obtained by solving the following system

3.3

$\{\begin{array}{c}{Z}_1^{\mathrm{\top }}{T}_{W}W(t_j)=\beta (F^{\mathrm{\top }}{T}_{W}W(t_j)+m{K}_1^{\mathrm{\top }}\widetilde{Z_1}P{T}_{W}W(t_j)+m{K}_2^{\mathrm{\top }}\widetilde{Z_2}{P}_S{T}_{W}W(t_j)),\hfill \\ Z_2^{\mathrm{\top }}{T}_{W}W(t_j)=\sigma (F^{\mathrm{\top }}{T}_{W}W(t_j)+m{K}_1^{\mathrm{\top }}\widetilde{Z_1}P{T}_{W}W(t_j)+m{K}_2^{\mathrm{\top }}\widetilde{Z_2}{P}_S{T}_{W}W(t_j)),\hfill \end{array}$ at $t_j=\frac{2j+1}{2m}$ for $j=0,1,\dots ,m-1$.

Error analysis

This section presents a rigorous analysis of the convergence properties of the numerical solution obtained by the proposed method. Specifically, we examine the convergence of the approximate solution to the exact solution of the given nonlinear stochastic integral equation (NLSIE). To quantify the accuracy of the approximation, we introduce the mean-square norm (or $L^2$-norm) for a stochastic process $X(t)$, defined as: ${\parallel X\parallel }_2={\left(\mathbb{E}[|X{|}^2]\right)}^{\frac{1}{2}},$ where $\mathbb{E}[\cdot ]$ denotes the expectation with respect to the underlying probability measure. This norm provides a natural measure for comparing the exact and approximate solutions in the mean-square sense. The analysis demonstrates that, under appropriate conditions, the numerical solution $X_m$ converges to the true solution X as the discretization parameter $h\to 0$.

Theorem 4.1

[18] If$f\in L^2[0,1)$satisfies the Lipschitz condition with Lipschitz constantC, then${\parallel e_m(t)\parallel }_2=\mathcal{O}(h)$, where$e_m(t)=|f(t)-{\sum }_{i=0}^{m-1}{c}_i{w}_i(t)|$and$c_i={\int }_0^{1}f(s)w_i(s)ds$.

Theorem 4.2

Suppose$x_m(t)$denote them-th order numerical approximation to the exact solution of the NLSIE (3.1). If

1. The mappings$f:[0,1)\to \mathbb{R}$and$k_i:[0,1)\to \mathbb{R}$ ($i=1,2$) belong to$L^2([0,1))$and exhibit uniform Lipschitz continuity with respective constantsC, $L_1$, and$L_2$respectively.

2. let$|k_1(t)|\le {\rho }_1$, $|k_2(t)|\le {\rho }_2$, $|\beta (x(t))|\le {\zeta }_1$and$|\sigma (x(t))|\le {\zeta }_2$, and

3. for${\eta }_1,{\eta }_2\ge 0$, $|\beta (x(t))-\beta (x_m(t))|\le {\eta }_1|x(t)-x_m(t)|$, $|\sigma (x(t))-\sigma (x_m(t))|\le {\eta }_2|x(t)-x_m(t)|$

then${\parallel x(t)-x_m(t)\parallel }_2^2=\mathcal{O}(h^2)$.

Proof

Consider $x_m(t)$, the Walsh function-based numerical approximation to the solution of equation (3.1). Then, we have $\begin{array}{rcl}x(t)-x_m(t)& =& f(t)-f_m(t)\\ & +& {\int }_0^t(k_1(s)\beta (x(s))-k_{1m}(s)\beta (x_m(s)))ds\\ & +& {\int }_0^t(k_2(s)\sigma (x(s))-k_{2m}(s)\sigma (x_m(s)))dB(s)\end{array}$ or, $\begin{array}{rcl}|x(t)-x_m(t)|& \le & |f(t)-f_m(t)|\\ & +& \left|{\int }_0^t\right(k_1(s)\beta (x(s))-k_{1m}(s)\beta (x_m(s))\left)ds\right|\\ & +& \left|{\int }_0^t\right(k_2(s)\sigma (x(s))-k_{2m}(s)\sigma (x_m(s)))dB(s)|.\end{array}$ We know that, ${(a+b+c)}^2\le 3a^2+3b^2+3c^2$, which implies that $\begin{array}{rcl}|x(t)-x_m(t){|}^2& \le & 3|f(t)-f_m(t){|}^2\\ & +& 3{\left|{\int }_0^t\right(k_1(s)\beta (x(s))-k_{1m}(s)\beta (x_m(s))\left)ds\right|}^2\\ & +& 3{\left|{\int }_0^t\right(k_2(s)\sigma (x(s))-k_{2m}(s)\sigma (x_m(s)))dB(s)|}^2.\end{array}$ It follows that

4.1

$\begin{array}{rcl}\mathbb{E}(|x(t)-x_m(t){|}^2)& \le & \mathbb{E}(3|f(t)-f_m(t){|}^2)\\ & +& \mathbb{E}\left(3{\left|{\int }_0^t\right(k_1(s)\beta (x(s))-k_{1m}(s)\beta (x_m(s))\left)ds\right|}^2\right)\\ & +& \mathbb{E}\left(3{\left|{\int }_0^t\right(k_2(s)\sigma (x(s))-k_{2m}(s)\sigma (x_m(s)))dB(s)|}^2\right).\end{array}$ Suppose, $I_1=\mathbb{E}\left(3{\left|{\int }_0^t\right(k_1(s)\beta (x(s))-k_{1m}(s)\beta (x_m(s))\left)ds\right|}^2\right),$ and $I_2=\mathbb{E}\left(3{\left|{\int }_0^t\right(k_2(s)\sigma (x(s))-k_{2m}(s)\sigma (x_m(s)))dB(s)|}^2\right).$ Using Theorem 4.1 in (4.1) we have

4.2

$\mathbb{E}(|x(t)-x_m(t){|}^2)\le C^2{h}^2+I_1+I_2.$

Now, we get $\begin{array}{rcl}|k_1(s)\beta (x(s))-k_{1m}(s)\beta (x_m(s))|=& & |k_1(s)\beta (x(s))-k_1(s)\beta (x_m(s))\\ & +& k_1(s)\beta (x_m(s))-k_{1m}(s)\beta (x_m(s))|\\ \le & & |k_1(s)||\beta (x(s))-\beta (x_m(s))|\\ & +& |k_1(s)-k_{1m}(s)||\beta (x_m(s))|\\ \le & & |k_1(s)||\beta (x(s))-\beta (x_m(s))|\\ & +& |k_1(s)-k_{1m}(s)||\beta (x(s))|\\ & +& |k_1(s)-k_{1m}(s)||\beta (x(s))-\beta (x_m(s))|.\end{array}$ Let $|k_1(s)|\le {\rho }_1$, $|\beta (x(t))|\le {\zeta }_1$, $|\beta (x(t))-\beta (x_m(t))|\le {\eta }_1|x(t)-x_m(t)|$. Then, we get $\begin{array}{rcl}|k_1(s)\beta (x(s))-k_{1m}(s)\beta (x_m(s))|& \le & {\rho }_1{\eta }_1|x(s)-x_m(s)|\\ & +& \sqrt{2}{L}_{1}h{\zeta }_1+\sqrt{2}{L}_{1}h{\eta }_1|x(s)-x_m(s)|,\end{array}$ which implies,

4.3

$|k_1(s)\beta (x(s))-k_{1m}(s)\beta (x_m(s))|\le \sqrt{2}{L}_{1}h{\zeta }_1+({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)|x(s)-x_m(s)|.$ This gives us $\begin{array}{rcl}{I}_1& \le & \mathbb{E}\left(3{\left({\int }_0^t|k_1(s)\beta (x(s))-k_{1m}(s)\beta (x_m(s))|ds\right)}^2\right)\\ & \le & \mathbb{E}\left(3{\left({\int }_0^t\right(\sqrt{2}{L}_{1}h{\zeta }_1+({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)|x(s)-x_m(s)|\left)ds\right)}^2\right).\end{array}$ For any $t>0$ and $f\in L^2[0,t)$, the Cauchy-Schwarz inequality yields: ${|{\int }_0^{t}f(s)ds|}^2\le t{\int }_0^t|f{|}^{2}ds.$ This implies that $\begin{array}{rcl}{I}_1& \le & \mathbb{E}\left(3{\int }_0^t\right(\sqrt{2}{L}_{1}h\zeta \\ & +& {({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)|x(s)-x_m(s)|)}^{2}ds)\\ & \le & \mathbb{E}\left(6{\int }_0^t\right({\left(\sqrt{2}{L}_{1}h{\zeta }_1\right)}^2\\ & +& {(({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)|x(s)-x_m(s)|)}^2\left)ds\right)\\ & \le & \mathbb{E}(6{\left(\sqrt{2}{L}_{1}h{\zeta }_1\right)}^2\\ & +& 6{({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)}^2{\int }_0^t|x(s)-x_m(s){|}^{2}ds).\end{array}$ Hence,

4.4

$\begin{array}{rcl}{I}_1& \le & 6{\left(\sqrt{2}{L}_{1}h{\zeta }_1\right)}^2\\ & +& 6{({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)}^2{\int }_0^t\mathbb{E}(|x(s)-x_m(s){|}^2)ds.\end{array}$

Now, $I_2=\mathbb{E}\left(3{\int }_0^t{|k_2(s)\sigma (x(s))-k_{2m}(s)\sigma (x_m(s))|}^{2}ds\right).$ Let $|k_2(s)|\le {\rho }_2$, $|\sigma (x(t))|\le {\zeta }_2$, $|\sigma (x(t))-\sigma (x_m(t))|\le {\eta }_2|x(t)-x_m(t)|$. Then, one gets $|k_2(s)\sigma (x(s))-k_{2m}(s)\sigma (x_m(s))|\le \sqrt{2}{L}_{2}h{\zeta }_2+({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)|x(s)-x_m(s)|.$ Therefore, we arrive at $\begin{array}{rcl}{I}_2\le & & \mathbb{E}\left(3{\int }_0^t{(\sqrt{2}{L}_{2}h{\zeta }_2+({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)|x(s)-x_m(s)|)}^{2}ds\right)\\ \le & & \mathbb{E}\left(6{\int }_0^t\right({(\sqrt{2}{L}_{2}h{\zeta }_2)}^2+{({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)}^2|x(s)-x_m(s){|}^2)ds)\\ & =& \mathbb{E}(6{(\sqrt{2}{L}_{2}h{\zeta }_2)}^2+6{({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)}^2{\int }_0^t|x(s)-x_m(s){|}^{2}ds)\end{array}$ Hence,

4.5

$I_2\le 6{(\sqrt{2}{L}_{2}h{\zeta }_2)}^2+6{({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)}^2{\int }_0^t\mathbb{E}(|x(s)-x_m(s){|}^2)ds.$

Using (4.4) and (4.5) in (4.2), we have $\begin{array}{rcl}\mathbb{E}(|x(t)-x_m(t){|}^2)& \le & C^2{h}^2+6{\left(\sqrt{2}{L}_{1}h{\zeta }_1\right)}^2\\ & +& 6{({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)}^2{\int }_0^t\mathbb{E}(|x(s)-x_m(s){|}^2)ds\\ & +& 6{(\sqrt{2}{L}_{2}h{\zeta }_2)}^2\\ & +& 6{({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)}^2{\int }_0^t\mathbb{E}(|x(s)-x_m(s){|}^2)ds,\end{array}$ which implies that

4.6

$\begin{array}{rcl}\mathbb{E}(|x(t)-x_m(t){|}^2)& \le & (C^2{h}^2+6{\left(\sqrt{2}{L}_{1}h{\zeta }_1\right)}^2+6{(\sqrt{2}{L}_{2}h{\zeta }_2)}^2)\\ & +& (6{({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)}^2+6{({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)}^2)\\ & \times & {\int }_0^t\mathbb{E}(|x(s)-x_m(s){|}^2)ds.\end{array}$ Let $\begin{array}{rcl}{R}_1& =& (C^2{h}^2+6{\left(\sqrt{2}{L}_{1}h{\zeta }_1\right)}^2+6{(\sqrt{2}{L}_{2}h{\zeta }_2)}^2),\\ R_2& =& (6{({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)}^2+6{({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)}^2).\end{array}$ Applying the Grönwall’s inequality yields the estimate

4.7

$\mathbb{E}(|x(t)-x_m(t){|}^2)\le R_{1}exp\left({\int }_0^t{R}_{2}ds\right).$ which implies that,

4.8

${\parallel x(t)-x_m(t)\parallel }_2^2=\mathbb{E}(|x(t)-x_m(t){|}^2)\le R_1{e}^{R_2}=\mathcal{O}(h^2).$ □

Numerical example

In this section, we present the numerical implementation of the proposed method for approximating solutions to stochastic integral equations. The convergence of the method is validated through comparison with exact solutions. To quantify the accuracy of the approximation, we define the maximum absolute error E between the exact solution X and the approximate solution Y as: ${\parallel E\parallel }_{\mathrm{\infty }}=\underset{1\le i\le m}{max}|X_i-Y_i|,$ where $X_i$ and $Y_i$ are the i-th Walsh coefficients of the exact and approximate solutions, respectively. Here, m is the number of basis functions used in the expansion.

In the following examples, iterations run for n steps, with error E having mean ${\overline{x}}_E$ and standard deviation $s_E$. Calculations were executed in Matlab 2022b.

Example 5.1

[3] Consider the NLSIE $x(t)=x_0+a^2{\int }_0^{t}cosh(x(s)){sinh}^3(x(s))ds-a{\int }_0^t{sinh}^2(x(s))dB(s),$ where $x_0=\frac{1}{20}$ and $a=\frac{1}{30}$. The term $B(t)$ refers to a standard Brownian motion, and $x(t)$ is a stochastic process evolving on the probability space $(\mathrm{\Omega },F_t,P)$. One can verify that the actual solution is $x(t)=arccoth(aB(t)+coth(x_0))$.

We can see from Fig. 1 that as the value of m increases from 16 to 32, the approximate solution and the exact solution come close to each other, justifying the feasibility and effectiveness of the method. Also, Fig. 2 gives the error trend of the method for $m=16$ and $m=32$ with 50 iterations. Tables 1 and 2 give a numerical data of the mean error, the standard deviation of the error, and the confidence interval 95% of the mean error for $m=16$ and $m=32$, respectively, at different values of t. The CPU time for Example 5.1 is recorded as 54 seconds for $m=16$ and 123 seconds for $m=32$.

[See PDF for image]

Figure 1

Approximate versus exact solutions in Example 5.1 with $m=16$ and $m=32$

[See PDF for image]

Figure 2

Error trend analysis for Example 5.1 with parameters $m=16$, $m=32$, and iteration count $n=50$

Table 1. Summary of error metrics-including mean error, standard deviation, and 95% confidence interval-for Example 5.1, evaluated across 50 iterations using $m=16$

t	${\overline{x}}_E$	$s_E$	95% confidence interval for mean error
			Lower	Upper
0.1	0.0130 × 10−6	0.2398 × 10−9	0.0129 × 10−6	0.0131 × 10−5
0.3	0.0391 × 10−6	0.3880 × 10−9	0.0390 × 10−6	0.0392 × 10−5
0.5	0.0696 × 10−6	0.4948 × 10−9	0.0694 × 10−6	0.0697 × 10−5
0.7	0.1001 × 10−6	0.4351 × 10−9	0.0999 × 10−6	0.1002 × 10−5
0.9	0.1262 × 10−6	0.5272 × 10−9	0.1261 × 10−6	0.1264 × 10−5

Table 2. Summary of error metrics-including mean error, standard deviation, and 95% confidence interval-for Example 5.1, evaluated across 50 iterations using $m=32$

t	${\overline{x}}_E$	$s_E$	95% confidence interval for mean error
			Lower	Upper
0.1	0.0153 × 10−6	0.4712 × 10−9	0.0151 × 10−6	0.0154 × 10−6
0.3	0.0412 × 10−6	0.4430 × 10−9	0.0411 × 10−6	0.0413 × 10−6
0.5	0.0696 × 10−6	0.4785 × 10−9	0.0695 × 10−6	0.0697 × 10−6
0.7	0.0979 × 10−6	0.4882 × 10−9	0.0977 × 10−6	0.0980 × 10−6
0.9	0.1240 × 10−6	0.6006 × 10−9	0.1239 × 10−6	0.1242 × 10−6

Example 5.2

[3] Consider the NLSIE $x(t)=x_0+\frac{a^2}{2}{\int }_0^{t}x(s)ds+a{\int }_0^t\sqrt{1+x^2(s)}dB(s),$ where $x_0=\frac{1}{4}$ and $a=\frac{1}{10}$$B(t)$. The term $B(t)$ refers to a standard Brownian motion, and $x(t)$ is a stochastic process evolving on the probability space $(\mathrm{\Omega },F_t,P)$. It is shown that the actual solution is given by $x(t)=sinh(aB(t)+arcsinh(x_0))$.

Figure 3 illustrates that when the value of m increases from 16 to 32, the approximate solution converges to the exact solution, thus validating the feasibility and efficacy of the method. Figure 4 illustrates the trend of error of the approach for $m=16$ and $m=32$ with 50 iterations. Tables 3 and 4 present numerical data on the mean error, standard deviation of error, and the confidence interval of mean error for $m=16$ and $m=32$, respectively, over different values of t. For $m=16$, the CPU time for Example 5.2 is recorded to be 58 seconds, while for $m=32$, it is said to be 130 seconds.

[See PDF for image]

Figure 3

Approximate versus exact solutions in Example 5.2 with $m=16$ and $m=32$

[See PDF for image]

Figure 4

Error trend analysis for Example 5.2 with parameters $m=16$, $m=32$, and iteration count $n=50$

Table 3. Summary of error metrics-including mean error, standard deviation, and 95% confidence interval-for Example 5.2, evaluated across 50 iterations using $m=16$

t	${\overline{x}}_E$	$s_E$	95% confidence interval for mean error
			Lower	Upper
0.1	0.11860 × 10−3	0.00445 × 10−3	0.11736 × 10−3	0.11983 × 10−3
0.3	0.36038 × 10−3	0.03798 × 10−3	0.34985 × 10−3	0.37090 × 10−3
0.5	0.64203 × 10−3	0.11773 × 10−3	0.60940 × 10−3	0.67467 × 10−3
0.7	0.92423 × 10−3	0.18942 × 10−3	0.87172 × 10−3	0.97674 × 10−3
0.9	1.16574 × 10−3	0.27573 × 10−3	1.08931 × 10−3	1.24217 × 10−3

Table 4. Summary of error metrics-including mean error, standard deviation, and 95% confidence interval-for Example 5.2, evaluated across 50 iterations using $m=32$

t	${\overline{x}}_E$	$s_E$	95% confidence interval for mean error
			Lower	Upper
0.1	0.13442 × 10−3	0.00877 × 10−3	0.13199 × 10−3	0.13685 × 10−3
0.3	0.35483 × 10−3	0.04688 × 10−3	0.34184 × 10−3	0.36783 × 10−3
0.5	0.59843 × 10−3	0.10231 × 10−3	0.57007 × 10−3	0.62679 × 10−3
0.7	0.83939 × 10−3	0.16003 × 10−3	0.79503 × 10−3	0.88375 × 10−3
0.9	1.06618 × 10−3	0.21692 × 10−3	1.00605 × 10−3	1.12630 × 10−3

Application

Stochastic logistic model

In 1838, Verhulst introduced a model to examine the population size $v(t)$ over time, known as the deterministic logistic equation [42], expressed as

6.1

$\frac{dv(t)}{dt}=rv(t)(1-\frac{v(t)}{k}),v(0)=v_0>0,$ where $k>0$ represents the external carrying capacity and $r(t)$ denotes the growth rate. The primary objective was to elucidate the exponential growth of the species described by Malthus in 1798, incorporating both the inherent self-limitation of the process and the external environmental conditions, which were subsequently examined by several authors, including Volterra. The aforementioned model assumes that the growth rate and carrying capacity remain constant and deterministic; nevertheless, fluctuations in the growth rate occur in real environments.

In 1958, Elton observed that one of the primary causes of variation is environmental instability. In many countries, climatic change, driven by phenomena such as global warming, disrupts growth rates and carrying capacity, referred to as environmental noise [43].

Assume the growth rate is influenced by ambient noise such that $r(t)\to r(t)+\alpha (t)\xi (t),$ where $\xi (t)$ denotes white noise and ${\alpha }^2(t)$ signifies the intensity of the white noise. Hence we obtain the generalisation of the integral form of (6.1) which is called the stochastic logistic equation, given as

6.2

$v(t)=v_0+{\int }_0^{t}r(s)v(s)(1-\frac{v(s)}{k})ds+{\int }_0^t\alpha (s)v(s)(1-\frac{v(s)}{k})dB(s),$ where $B(t)=\xi (t)dt$ is referred to as Brownian motion [17], as previously discussed. The stability of (6.2) is explored by Liu and Wang [44]. In this section, an approximation to the case with $x_0=0.3$, $r(t)=0.2$, $k=1$, and two different values of α, i.e., $\alpha (t)=0.3+0.3sin(t)$ and $\alpha (t)=0$ is offered for $m=16$ and $m=32$ [45] are illustrated in Figs. 5 and 6 respectively. Tables 5 and 6 give numerical data of the approximate solution of the model with $\alpha (t)=0.3+0.3sin(t)$ and $\alpha (t)=0$ for $m=16$ and $m=32$ respectively.

[See PDF for image]

Figure 5

Performance evaluation of the presented method applied to the logistic model (Sect. 6.1) using $m=16$

[See PDF for image]

Figure 6

Performance evaluation of the presented method applied to the logistic model (Sect. 6.1) using $m=32$

Table 5. Approximate solution of model (Sect. 6.1) with and without the intensity of white noise with the mean error ${\overline{x}}_E$ for $m=16$ with 50 iterations

t	α(t)=0.3 + 0.3sin(t)	α(t)=0	${\overline{x}}_E$
0.1	0.303863369	0.303953563	0.90194399 × 10−4
0.3	0.311662877	0.311941640	2.78762999 × 10−4
0.5	0.321567730	0.321537389	0.30341800 × 10−4
0.7	0.331301584	0.331027443	2.74140600 × 10−4
0.9	0.339154619	0.339373907	2.19287800 × 10−4

Table 6. Approximate solution of model (Sect. 6.1) with and without the intensity of white noise with the mean error ${\overline{x}}_E$ for $m=32$ with 50 iterations

t	α(t)=0.3 + 0.3sin(t)	α(t)=0	${\overline{x}}_E$
0.1	0.304759034	0.304514589	0.2444442 × 10−3
0.3	0.313540927	0.312522010	1.0189104 × 10−3
0.5	0.320844466	0.321438906	0.5944398 × 10−3
0.7	0.328109289	0.330237674	2.1283854 × 10−3
0.9	0.336252933	0.338550318	2.2973844 × 10−3

Conclusions

We presented a numerical method for solving nonlinear stochastic integral equations (NLSIEs) using the stochastic operational matrix derived from Itô integration of Walsh functions. To solve the integral equation, first the integral equation is converted to an algebraic system of equations, which has then been solved to determine an approximate solution to the problem. To demonstrate the applicability of the recommended method, a convergence study has been conducted. It was observed that the order of convergence of the method came to be $\mathcal{O}(h)$. A couple of examples are solved to support the authenticity of the method and the CPU time proves that the method takes much less time to solve the problem.

The study concludes with numerical solutions for the canonical stochastic logistic model in ecology, comparing scenarios with varying white noise intensities against the deterministic case. Based upon the data obtained, the method proved to be an important tool to solve such types of integral equations, like the non-linear pendulum problem, the Lotka-Volterra model, and many more. This method can be further extended to address multidimensional stochastic integral equations, including those driven by fractional Brownian motion. Based on the data obtained, the technique may prove to be an effective tool for solving several NLSIEs, including the Lotka-Volterra model, the nonlinear pendulum problem, and many more. This approach can be further refined to tackle multidimensional stochastic integral equations, including those governed by fractional Brownian motion.

Author contributions

All five authors have accepted responsibility for the entire content of this manuscript and approved its submission.

Data availability

All data generated or analyzed during this study are included in this article.

Declarations

Ethics approval and consent to participate

Not Applicable.

Consent for publication

Not Applicable.

Competing interests

The authors declare no competing interests.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Øksendal, B. Stochastic Differential Equations: An Introduction with Applications; 1998; 5 New York, Springer: [DOI: https://dx.doi.org/10.1007/978-3-662-03620-4]

2. Evans, L.C. An Introduction to Stochastic Differential Equations; 2012; Providence, Am. Math. Soc.:

3. Kloeden, P.E.; Platen, E. Numerical Solution of Stochastic Differential Equations; 1999; Berlin, Springer:

4. Zimmer, C.; Sahle, S. A termination criterion for parameter estimation in stochastic models in systems biology. Biosystems; 2015; 137, pp. 55-63. [DOI: https://dx.doi.org/10.1016/j.biosystems.2015.08.003]

5. Szekely, T., Jr; Burrage, K. Stochastic simulation in systems biology. Comput. Struct. Biotechnol. J.; 2014; 12, 20–21 pp. 14-25. [DOI: https://dx.doi.org/10.1016/j.csbj.2014.10.003]

6. Amar, P.; Paulevé, L. Hsim: a hybrid stochastic simulation system for systems biology. Electron. Notes Theor. Comput. Sci.; 2015; 313, pp. 3-21. [DOI: https://dx.doi.org/10.1016/j.entcs.2015.04.016]

7. Sabbar, Y.; Izadi, M.; Raezah, A.A.; Adel, W. Nonlinear dynamics of a general stochastic SIR model with behavioral and physical changes: analysis and application to zoonotic tuberculosis. Mathematics; 2024; 12, 13 [DOI: https://dx.doi.org/10.3390/math12131974] 1974.

8. Jenny, P.; Torrilhon, M.; Heinz, S. A solution algorithm for the fluid dynamic equations based on a stochastic model for molecular motion. J. Comput. Phys.; 2010; 229, 4 pp. 1077-1098.2576239 [DOI: https://dx.doi.org/10.1016/j.jcp.2009.10.008]

9. Peavoy, D.; Franzke, C.L.; Roberts, G.O. Systematic physics constrained parameter estimation of stochastic differential equations. Comput. Stat. Data Anal.; 2015; 83, pp. 182-199.3281805 [DOI: https://dx.doi.org/10.1016/j.csda.2014.10.011]

10. Frank, T. Kramers-Moyal expansion for stochastic differential equations with single and multiple delays: applications to financial physics and neuro-physics. Phys. Lett. A; 2007; 360, 4–5 pp. 552-562. [DOI: https://dx.doi.org/10.1016/j.physleta.2006.08.062]

11. Zimmer, C.; Sahle, S. Comparison of approaches for parameter estimation on stochastic models: generic least squares versus specialized approaches. Comput. Biol. Chem.; 2016; 61, pp. 75-85. [DOI: https://dx.doi.org/10.1016/j.compbiolchem.2015.10.003]

12. dos Reis, T.M.; Neto, A.C.; Campos, M.C.M. Gaussian basis sets for atomic and molecular calculations obtained from stochastic optimization. Comput. Theor. Chem.; 2017; 1099, pp. 133-139. [DOI: https://dx.doi.org/10.1016/j.comptc.2016.11.026]

13. Manninen, T.; Linne, M.L.; Ruohonen, K. Developing Itô stochastic differential equation models for neuronal signal transduction pathways. Comput. Biol. Chem.; 2006; 30, 4 pp. 280-291. [DOI: https://dx.doi.org/10.1016/j.compbiolchem.2006.04.002]

14. Laurini, M.P.; Caldeira, J.F. A macro-finance term structure model with multivariate stochastic volatility. Int. Rev. Econ. Finance; 2016; 44, pp. 68-90. [DOI: https://dx.doi.org/10.1016/j.iref.2016.03.008]

15. Platen, E.; Bruti-Liberati, N. Numerical Solution of Stochastic Differential Equations with Jumps in Finance; 2010; Berlin, Springer: [DOI: https://dx.doi.org/10.1007/978-3-642-13694-8]

16. Walsh, J.L. A closed set of normal orthogonal functions. Am. J. Math.; 1923; 45, 1 pp. 5-24.1506485 [DOI: https://dx.doi.org/10.2307/2387224]

17. Paikaray, P.P.; Parida, N.C.; Beuria, S. Numerical approximation of p-dimensional stochastic Volterra integral equation using Walsh function. J. Math. Comput. Sci.; 2023; 31, 4 pp. 448-464. [DOI: https://dx.doi.org/10.22436/jmcs.031.04.07]

18. Paikaray, P.P.; Beuria, S.; Parida, N.C. Numerical approximation of stochastic Volterra-Fredholm integral equation using Walsh function. Commun. Math. Appl.; 2023; 14, 5 pp. 1603-1613. [DOI: https://dx.doi.org/10.26713/cma.v14i5.2313]

19. Asgari, M.; Shekarabi, F.H. Numerical solution of nonlinear stochastic differential equations using the block pulse operational matrices. Math. Sci.; 2013; 7, 1 pp. 1-12.3877778 [DOI: https://dx.doi.org/10.1186/2251-7456-7-47]

20. Hashemi, S.A.S.; Saeedi, H. ADM-TF hybrid method for nonlinear Itô-Volterra integral equations. Math. Comput. Simul.; 2021; 185, pp. 783-798. [DOI: https://dx.doi.org/10.1016/j.matcom.2021.02.003]

21. Wu, J.; Jiang, G.; Sang, X. Numerical solution of nonlinear stochastic Itô-Volterra integral equations based on Haar wavelets. Adv. Differ. Equ.; 2019; 2019, 1 [DOI: https://dx.doi.org/10.1186/s13662-019-2440-6] 503.

22. Yadav, P.; Jahan, S.; Izadi, M. Taylor wavelet quasilinearization method for solving tumor growth model of fractional order. Results Control Optim.; 2024; 15, [DOI: https://dx.doi.org/10.1016/j.rico.2024.100437] 100437.

23. Zeghdane, R. Numerical approach for solving nonlinear stochastic Itô-Volterra integral equations using shifted Legendre polynomials. Int. J. Dyn. Syst. Differ. Equ.; 2021; 11, 1 pp. 69-83.4249576

24. Ahmadinia, M.; Afshariarjmand, H.; Salehi, M. Numerical solution of Itô-Volterra integral equations by the QR factorization method. J. Appl. Math. Comput.; 2023; 69, 4 pp. 3171-3188.4617132 [DOI: https://dx.doi.org/10.1007/s12190-023-01873-9]

25. Izadi, M.; Seifaddini, M.; Afshar, M. Approximate solutions of a SIR epidemiological model of computer viruses. Adv. Stud. Euro-Tbilisi Math. J.; 2021; 14, 4 pp. 203-219.4425168 [DOI: https://dx.doi.org/10.32513/asetmj/1932200822]

26. Singh, P.K.; Ray, S.S. A numerical approach based on Pell polynomial for solving stochastic fractional differential equations. Numer. Algorithms; 2024; 97, 4 pp. 1513-1534.4832578 [DOI: https://dx.doi.org/10.1007/s11075-024-01760-9]

27. Izadi, M.; Srivastava, H.M.; Mamehrashi, K. Numerical simulations of Rosenau-Burgers equations via Crank-Nicolson spectral Pell matrix algorithm. J. Appl. Math. Comput.; 2025; 71, 1 pp. 1009-1033.4854967 [DOI: https://dx.doi.org/10.1007/s12190-024-02273-3]

28. Barzegar, N.; Mirzaee, F.; Solhi, E. New technique based on Vieta-Lucas polynomials for solving nonlinear stochastic Itô-Volterra integral equation. Int. J. Numer. Model., Electron. Netw. Devices Fields; 2025; 38, 3 [DOI: https://dx.doi.org/10.1002/jnm.70044] e70044.

29. Mirzaee, F.; Naserifar, S.; Solhi, E. Meshless Barycentric rational interpolation method for solving nonlinear stochastic fractional integro-differential equations. Iran. J. Sci.; 2024; 48, 3 pp. 709-733.4753576 [DOI: https://dx.doi.org/10.1007/s40995-024-01621-z]

30. Mirzaee, F.; Solhi, E.; Naserifar, S. Approximate solution of stochastic Volterra integro-differential equations by using moving least squares scheme and spectral collocation method. Appl. Math. Comput.; 2021; 410, 4276454 126447.

31. Solhi, E.; Mirzaee, F.; Naserifar, S. Enhanced moving least squares method for solving the stochastic fractional Volterra integro-differential equations of Hammerstein type. Numer. Algorithms; 2024; 95, 4 pp. 1921-1951.4719437 [DOI: https://dx.doi.org/10.1007/s11075-023-01633-7]

32. Mirzaee, F.; Alipour, S. Cubic B-spline approximation for linear stochastic integro-differential equation of fractional order. J. Comput. Appl. Math.; 2020; 366, 4000739 [DOI: https://dx.doi.org/10.1016/j.cam.2019.112440] 112440.

33. Pearl, R.; Reed, L.J. On the rate of growth of the population of the United States since 1790 and its mathematical representation. Proc. Natl. Acad. Sci. USA; 1920; 6, 6 pp. 275-288. [DOI: https://dx.doi.org/10.1073/pnas.6.6.275]

34. Pielou, E.C. An Introduction to Mathematical Ecology; 1975; New York, Wiley-Interscience:

35. Volterra, V. Variations and fluctuations of the number of individuals in animal species living together. J. Cons. Int. Explor. Mer; 1928; 3, pp. 1-51. [DOI: https://dx.doi.org/10.1093/icesjms/3.1.3]

36. Izadi, M.; Srivastava, H.M. A discretization approach for the nonlinear fractional logistic equation. Entropy; 2020; 22, 114222131 [DOI: https://dx.doi.org/10.3390/e22111328] 1328.

37. Izadi, M. Comparison of various fractional basis functions for solving fractional-order logistic population model. Facta Univ., Ser. Math. Inform.; 2020; 35, 4 pp. 1181-1198.4246554

38. Izadi, M. A comparative study of two Legendre-collocation schemes applied to fractional logistic equation. Int. J. Appl. Comput. Math.; 2020; 6, 34094525 [DOI: https://dx.doi.org/10.1007/s40819-020-00823-4] 71.

39. Izadi, M.; Srivastava, H.M. Numerical approximations to the nonlinear fractional-order logistic population model with fractional-order Bessel and Legendre bases. Chaos Solitons Fractals; 2021; 145, 4228359 [DOI: https://dx.doi.org/10.1016/j.chaos.2021.110779] 110779.

40. Limaye, B.V. Functional Analysis; 1996; New Age International

41. Rao, G.P. Piecewise Constant Orthogonal Functions and Their Application to Systems and Control; 1983; Berlin, Springer: [DOI: https://dx.doi.org/10.1007/BFb0041228]

42. Verhulst, P.F. Notice sur la loi que la population suit dans son accroissement. Corresp. Math. Phys.; 1838; 10, pp. 113-121.

43. May, R.M. Stability and Complexity in Model Ecosystems; 2001; Princeton, Princeton University Press: [DOI: https://dx.doi.org/10.1515/9780691206912]

44. Liu, M.; Wang, K. A note on stability of stochastic logistic equation. Appl. Math. Lett.; 2013; 26, 6 pp. 601-606.3028061 [DOI: https://dx.doi.org/10.1016/j.aml.2012.12.015]

45. Heydari, M.H.; Cattani, C.; Razzaghi, M. An efficient computational method for solving nonlinear stochastic Itô integral equations: application for stochastic problems in physics. J. Comput. Phys.; 2014; 270, pp. 402-415.3209393 [DOI: https://dx.doi.org/10.1016/j.jcp.2014.03.064]