Operational matrix method for solving two-dimensional stochastic model with weakly singular kernel

Purpose The main objective of this study is to find out the numerical solution of the 2-D stochastic Ito-Volterra integral equation (SIVIE), characterized by kernels with singularities. Since solving these equations analytically is very complicated. Design/methodology/approach The study presents a novel operational matrix method for solving 2-D stochastic integral equations with weakly singular kernels. Operational matrices for product, weakly singular, and stochastic integrals are constructed using orthonormal Bernoulli polynomials. The method employs collocation at Newton–Cotes nodes, transforming the problem into a solvable system of algebraic equations. After determining these coefficients, we get the approximate solution. Error bounds and convergence analyses are also performed to ensure accuracy and stability. Findings The proposed method demonstrates high computational efficiency and precision in solving stochastic integral equations. Numerical experiments on two benchmark problems confirm its capability to handle singularities and randomness simultaneously. Comparative analysis with the bicubic B-spline method shows that the proposed approach achieves superior accuracy and reliability. Practical implications These equations are widely used in viscoelastic material modeling, financial mathematics (e.g. pricing path-dependent options), image processing (e.g. denoising and texture synthesis), and biological systems (e.g. drug diffusion or population dynamics). In engineering, they help model stress propagation and structural health monitoring. Social implications Two-dimensional stochastic Itô-Volterra integral equations with weakly singular kernels model complex systems influenced by memory effects and randomness, such as in climate dynamics, financial markets, and biomedical processes. Their solutions help predict and analyze behaviors under uncertainty, contributing to better risk assessment and decision-making. Originality/value The originality of this study lies in the development of a 2-D novel operational matrix method based on orthonormal Bernoulli polynomials to solve 2-D stochastic integral equations with weakly singular kernels. The approach uniquely constructs operational matrices for product, weakly singular, and stochastic integrals, enabling efficient numerical implementation through collocation at Newton–Cotes nodes. This integrated framework provides a new and unified strategy for addressing the combined challenges of singularity and randomness in stochastic systems.