Acta Polytechnica https://doi.org/10.14311/AP.2024.64.0128 Acta Polytechnica 64(2):128–141, 2024 © 2024 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague NUMERICAL SOLUTION FOR STOCHASTIC VOLTERRA-FREDHOLM INTEGRAL EQUATIONS WITH DELAY ARGUMENTS Kutorzi Edwin Yaoa,b, Yuxue Zhanga,b, Yufeng Shia,b,∗ a Shandong University, Institute for Financial Studies, 250100 Jinan, China b Shandong University, School of Mathematics, 250100 Jinan, China ∗ corresponding author: yfshi@sdu.edu.cn Abstract. We present a method for computing the stochastic operational matrix of integration to advance the study of stochastic Volterra-Fredholm integral equations (SVFIEs) based on delay arguments. First, the method evaluates the combined effects of the delay and its parameters on the accuracy improvement of the convergence rate. Our results can be applied to SVFIEs, with the operational delay matrices of the block pulse function simplified to algebraic ones. Numerical calculations were performed on a PC using Python 3 programs. Results also demonstrate the accuracy of approximate solutions; arithmetic operations are carried out without the need for derivation or integration. Keywords: Stochastic Volterra-Fredholm integral equations, block-pulse functions, Itô integral, delay operational matrix, error analysis. 1. Introduction Stochastic Volterra-Fredholm integral equations are an essential class of multi-dimensional integral equations that can rarely be solved exactly, and the computational complexity of mathematical operations is a critical obstacle in solving high-dimensional stochastic integral equations. Stochastic differential equations have various applications in various fields, such as medicine, economics, and social sciences, as well as engineering, biology, and financial mathematics. These equations play a crucial role in modelling population growth, where the stochastic Volterra-Fredholm integral equation is fundamental; see [1–7]. A stochastic Volterra-Fredholm integral equations can be modeled using several types of stochastic differential equations or, in more complicated cases, nonlinear stochastic differential equations of the Itô type [8–10]. There are some difficulties in finding exact solutions for SVIEs or SVFIEs, so the researchers have resorted to finding approximate solutions using numerical methods [11, 12]. SVFIEs with delay are used in applied sciences for modelling functions that contain time memory, such as mechanical systems, dynamical systems, and electric circuits, as well as in physical models, option pricing, and population growth [13]. On the one hand, in the theory of automatic systems, delay-differential equations are obtained [14–16]. On the other hand, some systems, such as the integral equations [17, 18], which were based on the operational matrices of integration, were estimated using polynomials. These included block pulse systems, the Fourier series, Legendre polynomials, Chebyshev polynomials, and Laguerre polynomials. Using numerical methods to approximate the solutions to such equations is often desirable since they cannot always be solved explicitly [19–27]. The Volterra integral equations with delay have received very little attention. We have developed approximation methods for SVFIEs with delay arguments. A stochastic operational matrix with time delay is presented to find an approximate solution of the Stochastic Volterra-Fredholm integral equations. Our focus is on the SVFIE: X(t) = f(t) + λ1 ∫ β α k1(t, s)X(s − τ)ds + λ2 ∫ t 0 k2(t, s)X(s − τ)ds + λ3 ∫ t 0 k3(t, s)X(s − τ)dB(s), where t ∈ [0, T ), τ ∈ [α, β], τ ∈ [0, t). In the above descriptions, X(t), f(t), k1(t, s), k2(t, s) and k3(t, s), for t, s ∈ [0, T ), are the stochastic processes defined on the same probability space (Ω, F ,P), and X(t) is unknown. B(t) is a one-dimensional standard Brownian motion process and ∫ t 0 k3(t, s)X(s − τ)dB(s) is the Itô integral. Both the j and k represent the Volterra kernel. The parameter in the variable of the function is calculated as τ = (q + λ)h with an integer q ≥ 0 and a fraction 0 ≤ λ < 1, chosen to approximate a function with a time delay. Following is an outline of the paper. A description of the fundamental properties of block-pulse functions is provided in Section 2, as well as the approximation of functions using block-pulse parts and an operational 128 https://doi.org/10.14311/AP.2024.64.0128 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en vol. 64 no. 2/2024 Stochastic Volterra-Fredholm with Delay integration matrix. A stochastic integration functional matrix is introduced in Section 3. The stochastic integration active matrix is used to solve stochastic delay Volterra integral equations in Section 4. We present the error estimation and rate of convergence in Section 5. The proposed scheme is accurate, which is proved by using numerical examples to demonstrate its effectiveness in Section 6. A brief conclusion is given in Section 7. 2. Block-pulse functions (BPFs) This section covers the notations, definitions, known results, and formulas related to BPFs, which are relevant to this paper. These details have been extensively discussed in [20, 21]. The block-pulse functions (BPF) Φi over the unit interval [0, 1) is defined as follows: for 0 ≤ i < m, and m ∈ {1, 2, . . .}: Φi(t) = { 1 (i − 1)h ≤ t < ih, 0 otherwise, (1) with t ∈ [0, T ), i = 1, 2, . . . , m, and h = T m . The block-pulse functions have the following properties: (1.) Disjointness: The BPFs are disjointed with each other in the interval t ∈ [0, T ): Φi(t)Φj(t) = δijΦi(t), (2) where i, j = 1, 2, . . . , m, and δij denotes the Kronecker delta. (2.) Orthogonality: The BPFs are orthogonal with each other in the interval t ∈ [0, T ):∫ T 0 Φi(t)Φj(t)dt = hδij , (3) where i, j = 1, 2, . . . , m. (3.) The third property is completeness: For every f ∈ L2[0, T ), when m → ∞ , Parseval’s identity holds, that is: ∫ T 0 f2(t)dt = ∞∑ i=1 f2 i ||Φi(t)||2, where fi = 1 h ∫ T 0 f(t)Φi(t)dt. The set of functions can be described by an m vector: Φ(t) = (Φ0(t), Φ1(t), · · · , Φm(t))T , where t ∈ [0, T ). Thus, we can write the relationship between BPFs and their integrals in the following matrix form. The above representation and disjointness property follows: Φ(t)ΦT (t) =  Φ1(t) 0 · · · 0 0 Φ2(t) · · · 0 ... ... . . . ... 0 0 · · · Φm(t)  m×m , (4) additionally, we deduce: ΦT (t)Φ(t) = 1, and: Φ(t)ΦT (t)F T = DF Φ(t), (5) the diagonal matrix DF corresponds to a constant vector F = (f1, f2, · · · , fm)T whose diagonal entries are related. 129 E. Y. Kutorzi, Y. Zhang, Y. Shi Acta Polytechnica 2.1. Functions approximation A real bounded function f(t), which f(t) ∈ L2[0, T ), can be expanded into a block pulse series as: f(t) ≃ f̂m(t) = m∑ i=1 fiΦi(t), (6) where fi is the block pulse coefficient with respect to the ith BPF Φi(t). The vector form is as follows: f(t) ≃ f̂m(t) = F T Φ(t) = ΦT (t)F, (7) where F = (f1, f2, · · · , fm)T . Let k(t, s) ∈ L2([0, T1) × [0, T2)). Similarly, it can be applied to BPFs such as: k(t, s) ≃ k̂m(t, s) = ΨT (s)KΦ(t) = ΦT (t)KT Ψ(s), (8) where Φ(t) and Ψ(s) are m1 and m2 dimensional BPFs vectors, respectively, and K = (kij), i = 1, 2, . . . , m1, j = 1, 2, . . . , m2 is the m1 × m2 block pulse coefficient matrix with: kij = 1 h1h2 ∫ T1 0 ∫ T2 0 k(t, s)Ψi(t)Φj(s)dsdt, where h1 = T1 m1 , h2 = T2 m2 . For convenience, we put m1 = m2 = m. 2.2. Integration operational matrix Computing ∫ t 0 Φi(s)ds follows: ∫ t 0 Φi(s)ds =  0 0 ≤ t < (i − 1)h, t − (i − 1)h (i − 1)h ≤ t < ih, h ih ≤ t < T. (9) Note that t − (i − 1)h, equals to h 2 at mid-point of [(i − 1)h, ih), thus we can approximate t − (i − 1)h, for (i − 1)h ≤ t < ih, by h 2 . From [20], we have: ∫ t 0 Φ(s)ds = PΦ(t). (10) As shown in the operational matrix of integration: P = h 2  1 2 2 · · · 2 0 1 2 · · · 2 0 0 1 · · · 2 ... ... ... . . . ... 0 0 0 · · · 1  m×m (11) Accordingly, each integral of f(t) can be approximated as follows:∫ t 0 f(s)ds ≃ ∫ t 0 F T Φ(s)ds ≃ F T PΦ(t). (12) 2.3. The operational matrix with time delay of BPFs The delay time is τ = (q + λ)h with an integer q ≥ 0, and a fraction 0 ≤ λ < 1, where the operational matrix of approximation is expressed as the time delay τ = qh, yields: ϕi(t − qh) = { ϕi+q(t) i ≤ m − q, 0 i > m − q, (13) and the function containing time delay, yields: ϕi(t − τ) =  ϕi+q(t) + ϕλ(t − (i + q)h) − ϕλ(t − (i + q − 1)h) i < m − q, ϕi+q(t) − ϕλ(t − (i + q − 1)h) i = m − q, 0 i > m − q, (14) 130 vol. 64 no. 2/2024 Stochastic Volterra-Fredholm with Delay alternatively, as vectors: ϕi(t − τ) = ∆T i HqΦ(t) − ∆T i HqΦλ(t) + ∆T i Hq+1Φλ(t). We expand the function ϕi(t − τ) into its block pulse series to avoid the expression Φλ(t) in the above equation: ϕi(t − τ) = (ci,1 ci,2 · · · ci,m)Φ(t) (15) While ci,j (j = 1, 2 · · · , m) are: ci,j = 1 h ∫ T 0 ϕi(t − τ)ϕj(t)dt = 1 h ∫ jh (j−1)h ϕi(t − τ)dt = 1 h ∆T i Hq (∫ jh (j−1)h Φ(t)dt − ∫ jh (j−1)h Φλ(t)dt + H ∫ jh (j−1)h Φλ(t)dt ) = ∆T i ((1 − λ)Hq + λHq+1)∆j . (16) We can develop the whole block pulse function vector containing time delay τ = (q + λ)h into its block pulse series in a vector form by noting that the expression ∆T i ((1 − λ)Hq + λHq+1)∆j is just one entry of the matrix ((1 − λ)Hq + λHq+1) with ith row and jth column: Φ(t − τ) = ((1 − λ)Hq + λHq+1)Φ(t). (17) Usually, the matrix (1 − λ)Hq + λHq+1 is referred to as the delay operational matrix. To put it in more concrete terms: (1 − λ)Hq + λHq+1 =  (q + 1)th︸ ︷︷ ︸ 0 · · · 0 1 − λ λ 0 · · · 0 0 · · · 0 0 1 − λ λ · · · 0 ... · · · ... ... ... ... . . . ... 0 · · · 0 0 0 0 · · · λ 0 · · · 0 0 0 0 · · · 1 − λ 0 · · · 0 0 0 0 · · · 0 ... · · · ... ... ... ... · · · ... 0 · · · 0 0 0 0 · · · 0  . (18) It is possible to obtain the block pulse series of a function with time delay τ = (q + λ)h by using Equation (17): f(t − τ) ≃ F T Φ(t − τ) = F T ((1 − λ)Hq + λHq+1)Φ(t). (19) 3. Stochastic integration operational matrix The integral of Itô of a single BPF ϕi(t) can be computed as follows: ∫ t 0 ϕi(s)dB(s) =  0 0 ≤ t < (i − 1)h, B(t) − B((i − 1)h), (i − 1)h ≤ t < ih, B(ih) − B((i − 1)h), ih ≤ t < T. (20) Now expressing ∫ t 0 ϕi(s)dB(s), in terms of the BPFs follows:∫ t 0 ϕi(s)dB(s) ≃ (B(ih/2) − B(i − 1)h/2) ϕi(t) + (B(ih) − B((i − 1)h)) m∑ j=i+1 ϕj(t). (21) Therefore: ∫ t 0 Φ(s)dB(s) ≃ PSΦ(t). (22) In this case, the stochastic operational matrix of integration can be expressed as follows: PS =  γ1 ρ1 ρ1 · · · ρ1 0 γ2 ρ2 · · · ρ2 0 0 γ3 · · · ρ3 ... ... ... . . . ... 0 0 0 · · · γm  m×m , (23) 131 E. Y. Kutorzi, Y. Zhang, Y. Shi Acta Polytechnica where ρi = B(ih) − B((i − 1)h), i = 1, 2, . . . , m − 1; γj = B(ih/2) − B((i − 1)h/2), j = 1, 2, . . . , m. This can be approximated by computing the Itô integral for every function f(t) as follows:∫ t 0 f(s)dB(s) ≃ F T Φ(s)dB(s) ≃ F T PSΦ(t). (24) 4. Solving stochastic Volterra-Fredholm integral equations with time delay The following linear stochastic Volterra-Fredholm integral equation is considered with a constant time delay τ > 0: X(t) = f(t) + λ1 ∫ β α k1(t, s)X(s − τ)ds + λ2 ∫ t 0 k2(t, s)X(s − τ)ds + λ3 ∫ t 0 k3(t, s)X(s − τ)dB(s), (25) where t ∈ [α, β], τ ∈ (0, β − α), t ∈ [0, T ). X(t) is a stochastic process whose coefficients are X(t), f(t), k1(t, s), k2(t, s) and k3(t, s), for α, β ∈ [α, β], t, s ∈ [0, T ), defined on the same probability space (Ω, F , P ). In addition, B(t) is a Brownian motion process, and ∫ t 0 k3(t, s)X(s − τ)dB(s) is integral for Itô. To facilitate block pulse functions, we typically set α = 0. Whenever α ̸= 0 we set s = t−α β−α T , where T = mh. Using BPFs to approximate functions X(t), f(t), K1(t, s), K2(t, s), and K3(t, s) by Equations (7), (8), and (19) gives the following result: X(t) ≃ XT Φ(t) = ΦT X, f(t) ≃ F T Φ(t) = ΦT F, k1(t, s) ≃ ΨT (t)K1Φ(s) = ΦT (s)KT 1 Ψ(t), k2(t, s) ≃ ΨT (t)K2Φ(s) = ΦT (s)KT 2 Ψ(t), k3(t, s) ≃ ΨT (t)K3Φ(s) = ΦT (s)KT 3 Ψ(t). According to Equation (19), X(s − τ) can be approximated as follows: X(s − τ) ≃ XT Ψ(s − τ) = XT ((1 − λ)Hq + λHq+1)Φ(t)Ψ(s), and by letting Z = ((1 − λ)Hq + λHq+1), we can write: X(s − τ) ≃ XT ZΨ(s). (26) The above approximates define X and F as stochastic block pulse coefficients vectors, respectively, and K1, K2, and K3 as stochastic block pulse coefficients matrices. Equation (25) is improved by substituting the above approximation: XT Φ(t) ≃ F T Φ(t) + XT Z(λ1 ∫ mh 0 Ψ(s)ΨT (s)ds)K1Φ(t) + XT Z(λ2 ∫ t 0 Ψ(s)ΨT (s)ds)K2Φ(t) + XT Z(λ3 ∫ t 0 Ψ(s)ΨT (s)dBs)K3Φ(t). (27) Let Ki j , j = 1, 2, 3, be the ith row of the constant matrix Kj , for j = 1, 2, 3. Ri be the ith row of the integration operational matrix P, Ri S be the ith row of the stochastic integration operational matrix PS , DKi j , be diagonal matrices with Ki j , for j = 1, 2, 3, as its diagonal entries. By the relation ∫mh 0 ΦT Φ(s)ds = hI, previous relations, and assuming m1 = m2, we have:(∫ mh 0 Ψ(s)ΨT (s)ds ) K1Φ(t) = hIK1Φ(t) = B1Φ(t), (28) 132 vol. 64 no. 2/2024 Stochastic Volterra-Fredholm with Delay where B1 = hK1, also:(∫ t 0 Ψ(s)ΨT (s)ds ) K1Φ(t) = (∫ t 0 Φ(s)ΦT (s)ds ) K2Φ(t) =  R1Φ(t)K1 2 Φ(t) R2Φ(t)K2 2 Φ(t) ... RmΦ(t)Km 2 Φ(t)  =  R1DK1 2 R2DK2 2 ... RmDKm 2 Φ(t) = B2Φ(t), (29) where: B2 = h 2  k2 11 2k2 12 2k2 13 · · · 2k2 1m 0 k2 22 2k2 23 · · · 2k2 2m 0 0 k2 33 · · · 2k2 3m ... ... ... . . . ... 0 0 0 · · · k2 mm  m×m , (30) also, we can consider the integral term Itô:(∫ t 0 Ψ(s)ΨT (s)dB(s) ) K3Φ(t) = (∫ t 0 Φ(s)ΦT (s)dB(s) ) K3Φ(t) =  R1Φ(t)K1 3 Φ(t) R2Φ(t)K2 3 Φ(t) ... Rm S Φ(t)Km 3 Φ(t)  =  R1 SDK1 3 R2 SDK2 3 ... Rm S DKm 3 Φ(t) = B3Φ(t), (31) where: B3 =  k3 11γ k3 12ρ k3 13ρ · · · k3 1mρ 0 k3 22γ k3 23ρ · · · k3 2mρ 0 0 k3 33γ · · · k3 3mρ(m − 2) ... ... ... . . . ... 0 0 0 · · · k3 mmγ(m − 1)  . (32) By substituting Equations (28), (29) and (31) in (27), we get: XT Φ(t) ≃ F T Φ(t) + XT Zλ1B1Φ(t) + Y T Zλ2B2Φ(t) + XT Zλ3B3Φ(t). Then: XT (I − Z(λ1B1 + λ2B2 + λ3B3)) ≃ F T . So, by setting M = (I − Z(λ1B1 + λ2B2 + λ3B3)) and replacing ≃ by =, we deduce: MT X = F. (33) It consists of a linear system of equations with lower triangular coefficients that yields the approximate block pulse coefficient of the stochastic process X(t). 5. Error estimation and rate of convergence The proposed method shows the fastest convergence rate for integral equations with time delay. There is a high-level agreement between the exact solution and numerical results. Theorem 1. Let f(t) be any arbitrary real bounded function, which is square integrable within the interval [0, 1), and e(t) = f(t) − f̂m(t), t ∈ I = [0, 1), where f̂m(t) = ∑m i=1 fiϕi(t) is the block pulse series of f(t). Then: ∥e(t)∥ ≤ h 2 √ 3 ∥f ′∥∞ , (34) in this case, ∥e(t)∥ = (∫ 1 0 |e(t)|2 dt ) 1 2 . 133 E. Y. Kutorzi, Y. Zhang, Y. Shi Acta Polytechnica Proof. See [23]. Theorem 2. Assume f(t, s) ∈ L2([0, 1) × [0, 1)) and e(t, s) = f(t, s) − f̂m(t, s), (t, s) ∈ A = [0, 1) × [0, 1), which f̂m(t, s) = ∑m i=1 ∑m j=1 fijΨi(t)Φj(s) is the block pulse series of f(t, s). Then: ∥e(t, s)∥ ≤ h 2 √ 3 ( ∥f ′ t∥ 2 ∞ + ∥f ′ t∥ 2 ∞ ) 1 2 , (35) where ∥e(t, s)∥ = (∫ 1 0 ∫ 1 0 |e(t, s)|2dsdt ) 1 2 . Proof. Let: eij(t, s) = { f(t, s) − fij (t, s) ∈ Aij , 0 (t, s) ∈ A − Aij , (36) where Aij = {(t, s) : (i − 1)h ≤ t < ih, (j − 1)h ≤ s < jh, h = 1 m }, and i, j = 1, 2, · · · , m. For i, j = 1, 2, · · · , m, thus, we get: eij(t, s) = f(t, s) − 1 h2 ∫ ih (i−1)h ∫ jh (j−1)h f(x, y)dydx = 1 h2 ∫ ih (i−1)h ∫ jh (j−1)h (f(t, s) − f(x, y)) dydx, now, by mean-value theorem, we deduce: eij(t, s) = 1 h2 ∫ ih (i−1)h ∫ jh (j−1)h ((t − x)f ′ t(ηi, ηj) + (s − y)f ′ s(ηi, ηj)) dydx = f ′ t(ηi, ηj) ( t + ( −i + 1 2 ) h ) + f ′ s(ηi, ηj) ( s + ( −j + 1 2 ) h ) , where (t, s), (ηi, ηj) ∈ Aij ; then: ∥eij(t, s)∥2 = ∫ ih (i−1)h ∫ jh (j−1)h |eij(t, s)|2 dsdt = h4 12(f ′2 t (ηi, ηj) + f ′2 s (ηi, ηj)), (37) where (ηi, ηj) ∈ Aij , i, j = 1, 2, . . . , m. Consequently, we have: ∥e(t, s)∥2 = ∫ 1 0 ∫ 1 0 |e(t, s)|2 dsdt = ∫ 1 0 ∫ 1 0  m∑ i=1 m∑ j=1 eij(t, s) 2 dsdt = m∑ i=1 m∑ j=1 ∫ 1 0 ∫ 1 0 e2 ij(t, s)dsdt = m∑ i=1 m∑ j=1 ∥eij(t, s)∥2 = h4 12 m∑ i=1 m∑ j=1 (f ′2 t (ηi, ηj) + f ′2 s (ηi, ηj)) ≤ h2 12 ( sup (x,y)∈A |f ′ t(x, y)|2 + sup (x,y)∈A |f ′ s(x, y)|2 ) , (38) or: ∥e(t, s)∥ ≤ h 2 √ 3 ( ∥f ′ t∥ 2 ∞ + ∥f ′ s∥2 ∞ ) 1 2 , hence, ∥e(s, t)∥ = O(h). □ Theorem 3. Let X(t) and X̂(t) be solutions of Equations (25) and (26), respectively, and let ∥X(t)∥ < C and ∥ki∥ < C for i = 1, 2, 3. Then: E( ∥∥∥X(t) − X̂(t) ∥∥∥2 ) ≤ O(h2), where t ∈ [0, T ), τ ∈ [0, 1); and: sup 0≤τ 0: X(t) = − t4 12 + t3 3 τ + (1 − τ2 2 )t2 + ∫ t 0 (t − s)X(s − τ)ds, (39) where s, t ∈ [0, T ], τ ∈ (0, T ); with the exact solution X(t) = t2, for 0 ≤ t ≤ T . In Tables 1–2, the numerical results are presented. The computations of mean, standard deviation, and mean confidence interval of error for n, χ̄E , and SE are provided in Tables 1–2. An approximate solution is depicted in Figure 1 as a trajectory based on the presented approach. The variation process of error is represented by curves in Figure 2. We observe a perfect agreement between the exact solution and the numerical results, achieving full convergence. Example 2 [16]. Consider the following Fredholm integral equation with (constant) time delay τ > 0: X(t) = t(T cos(T − τ) − sin(T − τ) − sin(τ)) + sin(t) + ∫ T 0 (ts)X(s − τ)ds, (40) where s, t ∈ [0, T ], τ ∈ (0, T ); with the exact solution X(t) = sin(t), for 0 ≤ t ≤ T . The numerical results are presented in Table 3. The trajectory of the approximate solution and exact solution are represented in Figures 3–7. 137 E. Y. Kutorzi, Y. Zhang, Y. Shi Acta Polytechnica Figure 2. Variation trend of error in Example 1 for m = 32, n = 50, n = 100, q = 0, λ = 0.5. λ = 0.1 λ = 0.3 λ = 0.5 λ = 0.7 λ = 0.9 m = 8 0.007788 0.007251 0.006711 0.00727 0.008367 m = 32 0.015303 0.015683 0.016064 0.016446 0.016827 m = 64 0.017717 0.017918 0.01812 0.018321 0.018523 Table 3. Error in Example 2 with q = 0. Figure 3. The trajectory of the approximate solution and exact solution of Example 2 for m = 32, m = 64, n = 50, q = 0, λ = 0.1. Figure 4. The trajectory of the approximate solution and exact solution of Example 2 for m = 32, m = 64, n = 50, q = 0, λ = 0.3. 138 vol. 64 no. 2/2024 Stochastic Volterra-Fredholm with Delay Figure 5. The trajectory of the approximate solution and exact solution of Example 2 for m = 32, m = 64, n = 50, q = 0, λ = 0.5. Figure 6. The trajectory of the approximate solution and exact solution of Example 2 for m = 32, m = 64, n = 50, q = 0, λ = 0.7. Figure 7. The trajectory of the approximate solution and exact solution of Example 2 for m = 32, m = 64, n = 50, q = 0, λ = 0.9. Example 3 [28]. Consider the stochastic Volterra Fredholm integral equation with time delay τ > 0: X(t) = −5τ − t + 12t2 − t3 − t4τ + ∫ t 0 (t − s)X(s − τ)ds + ∫ 1 0 (t + s)X(s − τ)ds + ∫ t 0 sX(s − τ)dB(s), (41) where s, t ∈ [0, T ], τ ∈ (0, T ); with the exact solution X(t) = exp ( (6t + 12t2)/2 + ∫ t 0 sdB(s) ) , {B(t) : 0 ≤ t ≤ T} is a Brownian motion process, and X(t) is an unknown stochastic process defined on the probability space (Ω,F,P). Tables 4–5 present numerical results for various values of m, for λ = 0.5 to compute the τ . 7. Conclusion It is possible to use a computational method based on the properties of BPFs with operational matrices to convert the problem into a system of linear algebraic equations. As a result, this technique transforms nonlinear 139 E. Y. Kutorzi, Y. Zhang, Y. Shi Acta Polytechnica n χ̄E SE 95 % confidence interval for mean of E Lower Upper 50 0.19136844 0.06088437 0.17406529 0.20867159 100 0.18904432 0.06630608 0.17588776 0.20220088 150 0.18060058 0.06734891 0.16973445 0.19146671 200 0.18457774 0.07436710 0.17420811 0.19494736 250 0.18513391 0.07368228 0.17595571 0.19431210 300 0.18425841 0.07426964 0.17582001 0.19269681 Table 4. Mean, standard deviation, and mean confidence interval for error in Example 3 with m = 32, q = 0, λ = 0.5. n χ̄E SE 95 % confidence interval for mean of E Lower Upper 50 0.18850030 0.07900637 0.16604694 0.21095366 100 0.19126831 0.07912381 0.17556843 0.20696819 150 0.19011608 0.08055573 0.17711915 0.20311301 200 0.19260138 0.08272020 0.18106700 0.20413575 250 0.18888323 0.08117339 0.17877191 0.19899455 300 0.18898994 0.08172767 0.17970417 0.19827572 Table 5. Mean, standard deviation, and mean confidence interval for error in Example 3 with m = 64, q = 0, λ = 0.5. stochastic Volterra-Fredholm integral equations into a system of linear algebraic equations whose coefficients represent BPFs that represent solutions to these equations. As well as error analysis, numerical examples provide a solid basis for combined effects and observe a perfect agreement between the exact solutions and the numerical results, achieving full convergence. These include stochastic integrals and ordinary differential equations; arithmetic operations are carried out without requiring derivatives or integration. A Python 3 environment was used to perform the computations associated with the examples. We observe the auspicious results and hope to extend the method to more general backward stochastic Volterra integral equations in sequels. Acknowledgements The authors thank the reviewers for their valuable comments and efforts to improve our article. This research was supported by National Key R&D Program of China (2023YFA1008903) and the Major Fundamental Research Project of Shandong Province of China (No. ZR2023DZ33). References [1] H. K. Dawood. Computational block-pulse functions method for solving Volterra integral equations with delay. Journal of University of Babylon for Pure and Applied Sciences 27(1):32–42, 2019. https://doi.org/10.29196/jubpas.v27i1.2063 [2] C. Kasumo. On the approximate solutions of linear Volterra integral equations of the first kind. Applied Mathematical Sciences 14(3):141–153, 2020. https://doi.org/10.12988/ams.2020.912176 [3] K. Maleknejad, P. Torabi, S. Sauter. Numerical solution of a non-linear Volterra integral equation. Vietnam Journal of Mathematics 44:5–28, 2016. https://doi.org/10.1007/s10013-015-0149-8 [4] E. Babolian, Z. Masouri. Direct method to solve Volterra integral equation of the first kind using operational matrix with block-pulse functions. Journal of Computational and Applied Mathematics 220(1–2):51–57, 2008. https://doi.org/10.1016/j.cam.2007.07.029 [5] T. S. Gutleb, S. Olver. A sparse spectral method for Volterra integral equations using orthogonal polynomials on the triangle. SIAM Journal on Numerical Analysis 58(3):1993–2018, 2020. https://doi.org/10.1137/19M1267441 [6] Y. Hamaguchi. On the maximum principle for optimal control problems of stochastic Volterra integral equations with delay. Applied Mathematics & Optimization 87:42, 2023. https://doi.org/10.1007/s00245-022-09958-w [7] K. Maleknejad, K. Mahdiani. Solving nonlinear mixed Volterra-Fredholm integral equations with two dimensional block-pulse functions using direct method. Communications in Nonlinear Science and Numerical Simulation 16(9):3512–3519, 2011. https://doi.org/10.1016/j.cnsns.2010.12.036 140 https://doi.org/10.29196/jubpas.v27i1.2063 https://doi.org/10.12988/ams.2020.912176 https://doi.org/10.1007/s10013-015-0149-8 https://doi.org/10.1016/j.cam.2007.07.029 https://doi.org/10.1137/19M1267441 https://doi.org/10.1007/s00245-022-09958-w https://doi.org/10.1016/j.cnsns.2010.12.036 vol. 64 no. 2/2024 Stochastic Volterra-Fredholm with Delay [8] M. Rabbani, K. Nouri. Solution of integral equations by using block-pulse functions. Mathematical Sciences Quarterly Journal 4(1):39–48, 2010. [9] S. H. Esmail Babolian, Zahra Masouri. New direct method to solve nonlinear Volterra-Fredholm integral and integro-differential equations using operational matrix with block-pulse functions. Progress In Electromagnetics Research B 8:59–76, 2008. https://doi.org/10.2528/PIERB08050505 [10] Y. Shi, T. Wang. Solvability of general backward stochastic Volterra integral equations. Journal of the Korean Mathematical Society 49(6):1301–1321, 2012. https://doi.org/10.4134/JKMS.2012.49.6.1301 [11] A. A. Khidir. A numerical technique for solving Volterra-Fredholm integral equations using Chebyshev spectral method. Ricerche di Matematica 2022. https://doi.org/10.1007/s11587-022-00692-7 [12] M. Samar, K. E. Yao, X. Zhu. Numerical solution of nonlinear backward stochastic Volterra integral equations. Axioms 12(9):888, 2023. https://doi.org/10.3390/axioms12090888 [13] A. Bellour, M. Bousselsal. A Taylor collocation method for solving delay integral equations. Numerical Algorithms 65(4):843–857, 2014. https://doi.org/10.1007/s11075-013-9717-8 [14] Q. Zhu, T. Huang. Stability analysis for a class of stochastic delay nonlinear systems driven by G-Brownian motion. Systems & Control Letters 140:104699, 2020. https://doi.org/10.1016/j.sysconle.2020.104699 [15] R. Song, B. Wang, Q. Zhu. Stabilization by variable-delay feedback control for highly nonlinear hybrid stochastic differential delay equations. Systems & Control Letters 157:105041, 2021. https://doi.org/10.1016/j.sysconle.2021.105041 [16] M. Nouri, K. Maleknejad. Numerical solution of delay integral equations by using block pulse functions arises in biological sciences. International Journal of Mathematical Modelling & Computations 6(3):221–232, 2016. [17] F. Toutounian, E. Tohidi, A. Kilicman. Fourier operational matrices of differentiation and transmission: Introduction and applications. Abstract and Applied Analysis 2013:198926, 2013. https://doi.org/10.1155/2013/198926 [18] J. Zhang, Y. Li, J. Xie. Numerical simulation of fractional control system using Chebyshev polynomials. Mathematical Problems in Engineering 2018:4270764, 2018. https://doi.org/10.1155/2018/4270764 [19] F. Stenger. Numerical Methods Based on Sinc and Analytic Functions, vol. 20. Springer New York, 1993. ISBN 978-1-4612-7637-1. https://doi.org/10.1007/978-1-4612-2706-9 [20] Z. Jiang, W. Schaufelberger. Block Pulse Functions and Their Applications in Control Systems. Springer-Verlag, Berlin, 1992. ISBN 978-3-540-55369-4. https://doi.org/10.1007/BFb0009162 [21] G. P. Rao. Piecewise Constant Orthogonal Functions and Their Application to Systems and Control. Springer-Verlag, Berlin, 1983. ISBN 978-3-540-12556-3. https://doi.org/10.1007/BFb0041228 [22] F. C. Klebaner. Introduction to Stochastic Calculus with Applications. Imperial College Press, 3rd edn., 2012. ISBN 978-1-84816-831-2. https://doi.org/10.1142/p821 [23] M. Khodabin, K. Maleknejad, M. Rostami, M. Nouri. Numerical approach for solving stochastic Volterra-Fredholm integral equations by stochastic operational matrix. Computers & Mathematics with Applications 64(6):1903–1913, 2012. https://doi.org/10.1016/j.camwa.2012.03.042 [24] M. A. Hussein, H. K. Jassim. Analysis of fractional differential equations with Antagana-Baleanu fractional operator. Progress in Fractional Differentiation and Applications 9(4):681–686, 2023. https://doi.org/10.18576/pfda/090411 [25] H. K. Jassim, M. A. Hussein, M. R. Ali. An efficient homotopy permutation technique for solving fractional differential equations using Atangana-Baleanu-Caputo operator. AIP Conference Proceedings 2845(1):060008, 2023. https://doi.org/10.1063/5.0157148 [26] V. Horvat. On collocation methods for Volterra integral equations with delay arguments. Mathematical Communications 4(1):93–109, 1999. [27] M. Mosleh, M. Otadi. Least squares approximation method for the solution of Hammerstein-Volterra delay integral equations. Applied Mathematics and Computation 258:105–110, 2015. https://doi.org/10.1016/j.amc.2015.01.100 [28] J.-H. He, M. H. Taha, M. A. Ramadan, G. M. Moatimid. Improved block-pulse functions for numerical solution of mixed Volterra-Fredholm integral equations. Axioms 10(3):200, 2021. https://doi.org/10.3390/axioms10030200 141 https://doi.org/10.2528/PIERB08050505 https://doi.org/10.4134/JKMS.2012.49.6.1301 https://doi.org/10.1007/s11587-022-00692-7 https://doi.org/10.3390/axioms12090888 https://doi.org/10.1007/s11075-013-9717-8 https://doi.org/10.1016/j.sysconle.2020.104699 https://doi.org/10.1016/j.sysconle.2021.105041 https://doi.org/10.1155/2013/198926 https://doi.org/10.1155/2018/4270764 https://doi.org/10.1007/978-1-4612-2706-9 https://doi.org/10.1007/BFb0009162 https://doi.org/10.1007/BFb0041228 https://doi.org/10.1142/p821 https://doi.org/10.1016/j.camwa.2012.03.042 https://doi.org/10.18576/pfda/090411 https://doi.org/10.1063/5.0157148 https://doi.org/10.1016/j.amc.2015.01.100 https://doi.org/10.3390/axioms10030200 Acta Polytechnica 64(2):128–141, 2024 1 Introduction 2 Block-pulse functions (BPFs) 2.1 Functions approximation 2.2 Integration operational matrix 2.3 The operational matrix with time delay of BPFs 3 Stochastic integration operational matrix 4 Solving stochastic Volterra-Fredholm integral equations with time delay 5 Error estimation and rate of convergence 6 Numerical examples 7 Conclusion Acknowledgements References