This work is focused on developing a hybrid numerical method that combines a higher-order finite difference method and multi-dimensional Hermite wavelets to address two-dimensional multi-term time-fractional integro-partial differential equations with multi-term weakly singular kernels having bounded and unbounded time derivatives at the initial time t=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0$$\end{document}. Specifically, the multi-term fractional operators are discretized using a higher-order approximation designed by employing different interpolation schemes based on linear, quadratic, and cubic interpolation leading to O(N-(4-alpha 1))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}(N<^>{-(4-\alpha _1)})$$\end{document} accuracy on a suitably chosen nonuniform mesh and O(N-alpha 1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}(N<^>{-\alpha _1})$$\end{document} accuracy on a uniformly distributed mesh. The weakly singular integral operators are approximated by a modified numerical quadrature, which is a combination of the composite trapezoidal approximation and the midpoint rule. The effects of the exponents of the weakly singular kernels over fractional orders are analyzed in terms of accuracy over uniform and nonuniform meshes for the solution having both bounded and unbounded time derivatives. The stability of the proposed semi-discrete scheme is derived based on L infinity\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L<^>\infty $$\end{document}-norm for uniformly distributed temporal mesh. Further, we employ the uniformly distributed collocation points in spatial directions to estimate the tensor-based wavelet coefficients. Moreover, the convergence analysis of the fully discrete scheme is carried out based on L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L<^>2$$\end{document}-norm leading to O(N-alpha 1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}(N<^>{-\alpha _1})$$\end{document} accuracy on a uniform mesh. It also highlights the higher-order accuracy over nonuniform mesh. Additionally, we discuss the convergence analysis of the proposed scheme in the context of the multi-term time-fractional diffusion equations involving time singularity demonstrating a O(N-(4-alpha 1))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}(N<^>{-(4-\alpha _1)})$$\end{document} accuracy on a nonuniform mesh with suitably chosen grading parameter. Note that the scheme reduces to O(N-alpha 1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}(N<^>{-\alpha _1})$$\end{document} accuracy on a uniform mesh. Several tests are performed on numerous examples in L infinity\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L<^>\infty $$\end{document}- and L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L<^>2$$\end{document}-norm to show the efficiency of the proposed method. Further, the solutions' nature and accuracy in terms of absolute point-wise error are illustrated through several isosurface plots for different regularities of the exact solution. These experiments confirm the theoretical accuracy and guarantee the convergence of approximations to the functions having time singularity, and the higher-order accuracy for a suitably chosen nonuniform mesh.
