EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6969 ISSN 1307-5543 – ejpam.com Published by New York Business Global A Comparative Study of Bernoulli Collocation and Hermite-Galerkin Methods for Solving Two-Dimensional Nonlinear Volterra Integral Equations of the Second Kind M. H. Ahmed1,∗, A. M. Aljabri1 1 Taibah University, College of Science, Department of Mathematics, P.O. Box 344, Madinah, 3002 Saudi Arabia Abstract. This article is devoted to the presentation of two numerical methods which give the solution of a two-dimensional nonlinear Volterra integral equation of the second kind. The first method, Bernoulli collocation, depend on approximating the unknown function using Bernoulli polynomials, while applying the collocation technique at shifted Chebyshev points over the interval [0,1]. The second method, Hermite-Galerkin method, relies on constructing an operational matrices and applying the Galerkin projection, which we have a system of nonlinear algebraic equations from Volterra integral equation. Discussion on the existence and uniqueness of the solution is provided. Finally, the effect of that two numerical methods is described. To illustrate the previously described methods, several numerical examples are provided. Numerical results show that the Bernoulli collocation method consistently provides more accurate and efficient results than the Hermite– Galerkin method for the same number of collocation points. Comparisons with previously published approaches further demonstrate the superiority of the proposed methods in terms of convergence and stability. 2020 Mathematics Subject Classifications: 11B39, 45B05, 45D05, 65R20 Key Words and Phrases: Non-linear two-dimensional Volterra integral equations, Bernoulli polynomials, shifted Chebyshev points, Hermite polynomials and Galerkin method 1. Introduction Volterra integral equations (VIEs), particularly in two dimensions, serve as an effective tool for modeling various phenomena in physical and engineering sciences. For instance, (VIEs) are used to describe viscoelastic rods and plates, where stress depends on the full strain history [1], and nonlinear viscoelastic solids [2], highlighting their relevance in modeling real-world materials and systems. These equations have attracted considerable ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6969 Email addresses: mramadan@taibahu.edu.sa (M. H. Ahmed), asmazaljabrii@gmail.com (A. M. Aljabri) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 2 of 16 attention in both theoretical analysis and the development of accurate and efficient nu- merical methods for their solution [3–5]. In the present study, we examine the following two-dimensional nonlinear Volterra integral equation of the second kind: u(χ, t) = f(χ, t) + λ ∫ t 0 ∫ χ 0 K(χ, t, s, y, u(s, y)) ds dy, χ, t ∈ [0, 1], (1) where u(χ, t) is the unknown function defined on the domain D = [0, 1]× [0, 1], and f(χ, t), K(χ, t, s, y, u(s, y)) are assumed to be known analytical functions. Many numerical methods have been developed to solve two-dimensional Volterra inte- gral equations due to the difficulty of obtaining analytical solutions. For example, several numerical techniques have been proposed for these equations, including the Chebysheve polynomials as the basis in the collocation method [6], differential transformation method [7], reproducing kernel function method [8], Legendre polynomials method [9], rationalized Haar functions [10] and two-dimensional block-pulse functions [11, 12]. In [13], a numerical method was presented for solving a class of nonlinear two-dimensional integral equations using Bernoulli polynomials. The approach relies on constructing an operational matrix of Bernoulli polynomials, resulting in an efficient and accurate approximation scheme. In [14], a numerical scheme based on shifted Jacobi operational matrices combined with the collo- cation method was developed to solve two-dimensional nonlinear fractional Volterra and Fredholm integral equations. This method effectively handles fractional orders while main- taining high accuracy. The authors of [15] introduced a Laguerre wavelet-based method for solving two-dimensional nonlinear integral equations. Their work included a detailed convergence analysis that confirmed the method’s accuracy and efficiency. Moreover, in [16], a numerical method using radial basis functions (RBFs) was developed to solve non- linear two-dimensional Volterra integral equations of the second kind on non-rectangular domains. The method exhibits strong flexibility in dealing with complex geometries while preserving computational accuracy. However, the field of integral equations is closely related to the solution of a wide class of differential models, including fractional ones, as they share both challenges and advanced numerical techniques. Several numerical [17], analytical [18], and semi-analytical techniques have been developed for solving differential equations, particularly focusing on accurate differentiation [19–21]. These studies provide a foundational background for extending such methods to more complex problems, including Volterra integral equations. In this work, we propose two numerical methods for solving a nonlinear two-dimensional Volterra integral equation. The first method, Bernoulli collocation, is based on approxi- mating the unknown function using Bernoulli polynomials, while applying the collocation technique at shifted Chebyshev points over the interval [0,1]. The integral equation is then transformed into a system of nonlinear algebraic equations. The second method, Hermite-Galerkin, relies on approximating the solution using Hermite functions and ap- plying the Galerkin projection, which also reduces the integral equation to a system of nonlinear algebraic equations. The outline of this paper is as follows: In Section 2, we review some fundamental formulations and properties of Bernoulli and Hermite polynomials. Section 3 is devoted M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 3 of 16 to the existence and uniqueness of the solution of the integral equation. In Section 4, we present computational methods for solving the two-dimensional Volterra integral equation (1) using both Bernoulli collocation method and Hermite-Galerkin method. Comparison of results are found in some numerical examples which provided in Section 5. The results are discussed in Section 6. Finally, the paper is concluded in Section 7. These previously proposed methods provide a foundation and motivation for the present study, which fo- cuses on Bernoulli collocation and Hermite–Galerkin approaches as efficient alternatives for solving two-dimensional nonlinear Volterra integral equations. 2. Some Definitions and Properties 2.1. Bernoulli polynomials Bernoulli polynomials were introduced by Jacob Bernoulli and later generalized by Euler to arbitrary values of the variable [22–24]. The generating function of Bernoulli polynomials BN (χ) is given by z1e χz1 ez1 − 1 = ∞∑ N=0 BN (χ) zN1 N ! , |z1| < 2π, (2) The Bernoulli numbers BN := BN (0) are defined via the generating function z1 ez1 − 1 = ∞∑ N=0 BN zN1 N ! , |z1| < 2π. (3) From this, the Bernoulli polynomials BN (χ) can be expressed as BN (χ) = N∑ k=0 ( N k ) Bkχ N−k. (4) They satisfy the identity BN (χ+ 1)−BN (χ) = NχN−1, N ∈ N0, (5) which implies that BN (0) = BN (1), N ∈ N \ {1}. (6) Substituting χ = 1 into (4), and using (6), yields BN = N∑ k=0 ( N k ) Bk. (7) These polynomials satisfy several important identities, such as B′ N (χ) = NBN−1(χ), N ≥ 1, M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 4 of 16∫ 1 0 BN (χ) dχ = 0, N ≥ 1, BN (χ+ 1)−BN (χ) = NχN−1, N ≥ 1, and BN (1− χ) = (−1)NBN (χ). The first five Bernoulli polynomials are B0(χ) = 1, B1(χ) = χ− 1 2 , B2(χ) = χ2 − χ+ 1 6 , B3(χ) = χ3 − 3 2 χ2 + 1 2 χ, B4(χ) = χ4 − 2χ3 + χ2 − 1 30 . 2.2. Hermite polynomials Hermite polynomials were first introduced in 1810 by Pierre-Simon Laplace. However, it was Charles Hermite who later defined a more generalized class of Hermite polynomials, which remain less known compared to the standard form, despite their significance. These polynomials have emerged as a fundamental tool in both pure and applied mathematics. Their relevance has recently increased due to their applications in quantum mechanics, engineering, physics, and other scientific fields. Hermite polynomials form a set of mutually orthogonal functions with respect to the weight function e−χ2 over the interval (−∞,∞) [23, 25, 26]. They can be generated using the Rodrigues formula: Hn(χ) = (−1)neχ2 dn dχn ( e−χ2 ) , n = 0, 1, 2, . . . The first few Hermite polynomials are given by H0(χ) = 1, H1(χ) = 2χ, H2(χ) = 4χ2 − 2, H3(χ) = 8χ3 − 12χ, H4(χ) = 16χ4 − 48χ2 + 12, H5(χ) = 32χ5 − 160χ3 + 120χ. These polynomials satisfy several important identities. The derivative of Hn(χ) is given by H ′ n(χ) = 2nHn−1(χ), n ≥ 1. M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 5 of 16 They also satisfy a recurrence relation Hn+1(χ)− 2χHn(χ) = −2nHn−1(χ), n = 1, 2, . . . Moreover, they exhibit symmetry properties based on the parity of n Hn(−χ) = (−1)nHn(χ). With respect to orthogonality, Hermite polynomials satisfy the following inner product identity over the entire real line∫ ∞ −∞ e−χ2 Hm(χ)Hn(χ) dχ = { 0, m ̸= n, 2n √ π n!, m = n. 3. Existence and uniqueness of the solution of the integral equation First, we need to prove the existence and uniqueness of the solution of (1). To this end, we employ the Banach fixed point theorem. Theorem 1. Let (C(D), ∥·∥) as the Banach space for continuous real-valued functions on D, with norm ∥u∥ = max (χ,t)∈D |u(χ, t)|. Assume that the function f(χ, t) is continuous on D, and the kernel k(χ, t, s, y, u) is continuous on D ×D × R, and satisfies the Lipschitz condition in the fifth argument |K(χ, t, s, y, u)−K(χ, t, s, y, v)| ≤ L|u− v| for all (χ, t, s, y) ∈ D ×D and all u, v ∈ R, with a constant L > 0. If α = |λ|L < 1, and 0 < α < 1 then the integral equation (1) has a unique solution u ∈ C(D). Proof. Define the operator T on C(D) by (Tu)(χ, t) = f(χ, t) + λ ∫ t 0 ∫ χ 0 K(χ, t, s, y, u(s, y)) ds dy. We aim to show that T is a contraction on C(D). For any u, v ∈ C(D), we have ∥Tu− Tv∥ = max (χ,t)∈D |(Tu)(χ, t)− (Tv)(χ, t)| = max (χ,t)∈D ∣∣∣∣λ ∫ t 0 ∫ χ 0 [K(χ, t, s, y, u(s, y))−K(χ, t, s, y, v(s, y))] ds dy ∣∣∣∣ ≤ |λ| max (χ,t)∈D ∫ t 0 ∫ χ 0 |K(χ, t, s, y, u(s, y))−K(χ, t, s, y, v(s, y))| ds dy ≤ |λ|L max (χ,t)∈D ∫ t 0 ∫ χ 0 |u(s, y)− v(s, y)| ds dy ≤ |λ|L∥u− v∥ max (χ,t)∈D ∫ t 0 ∫ χ 0 ds dy = |λ|L∥u− v∥ max (χ,t)∈D χ t = α∥u− v∥. M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 6 of 16 Since α = |λ|L < 1, the operator T is a contraction. Therefore, by the Banach Fixed Point Theorem, T has a unique fixed point u ∈ C(D), which is the unique solution of the integral equation. To illustrate the convergence, S.Bazm studied the Bernoulli polynomials for integral equations [27], while Mao and Shen analyzed the Hermite–Galerkin method for fractional PDEs [25]. 4. Description of the methods In this section, we solve Eq.(1) using Bernoulli collocation method and Hermite Galerkin method. 4.1. Bernoulli collocation method In this method, we approximate the unknown function ũ(χ, t) in equation (1) using a double series expansion of the form ũ(χ, t) = ∞∑ i=0 ∞∑ j=0 aijBi(χ)Bj(t), (8) where Bi(χ) and Bj(t) are Bernoulli polynomials and the unknown coefficients aij are be determined. This representation is then used to construct the approximate solution as follows. By truncating the infinite series in equation (8), we obtain the following finite approx- imation ũ(χ, t) ≈ N−1∑ i=0 N−1∑ j=0 aijBi(χ)Bj(t), (9) where N − 1 is the chosen degree of approximation in both variables. Substituting from (9) into (1) we get N−1∑ i=0 N−1∑ j=0 aijBi(χ)Bj(t) = f(χ, t) + λ ∫ t 0 ∫ χ 0 K χ, t, s, y, N−1∑ i=0 N−1∑ j=0 aijBi(s)Bj(y)  ds dy. (10) The unknown function ũ(χ, t) is approximated using Bernoulli polynomials as basis func- tions, while the collocation points χ l , t f are chosen as shifted Chebyshev points on the interval [0, 1], given by χ l = 1 2 ( 1 + cos ( (2l − 1)π 2N )) , t f = 1 2 ( 1 + cos ( (2f − 1)π 2N )) , l, f = 1, 2, . . . , N. (11) M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 7 of 16 Equation (10) can be represented as N−1∑ i=0 N−1∑ j=0 aijBi(χl )Bj(tf ) = f(χ l , t f )+λ ∫ t f 0 ∫ χ l 0 K χ l , t f , s, y, N−1∑ i=0 N−1∑ j=0 aijBi(s)Bj(y)  ds dy. (12) By substituting the collocation points defined in (11) into equation (10), we obtain a sys- tem of nonlinear algebraic equations containing N2 unknown coefficients aij , with indices i, j = 0, 1, ..., N − 1. In this study, N = 6 was chosen as it provides sufficient accuracy while keeping the computational cost reasonable. All double integrals, including the non- linear terms containing the kernel K, are numerically evaluated using MATLAB’s integral2 function (adaptive quadrature method). The resulting nonlinear system is solved itera- tively using the Picard iteration method with a convergence tolerance of 1× 10−6, with a maximum limit of 10 iterations to ensure numerical stability, leading to the approximate solution ũ(χ, t). 4.2. Hermite Galerkin method Assume that ũ(χ, t) is an approximate solution of the two-dimensional Volterra integral equation (1). The Galerkin method with Hermite polynomials is applied, yielding the following approximation ũ(χ, t) ≈ N−1∑ i=0 N−1∑ j=0 cijHi(χ)Hj(t), (13) where the Hermite polynomials are Hi(χ) and Hj(t) and the unknown Hermite coefficients cij are to be determined. From the right-hand side of equation (13) and substituting into equation (1), we obtain N−1∑ i=0 N−1∑ j=0 cijHi(χ)Hj(t) = f(χ, t) + η(χ, t), (14) where η(χ, t) = λ ∫ t 0 ∫ χ 0 K χ, t, s, y, N−1∑ i=0 N−1∑ j=0 cijHi(s)Hj(y)  ds dy. Multiplying equation (14) by Hb(χ)Hr(t), also integrating both sides of equation (14) with respect to χ and t over [0, 1], for b, r = 0, 1, . . . , N − 1, we have the following. N−1∑ i=0 N−1∑ j=0 cij ∫ 1 0 ∫ 1 0 Hi(χ)Hj(t)Hb(χ)Hr(t) dχ dt = ∫ 1 0 ∫ 1 0 f(χ, t)Hb(χ)Hr(t) dχ dt + ∫ 1 0 ∫ 1 0 η(χ, t)Hb(χ)Hr(t) dχ dt. (15) M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 8 of 16 By substituting all combinations of b, r = 0, 1, . . . , N − 1 into equation (15), we obtain a system of N2 nonlinear algebraic equations involving the unknown Hermite coefficients cij . In this work, N = 6 is chosen to balance computational efficiency and accuracy. All double integrals appearing in the formulation are numerically evaluated using MATLAB’s integral2 function. The resulting nonlinear system is solved iteratively using the Picard iteration method with a convergence tolerance of 1 × 10−6, leading to the approximate solution ũ(χ, t). 5. Numerical examples To illustrate the previously described methods, several numerical examples of the two- dimensional nonlinear Volterra integral equation (2D-NVIE) are provided. The results were obtained using MATLAB R2025a. Example 1. Consider the following 2D-NVIE [8, 28]: u(χ, t) = f(χ, t) + ∫ t 0 ∫ χ 0 (χs2 + cos y)u2(s, y) ds dy χ, t ∈ [0, 1], (16) where f(χ, t) = χ sin t(1− 1 9 χ2 sin2 t) + 1 10 χ6( 1 2 sin 2t− t), and exact solution is u(χ, t) = χ sin t. The comparison of absolute errors of equation (16) for various values of χ and t, with using Bernoulli collocation (BC), Hermite–Galerkin (HG) methods with N = 6, the reproducing kernel space [8] with N = 30 and the Taylor collocation method [28] with N = 64, is presented in Table (1). Figures (1) and (2) show the absolute error distributions obtained by these methods. Table 2 shows the absolute errors of equation (16) for (χ, t) = (0.5, 0.5) using BC and HG methods at different N . The errors decrease as N increases, with BC converging much faster than HG due to the latter’s higher computational cost from nested integrals. Figure 3 illustrates the same trend visually. Therefore, N = 6 was chosen as a suitable compromise between accuracy and efficiency. Figure 1: Absolute error of Example 1 by the BC method, N = 6. Figure 2: Absolute error for Example 1 by HG method, N = 6. M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 9 of 16 (χ, t) = ( 1 2p , 1 2p ) Bernoulli collo- cation method N = 6 Hermite- Galerkin method N = 6 Method of [8] with N = 30 Method of [28] with N = M = 64 p = 1 1.61× 10−7 3.70× 10−5 6.0× 10−5 3.45× 10−7 p = 2 7.54× 10−8 5.56× 10−5 1.29× 10−4 1.28× 10−10 p = 3 1.33× 10−8 8.06× 10−5 7.0× 10−5 2.22× 10−11 p = 4 1.77× 10−8 1.0× 10−4 5.33× 10−5 1.67× 10−13 p = 5 4.69× 10−9 3.54× 10−5 5.959× 10−5 3.23× 10−13 p = 6 2.93× 10−10 4.33× 10−5 7.4588× 10−4 9.42× 10−15 Table 1: Numerical results for Example 1. N BC method HG method 2 1.5726× 10−2 1.0118× 10−2 4 7.7782× 10−5 5.3089× 10−5 6 1.6191× 10−7 3.7095× 10−5 7 6.0440× 10−11 6.9195× 10−5 Table 2: Effect of N on the absolute error for (χ, t) = (0.5, 0.5). Figure 3: Effect of N on the absolute error for (χ, t) = (0.5, 0.5). Example 2. Let us present the following 2D-NVIE [10]: u(χ, t) = f(χ, t) + ∫ t 0 ∫ χ 0 u2(s, y) ds dy χ, t ∈ [0, 1], (17) M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 10 of 16 where f(χ, t) = χ2 + t2 − 1 45 χt(9χ4 + 10χ2t2 + 9t4), and the exact solution is u(χ, t) = χ2 + t2. Table (3) presents the absolute errors of equation (17) calculated using the Bernoulli collocation (BC) and Hermite–Galerkin (HG) methods for various values of χ and t. Fig- ures (4)–(5) illustrate the absolute error distributions obtained by these methods. Fur- thermore, a comparison with the rationalized Haar functions method [10] with N = 32 is also included. (χ, t) = ( 1 2p , 1 2p ) Bernoulli collo- cation method N = 6 Hermite- Galerkin method N = 6 Method of [10] with N = 32 p = 2 5.37× 10−12 6.92× 10−5 5.90× 10−5 p = 3 1.34× 10−12 5.06× 10−5 9.06× 10−7 p = 4 1.15× 10−12 9.52× 10−6 1.29× 10−8 p = 5 5.70× 10−13 2.49× 10−4 1.43× 10−10 p = 6 1.62× 10−13 4.67× 10−4 2.00× 10−12 Table 3: Numerical results for Example 2. Figure 4: Absolute error of Example 2 by the BC method with N = 6. Figure 5: Absolute error for Example 2 by the HG method with N = 6. Example 3. Let us present the following 2D-NVIE [29]: u(χ, t) = f(χ, t) + ∫ t 0 ∫ χ 0 (χ+ t)eu(s,y) ds dy χ, t ∈ [0, 1], (18) M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 11 of 16 where f(χ, t) = (χ+ t)(eχ + et − e(χ+t)). The exact solution is given by u(χ, t) = χ+ t. Table (4) shows the comparison of the absolute errors of equation (18) calculated using the Bernoulli collocation (BC) and Hermite–Galerkin (HG) methods for various values of χ and t. Figures (6)–(7) illustrate the absolute error distributions obtained by these methods. A comparison with the extrapolation method [29] for m = n = 26 is also provided. (χ, t) Bernoulli collo- cation method N = 6 Hermite- Galerkin method N = 6 Method of [29] with m = n = 26 (0.1, 0.1) 1.87× 10−11 5.24× 10−6 1.5× 10−7 (0.2, 0.2) 4.52× 10−12 4.36× 10−5 1.5× 10−6 (0.3, 0.3) 2.58× 10−11 8.14× 10−6 6.4× 10−6 (0.4, 0.4) 3.73× 10−12 1.91× 10−5 1.9× 10−5 (0.5, 0.5) 3.98× 10−11 1.67× 10−5 4.7× 10−5 (0.6, 0.6) 8.77× 10−11 3.09× 10−5 1.0× 10−4 (0.7, 0.7) 2.83× 10−11 5.62× 10−5 2.0× 10−4 (0.8, 0.8) 6.36× 10−12 5.79× 10−5 3.8× 10−4 (0.9, 0.9) 2.01× 10−9 1.5× 10−4 6.8× 10−4 Table 4: Numerical results for Example 3. Figure 6: Absolute error for Example 3 by the BC method with N = 6. Figure 7: Absolute error for Example 3 by HG method with N = 6. M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 12 of 16 6. Results and Discussion As shown in the numerical results for the three test examples (Tables (1)–(4) and Figures (1)–(7)), the Bernoulli collocation (BC) method provides smaller absolute errors compared to the Hermite–Galerkin (HG) method. Although N = 6 was generally selected as a suitable compromise between accuracy and computational cost, a higher value of N (e.g., N = 7 in Example (1)) was occasionally used to investigate the error behavior. The computational cost of the HG method increases significantly with N due to the evalu- ation of multiple nested integrals. The numerical results indicate that the BC method consistently outperforms the HG method. This superiority is mainly due to the better conditioning of the BC formulation and the suitability of the Bernoulli basis functions for smooth nonlinear kernels, which leads to faster convergence and reduced numerical oscillations. 7. Conclusions In this article, we employed Bernoulli collocation and Hermite–Galerkin methods to approximate the solution of two-dimensional nonlinear Volterra integral equations. Nu- merical results show that the Bernoulli collocation method consistently provides more accurate and efficient results than the Hermite–Galerkin method for the same number of points N . In our computations, we used N = 6, as increasing N for the Hermite–Galerkin method is difficult due to the complexity of its integrals, while the Bernoulli collocation method handles larger N more efficiently. Comparisons with previously published meth- ods, such as the reproducing kernel space [8], Taylor collocation [28], rationalized Haar functions [10], and extrapolation methods [29], further demonstrate the superiority of our approaches. The better performance of the Bernoulli collocation method is mainly due to the suitability of the Bernoulli basis functions for smooth nonlinear kernels and improved conditioning, leading to faster convergence and reduced numerical oscillations. Further- more, the present work can be extended to fractional-order versions of two-dimensional Volterra integral equations to capture memory and nonlocal effects, as illustrated in related studies [19, 30]. Acknowledgements The authors are deeply grateful to the editor and anonymous reviewers for their valu- able comments, insightful suggestions, and constructive recommendations that greatly improved the quality and organization of this manuscript. Competing Interests: The authors declare that they have no conflicts of interest re- garding the publication of this work. M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 13 of 16 References [1] Richard D Noren. A linear volterra integro-differential equation for viscoelastic rods and plates. Quarterly of applied mathematics, 45(3):503–514, 1987. [2] A25126821197 Wineman. Nonlinear viscoelastic solids—a review. Mathematics and mechanics of solids, 14(3):300–366, 2009. [3] Faheem Khan, Muhammad Omar, and Zafar Ullah. Discretization method for the numerical solution of 2d volterra integral equation based on two-dimensional bernstein polynomial. AIP Advances, 8(12), 2018. [4] Pouria Assari and Mehdi Dehghan. A meshless local discrete galerkin (mldg) scheme for numerically solving two-dimensional nonlinear volterra integral equations. Applied Mathematics and Computation, 350:249–265, 2019. [5] Manochehr Kazemi. An iterative method for solving two dimensional nonlinear volterra integral equations. Analytical and Numerical Solutions for Nonlinear Equa- tions, 6(1):41–57, 2021. [6] Zakieh Avazzadeh and Mohammad Heydari. Chebyshev polynomials for solving two dimensional linear and nonlinear integral equations of the second kind. Computational & Applied Mathematics, 31:127–142, 2012. [7] A Tari, MY Rahimi, S Shahmorad, and F Talati. Solving a class of two-dimensional linear and nonlinear volterra integral equations by the differential transform method. Journal of Computational and Applied Mathematics, 228(1):70–76, 2009. [8] A Fazli, T Allahviranloo, and Sh Javadi. Numerical solution of nonlinear two- dimensional volterra integral equation of the second kind in the reproducing kernel space. Mathematical Sciences, 11(2):139–144, 2017. [9] Somayeh Nemati, Pedro Miguel Lima, and Yadollah Ordokhani. Numerical solution of a class of two-dimensional nonlinear volterra integral equations using legendre polynomials. Journal of Computational and Applied Mathematics, 242:53–69, 2013. [10] Majid Erfanian and Hamed Zeidabadi. Solving two-dimensional nonlinear volterra in- tegral equations using rationalized haar functions. International Journal of Nonlinear Analysis and Applications, 14(8):95–105, 2023. [11] Esmail Babolian, Khosrow Maleknejad, Mohammad Mordad, and Bijan Rahimi. A numerical method for solving fredholm–volterra integral equations in two-dimensional spaces using block pulse functions and an operational matrix. Journal of Computa- tional and Applied Mathematics, 235(14):3965–3971, 2011. [12] Mohsen Jalalian, Kawa Wali Ali, Sarkawt Raouf Qadir, and Mohamad Reza Jalalian. A numerical method based on the radial basis functions for solving nonlinear two- dimensional volterra integral equations of the second kind on non-rectangular do- mains. Journal of Mathematical Modeling, 12(4):687–705, 2024. [13] Sohrab Bazm. Numerical solution of a class of nonlinear two-dimensional integral equations using bernoulli polynomials. Sahand Communications in Mathematical Analysis, 3(1):37–51, 2016. [14] Jalil Rashidinia, Tahereh Eftekhari, and Khosrow Maleknejad. Numerical solutions of two-dimensional nonlinear fractional volterra and fredholm integral equations using M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 14 of 16 shifted jacobi operational matrices via collocation method. Journal of King Saud University-Science, 33(1):101244, 2021. [15] K Maleknejad and M Soleiman Dehkordi. Numerical solutions of two-dimensional nonlinear integral equations via laguerre wavelet method with convergence analysis. Applied Mathematics-A Journal of Chinese Universities, 36(1):83–98, 2021. [16] K Maleknejad, S Sohrabi, and B Baranji. Application of 2d-bpfs to nonlinear in- tegral equations. Communications in Nonlinear Science and Numerical Simulation, 15(3):527–535, 2010. [17] Emad A Az-Zo’bi and Kamel Al-Khaled. A new convergence proof of the adomian decomposition method for a mixed hyperbolic elliptic system of conservation laws. Applied Mathematics and Computation, 217(8):4248–4256, 2010. [18] EA Az-Zo’bi. A reliable analytic study for higher-dimensional telegraph equation. J. Math. Comput. Sci, 18(4):423–429, 2018. [19] Emad Az-Zo’bi, Lanre Akinyemi, and Ahmed O Alleddawi. Construction of optical solitons for conformable generalized model in nonlinear media. Modern Physics Letters B, 35(24):2150409, 2021. [20] Jamil Abbas Haider, Shahbaz Ahmad, Hassan Ali Ghazwani, Mohamed Hussien, Musawa Yahya Almusawa, and Emad A Az-Zo’bi. Results validation by using finite volume method for the blood flow with magnetohydrodynamics and hybrid nanofluids. Modern Physics Letters B, 38(24):2450208, 2024. [21] Shumaila Kanwal, Syed Asif Ali Shah, Abdul Bariq, Bagh Ali, Adham E Ragab, and Emad A Az-Zo’bi. Insight into the dynamics of heat and mass transfer in nanofluid flow with linear/nonlinear mixed convection, thermal radiation, and activation energy effects over the rotating disk. Scientific Reports, 13(1):23031, 2023. [22] Da-Qian Lu and Qiu-Ming Luo. Some generalizations of 2d bernoulli polynomials. Journal of Inequalities and Applications, 2013(1):110, 2013. [23] DOAA SHOKRY MOHAMED and DINA MOHAMED ABDESSAMI. A comparison between bernoulli-collocation method and hermite-galerkin method for solving two- dimensional mixed volterra-fredholm singular integral equations. Transactions of A. Razmadze Mathematical Institute, 175(2), 2021. [24] Don Zagier. Curious and exotic identities for bernoulli numbers. Max Planck Insitute for Mathematics, Bonn, Germany, nd, 2014. [25] Zhiping Mao and Jie Shen. Hermite spectral methods for fractional pdes in unbounded domains. SIAM Journal on Scientific Computing, 39(5):A1928–A1950, 2017. [26] Changtao Sheng, Suna Ma, Huiyuan Li, Li-Lian Wang, and Lueling Jia. Generalised hermite functions and their applications in spectral approximations. In other words, 1(2):1–3, 2020. [27] Sohrab Bazm. Bernoulli polynomials for the numerical solution of some classes of linear and nonlinear integral equations. Journal of Computational and Applied Math- ematics, 275:44–60, 2015. [28] Hafida Laib, Aissa Boulmerka, Azzeddine Bellour, and Fouzia Birem. Numerical so- lution of two-dimensional linear and nonlinear volterra integral equations using taylor collocation method. Journal of Computational and Applied Mathematics, 417:114537, M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 15 of 16 2023. [29] Yubin Pan and Jin Huang. Extrapolation method for solving two-dimensional volter- ral integral equations of the second kind. Applied Mathematics and Computation, 367:124784, 2020. [30] Mohammad A Al Zubi, Kallekh Afef, and Emad A Az-Zo’bi. Assorted spatial optical dynamics of a generalized fractional quadruple nematic liquid crystal system in non- local media. Symmetry, 16(6):778, 2024. M. H. Ahmed, A. M. Aljabri / Eur. J. Pure Appl. Math, 18 (4) (2025), 6969 16 of 16 Appendix A. Implementation of the Bernoulli Collocation Method Algorithm 1 Bernoulli Collocation Maximum degree N ; Tolerance tol; Maximum iterations Imax; Nonlinear function g(u). 1. Initialize C(0) (coefficient vector) to zeros (or a suitable initial guess). 2. Define Collocation Points xi, tj using Chebyshev roots: xpoints = 1 2 [ cos ( (2i−1)π 2N ) + 1 ] . 3. Generate Bernoulli basis functions Bk(x) for k = 0, . . . , N − 1. for k = 1 to Imax do 4. Compute the approximate solution u(k−1)(s, y) from the previous iteration using C(k−1): u(k−1)(s, y) = ∑N−1 m=0 ∑N−1 n=0 C (k−1) m,n Bm(s)Bn(y). 5. Initialize the system matrices for the current iteration: A(k) ← zeros(N2, N2), F (k) ← zeros(N2, 1). 6. Assemble the system at each collocation point (xi, tj) for i, j = 1, . . . , N (Total N2 points): for i = 1 to N do for j = 1 to N do a. Compute the integral of the known non-linear term (from the previous approximation u(k−1)): J (k−1) i,j = λ ∫ tj 0 ∫ xi 0 K(xi, tj , s, y, u (k−1)(s, y)) ds dy b. Define the right-hand side vector F : F (k) (i−1)N+j = f(xi, tj) + J (k−1) i,j c. Define the system matrix A (which remains constant for the Bernoulli basis): A (k) (i−1)N+j,mN+n = Bm(xi)Bn(tj) 7. Solve the linear system for the new coefficients: C(k) ← (A(k))−1F (k). 8. Check for convergence: if ∥C(k) − C(k−1)∥∞ < tol then return C(k) and break. 9. Update the coefficients: C(k−1) ← C(k).