In this paper, an efficient matrix method based on 2D orthonormal Bernoulli polynomials are developed to obtain numerical solution of weakly singular fractional partial integro-differential equations (FPIDEs). Operational matrix of integration, almost operational matrix of integration, and product operational matrix are employed to transform the solution of Volterra singular FPIDEs to the system of linear or nonlinear algebraic equations. The obtained algebraic system can be easily solved by using an appropriate numerical method. Also, some useful theorems concerning the convergence and error analysis associated to the mentioned approach are proven in the current paper. Finally, some numerical examples are included to illustrate the accuracy, efficiency, and applicability of the proposed approach.
