A Bayesian Inverse Framework with Convolutional Regularization for Weakly Singular Volterra Integral Equations

1. Introduction

Volterra’s integral equations, named after Vito Volterra [1,2], constitute a fundamental class where the unknown function appears under an integral with variable limits. These equations model systems with memory effects across physics, biology, economics, and control theory [3,4].

The mathematical foundations were established by pioneering work of Volterra [3,4], Lalescu [5], Fredholm [6], and Hilbert [7], with later contributions from Wiener and Hopf [8,9], Taylor [10], and Kato [11].

Volterra equations are classified into first-kind and second-kind equations. While continuous kernels allow direct treatment, singular or weakly singular kernels pose considerable challenges due to their divergent behavior [12,13,14]. The theory was extensively developed by Tricomi [15], Muskhelishvili [16], and Erdélyi [17].

Existing solution methods include analytical approaches such as Laplace transforms [18,19], series solutions [20,21], and Picard iteration [22], alongside numerical methods including quadrature methods [23,24], collocation methods [13,25,26], and Galerkin methods [27,28]. Specialized techniques for singular kernels involve graded meshes [29,30] and product integration [31,32].

Traditional methods face significant limitations, especially regarding stability and uncertainty quantification near singularities. Conventional methods often fail near singularities, producing unreliable results due to inherent ill-conditioning and error propagation. Current techniques lack robust uncertainty management and adaptive strategies for handling singular behavior.

This study proposes an innovative hybrid approach combining Bayesian inference with Gaussian processes and adaptive convolutional regularization to solve Volterra integral equations with weak singularities. Our methodology addresses traditional limitations by exploiting Gaussian processes’ uncertainty quantification capabilities while utilizing regularization techniques to stabilize numerical behavior near singularities, enhanced by importance sampling [33] and logarithmic stabilization methods [34].

This paper is organized as follows: in Section 2, we provide a detailed description of the proposed method. In Section 3, we deal with the numerical examples and we conclude in Section 4.

2. Proposed Method

We consider the singular Volterra integral equation of the second kind defined by:

$$f\left(x\right)-\lambda {\int }_0^x{\displaystyle \frac{k(x,t)F\left(f\right(t\left)\right)}{{(x-t)}^{1-\alpha }}}dt=g\left(x\right),0\le x\le 1$$

(1)

where $f\left(x\right)$ represents the unknown function, while $k(x,t)$, $g\left(x\right)$, and F are given functions. The parameter $\lambda \in \mathbb{R}$ is a scalar multiplier, and $\alpha >0$ determines the type and strength of the singularity.

The Volterra operator is formally defined as:

$$\mathcal{L}\left[f\right]\left(x\right)=f\left(x\right)-\lambda {\int }_0^{x}K(x,t)dt$$

(2)

where ${\displaystyle K(x,t)=\frac{k(x,t)F\left(f\right(t\left)\right)}{{(x-t)}^{1-\alpha }}}$ represents the singular kernel. For $0<\alpha <1$, the exponent $(1-\alpha )\in (0,1)$ creates a weak singularity, and the integral

$${\int }_0^x{\displaystyle \frac{1}{{(x-t)}^{1-\alpha }}}dt={\displaystyle \frac{1}{\alpha }}{x}^{\alpha }<\infty$$

remains convergent, ensuring the well-posedness of the problem. Existence and uniqueness for similar singular kernels are discussed in [35,36].

2.1. Bayesian Formulation

2.1.1. Likelihood Function

Assuming that the observed data $g_{\mathrm{obs}}\left(x\right)=g_{\mathrm{true}}\left(x\right)$ represents noise-free observations, we define the likelihood function [37,38] as:

$$p\left(g_{\mathrm{obs}}\right|f)=\prod _{i=1}^n{\displaystyle \frac{1}{\sqrt{2\pi {\sigma }^2}}}exp\left(-{\displaystyle \frac{{[g_{\mathrm{obs}}\left(x_i\right)-\mathcal{L}\left[f\right]\left(x_i\right)]}^2}{2{\sigma }^2}}\right)$$

(3)

Let $A\left(x_i\right)=g_{\mathrm{obs}}\left(x_i\right)-\mathcal{L}\left[f\right]\left(x_i\right)$ denote the residual at point $x_i$. The likelihood becomes:

$$p\left(g_{\mathrm{obs}}\right|f)=\prod _{i=1}^n{\displaystyle \frac{1}{\sqrt{2\pi {\sigma }^2}}}exp\left(-{\displaystyle \frac{{\left[A\left(x_i\right)\right]}^2}{2{\sigma }^2}}\right)$$

(4)

The log-likelihood is derived as follows:

$$\begin{array}{cc}\hfill lnp\left(g_{\mathrm{obs}}\right|f)& =ln\left(\prod _{i=1}^n{\displaystyle \frac{1}{\sqrt{2\pi {\sigma }^2}}}exp\left(-{\displaystyle \frac{{\left[A\left(x_i\right)\right]}^2}{2{\sigma }^2}}\right)\right)\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & =\sum _{i=1}^{n}ln\left({\displaystyle \frac{1}{\sqrt{2\pi {\sigma }^2}}}exp\left(-{\displaystyle \frac{{\left[A\left(x_i\right)\right]}^2}{2{\sigma }^2}}\right)\right)\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & =\sum _{i=1}^{n}ln{\displaystyle \frac{1}{\sqrt{2\pi {\sigma }^2}}}-\sum _{i=1}^n{\displaystyle \frac{{\left[A\left(x_i\right)\right]}^2}{2{\sigma }^2}}\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & =-{\displaystyle \frac{n}{2}}ln\left(2\pi {\sigma }^2\right)-{\displaystyle \frac{1}{2{\sigma }^2}}\sum _{i=1}^n{\left[A\left(x_i\right)\right]}^2\hfill \end{array}$$

(8)

2.1.2. Prior Distribution

We impose a Gaussian process (GP) prior [38,39,40] on the unknown function f:

$$f\sim \mathcal{GP}(\mu \left(x\right),k(x,x^{\prime }))$$

(9)

where $\mu \left(x\right)$ is the mean function and $k(x,x^{\prime })$ is the covariance kernel. We employ a Radial Basis Function (RBF) kernel:

$$k_{\mathrm{RBF}}(x,x^{\prime })={\sigma }_f^{2}exp\left(-{\displaystyle \frac{|x-x^{\prime }{|}^2}{2{\ell }^2}}\right)$$

(10)

where ${\sigma }_f^2$ controls the variance and ℓ is the characteristic length scale.

2.1.3. Posterior Distribution

By Bayes’ theorem, the posterior distribution [37,38] is given by:

$$p\left(f\right|g_{\mathrm{obs}})={\displaystyle \frac{p\left(g_{\mathrm{obs}}\right|f)p\left(f\right)}{p\left(g_{\mathrm{obs}}\right)}}$$

(11)

Taking logarithms:

$$lnp\left(f\right|g_{\mathrm{obs}})=lnp\left(g_{\mathrm{obs}}\right|f)+lnp\left(f\right)-lnp\left(g_{\mathrm{obs}}\right)$$

(12)

2.2. Singularity Treatment

The primary computational challenge lies in handling the singularity [35,36] when ${(x-t)}^{-(1-\alpha )}\to \infty$ as $t\to x$. We propose a convolutional regularization approach:

$${(x-t)}^{-(1-\alpha )}\approx {[{(x-t)}^2+{ϵ}_h^2]}^{-(1-\alpha )/2}$$

(13)

where ${ϵ}_h$ is a small adaptive regularization parameter that prevents numerical overflow while preserving the singular behavior. The associated local regularization error is defined by:

$$E_{\mathrm{reg}}(x,t)={(x-t)}^{-(1-\alpha )}-{[{(x-t)}^2+{ϵ}_h^2]}^{-(1-\alpha )/2}$$

For $|x-t|\gg {ϵ}_h$, this error admits the expansion:

$$E_{\mathrm{reg}}(x,t)={\displaystyle \frac{(1-\alpha ){ϵ}_h^2}{2{(x-t)}^{3-\alpha }}}+{\displaystyle \frac{(1-\alpha )(3-\alpha ){ϵ}_h^4}{8{(x-t)}^{5-\alpha }}}+\mathcal{O}\left({\displaystyle \frac{{ϵ}_h^6}{{(x-t)}^{7-\alpha }}}\right)$$

By binomial expansion with $u={\displaystyle \frac{{ϵ}_h^2}{{(x-t)}^2}}\ll 1$, we have:

$${[{(x-t)}^2+{ϵ}_h^2]}^{-(1-\alpha )/2}={|x-t|}^{-(1-\alpha )}{\left(1+{\displaystyle \frac{{ϵ}_h^2}{{(x-t)}^2}}\right)}^{-(1-\alpha )/2}$$

Using the Taylor expansion:

$${(1+u)}^{-(1-\alpha )/2}=1-{\displaystyle \frac{1-\alpha }{2}}u+{\displaystyle \frac{(1-\alpha )(3-\alpha )}{8}}{u}^2+\mathcal{O}\left(u^3\right)$$

Substituting u and simplifying yields the desired expansion. We then define the relative regularization error as:

$${\displaystyle \frac{E_{\mathrm{reg}}(x,t)}{{(x-t)}^{-(1-\alpha )}}}={\displaystyle \frac{(1-\alpha ){ϵ}_h^2}{2{(x-t)}^2}}+\mathcal{O}\left({\displaystyle \frac{{ϵ}_h^4}{{(x-t)}^4}}\right)$$

Furthermore, if the error in the data is $\delta g$, then the propagated error in the solution satisfies the bound:

$${\parallel \delta f\parallel }_{L^2}\le {\displaystyle \frac{1}{{ϵ}_h^{(1-\alpha )/2}}}{\parallel \delta g\parallel }_{L^2}$$

This relationship demonstrates that values of ${ϵ}_h$ that are too small amplify data errors.

To concentrate quadrature points near the singularity, we employ a graded mesh transformation:

$$t_j=x·{\left({\displaystyle \frac{j}{n}}\right)}^{\beta },j=0,1,\dots ,n$$

(14)

with $\beta >1$ to increase point density near $t=0$. The adaptive weights [35,41] are defined as:

$$w_j(x,t_j)=1+Cexp\left(-\gamma |x-t_j|\right)$$

(15)

where C and $\gamma$ are positive parameters that enhance accuracy near the singular point. The local density of quadrature points is characterized by the spacing between consecutive points:

$$\Delta t_j=t_{j+1}-t_j=x\left[{\left({\displaystyle \frac{j+1}{n}}\right)}^{\beta }-{\left({\displaystyle \frac{j}{n}}\right)}^{\beta }\right]$$

Using the mean value theorem, there exists ${\xi }_j\in (j,j+1)$ such that:

$$\Delta t_j=x·{\displaystyle \frac{\beta }{n}}·{\left({\displaystyle \frac{{\xi }_j}{n}}\right)}^{\beta -1}·{\displaystyle \frac{1}{n}}={\displaystyle \frac{x\beta }{n^{\beta +1}}}{\xi }_j^{\beta -1}$$

We define the relative density of points near the singularity as:

$${\displaystyle \frac{\Delta t_1}{\Delta t_n}}={\left({\displaystyle \frac{1}{n}}\right)}^{\beta -1}$$

Indeed, for $j=1$, ${\displaystyle \Delta t_1\approx \frac{x\beta }{n^{\beta +1}}}$, and for $j=n$, ${\displaystyle \Delta t_n\approx \frac{x\beta }{n^2}}$. Therefore, ${\displaystyle \frac{\Delta t_1}{\Delta t_n}\approx \frac{1}{n^{\beta -1}}}$. This relationship demonstrates that the larger $\beta$ becomes, the greater the concentration of points near the singularity. For the singular integral

$$I={\int }_0^x{\displaystyle \frac{k(x,t)F\left(f\right(t\left)\right)}{{(x-t)}^{1-\alpha }}}dt$$

the quadrature error with the graded mesh transformation is given by:

$$E_{\beta }^{\left(n\right)}=I-\sum _{j=0}^{n-1}{w}_j{\displaystyle \frac{k(x,t_j)F\left(f\left(t_j\right)\right)}{{(x-t_j)}^{1-\alpha }}}$$

This error satisfies the bound [13,30]

$$|E_{\beta }^{\left(n\right)}|\le C_0·n^{-min(2,{\displaystyle \frac{\beta (1+\alpha )}{1+\alpha }})}·{\parallel g^{\left(p\right)}\parallel }_{\infty }$$

(16)

where $p=\lfloor min(2,{\displaystyle \frac{\beta (1+\alpha )}{1+\alpha }})\rfloor$, $C_0$ is a constant independent of n, and g is the function defined by ${\displaystyle g\left(t\right)=\frac{k(x,t)F\left(f\right(t\left)\right)}{{(x-t)}^{1-\alpha }}}$ with $0<\alpha <1$.

When we perform the transformation $s={\left({\displaystyle \frac{t}{x}}\right)}^{1/\beta }$, we obtain $t=x{s}^{\beta }$ and $dt=\beta x{s}^{\beta -1}ds$. The transformed integral becomes:

$$I=x^{\alpha }{\int }_0^1{\displaystyle \frac{k(x,x{s}^{\beta })F\left(f\left(x{s}^{\beta }\right)\right)}{{(1-s^{\beta })}^{1-\alpha }}}\beta s^{\beta -1}ds$$

The singularity at $s=1$ has a modified order of ${(1-s^{\beta })}^{-(1-\alpha )}$, and the error depends on the regularity of the transformed integrand. We apply a uniform composite quadrature formula in s (trapezoid in this case), of order p (typically $p\le 2$). The error is proportional to:

$$E_n\sim {\displaystyle \frac{1}{n^p}}\underset{0\le s\le 1}{max}\left|h^{\left(p\right)}\left(s\right)\right|.$$

By repeated derivation of $h\left(s\right)=\beta s^{\beta -1}{(1-s^{\beta })}^{-(1-\alpha )}·k(·)·F\left(f(·)\right),$ we obtain that each derivative $h^{\left(p\right)}\left(s\right)$ generates dominant terms of the type $C{(1-s^{\beta })}^{-(1-\alpha )-m}·s^{\beta -1-p},m\le p.$ For these terms to remain bounded near $s=1$, we require $(1-\alpha )+m\le \beta m,\forall m\le p,$ the most restrictive condition being for $m=p$. Thus we need:

$$\beta \ge {\displaystyle \frac{(1-\alpha )+p}{p}}.$$

In the context $p\le 2$, we see that the real order of convergence is:

$$min\left(p,{\displaystyle \frac{\beta (1+\alpha )}{1+\alpha }}\right)=min\left(2,{\displaystyle \frac{\beta (1+\alpha )}{1+\alpha }}\right).$$

Hence the origin of the bound ${\displaystyle |E_{\beta }^{\left(n\right)}|\le C_0{n}^{-min(2,\frac{\beta (1+\alpha )}{1+\alpha })},{\parallel g^{\left(p\right)}\parallel }_{\infty }.}$ The optimal convergence order is achieved for ${\displaystyle {\beta }_{\mathrm{opt}}=\frac{2-\alpha }{1-\alpha }}$. This optimal value is obtained by maximizing the convergence exponent in Equation (16). The aim is to maximize the exponent:

$$\gamma \left(\beta \right)=min\left(2,{\displaystyle \frac{\beta (1+\alpha )}{1+\alpha }}\right).$$

If ${\displaystyle \beta \le \frac{2(1+\alpha )}{1+\alpha }=2}$, then ${\displaystyle \gamma \left(\beta \right)=\frac{\beta (1+\alpha )}{1+\alpha }}$ and if $\beta \ge 2$, then $\gamma \left(\beta \right)=2$. However, the derivability bound restricts $\beta$ to $\beta \le {\displaystyle \frac{(1-\alpha )+p}{p}}$. For $p=2$, this gives:

$$\beta \le {\displaystyle \frac{2-\alpha }{2}}+{\displaystyle \frac{1-\alpha }{2}}={\displaystyle \frac{2-\alpha }{1-\alpha }}.$$

Therefore, the tipping point where ${\displaystyle \frac{\beta (1+\alpha )}{1+\alpha }}=2$ and the maximum derivability are equal is indeed

$${\beta }_{\mathrm{opt}}={\displaystyle \frac{2-\alpha }{1-\alpha }}.$$

(17)

We then employ importance sampling by drawing samples from the prior $p\left(f\right)$ and weighting them by the likelihood. The importance weight for sample i is:

$$w_i={\displaystyle \frac{p\left(g_{\mathrm{obs}}\right|f_i)}{p\left(f_i\right)}}\propto p\left(g_{\mathrm{obs}}\right|f_i)$$

(18)

To prevent numerical overflow in the exponential function, we apply log-space stabilization:

$$ln{w}_i=lnp\left(g_{\mathrm{obs}}\right|f_i)-\underset{j}{max}lnp\left(g_{\mathrm{obs}}\right|f_j)$$

(19)

2.3. Posterior Approximation and Reconstruction

After computing the importance weights ${\left\{w_i\right\}}_{i=1}^N$ for all N prior samples ${\left\{f_i\right\}}_{i=1}^N$, we normalize them:

$${\tilde{w}}_i={\displaystyle \frac{w_i}{{\sum }_{j=1}^N{w}_j}}$$

(20)

The posterior distribution is then approximated by the weighted empirical distribution:

$$p\left(f\right|g_{\mathrm{obs}})\approx \sum _{i=1}^N{\tilde{w}}_i{\delta }_{f_i}\left(f\right)$$

(21)

where ${\delta }_{f_i}$ denotes the Dirac delta function centered at $f_i$.

The posterior mean function is computed as:

$$\mathbb{E}\left[f\left(x\right)\right|g_{\mathrm{obs}}]=\sum _{i=1}^N{\tilde{w}}_i{f}_i\left(x\right)$$

(22)

The posterior variance is given by:

$$\mathrm{Var}\left[f\left(x\right)\right|g_{\mathrm{obs}}]=\sum _{i=1}^N{\tilde{w}}_i{\left[f_i\left(x\right)\right]}^2-{\left[\sum _{i=1}^N{\tilde{w}}_i{f}_i\left(x\right)\right]}^2$$

(23)

For uncertainty quantification, we construct pointwise credible intervals. At each point x, we sort the weighted samples ${\left\{(f_i\left(x\right),{\tilde{w}}_i)\right\}}_{i=1}^N$ and compute the $100(1-{\alpha }_0)\%$ credible interval as:

$${\mathrm{CI}}_{1-{\alpha }_0}\left(x\right)=[Q_{{\alpha }_0/2}\left(x\right),Q_{1-{\alpha }_0/2}\left(x\right)]$$

(24)

where $Q_p\left(x\right)$ is the p-th weighted quantile of the posterior samples at point x and ${\alpha }_0$ is the confidence level.

2.3.1. Convergence Assessment

The quality of the approximation can be assessed through several complementary metrics that provide insights into different aspects of the reconstruction performance. The Effective Sample Size (ESS) [37,38], defined as ${\displaystyle \mathrm{ESS}=\frac{1}{{\sum }_{i=1}^N{\tilde{w}}_i^2}}$, serves as a crucial indicator of sampling efficiency. When the ESS is low, it suggests that the importance weights are poorly distributed, with only a few samples carrying most of the weight, thereby indicating poor sampling efficiency and potential numerical instability.

Additionally, residual analysis provides a direct measure of how well the reconstructed solution satisfies the original integral equation. The residual function $R\left(x\right)=g_{\mathrm{obs}}\left(x\right)-\mathcal{L}\left[\mathbb{E}\left[f\right|g_{\mathrm{obs}}]\right]\left(x\right)$ quantifies the discrepancy between the observed data and the Volterra operator applied to the posterior mean. Small residuals across the domain indicate that the reconstructed solution provides a good fit to the integral equation, while large or systematic residuals may suggest inadequate prior specification or insufficient sampling.

Furthermore, the Monte Carlo Standard Error (MCSE) quantifies the uncertainty in our posterior estimates due to the finite sample approximation. Given by ${\displaystyle \mathrm{MCSE}\left(x\right)=\sqrt{\frac{\mathrm{Var}\left[f\left(x\right)\right|g_{\mathrm{obs}}]}{\mathrm{ESS}}}}$, this metric provides a measure of the precision of our posterior mean estimates. A decreasing MCSE with increasing sample size indicates convergence of the Monte Carlo approximation, while persistently large MCSE values suggest the need for additional samples or improved sampling strategies.

2.3.2. Proposed Algorithm

The complete Bayesian reconstruction algorithm [38,39] proceeds as follows:

• Initialization: Set up the GP prior with chosen hyperparameters $({\sigma }_f,\ell )$ and discretization grid ${\left\{x_i\right\}}_{i=1}^n$.
• Prior Sampling: Generate N function samples ${\left\{f_i\right\}}_{i=1}^N$ from the GP prior.
• Likelihood Evaluation: For each sample $f_i$:
  • Compute $\mathcal{L}\left[f_i\right]\left(x_j\right)$ using adaptive quadrature with graded mesh transformation
  • Apply singularity regularization: ${(x-t)}^{-(1-\alpha )}\approx {[{(x-t)}^2+{ϵ}_h^2]}^{-(1-\alpha )/2}$
  • Evaluate $lnp\left(g_{\mathrm{obs}}\right|f_i)$ using the likelihood formula
• Weight Computation: Calculate importance weights $w_i$ with log-space stabilization:

  $$ln{w}_i=lnp\left(g_{\mathrm{obs}}\right|f_i)-\underset{j}{max}lnp\left(g_{\mathrm{obs}}\right|f_j)$$
• Posterior Approximation: Normalize weights ${\tilde{w}}_i={\displaystyle \frac{w_i}{{\sum }_{j=1}^N{w}_j}}$ and compute posterior statistics:
  • Posterior mean: $\mathbb{E}\left[f\left(x\right)\right|g_{\mathrm{obs}}]={\sum }_{i=1}^N{\tilde{w}}_i{f}_i\left(x\right)$
  • Posterior variance: $\mathrm{Var}\left[f\left(x\right)\right|g_{\mathrm{obs}}]={\sum }_{i=1}^N{\tilde{w}}_i{\left[f_i\left(x\right)\right]}^2-{\left[{\sum }_{i=1}^N{\tilde{w}}_i{f}_i\left(x\right)\right]}^2$
• Reconstruction: Output the posterior mean as the solution estimate along with pointwise credible intervals ${\mathrm{CI}}_{1-{\alpha }_0}\left(x\right)=[Q_{{\alpha }_0/2}\left(x\right),Q_{1-{\alpha }_0/2}\left(x\right)]$ for uncertainty quantification.

2.4. Convergence and Stability

Let $C[0,1]$ be the space of continuous functions on $[0,1]$ equipped with the uniform norm ${\parallel f\parallel }_{\infty }={max}_{x\in [0,1]}\left|f\left(x\right)\right|$. We define the weighted Sobolev space by:

$$H_{\alpha }^s[0,1]=\{f\in L^2{[0,1]:\parallel f\parallel }_{H_{\alpha }^s}<\infty \}$$

where

$${\parallel f\parallel }_{H_{\alpha }^s}^2={\int }_0^1{\left|f\left(x\right)\right|}^2{x}^{2\alpha }dx+{\int }_0^1{\left|f^{\left(s\right)}\left(x\right)\right|}^2{x}^{2\alpha }dx$$

We assume kernel regularity conditions, namely that the kernel $k(x,t)$ satisfies: $k\in C^2([0,1]\times [0,1])$, $\left|k(x,t)\right|\le M_k$ for some constant $M_k>0$, and ${\displaystyle |\frac{\partial k}{\partial x}(x,t)|+|\frac{\partial k}{\partial t}(x,t)|\le M_k^{\prime }}$.

For ${ϵ}_h>0$ and the regularization ${(x-t)}^{-(1-\alpha )}\approx {[{(x-t)}^2+{ϵ}_h^2]}^{-(1-\alpha )/2}$, we have:

$$\left|{(x-t)}^{-(1-\alpha )}-{[{(x-t)}^2+{ϵ}_h^2]}^{-(1-\alpha )/2}\right|\le C_{\alpha }{ϵ}_h^{1-\alpha }{|x-t|}^{-(2-\alpha )}$$

for $|x-t|\ge {ϵ}_h$, where $C_{\alpha }$ depends only on $\alpha$.

Using Taylor expansion and the mean value theorem, for $u={(x-t)}^2$ and $v={(x-t)}^2+{ϵ}_h^2$, we have:

$$u^{-(1-\alpha )/2}-v^{-(1-\alpha )/2}=-{\displaystyle \frac{1-\alpha }{2}}{\xi }^{-(3-\alpha )/2}{ϵ}_h^2$$

where $\xi \in [u,v]$. The result follows by estimating $\xi$.

2.4.1. Monte Carlo Estimation and Stability Analysis

The Monte Carlo estimator of the posterior mean satisfies:

$$\mathbb{E}\left[{\left|{\displaystyle \frac{1}{N}}\sum _{i=1}^N{\tilde{w}}_i{f}_i\left(x\right)-\mathbb{E}\left[f\left(x\right)\right|g^{\mathrm{obs}}]\right|}^2\right]\le {\displaystyle \frac{C}{N}}$$

where C depends on ${\sigma }_f$, ℓ, and the problem parameters, but is independent of N. Using the variance decomposition:

$$\mathrm{Var}\left[\sum _{i=1}^N{\tilde{w}}_i{f}_i\left(x\right)\right]=\mathrm{Var}\left[{\displaystyle \frac{{\sum }_{i=1}^N{w}_i{f}_i\left(x\right)}{{\sum }_{i=1}^N{w}_i}}\right]$$

By the Cauchy-Schwarz inequality and Gaussian process properties:

$$\mathbb{E}\left[w_i^2\right]\le \mathbb{E}{\left[w_i\right]}^2·exp\left(C{\sigma }_f^2\right)$$

The strong law of large numbers applies since $\mathbb{E}\left[w_i{f}_i\left(x\right)\right]<\infty$ under our assumptions.

2.4.2. Stability with Respect to Data Perturbations

Let $g_1,g_2$ be two data functions with $\parallel g_1-g_2{\parallel }_{\infty }\le \delta$. Then the corresponding posterior means ${\mu }_1,{\mu }_2$ satisfy:

$$\parallel {\mu }_1-{\mu }_2{\parallel }_{\infty }\le C\delta$$

where C depends on the hyperparameters but not on $\delta$. The log-likelihoods satisfy:

$$|logp\left(g_1\right|f)-logp\left(g_2\right|f)|\le {\displaystyle \frac{\parallel g_1-g_2{\parallel }_{\infty }}{2{\sigma }^2}}\left(\parallel g_1{\parallel }_{\infty }+\parallel g_2{\parallel }_{\infty }+2{\parallel L\left[f\right]\parallel }_{\infty }\right)$$

By continuity of the operator L on $C[0,1]$ and compactness of the Gaussian process support, the result follows.

2.4.3. Minimax Risk Analysis

For the function class $\Theta =\{f\in H_{\alpha }^s[0,1]:\parallel f{\parallel }_{H_{\alpha }^s}\le M\}$, the minimax risk satisfies:

$$\underset{\widehat{f}}{inf}\underset{f\in \Theta }{sup}\mathbb{E}[\parallel \widehat{f}-f{\parallel }_{\infty }^2]\ge c{n}^{-{\displaystyle \frac{2s}{2s+1+2\alpha }}}$$

where $c>0$ is a constant. Under the optimal choice of hyperparameters $\ell \sim n^{-{\displaystyle \frac{1}{2s+1+2\alpha }}}$ and with $N=O\left(n^2\right)$ samples, our Bayesian estimator ${\widehat{f}}_n$ satisfies:

$$\underset{f\in \Theta }{sup}\mathbb{E}[\parallel {\widehat{f}}_n-f{\parallel }_{\infty }^2]\le C{n}^{-{\displaystyle \frac{2s}{2s+1+2\alpha }}}$$

3. Numerical Examples

To validate the effectiveness of the proposed Bayesian approach and demonstrate its capability in handling weakly singular Volterra integral equations, we present a series of numerical experiments. These simulations are designed to assess the accuracy, convergence properties, and robustness of our method under various conditions. For the quadrature points, we calculate $\beta$, using his optimal equation form in (17).

3.1. Example 1: Abel Integral Equation

We begin our numerical investigation with a classical Abel integral equation, which serves as a fundamental benchmark for methods dealing with weakly singular kernels. This example allows us to evaluate the performance of our Bayesian inverse approach against well-established analytical and numerical solutions.

We consider the following Abel integral equation [42] of the second kind:

$$f\left(x\right)=1-{\int }_0^x{\displaystyle \frac{f\left(x\right)}{{(x-t)}^{0.5}}}dt,0\le x\le 1$$

(25)

This equation represents a specific instance of our general formulation with the parameter values $\lambda =1$, $\alpha =0.5$, and the right-hand side function $g\left(x\right)=1$. The kernel function $k(x,t)=1$ and the nonlinear function $F\left(f\right(t\left)\right)=f\left(t\right)$, making this a linear Abel equation with a characteristic weak singularity of order $1/2$ at $t=x$.

The choice of $\alpha =0.5$ corresponds to the classical Abel kernel, which exhibits the typical square-root singularity that appears in numerous physical applications, including particle physics, plasma dynamics, and image reconstruction problems. This test case serves as an ideal benchmark for validating our numerical approach due to its analytical tractability. The exact solution of this equation is ${\displaystyle f\left(x\right)=\frac{1}{1+2\sqrt{x}}}$.

Figure 1 and Figure 2 showcases the compelling results obtained from our comprehensive numerical experiments for Example 1. The striking visual evidence presented in Figure 1 reveals an exceptional alignment between the approximate and exact solutions, providing unequivocal validation of our method’s remarkable effectiveness. This outstanding performance is further substantiated by the quantitative analysis in Table 1, where our approach demonstrates extraordinary precision with an average absolute error of merely ${10}^{-6}$ a testament to the method’s superior computational accuracy.

To comprehensively evaluate the influence of the critical hyperparameters $N_{\mathrm{points}}$, C, $\gamma$, and ${\sigma }_f$, we conducted an extensive parametric study across diverse computational scenarios. Figure 2 provides crucial insights into the evolution of both mean squared error and residual norm as functions of $N_{\mathrm{points}}$. The results demonstrate a clear and consistent pattern: as the number of quadrature points increases, we observe a remarkable dual improvement in both solution accuracy and convergence behavior, highlighting the robust scalability of our approach.

The hyperparameter sensitivity analysis reveals fascinating computational dynamics. Table 2 demonstrates that increasing ${\sigma }_f$ values lead to substantial MSE reduction, thereby significantly enhancing convergence characteristics. Similarly, Table 3 corroborates this trend for the $\gamma$ parameter, indicating optimal performance tuning capabilities. However, Table 4 reveals the expected computational trade-off: while larger C values may offer certain advantages, they inevitably result in increased computational overhead a critical consideration for practical implementation.

Most importantly, the comparative analysis presented in Table 1 illustrates the comparative performance metrics, indicating improved accuracy of the proposed method across all test points. Our proposed technique consistently outperforms both the collocation method and the Legendre wavelet approach, achieving significantly lower absolute errors across all test cases. This decisive advantage unambiguously establishes the superior performance and enhanced reliability of our methodology compared to these well-established competing methods. The substantial error reduction achieved by our approach not only validates its theoretical foundations but also demonstrates its practical value for high-precision computational applications.

3.2. Example 2: Nonlinear Volterra Equation with Polynomial Right-Hand Side

To further evaluate the robustness and versatility of our Bayesian approach, we consider a more challenging nonlinear Volterra integral equation that incorporates both the characteristic weak singularity and a complex polynomial forcing term. This test case allows us to assess the method’s performance in handling nonlinear kernels while maintaining computational efficiency.

We examine the following nonlinear Volterra integral equation [42]:

$$f\left(x\right)-{\int }_0^x{\displaystyle \frac{xt}{{(x-t)}^{0.5}}}dt=x^3-4096x^{8.5}$$

(26)

The exact solution to this equation is known to be $f\left(x\right)=x^3$, which provides a clear benchmark for accuracy assessment. We implemented our Bayesian algorithm to solve this problem and conducted extensive numerical experiments to evaluate its performance across different parameter configurations.

The results showcased in Figure 3 and Figure 4 for Example 2 present an even more impressive demonstration of our method’s exceptional capabilities. The comprehensive analysis includes both the direct comparison between exact and approximate solutions (Figure 3) and the detailed performance evolution as a function of computational points (Figure 4). The visual inspection of Figure 3 reveals a virtually flawless correspondence between the exact solution and our method’s approximation a remarkable achievement that underscores the method’s extraordinary precision and reliability.

This exceptional visual agreement is powerfully reinforced by the quantitative metrics presented in Table 5, which demonstrates our approach’s outstanding computational excellence with an average absolute error reaching an impressive ${10}^{-7}$—representing a full order of magnitude improvement compared to Example 1. This exceptional accuracy level positions our method at the forefront of high-precision numerical techniques.

Consistent with our findings from Example 1, the comprehensive hyperparameter analysis reveals intriguing computational characteristics. Figure 4 corroborates our previous observations, confirming that convergence performance exhibits consistent enhancement as $N_{\mathrm{points}}$ increases, thereby validating the method’s robust scalability across different problem domains.

Remarkably, the hyperparameter sensitivity analysis for Example 2 unveils a fascinating aspect of our method’s robustness. Table 6, Table 7 and Table 8 demonstrate that the MSE and maximum absolute error remain virtually unaffected by variations in the hyperparameters ${\sigma }_f$, C, and $\gamma$. This exceptional stability indicates that our method possesses inherent robustness against parameter fluctuations, with only marginal variations in computational time observed across different parameter configurations due to intermediate calculations that influence the overall computational overhead.

The definitive validation of our method’s supremacy emerges from the comprehensive comparative analysis presented in Table 5. Our proposed approach decisively outperforms all competing methodologies, consistently achieving the lowest absolute errors across all evaluation metrics. This comprehensive superiority not only validates the theoretical foundations of our approach but also establishes its undisputed dominance in the field of high-precision numerical computations, making it the method of choice for demanding computational applications.

3.3. Limitations

Despite the exceptional performance demonstrated by our approach, it is important to acknowledge certain limitations that define its current scope of application. First, the proposed method is specifically designed to handle Volterra integral equations with weak singularities ($0<\alpha <1$) and has not been validated for strong singularity cases ($\alpha \ge 1$), thereby restricting its applicability to this particular class of problems. Second, scalability considerations arise as increasing the number of samples N required for adequate convergence leads to substantial growth in computational cost, which becomes particularly challenging for high-dimensional problems. Finally, the importance sampling algorithm may experience weight degeneracy in ill-conditioned problems, where only a few samples carry the majority of the statistical weight, potentially compromising numerical stability and reliability of posterior estimates.

4. Conclusions

In this study, we have proposed a novel numerical approach for solving Volterra’s integral equations with weak singularities by approximating the singular kernel through convolutional regularization and employing RBF kernels within a Gaussian process prior framework. A graded mesh transformation technique was implemented to achieve strategic sampling that concentrates the quadrature points around the singular regions of the equation, thereby enhancing computational efficiency and accuracy.

The comprehensive numerical experiments conducted across multiple test cases have demonstrated the robustness and effectiveness of our proposed methodology. The results consistently show strong convergence properties, achieving remarkable accuracy with absolute errors ranging from ${10}^{-6}$ to ${10}^{-7}$. These performance metrics validate the theoretical foundations of our approach and confirm its practical applicability for solving weakly singular Volterra’s integral equations.

Furthermore, our parametric studies revealed that the method exhibits stable behavior across various hyperparameter configurations, with the number of quadrature points $N_{\mathrm{points}}$ being the most influential factor for convergence enhancement. The computational efficiency of the algorithm, combined with its high accuracy, makes it a valuable tool for practitioners dealing with singular integral equations in various scientific and engineering applications.

Future work will focus on extending the method to strongly singular Volterra equations, as well as exploring adaptive mesh refinement and automated hyperparameter tuning using machine learning.