
 
%%Manipulating the above equation system, we determine the hybrid coefficients c_{i0 }.

%%
%% Copyright 2007, 2008, 2009 Elsevier Ltd
%%
%% This file is part of the 'Elsarticle Bundle'.
%% ---------------------------------------------
%%
%% It may be distributed under the conditions of the LaTeX Project Public
%% License, either version 1.2 of this license or (at your option) any
%% later version.  The latest version of this license is in
%%    http://www.latex-project.org/lppl.txt
%% and version 1.2 or later is part of all distributions of LaTeX
%% version 1999/12/01 or later.
%%
%% The list of all files belonging to the 'Elsarticle Bundle' is
%% given in the file `manifest.txt'.
%%

%% Template article for Elsevier's document class `elsarticle'
%% with numbered style bibliographic references
%% SP 2008/03/01
%%
%%
%%
%% $Id: elsarticle-template-num.tex 4 2009-10-24 08:22:58Z rishi $
%%
%%

%\begin{filecontents*}{fredholm_reference_patch.bib}
%@article{sathar2018numerical,
%  author  = {Sathar, Mohammad Hasan Abdul and Rasedee, Ahmad Fadly Nurullah and Ahmedov, Anvarjon A. and Asbullah, Muhammad Asyraf},
%  title   = {Numerical Integration Based on Linear Legendre Multi Wavelets},
%  journal = {Journal of Physics: Conference Series},
%  volume  = {1132},
%  pages   = {012008},
%  year    = {2018},
%  doi     = {10.1088/1742-6596/1132/1/012008}
%}
%\end{filecontents*}

\documentclass[final,3p,times]{elsarticle}

%% Use the option review to obtain double line spacing
%% \documentclass[preprint,review,12pt]{elsarticle}

%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%% \documentclass[final,1p,times]{elsarticle}
%% \documentclass[final,1p,times,twocolumn]{elsarticle}
%% \documentclass[final,3p,times]{elsarticle}
%% \documentclass[final,3p,times,twocolumn]{elsarticle}
%% \documentclass[final,5p,times]{elsarticle}
%% \documentclass[final,5p,times,twocolumn]{elsarticle}

%% if you use PostScript figures in your article
%% use the graphics package for simple commands
%%\usepackage{graphics}
%% or use the graphicx package for more complicated commands
\usepackage{graphicx}
%% or use the epsfig package if you prefer to use the old commands
%% \usepackage{epsfig}

%% The amssymb package provides various useful mathematical symbols
\usepackage{amssymb}

%% The amsthm package provides extended theorem environments
%% \usepackage{amsthm}

%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers after \end{frontmatter}.
%% \usepackage{lineno}

%% natbib.sty is loaded by default. However, natbib options can be
%% provided with \biboptions{...} command. Following options are
%% valid:

%%   round  -  round parentheses are used (default)
%%   square -  square brackets are used   [option]
%%   curly  -  curly braces are used      {option}
%%   angle  -  angle brackets are used    <option>
%%   semicolon  -  multiple citations separated by semi-colon
%%   colon  - same as semicolon, an earlier confusion
%%   comma  -  separated by comma
%%   numbers-  selects numerical citations
%%   super  -  numerical citations as superscripts
%%   sort   -  sorts multiple citations according to order in ref. list
%%   sort&compress   -  like sort, but also compresses numerical citations
%%   compress - compresses without sorting
%%
%% \biboptions{comma,round}

% \biboptions{}


  \usepackage{newtxtext,newtxmath}
\usepackage{bm}


 \usepackage{amsmath}
 \usepackage{verbatim}
 \usepackage{booktabs}
 \usepackage{tabularx}
 \usepackage{multirow}
 \usepackage{siunitx}
 \usepackage{mathtools}
 \usepackage{hyperref}

 \hypersetup{colorlinks=true,linkcolor=black,citecolor=black,urlcolor=black}
\pdfstringdefDisableCommands{%
  \def\corref#1{}%
  \def\cnotenum#1{}%
  \def\notenum#1{}%
  \def\tnoteref#1{}%
  \def\thanks#1{}%
}


\journal{Applied Numerical Mathematics}

\begin{document}




\begin{frontmatter}

\title{Regularized Multiwavelet Collocation Methods for Ill-Posed Fredholm Integral Equations of the First Kind}
%\tnotetext[label0]{This is only an example}


\author[label1]{Ahmed Kaffel \corref{cor1}}
\address[label1]{Department of Mathematical Sciences\\
University of Wisconsin--Milwaukee\\
Milwaukee, WI 53211, USA}

\cortext[cor1]{Corresponding author}

\ead{kaffel@uwm.edu}




\begin{abstract}
Fredholm integral equations of the first kind arise naturally in inverse problems where the available measurements are indirect, smoothed, or contaminated by noise. Because the underlying integral operators are compact, direct discretization often produces unstable reconstructions unless an appropriate regularization strategy is introduced. This paper develops a unified regularized multiwavelet-collocation framework for such ill-posed problems and uses it to compare Linear Legendre multiwavelets, Legendre wavelets, Haar wavelets, B-spline wavelets, and hybrid wavelet bases under the same computational protocol. The discrete first-kind system is stabilized by Tikhonov-type regularization, with the regularization parameter selected a posteriori through Morozov's discrepancy principle. The resulting methods are tested on benchmark integration problems and one- and two-dimensional Fredholm equations, with attention not only to approximation error but also to residual behavior, conditioning, noise sensitivity, Morozov-selected regularization parameters, and computational cost. The results show that basis selection is not merely an approximation issue: it directly affects the stability and regularization behavior of the inverse recovery. Higher-order local bases, especially Legendre, B-spline, and Linear Legendre multiwavelet discretizations, are most effective for smooth reconstructions, whereas Haar-type bases remain useful when simplicity, low cost, or nonsmooth features are important. The results demonstrate that basis selection should be viewed as part of the regularization strategy itself rather than solely as an approximation choice. Beyond identifying effective bases for different smoothness regimes, the study provides practical guidelines for selecting approximation spaces in regularized inverse problems by balancing accuracy, stability, conditioning, robustness to noise, and computational cost.
\end{abstract}

\begin{keyword}
Fredholm integral equations of the first kind \sep
Ill-posed inverse problems \sep
Tikhonov regularization \sep
Morozov discrepancy principle \sep
Multiwavelet collocation \sep
Inverse problem stabilization \sep
Conditioning analysis
\end{keyword}

\end{frontmatter}

%%
%% Start line numbering here if you want
%%
% \linenumbers

%% main text

\section{Introduction}
\label{sec1}


Integral equations play an important role in many scientific and engineering problems \cite{lonseth1977sources,kress2000numerical}. In inverse problems, Fredholm equations of the first kind arise when the measured data are indirect or smoothed observations of an unknown quantity. Typical examples include computerized tomography \cite{natterer2001mathematics}, image restoration \cite{banham1997digital}, inverse scattering and electromagnetic imaging \cite{colton1996simple,kress1995integral}, and other ill-posed recovery problems described in classical inverse-problem references \cite{engl2014inverse,hansen2010discrete,vogel2002computational}. These examples motivate stable discretizations that account explicitly for the smoothing effect of the integral operator and the presence of data perturbations.

In this paper, we consider the Fredholm integral equation of the first kind of the following generic form:
\begin{equation} \label{eq1}
    \int_{a}^{b} \mathcal{K}(x,y)f(y)dy= g(x), \hspace{5mm} c \leq x \leq d, 
\end{equation}

where the functions $g(x)$ and $\mathcal{K}(x,y)$ represent the data and the kernel, respectively, and the unknown function $f(y)$ is the solution to be determined. Throughout the paper, $\mathcal{K}$ denotes the continuous kernel, $K$ denotes the induced operator or a discrete coefficient matrix when the context is algebraic, and $\alpha$ denotes the Tikhonov regularization parameter. Note that Eq. (\ref{eq1}) can be expressed as 
\begin{equation} \label{eq1a}
Kf= g, 
\end{equation}
where the operator 
$$Kf= \int_{a}^{b} \mathcal{K}(x,y)f(y)dy$$
is linear; here the kernel is assumed to be smooth.

According to Hadamard, a problem is well-posed if it has a solution (existence), not more than one solution (uniqueness), and the solution depends continuously on the data (stability). Problems are called ill-posed if they are not well-posed, i.e. small changes in the data cause large changes in the unknown function $f$ \cite{groetsch1984theory}. This loss of continuous dependence is the main difficulty in first-kind Fredholm equations, and it explains why a direct discretization of Eq. (\ref{eq1}) can be unstable. For a discussion of existence, uniqueness, and stability issues for integral equations and inverse problems, see \cite{groetsch1984theory,wing1991primer,kress2000numerical}.

Wavelet and polynomial bases have been used extensively to discretize integral equations. Relevant examples include Haar and block-pulse bases, Legendre and Chebyshev wavelets, B-spline wavelets, rationalized Haar functions, hybrid block-pulse--Legendre bases, and operational-matrix formulations \cite{hsiao2009hybrid,sohrabi2013study,maleknejad2004numerical,reihani2007rationalized,maleknejad2010method,adibi2010chebyshev,maleknejad2007numerical,lin2009haar,aziz2010comparative,wang2008numerical}. These references are retained because they are directly connected with the wavelet discretizations and quadrature formulas used below. In contrast, broad nonlinear solvers such as homotopy, variational-iteration, and decomposition methods are not emphasized here, since the present paper is focused on regularized wavelet collocation for noisy first-kind inverse recovery.

Numerical treatment of first-kind integral equations is more delicate than that of second-kind equations because compact integral operators typically produce ill-conditioned algebraic systems. If the data are perturbed by a noise level $\delta$, small changes in the observation data may cause large changes in the recovered solution. Thus a feasible numerical method must combine approximation with regularization and with a reliable strategy for choosing the regularization parameter. Classical regularization ideas go back to Phillips and Tikhonov \cite{phillips1962technique,tikhonov1963regularization,tikhonov1963solution,tikhonov1987ill,tikhonov2013numerical}, and modern computational treatments are described in foundational references such as \cite{engl1996regularizatione,hansen1992analysis,hansen2010discrete,vogel2002computational,morozov2012methods}. Parameter-choice rules, including Morozov-type discrepancy principles and related Newton or quasi-Newton procedures, have been studied in \cite{baumeister1987stable,cheng2000one,kunisch1998iterative,dennis1996numerical,ito1992choice,engl2000new,kunisch1993class,groetsch1993inverse,wang2001fast,zou2009cubically}. These older references are kept because they support the regularization and parameter-selection mechanism used in the paper.

The present work applies a Tikhonov-type regularization strategy to first-kind Fredholm equations and combines it with several wavelet-collocation discretizations. The regularization transforms the unstable first-kind formulation into a stabilized algebraic problem whose solution can be computed reliably even when the right-hand side is contaminated by measurement noise. The regularization parameter is selected using an a posteriori strategy based on Morozov's discrepancy principle, implemented through the cubic-convergence procedure of \cite{wang2001fast,zou2009cubically}.

Recent work confirms that first-kind Fredholm integral equations remain an active area of numerical analysis, particularly when regularization is combined with collocation, projection, or wavelet-based discretization. Regularization--wavelet collocation methods have recently been used to stabilize ill-posed first-kind equations \cite{pooja2024waveletcollocation}, while regularized Legendre-collocation, weakly singular projection methods, and iterated Tikhonov approaches have also been proposed for this class of problems \cite{bilel2025regularizedlegendre,bechouat2024weaklysingular,pang2025iterated}. Comparative studies of wavelet and orthogonal-polynomial bases further show that the choice of approximation space can significantly affect accuracy and efficiency \cite{temirbekov2025comparative}. Modern kernel and Gaussian-process approaches also emphasize conditioning, minimum-norm recovery, and the role of sampling design in first-kind Fredholm problems \cite{qiu2025gaussian}. These developments motivate a systematic comparison in which several compactly supported bases are tested under the same regularized collocation framework.

All these recent studies have advanced several important directions for first-kind Fredholm integral equations, including regularized wavelet collocation, regularized projection or collocation schemes, and comparative studies of wavelet and orthogonal-polynomial bases. These contributions establish that both regularization and basis construction are important. The present work addresses a different and more diagnostic question: when the same ill-posed Fredholm problem is solved under one fixed Tikhonov--Morozov collocation framework, how does the choice of basis affect the behavior of the regularized inverse problem itself?

This distinction is important because, for first-kind equations, a basis is not only an approximation device. Once noisy data and regularization are introduced, the basis also influences the spectrum of the discrete operator, the conditioning of the normal equations, the residual level enforced by the discrepancy principle, the selected value of the regularization parameter, and the cost of the resulting computation. Thus two bases that appear comparable in a noise-free error table may behave quite differently when the data are perturbed and the inverse problem must be stabilized.

The main contribution of this paper is therefore a unified and reproducible diagnostic comparison of Linear Legendre multiwavelets, Legendre wavelets, Haar wavelets, B-spline wavelets, and hybrid wavelets for ill-posed Fredholm integral equations of the first kind. All bases are tested under the same discretization, regularization, parameter-selection, and error-measurement protocol. The comparison goes beyond approximation error by reporting the solution error $E_{\rm sol}$, residual error $E_{\rm res}$, Morozov-selected parameter $\alpha(\delta)$, conditioning of the discrete operator $\kappa(A_N)$, conditioning of the regularized normal equations $\kappa(A_N^T A_N+\alpha I)$, L-curve behavior, noise robustness, and computational cost. The novelty is not the isolated use of wavelets or Tikhonov regularization; rather, it is the demonstration that basis selection materially changes the stabilization behavior of noisy and ill-conditioned first-kind Fredholm discretizations. This provides practical guidance on which basis is preferable for smooth recovery, nonsmooth recovery, or low-cost computation.

To the best of the author's knowledge, no previous study has reported a systematic comparison of Linear Legendre multiwavelets, Legendre wavelets, Haar wavelets, B-spline wavelets, and hybrid wavelets under a common Tikhonov--Morozov regularization framework while simultaneously examining conditioning, discrepancy-selected regularization parameters, residual behavior, noise robustness, and computational cost. This unified diagnostic perspective constitutes the principal contribution of the present work.

The practical significance of this study is that inverse problems are rarely solved in the absence of measurement uncertainty. In such situations, the approximation basis and the regularization strategy cannot be viewed as independent design choices. The numerical results demonstrate that basis selection directly affects conditioning, regularization strength, noise propagation, and computational cost. Consequently, the choice of approximation space becomes part of the stabilization strategy itself. This observation provides a more informative criterion for basis selection than approximation accuracy alone.

This article is organized as follows. Section~\ref{sec:mra} introduces the multiwavelet systems used in the discretization. Section~\ref{sec:integration} describes the corresponding quadrature and numerical integration formulas. Section~\ref{sec:method} presents the regularized collocation formulation and the Tikhonov--Morozov parameter-selection procedure. Section~\ref{sec:examples} reports the numerical examples, including integration tests, Fredholm benchmark problems, and noise-sensitivity diagnostics. Section~\ref{sec:conclusion} summarizes the conclusions.

\section{Multiwavelet systems}\label{sec:mra}
For completeness, this section briefly summarizes the wavelet notation used in the collocation discretization. The purpose is not to reproduce the full theory of multiresolution analysis, but to fix the scaling, translation, support, and coefficient conventions needed later. Detailed theoretical developments can be found in the standard wavelet references cited below.

A wavelet basis is generated from a scaling function $\phi$ and a mother wavelet $\psi$ by dyadic dilation and translation. The scaling function satisfies a two-scale relation of the form
\[
\phi(x)=\sum_k p_k\phi(2x-k),
\]
where the sequence $\{p_k\}$ is finite for the compactly supported bases considered here. The associated finite-dimensional approximation spaces are generated by translated and dilated functions of the form
\[
\phi_{j,k}^{m}(x)=\phi^m(2^j x-k), \qquad
\psi_{j,k}^{m}(x)=\psi^m(2^j x-k),
\]
where $j$ denotes the resolution level, $k$ the translation index, and $m$ the local polynomial or multiwavelet component.

The construction is based on a nested sequence of approximation spaces
\[
\cdots \subset V_{-1}\subset V_0\subset V_1\subset \cdots,
\]
with complementary detail spaces $W_j$ satisfying $V_{j+1}=V_j\oplus W_j$. Thus functions are approximated by finite sums of scaling and wavelet components. In the present work, this general framework is specialized to Linear Legendre multiwavelets, Legendre wavelets, Haar wavelets, hybrid wavelets, and B-spline wavelets. These bases are then used to approximate the unknown function, the right-hand side, and the kernel in the regularized collocation discretization of the first-kind Fredholm equation.

\section{Numerical integration}\label{sec:integration}

\subsection{Linear Legendre Multi-Wavelets (LLMW)}
\label{subsec:llmw-integration}

In this subsection we consider the problem of integral approximation by the Linear Legendre Multi Wavelets (LLMW). These wavelets consist of a family of functions constructed from dilation and translation of a single function called the mother wavelet. Here we start by introducing 
the linear Legendre mother wavelet  functions $\psi^0 (x)$, and $\psi^1 (x)$,  which are given by : 

$$
\boldsymbol{ \psi^{0} (x)}= \left\{
    \begin{array}{ll}
        -\sqrt{3}(4x-1)  & \mbox{ [0,$\frac{1}{2}),$} \\
        \sqrt{3}(4x-3)  & \mbox{ [$\frac{1}{2},1),$} \\
    \end{array}
\right.
\boldsymbol{\psi^1 (x)}= \left\{
    \begin{array}{ll}
        6x-1  & \mbox{ [0,$\frac{1}{2}),$} \\
        6x-5  & \mbox{ [$\frac{1}{2},1),$} \\
    \end{array}
\right.
$$
and by performing translations and dilations of the mother functions, the linear Legendre multiwavelets are generated as follows :
$$
\boldsymbol{\psi^j_{kn} (x)}= \left\{
    \begin{array}{ll}
        2^{\frac{k}{2}}\psi^j(2^{k} x -n), & \mbox{  $  \frac{n}{ 2^{k}} \leq x  <   \frac{n+1}{ 2^{k}},$} \\
        0  & \mbox{ otherwise,} \\
    \end{array}
\right.
$$

where $n = 0,1,\ldots ,2^k-1$, $k \in Z^+$, and $j = 0,1$. The family $\{\boldsymbol{\psi^j_{kn} (x)}\} $ forms an orthonormal basis for $L^2({\rm I\!R})$. 
%These functions are defined on the interval $[0,1)$. 
We shall continue further the description of LLMW properties with introducing the scaling functions 
($\phi_0 (x)$,$\phi_1 (x)$) as follows:

$$
\phi_0 (x) = 1 , \phi_1 (x) =  \sqrt{3}(2x-1), \hspace{5 mm} 0 \leq x <1.
$$

A function $f(x)$ defined over the interval $[0, 1)$  can be approximated using the LLMW method as follows:

$$
 f_M(x) \approx c_0 \phi_0 (x)+c_1 \phi_1 (x) + \sum_{k=0 }^{M } \sum_{j=0 }^{1}\sum_{n=0 }^{2^k-1} c_{kn}^{j}\psi_{kn}^{j}(x)= C^{T}\Psi(x),
$$
where $M$ is the dilation for the LLMW and 
$c_{kn}^{j}=<f(x), \psi_{kn}^{j}(x)>$, in which $ <.,. >$ denotes the inner product. $C$ and $\psi(x)$ are given by:
$$
 C =[c_0,c_1,c_{00}^{0},...,c_{M0}^{0},c_{M1}^{0},...,c_{M(2^{M}-1)}^{0},...,c_{M0}^{1},c_{M1}^{1},...,c_{M(2^{M}-1)}^{1} ]^{T},
 $$

 $$
 \Psi =[\psi_0,\psi_1,\psi_{00}^{0},...,\psi_{M0}^{0},\psi_{M1}^{0},...,\psi_{M(2^{M}-1)}^{0},...,\psi_{M0}^{1},\psi_{M1}^{1},...,\psi_{M(2^{M}-1)}^{1} ]^{T}.
$$
 Similarly, using Linear Legendre Multiwavelet functions we can approximate the function $f(x,y)$ as:
$$
f(x,y)  \approx \Psi^T(x) K \Psi(y), 
$$
where
\begin{equation*}
K = 
\begin{pmatrix}
c_{0}c_{0} & c_{0}c_{1} & \cdots & c_{0}c_{M(2^M-1)}^{1} \\
c_{1}c_{0} & c_{1}c_{1} & \cdots & c_{1}c_{M(2^M-1)}^{1} \\
\vdots  & \vdots  & \ddots & \vdots  \\
c_{M(2^M-1)}^{1}c_{0} & c_{M(2^M-1)}^{1}c_{1} & \cdots & c_{M(2^M-1)}^{1}c_{M(2^M-1)}^{1}
\end{pmatrix}.
\end{equation*}

The Linear Legendre multiwavelet quadrature procedure is used below to approximate single, double, triple, and improper integrals, and later as part of the Fredholm collocation discretization.  

In one-dimensional space, the integral of a function $f(x) \in L^2({\rm I\!R})$ defined over an interval $[a,b]$ can be approximated using the Linear Legendre Multi-Wavelet method as follows : 
$$
\int_a^{b} f(x) dx \approx \frac{1}{ 2^{k}-1}\sum_{i=0 }^{2^{k+2} -1}f(x_i) = \frac{(b-a)}{ 2^{k+2}}\sum_{i=0 }^{2^{k+2} -1}f(a+(b-a)\frac{i+0.5}{ 2^{k+2}})
$$

This formula can be extended to double integral for a function $f(x,y) $  defined over an interval $[a,b] \times [c,d]$, as follows: 
$$
\int_c^{d}\int_a^{b} f(x,y) dx dy \approx \frac{(d-c)(b-a)}{ 2^{2k+4}}\sum_{i=0 }^{2^{k+2} -1}\sum_{j=0 }^{2^{k+2} -1}f(x_i,y_j),
$$
where, 

$$
x_i=a+(b-a)\frac{i+0.5}{ 2^{k+2}},   i = 0,1..2^{k+2}-1
$$
$$
y_j=c+(d-c)\frac{j+0.5}{ 2^{k+2}}, j = 0,1..2^{k+2}-1.
$$

Similarly, the triple integral of a function $f(x,y,z) $, defined over  $[a,b]\times[c,d]\times[e,f]$, can be approximated as:
$$
\int_e^{f}\int_c^{d}\int_a^{b} f(x,y,z) dx dy dz \approx \frac{(f-e)(d-c)(b-a)}{ 2^{3k+4}}\sum_{i=0 }^{2^{k+2} -1}\sum_{j=0 }^{2^{k+2} } 
\sum_{l=0 }^{2^{k+2} -1} f(x_i,y_j,z_l),
$$
where, 

$$
x_i=a+(b-a)\frac{i+0.5}{ 2^{k+2}}, \hspace{4 mm}  i = 0,1..2^{k+2}-1,
$$
$$
y_j=c+(d-c)\frac{j+0.5}{ 2^{k+2}},\hspace{4 mm} j = 0,1..2^{k+2}-1,
$$
$$
z_l=e+(f-e)\frac{l+0.5}{ 2^{k+2}},\hspace{4 mm} l = 0,1..2^{k+2}-1.
$$
\subsection{Legendre Wavelets }
Legendre wavelets constitute a family of functions constructed from dilation and translation of a single function called mother wavelet. For any positive integer $k$, the one-dimensional Legendre wavelets are defined on the interval $[0,1)$ as follows:  

\begin{equation}
\begin{aligned}
\boldsymbol{ \psi_{nm} (t)}= \left\{
    \begin{array}{ll}
        \sqrt{m + \frac{1}{2}} 2^{\frac{k}{2}} L_m(2^k t-2n+1), & t \in [ \frac{n-1}{2^{k-1}}, \frac{n}{2^{k-1}} ) , \\
        0,  & \mbox{otherwise},\\
    \end{array}
\right.
\end{aligned}
\end{equation}


where $n = 1,2.. ,2^{k-1}$  and $m = 0,1..M-1$. Here, $L_m(t)$ are the Legendre polynomials of order $m$, which can be defined  as follows:

\begin{equation}
\begin{aligned}
L_0(t)= 1, \, L_1(t)= t,
\end{aligned}
\end{equation}

\begin{equation}
\begin{aligned}
L_{m+1}(t)= \frac{2m+1}{m+1}tL_{m}(t)- \frac{m}{m+1}L_{m-1}(t),\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, m=1,2,3...
\end{aligned}
\end{equation}

The set of $\{ L_m \} $ in the Hilbert space $L^2[-1,1]$ is a complete orthogonal set. Orthogonality of Legendre Polynomials on the interval $[-1,1]$ implies that

\begin{equation}
\begin{aligned}
<L_{m}(t),L_p(t)>=\int_{-1}^{1}L_{m}(t)L_p(t)dt= \left\{
    \begin{array}{ll}
       \frac{2}{2m+1}, & m = p, \\
        0,  &  m \neq p ,\\
    \end{array}
\right.
\end{aligned}
\end{equation}

for $m,p=0,1,2,....$.

Furthermore, the set of wavelets  $\psi_{nm}$ makes an orthonormal basis in $L^2[0, 1]$, that is

\begin{equation}
\int_{0}^{1}\psi_{nm}(t)\psi_{kl}dt= \delta_{nk} \delta_{ml},
\end{equation}
 in which $\delta$ denotes the Kronecker delta function.

A function $f(t)$ in $L_2[0; 1)$ can be approximated in terms of Legendre wavelets as follows:

\begin{equation}
\begin{aligned}
f(t) \approx \sum_{n=1 }^{2^{k-1}}\sum_{m=0}^{M-1} c_{nm}\psi_{nm}(t) = C^{T}\Psi(t)
\end{aligned}
\end{equation}
 in which $ C$ and $\Psi(t)$ are $2^{k-1}M$ of the form
\begin{equation}
\begin{aligned}
 C =[c_{10},c_{11},...,c_{1(M-1)},c_{20},..,c_{2(M-1)},.....,c_{2^{k-1}0},..,c_{2^{k-1}(M-1)} ]^{T},
 \end{aligned}
\end{equation}
and
\begin{equation}
\begin{aligned}
 \Psi(t) =[\psi_{10}(t),...,\psi_{1(M-1)}(t),\psi_{20}(t),...,\psi_{2(M-1)}(t),...,\psi_{2^{k-1}0}(t),...,\psi_{2^{k-1}(M-1)}(t)]^{T}.
 \end{aligned}
\end{equation}

The wavelet coefficients, $c_{nm}$, are determined by:
\begin{equation}
\begin{aligned}
c_{nm}=<f(t),\Psi_{mn}(t)>, \hspace{5mm} \text{for} \hspace{2mm} n=1,2,..,2^{k-1}, \hspace{3 mm} m=0,1,2, M-1.
\end{aligned}
\end{equation}

Similarly, the function $\mathcal{K}(x,y) \in L_2([0, 1) \times [0; 1) )$  can be approximated as:

\begin{equation}
\begin{aligned}
\mathcal{K}(x,y)  \approx \sum_{i=1 }^{2^{k-1}}\sum_{j=0}^{M-1}\sum_{k=1 }^{2^{k-1}}\sum_{l=0}^{M-1}\psi_{ij}(x)c_{ijkl}\psi_{kl}(y) = \Psi^{T}(x)K\Psi(y), 
\end{aligned}
\end{equation}

where $K$ is a $2^{k-1}M \times 2^{k-1}M$ matrix which satisfies:
\begin{equation}
K_{ij} =<\psi_{i}(x),<\mathcal{K}(x,y),\psi_{j}(y)>>, \hspace{3 mm} i,j=1,2,..., 2^{k-1}M.
\end{equation}
\subsection{Haar-Wavelets }
\label{subsec:legendre-wavelets}

Used in various numerical approximations problems such as integral equations, partial differential equations and ordinary differential equations, Haar wavelet method is a powerful numerical tool,  characterized by a simple orthonormal wavelets basis and a compact support. The scaling function $h_1(x)$ for the Haar wavelets defined over a sub-interval $[a,b)$ is given by: 

$$
\boldsymbol{h_1(x)} = \left\{
    \begin{array}{ll}
        1  & \mbox{for  x $\in [a,b)$ } \\
        \ 0     & \mbox{     elsewhere}
    \end{array}
\right.
$$

The mother wavelet function $h_2(x)$, which is also defined over an interval $[a,b)$ is denoted by:

$$
\boldsymbol{h_2(x)}= \left\{
    \begin{array}{ll}
        1  & \mbox{ for x $\in [a,\frac{a+b}{2})$ } \\
        \ -1     & \mbox{ for  x $\in [\frac{a+b}{2},b)$ }\\
        
        \ 0     & \mbox{     elsewhere.}
    \end{array}
\right.
$$

The Haar wavelets family $\{h_i(x)\}$, for $ i=3,4,..,2M$, is generated by translations and dilations of the mother function as follows:
$$
\text{For} \hspace{3mm}  i=3,4,..,2M \hspace{3mm} \text{we have}\\
\hspace{0.4cm}
\boldsymbol{h_i(x)} = \left\{
    \begin{array}{ll}
        1  & \mbox{ x $\in [\alpha,\beta)$ } \\
        \ -1     & \mbox{   x $\in [\beta,\gamma)$ } \\
        
        \ 0     & \mbox{     elsewhere,}
    \end{array}
\right.
$$
where $\alpha = a +k(b-a)/m$, $\beta = a +(b-a)(2k+1)/2m$ and $ \gamma = a +(k+1)(b-a)/m $.

The integer $m = 2^j$, where $j = 0,1.. 2^M$ and integer $k = 0,1..,m-1$. The integer $j$ indicates the level of the wavelet and $k$ is the translation parameter. The maximum level of the resolution is the integer $J$. The relation between $i$, $m$ and $k$ is given by $i=m+k+1$.

Any function f(x) which is square integrable in the interval $[a,b)$ can be represented as an infinite sum of Haar wavelets, as follows:

$$
f(x)=\sum_{i=1 }^{\infty}\alpha_i h_i(x)=C^T H(x),
$$
 where $H$ represents the Haar-wavelet family and $C$ is the Haar-wavelet coefficient vector, with components $\alpha_i$ given by:
 $$\alpha_i=2^j<h_i(x),f(x)>,$$
 where $<.,.>$ denotes the inner product and $j$ denotes the corresponding level of resolution.
 
 Similarly, using Haar wavelet functions we can approximate the function $f(x,y)$ as:
$$
f(x,y)  \approx \sum_{i=1 }^{2M}\sum_{k=0}^{2M}  h_{i}(x)c_{ij} h_{k}(y)=H^T(x) C H(y). 
$$
  For $i=1$ and $k=2,.., 2M$, the wavelet coefficient $c_{1k}$ are computed as follows:
\begin{equation}
\begin{aligned}
c_{1k}= 2^{j1}<h_i(x),<f(x,y),h_k(y)>>. 
\end{aligned}
\end{equation}
For $k=1$ and $i=2,.., 2M$, the wavelet coefficient $c_{i1}$ are computed as follows:
\begin{equation}
\begin{aligned}
c_{i1}= 2^{j2}<h_i(x),<f(x,y),h_k(y)>>.
\end{aligned}
\end{equation}
For $ i=2,.., 2M$ and $k=2,.., 2M$, the wavelet coefficients $c_{ik}$ are computed as follows:
\begin{equation}
\begin{aligned}
c_{ik}= 2^{j1+j2}<h_i(x),<f(x,y),h_k(y)>>. 
\end{aligned}
\end{equation}

The same Haar-wavelet quadrature framework is then used to approximate single, double, and triple improper integrals. Following the numerical-integration procedure based on Haar and Linear Legendre multiwavelets, we obtain the following quadrature formula for a function $f(x)$ defined on the interval $[a,b]$:
$$
\int_a^{b} f(x) dx \approx \frac{b-a}{ 2M}\sum_{k=1 }^{2M} f(x_k) = \frac{(b-a)}{ 2M}\sum_{k=1 }^{2M }f(a+\frac{(b-a)(k-0.5)}{ 2M}).
$$

This formula can be extended to approximate a double integral for a function $f(x,y)$  defined over $[a,b] \times [c,d]$, as follows: 

$$
\int_c^{d} \int_a^{b} f(x,y) dx dy\approx \frac{(b-a)(d-c)}{ 2M^2}\sum_{l=1 }^{2M}\sum_{k=1 }^{2M} f(a+(b-a)\frac{(k-0.5)}{ 2M},c+(d-c)\frac{(l-0.5)}{ 2M}).
$$

Similarly, this formula can be extended to approximate a triple integral for a function $f(x,y,z)$  defined over  $[a,b] \times [c,d] \times [e,f]$, as follows: 

$$
\int_e^{f} \int_c^{d} \int_a^{b} f(x,y,z) dx dy dz \approx \frac{(b-a)(d-c)(f-e)}{ 8M^3}\sum_{i=1 }^{2M}\sum_{j=1 }^{2M}\sum_{l=1 }^{2M} f(x_i,y_j,z_l).
$$


\subsection{Hybrid-Wavelets  }
\label{subsec:haar-wavelets}

The orthogonal set of Hybrid wavelet functions $\psi_{ij}(x)$, $i = 1,..,n$ and $ j = 0,..,m-1$ is defined over the interval $[0,1)$ as:

$$
\boldsymbol{\psi_{ij} (x)} = \left\{
    \begin{array}{ll}
        L_j(2nx-2i +1), & \mbox{ for $ x \in [\frac{i-1}{n},\frac{i}{n})$} \\
       0,  & \mbox{ otherwise,} \\
    \end{array}
\right.    
$$
where $n$ is the order of the block pulse function and $m$ is the order of the Legendre polynomials. The Legendre polynomials, $L_{k}(x)$, of order $k$ can be calculated recursively as:

$$
L_0(x)= 1, \, L_1(x)= x,
$$

$$
L_{k+1}(x)= \frac{2k+1}{k+1}xL_{k}(x)- \frac{k}{k+1}L_{k-1}(x)\,\, \hspace{8 mm}  k=1,2,3...
$$

Any function $ f(x)$ which is square integrable in the interval $[0,1) $ can be represented as:
$$
f(x)= \sum_{i=1 }^{\infty}\sum_{j=0}^{\infty} c_{ij}\psi_{ij} \,\, \hspace{6mm} i=0,1,..\infty \,\, ,\,\,  j=0,1,..\infty \,\, \hspace{6mm}  x \in [0,1).
$$

If the function $ f(x)$ is piecewise constant or may be approximated by a piecewise constant function, $f(x)$ can be approximated as follows:

$$
f(x)  \approx \sum_{i=1 }^{n}\sum_{j=0}^{m-1} c_{ij}\psi_{ij}=C^T \Psi, 
$$
where $\Psi$ represents the wavelet family and $C$ denotes the wavelet coefficient vector defined by:

\begin{equation}
\begin{aligned}
c_{ij}=\frac{<\Psi_{ij},f(x)>}{<\Psi_{ij},\Psi_{ij}>}.
\end{aligned}
\end{equation}
 A single integration of the function $f(x)$ over the interval $[0,1]$ leads to the following formula:

$$
\int_0^{1} f(x) dx \approx \int_0^{1}\sum_{i=1 }^{n}\sum_{j=0}^{m-1} c_{ij}\psi_{ij} .
$$

In order to compute the coefficients $c_{ij}$ of the hybrid functions we consider the nodal points:
$$
x_k=\frac{2k-1}{2mn},\,\, \hspace{6mm} k=1,2,..mn, 
$$

where $f(x_k) $ can be written as:
$$
f(x_k)=\sum_{i=1 }^{n}\sum_{j=0}^{m-1} c_{ij}\psi_{ij}(x_k),\,\, \hspace{6mm} k=1,2,..mn. 
$$

Hence, the approximate value of the single integral, for different values of the of the Legendre polynomials $m$, are given by:
$$$$
\hspace{5 mm} For $m = 1$,
$$
\int_0^{1} f(x) dx  \approx  \frac{1}{n}\sum_{i=1}^{n} f(\frac{2i-1}{2n}).
$$

For $m = 2$,
$$
\int_0^{1} f(x) dx  \approx  \frac{1}{2n}\sum_{i=1}^{2n} f(\frac{2i-1}{4n}).
$$

For $m = 3$,
$$
\int_0^{1} f(x) dx  \approx  \frac{1}{8n}\sum_{i=1}^{n} 3f(\frac{6i-5}{6n})+2f(\frac{6i-3}{6n})+3f(\frac{6i-1}{6n}).
$$

For $m = 4$,
$$
\int_0^{1} f(x) dx  \approx  \frac{1}{48n}\sum_{i=1}^{n} 13f(\frac{8i-7}{8n})+11f(\frac{8i-5}{8n})+11f(\frac{8i-3}{8n})+13f(\frac{8i-1}{8n}).
$$

For $m = 5$,
$$
\int_0^{1} f(x) dx  \approx  \frac{1}{1152n}\sum_{i=1}^{n} 275f(\frac{10i-9}{10n})+100f(\frac{10i-7}{10n})+ 
$$ 
$$
402f(\frac{10i-5}{10n})+100f(\frac{10i-3}{10n})+275f(\frac{10i-1}{10n}).
$$

For $m = 6$,
$$
\int_0^{1} f(x) dx  \approx  \frac{1}{1280n}\sum_{i=1}^{n} 247f(\frac{12i-11}{12n})+139f(\frac{12i-9}{12n})+
$$ 
$$
254f(\frac{12i-7}{12n})+254f(\frac{12i-5}{12n})+139f(\frac{12i-3}{12n})+247f(\frac{12i-1}{10n}).
$$



Similarly, using hybrid functions we can approximate the function $f(x,y)$ as:
$$
f(x,y)  \approx \sum_{i=1 }^{n}\sum_{j=0}^{m-1} \sum_{k=1 }^{n}\sum_{l=0}^{m-1} c_{ijkl}\psi_{ij}(x) \psi_{kl}(y)=\Psi^T C \Psi, 
$$
where the wavelet coefficients $c_{ijkl}$ are given by:
\begin{equation}
\begin{aligned}
c_{ijkl}=\frac{<\Psi_{ij},<f(x,y),\Psi_{kl}>>}{<\Psi_{ij},\Psi_{ij}><\Psi_{kl},\Psi_{kl}>}
\end{aligned}
\end{equation}
The double integral approximations are given below for the first few values of $m$.  
$$$$
\hspace{5 mm} For $m = 1$,
$$
\int_0^{1} \int_0^{1} f(x,y) dx dy  \approx  \frac{1}{n^2}\sum_{i=1}^{n}\sum_{k=1}^{n} f(\frac{2i-1}{2n},\frac{2k-1}{2n}).
$$

For $m = 2$,
$$
\int_0^{1} \int_0^{1} f(x,y) dx dy \approx  \frac{1}{4n^2}\sum_{i=1}^{2n}\sum_{k=1}^{2n} f(\frac{2i-1}{4n},\frac{2k-1}{4n}).
$$

For $m = 3$,
$$
\int_0^{1} \int_0^{1} f(x,y) dx dy   \approx  \frac{1}{64n^2}\sum_{i=1}^{n}\sum_{k=1}^{n} 9f(\frac{6i-5}{6n},\frac{6k-5}{6n})+
6f(\frac{6i-5}{6n},\frac{6i-3}{6n})+9f(\frac{6i-5}{6n},\frac{6i-1}{6n})
$$
$$
6f(\frac{6i-3}{6n},\frac{6k-5}{6n})+
4f(\frac{6i-3}{6n},\frac{6i-3}{6n})+6f(\frac{6i-3}{6n},\frac{6i-1}{6n})
$$
$$
9f(\frac{6i-1}{6n},\frac{6k-5}{6n})+
6f(\frac{6i-1}{6n},\frac{6i-3}{6n})+9f(\frac{6i-1}{6n},\frac{6i-1}{6n}).
$$
\begin{comment}

\section{B-spline Wavelets}
\label{appendix-sec1}
Recently used in solving integral and integro-differential equations [6], B-spline Wavelets has a  semi-orthogonal basis of scaling functions with vanishing moments, characterized by a finite support with even and odd symmetry and generating easily implemented analytical expressions. Some of the important properties relevant to B-spline Wavelets are given in this subsection. We shall start from introducing the cardinal B-spline functions of order m $\leq$ the k not sequence X:

$$
B_{m,j,X}(x) = \frac{x-x_j}{x_{j+m-1}-x_j}B_{m-1,j,X}(x) + 
$$
$$
\frac{x_{j+m}-x}{x_{j+m}-x_j}B_{m-1,j+1,X}(x)
$$
$$
\boldsymbol{B_{1,j,X}(x)} = \left\{
    \begin{array}{ll}
       1 & \mbox{ for $ x \in [x_j,x_{j+1})$} \\
       0  & \mbox{ else} \\
    \end{array}
\right.    
$$
Where the knot sequence X is denoted by:

$$
X=\{x_k^j\}^{2^{j}+m-1}_{k=1-m}
$$

We can consider for example the knot sequence defined over the interval [0,1], at any level
j ∈ $Z^+$, the discrete step is $2^{−j}$, and this generates n = $2^{j}$ segments in [0,1]: 

$$
\boldsymbol{X^j} = \left\{
    \begin{array}{ll}
        &x_{1-m}^{j}=x_{2-m}^{j}=...=x_{0}^{j}=0 \mbox{ }\\
        &x_{k}^{j}=k 2^{-j}, \hspace{3 mm} \mbox{$ k=1,..,2^{-j}-1$}\\
        &x_{2^{-j}}^{j}=x_{2^{-j+1}}^{j}=...=x_{2^{-j+m-1}}^{j}=0 \mbox{ }\\
    \end{array}
\right.    
$$

m denotes here the order of the B-Spline wavelets, which are constructed starting from the lowest octave level $j_0$ evaluated from:

$$2^{j_0}\geq 2m-1$$
For each level $j \geq j_0 $, the wavelet function of order $m$ can be determined as follows:
$$
\boldsymbol{ \psi_{m,j,k}} = \left\{
    \begin{array}{ll}
        &B_{m,j_0,k}(2^{j-j_0}x) \mbox{\, $ k=1-m,..,-1$}\\
        &B_{m,j_0,2^{j}-m-k}(1-2^{j-j_0}x) \mbox{\, $k = 2^{j}-m+1,..,2^{j}-1$}\\
        &B_{m,j_0,0}(2^{j-j_0}x-2^{-j_0}k) \mbox{\, $k = 0,..,2^{j}-m$}\\
    \end{array}
\right.    
$$
The two-scale relation for the m-order semi-orthogonal B-spline wavelet are denoted by:
$$
\psi_{m,j,k-m}(x) = \sum_{j=k}^{2k+2m-1} q_{k,j} B_{m,j,k-1}(x) \,\, k=1-m,..,-1
$$

$$
\psi_{m,j,k-m}(x) = \sum_{j=2k-m}^{2k+2m-1} q_{k,j} B_{m,j,k-1}(x) \,\, k=m,..,n-m+1
$$

$$
\psi_{m,j,k-m}(x) = \sum_{j=k}^{n+km-1} q_{k,j} B_{m,j,k-1}(x) \,\, k=n-m+2,..,n
$$
Where, $$q_{k,j}=q_{j-2k}$$
So, for each level $j \geq j_0 $, the B-Spline wavelet of order m is given by:

$$
\boldsymbol{\Psi_{m,j,k}} = \left\{
    \begin{array}{ll}
        & \psi_{m,j_0,k}(2^{j-j_0}x) \mbox{\, $ k=-m+1,..,-1$}\\
        &\psi_{m,j_0,2^{j}-m-k}(1-2^{j-j_0}x) \mbox{\, $k = 2^{j}-2m+2,..,2^{j}-m$}\\
        &\psi_{m,j_0,0}(2^{j-j_0}x-2^{-j_0}k) \mbox{\, $k = 0,..,2^{j}-2+1m$}\\
    \end{array}
\right.    
$$
Recursively, The B-Spline wavelet can be computed applying the convolution:
$$
\rho_m(x)=\int_{-\infty}^{+\infty}\rho_{m-1}(x-t)\rho_{1}(t)dt=\int_{0}^{1}\rho_{m-1}(x-t)
$$
Where,

$$
\boldsymbol{\rho_{1}(x)} = \left\{
    \begin{array}{ll}
        & 0 \mbox{\, $0\leq x <1$}\\
        &1 \mbox{\, else}\\
    \end{array}
\right.    
$$

The construction of the scaling function of m-th order B-spline wavelet is based on the two scale relation:
$$
\rho(x) = \sum_{j=0}^{m} p_{k} \rho_{m}(2x-k) \,\, k=n-m+2,..,n.
$$
Where, $p_k$ denotes the two scale sequence :
$$ p_k = 2^{1-m}(\begin{array}{ll}
        & m\\
        & k\\
    \end{array}),\,\, 0 \leq k < m$$

The two scale relation for the m-th order B-spline is denoted by: 
$$
\psi_{m} = \sum_{k=0}^{2m-2} q_k\rho_{m}(2x-k)  
$$
Where,
Where, $p_k$ denotes the two scale sequence:
$$ p_k = (-1)^k 2^{1-m} \sum_{l=0}^{m} (\begin{array}{ll}
        & m\\
        & l\\
    \end{array}) \rho_{2m}(k-l+m)   $$
    
A function f(x), which is defined over [0,1], can be approximated by B-spline scaling functions as with admitting (M = $j_0$):
$$
f(x)= \sum_{k=-4}^{2^{m}-1}s_k\rho_{M,k}= S^{T}\phi_m
$$
Where,
$$ S=[s_{-4},s_{-3},..,s_{2^{m-1}}]$$
$$ \Phi=[\phi_{-4},\phi_{-3},..,\phi_{2^{m-1}}]$$
$$
s_k=\int_{0}^{1}f(x)\tilde{\rho}_{m,k}(x)dx \, \, k=-4,-3,..,2^{m}-1
$$
The last equation that denotes the approximation of a function f(x) defined over [0,1] can be expressed in a simpler manner as follows:

$$
f(x)= \sum_{k=-4}^{2^{j}-1}c_{m,k}\rho_{5,k}(x)+\sum_{j=5}^{\infty}\sum_{k=-4}^{2^{j}-5}d_{j,k}\psi_{j,k}(x)
$$

If the infinite series is truncated at M, then the last equation becomes simpler:

$$
f(x) \approx \sum_{k=-4}^{15}c_{m,k}\rho_{5,k}(x)+\sum_{j=5}^{M}\sum_{k=-4}^{2^{j}-5}d_{j,k}\psi_{j,k}(x)=C^{T}\Psi(x)
$$

Where,

$$ C=[c_{-4},c_{-3},..,c_{15},d_{5,-4},d_{5,-3},..,d_{5,7},..,d_{M,-M+1},..,d_{2^{M},-5}]
$$

$$
\Psi=[\rho_{5,-4},\rho_{5,-3},..,\rho_{5,15},\phi_{5,-4},\phi_{5,-3},..,\phi_{5,7},..,\phi_{M,-M+1},..,\phi_{2^{M},-5}]
$$

With
$$
c_k=\int_{0}^{1}f(x)\tilde{\rho}_{5,k}(x)dx \, \, k=-4,-3,..,15
$$
$$
d_{jk}=\int_{0}^{1}f(x)\tilde{\psi}_{j,k}(x)dx \, \, \, \,k=5,-4,..M, \, \, k=-4,-3,..,2^{j}-5
$$

This is my comment.
Note that it can span multiple lines.
This is very useful.
\end{comment}

\section{Numerical method}\label{sec:method}
In this section we suggest an efficient method for solving $1D$ and $2D$ Fredholm integral equations of the first kind. The wavelets in the previous section are first used as a basis to reduce the solution of the integral equations to a system of algebraic equations. Then an algorithm based on the regularization to handle the ill-posed Fredholm equations will be introduced.

\subsection{Solution of Fredholm integral equations of the first kind}
We consider the $1D$ Fredholm integral equation:
\begin{equation} \label{eqe1}
\int_{0}^{1} \mathcal{K}(x,y)f(y)dy=g(x). 
\end{equation}

Using the wavelets described in the previous section, we approximate $f(x)$, $g(x)$, and the kernel $\mathcal{K}(x,y)$ as follows:

\begin{equation} \label{eqd1}
f(x)  \approx C^T \Psi(x), 
\end{equation}
\begin{equation} \label{eqd2}
g(x)  \approx G^T \Psi(x), 
\end{equation}
\begin{equation} \label{eqd3}
\mathcal{K}(x,y)  \approx \Psi^T (x) K  \Psi(y), 
\end{equation}

where $C$ is the unknown vector associated to the unknown function $f(x)$, $G$ is a known vector associated with the data function $g(x)$ and $K$ is the known wavelet coefficient matrix associated with the kernel $\mathcal{K}(x,y)$.

By substituting Eqs. (\ref{eqd1})-(\ref{eqd3}) into Eq. (\ref{eq1}), we get:
\begin{equation}
 \int_{0}^1 \Psi^T(x)K\Psi(y)\Psi^T(y)C dy = \Psi^T(x)G,
\end{equation}
or
\begin{equation}\label{eqd4}
\begin{aligned}
 \Psi^T(x)K \left(\int_{0}^1 \Psi(y)\Psi^T(y) dy \right ) C= \Psi^T(x)G.
\end{aligned}
\end{equation}

Let 
\begin{equation}\label{eqd5}
P=\int_{0}^1 \Psi(y)\Psi^T(y) dy, 
\end{equation}
so $P$ is a matrix. By substituting Eq.(\ref{eqd5}) into Eq.(\ref{eqd4}), we get:

\begin{equation}\label{eqd6}
 \Psi^T(x)K P C= \Psi^T(x)G.
\end{equation}

Therefore, we have the following linear system of equations:
\begin{equation}\label{eqd7}
K_P C= G, \hspace{1cm} with \hspace{1cm} K_P=KP
\end{equation}
which can be solved for the unknown vector $C$.

Note that for the Haar wavelets and the Legendre wavelets, we have $P=I$ and Eq.(\ref{eqd7}) is reduced to:
\begin{equation}\label{eqd8}
K C= G.
\end{equation}

For the hybrid wavelets, the integration of the cross product of two wavelet vectors shows that the matrix $P$ can be written as an $m \times m$ diagonal matrix $D$ in the form: 
\begin{equation*}
D = \frac{1}{n}
\begin{pmatrix}
1 & 0 & \cdots & 0 \\
0 & \frac{1}{3} & \cdots & 0 \\
\vdots  & \vdots  & \ddots & \vdots  \\
0 & 0 & \cdots & \frac{1}{2m-1}
\end{pmatrix}.
\end{equation*}
Therefore, the 1D- Fredholm integral equation of the first kind is reduced to the following linear system of equations:
\begin{equation}\label{eqld7}
K_D C= G, \hspace{5 mm} \text{with} \hspace{5 mm} K_D=KD, 
\end{equation}
which can be solved for the unknown vector $C$.

 Similarly, if we consider the $2D$ Fredholm integral equation:

\begin{equation} \label{eqe2}
\int_{0}^{1} \int_{0}^{1} k(s,t,x,y)f(x,y)dy=g(x,y), 
\end{equation}

and we follow the same steps, we will approximate the functions $f(x,y)$, $g(x,y)$ and
$k(s,t,x,y)$ with respect to the $2D$ wavelets as:
\begin{equation} \label{eqdd1}
f(x,y)  \approx  \Psi^T(x,y).C, 
\end{equation}
\begin{equation} \label{eqdd2}
g(x,y)  \approx \Psi^T(x,y).G, 
\end{equation}
\begin{equation} \label{eqdd3}
k(s,t,x,y)  \approx \Psi^T (s,t) K  \Psi(x,y). 
\end{equation}
Therefore, we obtain a similar linear system of equations in the form:
\begin{equation}\label{eqdd7}
K_{\tilde{D}} C= G, \hspace{3mm} \text{with} \hspace{3mm} K_{\tilde{D}}=KP
\end{equation}
which can be solved for the unknown vector $C$, with $P$ is given by:
\begin{equation}\label{eqdd5}
P=\int_{0}^1 \int_{0}^1 \Psi(x,y)\Psi^T(x,y) dydx. 
\end{equation}
Note that for the Haar wavelets and the Legendre wavelets, we have $P=I$. For the hybrid wavelets, the integration of the cross product of two wavelet vectors shows that the matrix $P=\tilde{D}=D^2$.



\subsection{ The Tikhonov regularization method}

Most of the standard numerical procedures used to solve such integral equation or the equivalent linear system given by Eq.(\ref{eqdd7}) are unstable and pose serious problems since the first-kind Fredholm problem is ill-posed, i.e. when the data $G$ are perturbed by a noise $\delta$,  small changes in the observation data can cause arbitrarily large changes in the results. Thus, some kind of regularization methods such as Tikhonov method  are needed to obtain stable solutions for ill-posed problems. These methods depend strongly on the
choice of the regularization parameter. Hence, it is necessary to consider appropriate strategies for selecting the regularization parameters.

In this subsection, we introduce a special regularization method which  is applied  to handle the ill-posed Fredholm problems. 
This regularization method will produce stable approximations even if the solution $f$ does not depend continuously on the data $g$, which means
 that the inverse of the operator is unbounded, $||K^{-1}||=\infty$. 

The idea of the regularization method is to transform an ill-posed Fredholm integral equation of the first kind
\begin{equation} \label{eqq2}
    \int_{0}^{1} \mathcal{K}(x,y)f(y)dy= g(x), 
\end{equation}
to a nearly well-posed Fredholm integral equation of the second kind
\begin{equation} \label{eq:tikhonov-second-kind}
  f_{\alpha}(x)=\frac{g(x)}{\alpha} - \frac{1}{\alpha} \int_{0}^{1} \mathcal{K}(x,y) f_{\alpha}(y)dy, 
\end{equation}
where $\alpha$ is a small positive parameter called the regularization parameter. Moreover, it was proved by  Tikhonov et al. \cite{tikhonov1963solution, tikhonov1987ill, tikhonov2013numerical, tikhonov1963regularization}, Philips \cite{phillips1962technique} and \cite{delves1974numerical}, Groetsch \cite{groetsch1984theory} that the solution $f_{\alpha}(x)$ of Eq. (\ref{eq:tikhonov-second-kind}) converges to the solution $f(x)$ of Eq. (\ref{eqq2}). In other words, it is shown that 
\begin{equation} \label{eq:tikhonov-convergence}
 f(x)=\lim\limits_{\alpha 	\rightarrow 0^+} f_{\alpha}(x), 
\end{equation}

By converting the first kind to a second kind, we can use the existing techniques of the second kind to the transformed equation. 

 We want to solve Eq. (\ref{eqq2}) with noisy data $g^{\delta}$. A common way
to obtain a solution is to determine the best fit in the sense that we minimize the functional
\begin{equation} \label{eq4}
J(f)= ||Kf-g^{\delta}||^2, 
\end{equation}
where the operator $K$ is given by $Kf= \int_{0}^{1} \mathcal{K}(x,y)f(y)dy$. But again, this solution does not depend continuously on the data. To ensure stability, a penalty term is added to this functional. Then we get the Tikhonov functional:
\begin{equation} \label{eq5}
J_{\alpha}(f)= ||Kf-g^{\delta}||^2 + \alpha ||f||^2 , 
\end{equation}
 with a regularization parameter $\alpha$.

The Tikhonov functional $J_{\alpha}(f)$ has a unique minimum  $f_{\alpha}^{ \delta}$ for all $g^{\delta}$ and $\alpha > 0$. This minimum is also the unique solution of the normal equation

\begin{equation} \label{eq6K}
(K^{*}K+\alpha I)f_{\alpha}^{ \delta} =K^{*} g^{\delta}, 
\end{equation}

where $K^{*}$ denotes the adjoint operator of $K$ and the operator $K^{*}K+\alpha I$ is bounded and invertible, i.e. $f_{\alpha}^{ \delta}$ depends continuously on $g^{\delta}$. 
Note that $f_{\alpha}^{ \delta}$ defined by $f_{\alpha}^{ \delta}=(K^{*}K+\alpha I)^{-1}K^{*} g^{\delta}$ minimizes the Tikhonov functional $J_{\alpha}(f)$.

In order to solve the equivalent linear system given by the matrix form Eq.(\ref{eqd7}) for the unknown vector $C$, we consider solving an equivalent linear system to Eq. (\ref{eq6K}), given by

\begin{equation} \label{eqsys}
(K_P^{T}K_P+\alpha I)C_{\alpha}^{ \delta} =K_P^{T} G^{\delta}, 
\end{equation}
where $G^{\delta}=G+\delta rand(size(G))$, and $C_{\alpha}^{ \delta}$ is the regularized form of $C$. Finally, the unknown function $f$ will be approximated as
\begin{equation} \label{eqds1}
f(x)  \approx f_{\alpha}^{ \delta}= C_{\alpha}^{ \delta} {}^{T} \Psi(x), 
\end{equation}

In general, if the regularization parameter $\alpha$ is chosen to be large, then the error in $g^{\delta}$ will be damped, and if  $\alpha$ is chosen to be small, then 
$J_{\alpha}$ will be a good approximation of $J$ but the error will affect the solution. The well-known Morozov discrepancy principle has been used for linear ill-posed problems to choose the good regularization parameter. The principle is devoted to a posteriori choice of the regularization parameter $\alpha(\delta)>0$ . The idea of the strategy is to
choose the regularization parameter so that the residual $||K f_{\alpha}^{ \delta}-g^{\delta}|| $ or $||K_P C_{\alpha}^{ \delta}-G^{\delta}|| $ has the same error level as the observation data.
That is, we require that $\alpha(\delta)>0$ be determined by the equation
\begin{equation} \label{eq6}
||K f_{\alpha}^{ \delta}-g^{\delta}||  ={\delta}. 
\end{equation}

In some applications, the Morozov principle may not be so satisfactory and the optimal convergence
of regularized solutions is not obtained \cite{engl1996regularizatione, morozov2012methods}. We therefore consider a more general class of
the damped Morozov discrepancy principle \cite{ito1992choice,engl2000new,kunisch1993class} given by
\begin{equation} \label{eq7}
||K f_{\alpha}^{ \delta}-g^{\delta}|| + \alpha^{\gamma} ||f_{\alpha}^{ \delta}||=\delta^2, 
\end{equation}
where $\gamma \in [1,\infty]$. Obviously, the exact Morozov principle (\ref{eq6}) is a special case of the damped case with $\gamma = \infty$.

In this study, this special Tikhonov regularization method is applied to handle the ill-posed Fredholm problems and the cubic convergence algorithm 
in \cite{wang2001fast, zou2009cubically}, based on the Morozov’s discrepancy principle, is incorporated. Like the Newton
method and the quasi-Newton method, this special method is used in our study for most of the posteriori regularization parameter choice strategies.



\subsection{Numerical diagnostics and reproducibility}
\label{subsec:diagnostics}

The goal of the numerical study is not only to compare pointwise errors for noise-free manufactured examples, but also to evaluate the behavior of each basis when the Fredholm equation is treated as an ill-posed inverse problem. For this reason, each wavelet discretization is assessed through a common set of stabilization diagnostics. The coefficient vector is denoted by $C_{\alpha,N}^{\delta}$, where $N$ is the number of retained basis functions, $\alpha$ is the Tikhonov parameter, and $\delta$ is the prescribed relative perturbation level in the data. The discrete system is written in the form
\[
A_N C \simeq G_N,
\]
where $A_N$ is the collocation matrix obtained from the corresponding wavelet basis.

The two primary accuracy measures are the relative reconstruction error and the relative residual,
\begin{equation}
E_{\rm sol}=
\frac{\|f_{\alpha,N}^{\delta}-f_{\rm exact}\|_2}{\|f_{\rm exact}\|_2},
\qquad
E_{\rm res}=
\frac{\|A_N C_{\alpha,N}^{\delta}-G_N^{\delta}\|_2}{\|G_N^{\delta}\|_2}.
\label{eq:diagnostic-error-definitions}
\end{equation}
The reconstruction error measures accuracy of the unknown function, while the residual measures how closely the regularized solution satisfies the noisy algebraic system. These two quantities need not rank the bases in the same way, which is why both are reported.

To quantify the stabilizing effect of Tikhonov regularization, we also report
\begin{equation}
\kappa(A_N), \qquad \kappa(A_N^T A_N+\alpha I).
\label{eq:condition-number-diagnostics}
\end{equation}
The first quantity measures the ill-conditioning inherited from the compact smoothing operator, whereas the second measures the conditioning of the regularized normal equations actually solved in the computation. This distinction is central to the novelty of the paper: the basis is not judged only by approximation order, but by how it behaves after regularization.

When noisy data are used, the perturbed right-hand side is generated at the collocation points by
\begin{equation}
g^{\delta}(x_i)=g(x_i)+\delta\|g\|_{\infty}\eta_i,
\qquad -1\leq \eta_i\leq 1,
\label{eq:noise-model-reproducibility}
\end{equation}
where $\eta_i$ are random perturbations. The noise levels considered in the diagnostic experiments are
\[
\delta=10^{-4},\;10^{-3},\;10^{-2},\;5\times 10^{-2}.
\]
The reported computations use the stated basis sizes, noise levels, Gaussian kernel width, Morozov tolerance, and a fixed random seed for the noisy-data experiments. Unless otherwise stated, the Gaussian-smoothing diagnostics use the same number of unknowns for all bases, $N=32$, the same collocation grid, and the same random perturbation at each noise level. The random seed used to generate the noisy data reported is fixed at 20250818. 

The numerical framework developed in this work serves not only as a solver for first-kind Fredholm equations, but also as a diagnostic environment for studying the interaction between discretization and regularization. All reported experiments are generated using the same implementation and computational protocol. The basis dimensions, noise levels, Gaussian kernel widths, Morozov tolerances, quadrature rules, random seeds, and post-processing scripts are fixed and recorded to ensure reproducibility. This common protocol allows the observed differences among the tested bases to be attributed to the discretization and stabilization behavior rather than to implementation-dependent choices.

The numerical investigations are designed to address a central question that extends beyond approximation accuracy alone. For ill-posed first-kind Fredholm equations, the effectiveness of a wavelet basis is governed not only by its approximation properties, but also by its influence on the stability of the inverse recovery process. In particular, the study examines how basis selection affects the conditioning of the discrete operator, the conditioning of the regularized normal equations, the discrepancy-selected regularization parameter, residual reduction, robustness to data perturbations, and computational cost. The diagnostics reported below therefore provide a unified framework for understanding the interplay between approximation quality, regularization, conditioning, and noise propagation. This perspective forms the principal theme of the paper and provides practical guidance for selecting appropriate wavelet bases in the numerical solution of ill-posed Fredholm inverse problems.


\section{Numerical examples}\label{sec:examples}
The numerical experiments are designed not only to compare approximation accuracy, but also to investigate how basis selection affects conditioning, regularization strength, residual behavior, robustness to noise, and computational efficiency under a common Tikhonov--Morozov framework. We first compare numerical integration formulas based on Haar wavelets, hybrid functions, and Linear Legendre multiwavelets. We then study one- and two-dimensional Fredholm benchmark problems and finally report stabilization diagnostics for noisy inverse recovery.
\subsection{Test problems for a comparative study of numerical integration}

Example 1: 
\\
$$
\int_{0}^{1}sin(x^2)dx=S(1)=\sum_{n=0}^{\infty} \frac{(-1)^n}{(2n+1)!(4n+3)} ,
$$

 Example 2: \\
%\approx 0.3102683017233811
$$
\int_{0}^{5} \sqrt{x^2-5x+31}dx = \frac{99}{44}+sinh^{-1}(\frac{5}{3\sqrt{11}}),
$$
\hspace{0.6 cm} Example 3: \\
$$
\int_{0}^{1}\frac{e^{\frac{-1}{x}}}{x^2}dx =e^{-1},
$$
\hspace{0.6cm} Example 4: \\
$$
\int_{0}^{\pi}\int_{0}^{\pi}sin(x+y)dxdy = 2,
$$
\hspace{0.6cm} Example 5: \\
$$
\int_{0}^{1}\int_{0}^{1}\frac{1}{\sqrt{x^2+y^2}}dxdy =ln(1+\sqrt2)+tanh^{-1}(\frac{\sqrt{2}}{2}),
$$
\hspace{0.6cm} Example 6: \\
$$
\int_{1}^{2}\int_{1}^{2}\int_{1}^{2}\frac{1}{\sqrt{x+y+z}}dxdydz =4log(\frac{19943936}{9765625})+log(\frac{2500\sqrt{\frac{5}{3}}}{59049}). 
$$

\subsection{Discussion of integration tests}
%\textbf{}\\

To illustrate the accuracy of the different methods we reproduced all the numerical results of the six examples in \citep{aziz2010comparative,sathar2018numerical} by applying Haar wavelets and the Hybrid wavelets.
The comparative analysis between Haar, hybrid and LLMW functions in terms of relative error is performed to find numerical approximations of different types of integrals. As shown in Tables~\ref{tab:int-ex1}--\ref{tab:int-ex6}, the comparison of the numerical results, obtained by LLMW, hybrid functions and Haar wavelets, showed the advantage of the LLMW method over the Haar wavelet method.


\begin{table}[p]
\caption{Relative errors of example 1}\label{tab:int-ex1}
\begin{center}
\resizebox{0.7\textwidth}{!}{ 
\begin{tabular}{|l|l|l|l|l|l|l|r|}
  \hline \multicolumn{2}{c|}{\textit{\textbf{Haar}}} & \multicolumn{2}{c|}{\textit{\textbf{Hybrid}}} & \multicolumn{2}{c|}{\textit{\textbf{LLMW}}}\\  
  
  \hline \multicolumn{2}{c|}{j=4 \;\; 1.4771E-04 } &\multicolumn{2}{c|}{m=3\;n=5\;\; 1.0110E-05 }&\multicolumn{2}{c|}{k=4 \;\; 3.54131E-05 } \\
  
   \hline \multicolumn{2}{c|}{j=5 \;\; 3.5432E-05 } &\multicolumn{2}{c|}{m=4\;n=8\;\; 7.9460E-07 }&\multicolumn{2}{c|}{k=5 \;\; 8.8562E-06 } \\
   
    \hline \multicolumn{2}{c|}{j=6 \;\; 8.8574E-06 } &\multicolumn{2}{c|}{m=5\;n=12\;\; 3.4471E-11 }&\multicolumn{2}{c|}{k=6 \;\; 2.2122E-06 } \\
    
    \hline \multicolumn{2}{c|}{j=7 \;\; 2.2143E-06 } &\multicolumn{2}{c|}{m=6\;n=20\;\; 9.6974E-07 }&\multicolumn{2}{c|}{k=7 \;\; 5.5307E-07 } \\
  
  \hline
\end{tabular}
}
\end{center}
\end{table}

\begin{table}[p]
\caption{Relative errors of example 2}\label{tab:int-ex2}
\begin{center}
\resizebox{0.7\textwidth}{!}{ 
\begin{tabular}{|l|l|l|l|l|l|l|r|}
  \hline \multicolumn{2}{c|}{\textit{\textbf{Haar}}} & \multicolumn{2}{c|}{\textit{\textbf{Hybrid}}} & \multicolumn{2}{c|}{\textit{\textbf{LLMW}}}\\  
  
  \hline \multicolumn{2}{c|}{j=4 \;\; 3,5293E-05 } &\multicolumn{2}{c|}{m=3\;n=5\;\; 3.6467E-07 }&\multicolumn{2}{c|}{k=4 \;\; 8,8228E-06 } \\
  
   \hline \multicolumn{2}{c|}{j=5 \;\; 8.8229E-06 } &\multicolumn{2}{c|}{m=4\;n=8\;\; 2.8648E-08 }&\multicolumn{2}{c|}{k=5 \;\; 2.2063E-06 } \\
   
    \hline \multicolumn{2}{c|}{j=6 \;\; 2.2057E-06 } &\multicolumn{2}{c|}{m=5\;n=12\;\; 1.3262E-12 }&\multicolumn{2}{c|}{k=6 \;\; 5.5321E-07 } \\
    
    \hline \multicolumn{2}{c|}{j=7 \;\; 5.5143E-06 } &\multicolumn{2}{c|}{m=6\;n=20\;\; 3.7148E-14 }&\multicolumn{2}{c|}{k=7 \;\; 1.3946E-07 } \\
  
  \hline
\end{tabular}
}
\end{center}
\end{table}


\begin{table}[p]
\caption{Relative errors of example 3}\label{tab:int-ex3}
\begin{center}
\resizebox{0.7\textwidth}{!}{ 
\begin{tabular}{|l|l|l|l|l|l|l|r|}
  \hline \multicolumn{2}{c|}{\textit{\textbf{Haar}}} & \multicolumn{2}{c|}{\textit{\textbf{Hybrid}}} & \multicolumn{2}{c|}{\textit{\textbf{LLMW}}}\\  
  
  \hline \multicolumn{2}{c|}{j=4 \;\; 4.0642E-05 } &\multicolumn{2}{c|}{m=3\;n=5\;\; 2.9376E-04 }&\multicolumn{2}{c|}{k=4 \;\; 1.0172E-05 } \\
  
   \hline \multicolumn{2}{c|}{j=5 \;\; 1.0173E-05 } &\multicolumn{2}{c|}{m=4\;n=8\;\; 6.6406E-05 }&\multicolumn{2}{c|}{k=5 \;\; 2.5427E-06 } \\
   
    \hline \multicolumn{2}{c|}{j=6 \;\; 2.5431E-06 } &\multicolumn{2}{c|}{m=4\;n=8\;\; 1.8220E-06}&\multicolumn{2}{c|}{k=6 \;\; 6.3526E-07 } \\
    
    \hline \multicolumn{2}{c|}{j=7 \;\; 6.3578E-07 } &\multicolumn{2}{c|}{m=6\;n=20\;\; 3.7947E-08 }&\multicolumn{2}{c|}{k=7 \;\; 1.5602E-07 } \\
  
  \hline
\end{tabular}
}
\end{center}
\end{table}

\begin{table}[p]
\caption{Relative errors of example 4}\label{tab:int-ex4}
\begin{center}
\resizebox{0.7\textwidth}{!}{ 
\begin{tabular}{|l|l|l|l|l|l|l|r|}
  \hline \multicolumn{2}{c|}{\textit{\textbf{Haar}}} & \multicolumn{2}{c|}{\textit{\textbf{Hybrid}}} & \multicolumn{2}{c|}{\textit{\textbf{LLMW}}}\\  
  
  \hline \multicolumn{2}{c|}{j=4 \;\; 5.0215E-05 } &\multicolumn{2}{c|}{m=3\;n=5\;\; 5.3261E-07 }&\multicolumn{2}{c|}{k=4 \;\; 1.2550E-04 } \\
  
   \hline \multicolumn{2}{c|}{j=5 \;\; 1.2551E-05 } &\multicolumn{2}{c|}{m=4\;n=8\;\; 3.4122E-08 }&\multicolumn{2}{c|}{k=5 \;\; 3.1356E-05 } \\
   
    \hline \multicolumn{2}{c|}{j=6 \;\; 3.1375E-05 } &\multicolumn{2}{c|}{m=4\;n=8\;\; 2.4702E-13}&\multicolumn{2}{c|}{k=6 \;\;7.8390E-06 } \\
    
    \hline \multicolumn{2}{c|}{j=7 \;\; 7.8437E-06 } &\multicolumn{2}{c|}{- \;\;\;\; - }&\multicolumn{2}{c|}{k=7 \;\; 1.92752E-06 } \\
  
  \hline
\end{tabular}
}
\end{center}
\end{table}

\begin{table}[p]
\caption{Relative errors of example 5}\label{tab:int-ex5}
\begin{center}
\resizebox{0.7\textwidth}{!}{ 
\begin{tabular}{|l|l|l|l|l|l|l|r|}
  \hline \multicolumn{2}{c|}{\textit{\textbf{Haar}}} & \multicolumn{2}{c|}{\textit{\textbf{Hybrid}}} & \multicolumn{2}{c|}{\textit{\textbf{LLMW}}}\\  
  
  \hline \multicolumn{2}{c|}{j=4 \;\; 7.1275E-03 } &\multicolumn{2}{c|}{m=3\;n=5\;\; 8.98636E-03 }&\multicolumn{2}{c|}{k=4 \;\; 3.5718E-03 } \\
  
   \hline \multicolumn{2}{c|}{j=5 \;\; 3.5719E-03 } &\multicolumn{2}{c|}{m=4\;n=8\;\; 4.6532E-03 }&\multicolumn{2}{c|}{k=5 \;\; 1,7780E-03 } \\
   
    \hline \multicolumn{2}{c|}{j=6 \;\; 1.7880E-03 } &\multicolumn{2}{c|}{m=4\;n=8\;\; 1.8550E-03}&\multicolumn{2}{c|}{k=6 \;\; 8,9449E-04 } \\
    
    \hline \multicolumn{2}{c|}{j=7 \;\; 8.9450E-04 } &\multicolumn{2}{c|}{- \;\;\;\; - }&\multicolumn{2}{c|}{k=7 \;\; 4.4740E-04 } \\
  
  \hline
\end{tabular}
}
\end{center}
\end{table}

\begin{table}[p]
\caption{Relative errors of example 6}\label{tab:int-ex6}
\begin{center}
\resizebox{0.4\textwidth}{!}{ 
\begin{tabular}{|l|l|l|l|l|l|l|r|}
  \hline \multicolumn{2}{c|}{\textit{\textbf{Haar}}} & \multicolumn{2}{|}{\textit{\textbf{}}} & \multicolumn{2}{c|}{\textit{\textbf{LLMW}}}\\  
  
  \hline \multicolumn{2}{c|}{j=4 \;\; 1.2863E-05 } &\multicolumn{2}{|}{}&\multicolumn{2}{c|}{k=4 \;\; 3.2103E-06 } \\
  
   \hline \multicolumn{2}{c|}{j=5 \;\; 3.2159E-06 } &\multicolumn{2}{|}{}&\multicolumn{2}{c|}{k=5 \;\; 1.7970E-07 } \\
   
    \hline \multicolumn{2}{c|}{j=6 \;\; 8.0398E-07 } &\multicolumn{2}{|}{}&\multicolumn{2}{c|}{k=6 \;\; 1.9506E-07 } \\
    

  
  \hline
\end{tabular}
}
\end{center}
\end{table}

Another important feature we found is that the hybrid function approach is slightly better than the Haar Wavelet approach for all cases excluding the improper integral cases ( example 3 and example 5 ) due to the effect of singularity. A better accuracy is achieved as we increase the number of nodes. For example, the relative error for the example 2  is in the order of $10^{-7}$ for a small number of nodes $(m=3, n=5)$; and it is in the order of $10^{-14}$ for a larger number of nodes $(m=6, n=20)$, which clearly indicates the efficiency and the superiority of this hybrid approach.
 

It is also noted that LLMW functions have slightly more accurate results than the Haar wavelet functions as shown in Tables~\ref{tab:int-ex1}--\ref{tab:int-ex6}.

 In terms of CPU performance, the CPU time is computed for different values of dilation factor k in the cases of the Haar wavelets and LLMW functions and for different values of m in the case of the hybrid functions.
 As shown in Figure~\ref{fig:cpu-wavelet-methods}, the numerical results show that the Haar wavelet approach is $10$ times faster than the hybrid method, and the LLMW  method is much faster than the Haar wavelet method. The simple applicability of Haar wavelets and the fast convergence of LLMW method provide a solid foundation for using these functions in the context of numerical approximations of integral equations and differential equations. In what follows, we will use these tools to solve a few examples of ill-posed Fredholm integral equations of the first kind.
 

\begin{figure}[p]
    \centering
 \includegraphics[width=11cm]{courbe.jpg}
    \caption{CPU time comparison of different wavelet methods}
    \label{fig:cpu-wavelet-methods}
   \end{figure}

\subsection{Fredholm integral equations of the first kind }
In this section numerical results are presented to illustrate the efficiency and accuracy of the numerical methods that were described in the previous section. We consider the following Fredholm integral equations:\\


{\textbf{Problem 1}}\\

\begin{equation}
\begin{aligned}
\int_0^1 \sin(st) f(s) ds=\frac{\sin(t)-t \cos(t)}{t^2} .
\end{aligned}
\end{equation}

with the exact solution  $f(t)=t$. \\

{\textbf{Problem 2}}\\

\begin{equation}
\begin{aligned}
\int_0^1 e^{s^2t}f(t)dt=g(s)=\frac{e^{s^2+1}-1}{s^2+1} 
\end{aligned}
\end{equation}

with the exact solution  $f(t)=e^t$. \\

%{\textbf{Problem 3}}\\
%\begin{equation} \label{eq11}
%\begin{aligned}
%\int_{0}^{1} \sqrt{se^{t}} f(s)ds= \frac{e^{t/2}}{2},
%\end{aligned}
%\end{equation}

%with the exact solution $f(t)=t^{1/2}$. \\

%{\textbf{Problem 4}}\\
%\begin{equation} \label{eq12}
%\begin{aligned}
%\int_{0}^{\pi/2} \cos(t^2+3t+1) f(s)ds= \frac{1}{6} %\cos(t^2+2)-\frac{\pi}{4} \sin(t^2),
%\end{aligned}
%\end{equation}
%with the exact solution $f(t)=sin(3t+1)$. \\

%{\textbf{Problem 5}}\\
%\begin{equation}
%\begin{aligned}
%\frac{1}{2}(e^2-1)e^{(x+y)} = \int_0^1\int_0^1e^{(x+y+s+t)}g(s,t)dtds,
%\end{aligned}
%\end{equation}
%with 
%\begin{equation*}
%\begin{aligned}
%g(x,y)= \frac{2e^{(x+y)}}{e^2-1}.
%\end{aligned}
%\end{equation*} \\



\subsection{Discussion of Fredholm benchmark problems}

We used the Linear Legendre Multiwavelet approach and the Tikhonov regularization and applied the cubic discrepancy principle, which were discussed in Section~\ref{sec:method}, to solve the ill-posed integral Fredholm equations of the problems $1-3$.

The numerical solution of the integral equation for the two problems is obtained and excellent agreement with the analytical solution is shown.

In Table~\ref{tab:abs-prob1} the absolute errors of the proposed numerical method for problem 1 are compared with the numerical results obtained by the Legendre wavelet method in \cite{maleknejad2007numerical} for $M=1$ and with the Chebyshev wavelet method in \cite{adibi2010chebyshev} for $M=2, k=2$. Furthermore, another comparison between the absolute errors of our method for problem 2  with the Legendre wavelets method \cite{maleknejad2007numerical} for $M=1$ is given in Table~\ref{tab:abs-prob2}. This also demonstrates the superiority of the proposed method for the solution of the problems 1 and 2.

\begin{table}[ht]
\caption{Absolute errors for Problem 1}\label{tab:abs-prob1} % title of Table
\begin{center}
\resizebox{0.9\textwidth}{!}{ 
\begin{tabular}{c c c c c c} % centered columns (4 columns)
\hline %inserts double horizontal lines
%\hspace{5 mm}
t & Present Method   & Method in \cite{maleknejad2007numerical} & Present Method  & Present Method & Method in \cite{adibi2010chebyshev}\\
 & with M=1   & with M=1 & with M=2, k=2  & with M=3, k=2 & with M=2, k=2 \\ [0.4ex] % inserts table
%heading
\hline % inserts single horizontal line
0.1 & 8.4554E-06 & 3.5427E-05 & 1.93E-07 & 5.37E-11 & 5.32E-05\\ % inserting body of the table
0.2 & 1.5395 E-05 & 1.1933E-05 & 9.22E-08 & 1.81 E-11 & 1.37E-04 \\
0.3 & 9.5703E-06 & 2.3777E-05 & 8.48E-09 & 3.61E-11 & 2.21E-04 \\
0.4 & 5.9440 E-07 & 1.0392E-07 & 1.09E-07 & 1.57E-13 & 3.04E-04 \\
0.5 & 1.9835E-06 & 1.1210E-05 &  5.23E-08 & 1.70E-11 & 5.98E-04 \\
0.6 & 8.7416E-07 & 5.3623E-06 & 3.85E-08 & 8.14E-12 & 3.80E-04 \\
0.7 & 2.3518E-07 & 1.4280E-06 & 2.47E-08 & 2.16E-12 & 5.51E-04 \\
0.8 & 4.4564E-08 & 5.9249E-07 & 1.06E-08 & 8.99E-13 & 7.21E-04 \\
0.9 & 1.5619E-08 & 6.9940E-07 & 2.82E-09 & 1.06E-12 & 8.92E-04 \\ [1ex] % [1ex] adds vertical space
 \hline
\end{tabular}
%\label{table:nonlin} % is used to refer this 
 
%\end{tabular}
}
\end{center}
\end{table}



\begin{table}[p]
\caption{Absolute errors for Problem 2}\label{tab:abs-prob2} % title of Table
\begin{center}
\resizebox{0.4\textwidth}{!}{ 
\begin{tabular}{c c c } % centered columns (4 columns)

 \hline %inserts double horizontal lines

t & Present Method     & Method in \cite{maleknejad2007numerical} \\
 & with M=1     &  with M=1 \\ [0.4ex] % inserts table
%heading
\hline % inserts single horizontal line
0.1 & 2.098E-04  & 2.469E-04\\ % inserting body of the table
0.2 & 2.229E-04 &   2.802E-04 \\
0.3 & 4.640E-04  &  3.790E-04 \\
0.4 & 7.721E-04 &   8.005E-04 \\
0.5 & 6.144E-03 &   3.136E-02 \\
0.6 & 5.760E-04 &  3.668E-04 \\
0.7 & 9.013E-04 &   1.386E-03 \\
0.8 & 2.060E-04 & 1.051E-03 \\
0.9 & 5.423E-04 &  4.115E-03 \\ %[1ex] % [1ex] adds vertical space
 \hline
\end{tabular}
%\label{table:nonlin} % is used to refer this 

}
\end{center}
\end{table}

%Comparison between the results of this method with the results of method given in \cite{bahmanpour2019solving}.
%and with other wavelets methods in the literature is given in Table2 and in Figures...
%We used the Tikhonov regularization and applied the cubic discrepancy principle to solve this ill-integral Fredholm equation. Comparison between the results of this method with the method presented in [] and with other wavelets methods in the literature is given in Table2 and in Figures...

%It is observed that the current method can provide slightly better results for certain cases, like in example 3, as shown in Figure 4 and Figure 5. This shows that for this example, when the power of the solution is fractional, the method given in \cite{bahmanpour2019solving}, which is based on Muntz wavelets provided very accurate results but not better than our method. This also demonstrates the superiority of the proposed method for solutions containing fractional powers.

\begin{comment}


%\begin{figure}[p]
 %   \centering
 \includegraphics[width=9cm]{fg4.png}
 %\vspace{-1 cm}
    \caption{Absolute error for Example 2 using the current method with $N=3$,  $\delta=10^{-5}$.}
    \label{error-ex1 }
   \end{figure}
 % \vspace{-2cm} 
\begin{figure}[p]
    \centering
 \includegraphics[width=10cm]{fg1.pdf}
 \vspace{-1 cm}
    \caption{Absolute error for Example 2  using the Muntz wavelets with $k=1, M=5$ \cite{bahmanpour2019solving}}, 
    \label{error-ex1 }
   \end{figure}
   
  \begin{figure}[p]
    \centering
 \includegraphics[width=9cm]{fg5.png}
 %\vspace{-1 cm}
    \caption{Absolute error for Example 3 using the current method with $N=3$,  $\delta=10^{-5}$.}
    \label{error-ex1 }
   \end{figure}
   
\begin{figure}[p]
    \centering
 \includegraphics[width=9.7cm]{fg2.pdf}
 \vspace{-1 cm}
    \caption{Absolute error for Example 3  using the Muntz wavelets with $k=1, M=6$ \cite{bahmanpour2019solving}.}
    \label{error-ex1 }
   \end{figure}   
   
    \begin{figure}[p]
    \centering
 \includegraphics[width=9cm]{fg6.png}
 %\vspace{-1 cm}
    \caption{Absolute error for Example 2 using the current method with $N=3$,  $\delta=10^{-4}$.}
   \label{error-ex1 }
  \end{figure} 
 \end{comment}  
   

%We report in figure (5.1) the obtained results of the relative error and the approximated functions obtained applying the different approximating methods. The noise level admitted here is ($\delta = 1e-3 $). For each method, the choice of the regularization parameter $\alpha$ was effectuated applying the Damped Morozov's discrepancy principle.    

%{%
  % \begin{figure}[p]
   %  \begin{center}
     % \includegraphics[width=15cm]{fig51.png}
               % \caption{Relative error and Approximated functions applying % Haar-wavelet, Legendre-wavelet, Hybrid-wavelet and the  % Projected-Tikhonov method }
 %    \label{Absolute error corresponding to example (1,2,3) on the interval % [0,1] }
%     \end{center}
%  \end{figure}
%  }%
%As it is shown in figure (5.1), the best results of relative error are obtained applying the Legendre wavelet method and the Projected Tikhonov method. For Legendre wavelet, the relative error is of the order (1e-3/1e-4). For Projected method, the relative error is of the order (1e-2/1e-4). For Hybrid wavelet, the relative error is of the order (1e-2) and for Haar wavelet, the relative error is of the order (1e-1/1e-2). In terms of accuracy, the orthonormal Legender wavelet method and the projected method are the best. For Hybrid wavelet, which is not an orthonormal wavelet method, the solution is not very accurate. For Haar wavelet, it is an orthonormal method, it is not accurate but it can be considered as the simplest and the fastest wavelet method. Furthermore, one strong point of the wavelet approach is that,as it is meshless method, it is economic in terms of CPU time, comparing to trapezoidal method implemented in the Projected approach. Furthermore, another criteria can be adopted in this discussion, which is the code complexity. As a matter of fact, the wavelet methods are easily implemented, characterized by its multi-resolution capability and can be conveniently extended to 2D/3D cases.  And as they are a mesh-less method, it is very economic in terms of CPU allocated memory, comparing to finite elements approaches.\\

%In order to highlight the stability of these approximating methods, we represent in figure (5.2) the average of the relative error over [0,1] of the approximate solution against the noise level $\delta$. Comparing the obtained results, we can notice that all the approximated method are relatively stable against the noise level, but the most stable one remains the hybrid method.\\

% \begin{figure}[p]
 %    \centering
  %   \includegraphics[width=15cm]{fig52.png}
  %             \caption{Average value of relative error obtained applying Haar-wavelet, Legendre-wavelet, Hybrid-wavelet and the Projected-Tikhonov method against the perturbation noise level admitted $\delta$}
 % \end{figure}




%






% order to provide a very good test for the 2D Legendre wavelet method (k=2,k'=1,M= M'=2) to perform the inversion, we consider the 2D-Fredholm integral equation of the second kind given as follows:
%begin{equation}
%\begin{aligned}
%\frac{1}{2}(e^2-1)e^{(x+y)} = %\int_0^1\int_0^1e^{(x+y+s+t)}g(s,t)dtds
%\end{aligned}
%end{equation}
%with : 
%begin{equation*}
%begin{aligned}
%g(x,y)= \frac{2e^{(x+y)}}{e^2-1}
%end{aligned}
%\end{equation*}
% a simpler way, we are going to use the 2D- Legendre wavelet transformer to convert the Fredholm integral equation of the first kind to a system of linear equations. The solution to this ill-posed problem is determined applying the Tikhonov regularization method. Finaly, the inverse 2D wavelet transformer is used to evaluate the unknown function solution to the Fredholm problem. One important feature of this approximation method, is that it takes into consideration the error due to measurements. We report in figure (5.3) the results of the relative error and the approximated functions obtained applying the proposed method. The noise level admitted here is ($\delta = 1e-3 $).

% \begin{figure}[p]
 %    \centering
 %    \includegraphics[width=15cm]{fig53.png}
 %               \caption{Relative error and Exact function versus approximated function obtained applying the 2D Legendre-wavelet approach}
%  \end{figure}
%As it is shown in figure (5.3), the Legendre wavelet method generates a good range of accuracy. The extension of the proposed method to 2D space dimension is very flexible. Furthermore, this method is not cost in terms of CPU time and allocated memory, as it is a mesh-less method. 

%Similarly to the 1D case, we represent in figure (5.4) the average value of the relative error over [0,1] of the approximate solution against the noise level $\delta$. We can obviously notice that the stability can be considered as a strong point of the proposed method.
% \begin{figure}[p]
%     \centering
 %    \includegraphics[width=15cm]{fig54.png}
 %               \caption{Relative error and Exact function versus approximated function obtained applying the 2D Legendre-wavelet approach}
 % \end{figure}







\subsection{Basis-dependent stabilization diagnostics for a Gaussian smoothing problem}
\label{subsec:gaussian-diagnostics}

The preceding examples confirm that the proposed wavelet discretizations reproduce standard benchmark solutions. We now add a more targeted experiment whose purpose is to expose the part of the method that is most relevant for ill-posed inverse recovery: the interaction between basis choice and regularization. We consider the compact smoothing operator
\begin{equation}
\int_0^1 \exp\left(-\frac{(x-y)^2}{\sigma^2}\right) f(y)\,dy=g(x),
\qquad 0\le x\le 1,
\label{eq:gaussian-smoothing}
\end{equation}
with $\sigma=0.1$. This kernel strongly smooths the unknown function, so the inverse problem is sensitive to perturbations in $g$. The right-hand side is generated by applying the operator to a known exact solution. Two targets are used:
\begin{equation}
f_1(y)=\sin(\pi y)+y(1-y),
\qquad
f_2(y)=|y-1/2|.
\label{eq:smooth-nonsmooth-targets}
\end{equation}
The first target is smooth and favors higher-order local approximation. The second target has a corner and tests whether the basis remains reliable when the recovered solution has reduced regularity.

All bases are compared at the same algebraic size, $N=32$, with the same collocation points, the same noise realizations, and the same discrepancy-based selection of $\alpha$. This equal-degree-of-freedom protocol is important because it separates the effect of the approximation basis from the effect of changing the size of the discrete inverse problem. The data perturbation model is the one in Eq.~\eqref{eq:noise-model-reproducibility}, and the reported diagnostics are those defined in Eqs.~\eqref{eq:diagnostic-error-definitions} and \eqref{eq:condition-number-diagnostics}.


\begin{table}[htbp]
\centering
\scriptsize
\caption{Noise-sensitivity and stabilization behavior for the Gaussian smoothing Fredholm problem with a smooth exact solution. For each basis and noise level, the table reports the Morozov-selected parameter $\alpha$, the relative solution error, the relative residual error, and the conditioning of the regularized normal matrix.}
\label{tab:noise-results}
\begin{tabular}{lccccc}
\toprule
Basis & $\delta$ & $\alpha$ & $E_{\rm sol}$ & $E_{\rm res}$ & $\kappa(A^T A+\alpha I)$ \\
\midrule
Haar & $1\times10^{-4}$ & $1.36\times10^{-4}$ & $5.12\times10^{-2}$ & $9.58\times10^{-5}$ & $3.81\times10^{4}$ \\
Haar & $1\times10^{-3}$ & $1.35\times10^{-3}$ & $5.19\times10^{-2}$ & $9.64\times10^{-4}$ & $3.84\times10^{3}$ \\
Haar & $1\times10^{-2}$ & $1.53\times10^{-2}$ & $5.70\times10^{-2}$ & $9.68\times10^{-3}$ & $3.40\times10^{2}$ \\
Haar & $5\times10^{-2}$ & $7.71\times10^{-2}$ & $6.06\times10^{-2}$ & $5.05\times10^{-2}$ & $6.82\times10^{1}$ \\
Hybrid & $1\times10^{-4}$ & $1.36\times10^{-4}$ & $3.30\times10^{-2}$ & $8.84\times10^{-5}$ & $3.81\times10^{4}$ \\
Hybrid & $1\times10^{-3}$ & $1.77\times10^{-3}$ & $3.41\times10^{-2}$ & $9.82\times10^{-4}$ & $2.94\times10^{3}$ \\
Hybrid & $1\times10^{-2}$ & $1.75\times10^{-2}$ & $5.02\times10^{-2}$ & $9.87\times10^{-3}$ & $2.97\times10^{2}$ \\
Hybrid & $5\times10^{-2}$ & $7.71\times10^{-2}$ & $6.42\times10^{-2}$ & $4.67\times10^{-2}$ & $6.82\times10^{1}$ \\
Legendre & $1\times10^{-4}$ & $1.36\times10^{-4}$ & $3.71\times10^{-2}$ & $1.00\times10^{-4}$ & $3.81\times10^{4}$ \\
Legendre & $1\times10^{-3}$ & $1.54\times10^{-3}$ & $3.80\times10^{-2}$ & $1.01\times10^{-3}$ & $3.36\times10^{3}$ \\
Legendre & $1\times10^{-2}$ & $1.34\times10^{-2}$ & $3.90\times10^{-2}$ & $9.31\times10^{-3}$ & $3.89\times10^{2}$ \\
Legendre & $5\times10^{-2}$ & $1.01\times10^{-1}$ & $7.61\times10^{-2}$ & $4.64\times10^{-2}$ & $5.24\times10^{1}$ \\
B-spline & $1\times10^{-4}$ & $3.15\times10^{-7}$ & $3.35\times10^{-2}$ & $8.64\times10^{-5}$ & $5.64\times10^{5}$ \\
B-spline & $1\times10^{-3}$ & $2.36\times10^{-5}$ & $4.98\times10^{-2}$ & $9.94\times10^{-4}$ & $7.52\times10^{3}$ \\
B-spline & $1\times10^{-2}$ & $4.00\times10^{-4}$ & $6.24\times10^{-2}$ & $9.40\times10^{-3}$ & $4.44\times10^{2}$ \\
B-spline & $5\times10^{-2}$ & $2.31\times10^{-3}$ & $7.76\times10^{-2}$ & $4.76\times10^{-2}$ & $7.77\times10^{1}$ \\
LLMW & $1\times10^{-4}$ & $1.36\times10^{-4}$ & $5.26\times10^{-2}$ & $9.92\times10^{-5}$ & $3.81\times10^{4}$ \\
LLMW & $1\times10^{-3}$ & $1.77\times10^{-3}$ & $5.42\times10^{-2}$ & $9.46\times10^{-4}$ & $2.94\times10^{3}$ \\
LLMW & $1\times10^{-2}$ & $1.75\times10^{-2}$ & $6.01\times10^{-2}$ & $9.27\times10^{-3}$ & $2.97\times10^{2}$ \\
LLMW & $5\times10^{-2}$ & $7.71\times10^{-2}$ & $7.85\times10^{-2}$ & $4.57\times10^{-2}$ & $6.82\times10^{1}$ \\
\bottomrule
\end{tabular}
\end{table}

\begin{table}[htbp]
\centering
\scriptsize
\caption{Conditioning diagnostic at noise level $\delta=10^{-3}$. The table compares the ill-conditioned collocation matrix with the corresponding Tikhonov-regularized normal matrix.}
\label{tab:conditioning-results}
\begin{tabular}{lccccc}
\toprule
Basis & $N$ & $\kappa(A)$ & $\alpha$ & $\kappa(A^T A+\alpha I)$ & $E_{\rm sol}$ \\
\midrule
Haar & 32 & $3.09\times10^{16}$ & $1.35\times10^{-3}$ & $3.84\times10^{3}$ & $5.19\times10^{-2}$ \\
Hybrid & 32 & $3.28\times10^{16}$ & $1.77\times10^{-3}$ & $2.94\times10^{3}$ & $3.41\times10^{-2}$ \\
Legendre & 32 & $2.09\times10^{16}$ & $1.54\times10^{-3}$ & $3.36\times10^{3}$ & $3.80\times10^{-2}$ \\
B-spline & 32 & $1.15\times10^{16}$ & $2.36\times10^{-5}$ & $7.52\times10^{3}$ & $4.98\times10^{-2}$ \\
LLMW & 32 & $1.79\times10^{16}$ & $1.77\times10^{-3}$ & $2.94\times10^{3}$ & $5.42\times10^{-2}$ \\
\bottomrule
\end{tabular}
\end{table}

\begin{table}[htbp]
\centering
\scriptsize
\caption{Computational cost for the clean Gaussian smoothing benchmark. The table reports assembly time, solve time, total time, and the corresponding solution error.}
\label{tab:cost-results}
\begin{tabular}{lccccc}
\toprule
Basis & $N$ & Assembly & Solve & Total & $E_{\rm sol}$ \\
\midrule
Haar & 32 & $3.499\times10^{0}$ & $1.694\times10^{-4}$ & $3.499\times10^{0}$ & $4.98\times10^{-2}$ \\
Hybrid & 32 & $6.838\times10^{-1}$ & $1.045\times10^{-4}$ & $6.839\times10^{-1}$ & $3.18\times10^{-2}$ \\
Legendre & 32 & $5.733\times10^{-1}$ & $1.613\times10^{-4}$ & $5.734\times10^{-1}$ & $3.64\times10^{-2}$ \\
B-spline & 32 & $1.474\times10^{0}$ & $1.137\times10^{-4}$ & $1.474\times10^{0}$ & $9.59\times10^{-3}$ \\
LLMW & 32 & $1.274\times10^{0}$ & $1.027\times10^{-4}$ & $1.274\times10^{0}$ & $5.17\times10^{-2}$ \\
\bottomrule
\end{tabular}
\end{table}

\begin{table}[htbp]
\centering
\scriptsize
\caption{Nonsmooth reconstruction test for $f(y)=|y-1/2|$ at $\delta=10^{-3}$. The table reports error, residual, selected $\alpha$, and stabilization behavior.}
\label{tab:nonsmooth-results}
\begin{tabular}{lccccc}
\toprule
Basis & $N$ & $\alpha$ & $E_{\rm sol}$ & $E_{\rm res}$ & $\kappa(A^T A+\alpha I)$ \\
\midrule
Haar & 32 & $4.58\times10^{-4}$ & $5.13\times10^{-2}$ & $8.17\times10^{-4}$ & $1.13\times10^{4}$ \\
Hybrid & 32 & $6.87\times10^{-4}$ & $4.59\times10^{-2}$ & $8.17\times10^{-4}$ & $7.55\times10^{3}$ \\
Legendre & 32 & $4.58\times10^{-4}$ & $3.84\times10^{-2}$ & $8.40\times10^{-4}$ & $1.13\times10^{4}$ \\
B-spline & 32 & $1.59\times10^{-6}$ & $1.32\times10^{-1}$ & $9.31\times10^{-4}$ & $1.12\times10^{5}$ \\
LLMW & 32 & $5.24\times10^{-4}$ & $4.91\times10^{-2}$ & $8.42\times10^{-4}$ & $9.89\times10^{3}$ \\
\bottomrule
\end{tabular}
\end{table}


Tables~\ref{tab:noise-results}--\ref{tab:nonsmooth-results} summarize the new diagnostic evidence. Table~\ref{tab:noise-results} shows how the reconstruction error, residual error, selected value of $\alpha$, and regularized condition number vary with the noise level. Table~\ref{tab:conditioning-results} isolates the stabilization effect at $\delta=10^{-3}$ by comparing $\kappa(A_N)$ with $\kappa(A_N^TA_N+\alpha I)$. Table~\ref{tab:cost-results} reports wall-clock assembly and solve times for the clean Gaussian benchmark, and Table~\ref{tab:nonsmooth-results} repeats the diagnostic comparison for the nonsmooth target. Taken together, these tables are the main additional contribution of this section: they show that basis selection is not a purely approximation-theoretic question once the equation is noisy and ill-conditioned.

Table~\ref{tab:noise-results} shows that increasing the perturbation level generally increases the discrepancy-selected regularization parameter and the residual error. This behavior is desirable: the method does not use a fixed artificial damping parameter, but adapts the amount of stabilization to the amount of noise in the data. The table also shows that the ordering of the bases is not determined by clean-data accuracy alone; the same basis can behave differently as the noise level changes.

Table~\ref{tab:conditioning-results} gives the clearest evidence that the new diagnostics are measuring stabilization, not only approximation. The unregularized collocation matrices have condition numbers of order $10^{16}$, while the regularized normal equations are reduced by many orders of magnitude. This explains why a regularization-centered comparison is needed for first-kind equations: without this stabilization information, a small approximation error may hide a numerically unstable inverse recovery.

Table~\ref{tab:cost-results} adds the computational dimension to the comparison. Since all bases are tested with the same number of unknowns, the differences in cost are due mainly to basis evaluation, support structure, and assembly complexity. These timings make the conclusions more useful in practice, because the preferred basis depends on the combined balance of accuracy, conditioning, robustness, and cost. Figure~\ref{fig:cpu-cost-diagnostic} visualizes the same cost comparison and highlights that accuracy and time-to-solution do not rank the bases in the same order.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.78\textwidth]{figures/fig_cpu_cost.png}
\caption{Computational cost for the Gaussian smoothing diagnostic test. The bars show total wall-clock time for assembly and one regularized solve at the same basis size $N=32$; the annotations report the corresponding clean-data solution error. This figure shows that the basis with the smallest error is not necessarily the least expensive basis.}
\label{fig:cpu-cost-diagnostic}
\end{figure}

Figure~\ref{fig:cpu-cost-diagnostic} gives a practical interpretation of the cost table. Haar is simple but not always the fastest in this implementation because matrix assembly depends on the details of basis evaluation and support handling. B-spline and LLMW bases require more work than Legendre or hybrid bases in the clean Gaussian test, but they may be preferable when their approximation properties or local smoothness are important. This is why the paper reports cost together with error and conditioning, instead of reporting accuracy alone.

Table~\ref{tab:nonsmooth-results} shows that reduced regularity changes the basis ranking. Polynomial-type and spline-type bases are highly effective for smooth targets, but a nonsmooth point limits the advantage of high local order and makes compact local behavior more important. This supports the practical conclusion that basis choice should be matched to the expected smoothness of the unknown function as well as to the noise level.

Figure~\ref{fig:reconstruction-smooth} shows representative reconstructions for the smooth target at $\delta=10^{-3}$. The higher-order local bases capture the smooth profile well, while the low-order Haar approximation remains robust but less smooth. Figure~\ref{fig:reconstruction-nonsmooth} repeats the comparison for the nonsmooth target. In this case, the relative ranking changes: the best basis is no longer determined only by polynomial order, because local behavior near the corner becomes important.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.82\textwidth]{figures/fig_reconstruction_smooth.png}
\caption{Gaussian smoothing problem with smooth target $f_1(y)=\sin(\pi y)+y(1-y)$ at $\delta=10^{-3}$. The plot compares the exact function with reconstructions obtained from the regularized wavelet-collocation systems.}
\label{fig:reconstruction-smooth}
\end{figure}

The smooth reconstruction in Figure~\ref{fig:reconstruction-smooth} illustrates the advantage of bases that can represent smooth local variation efficiently. The differences between the curves should not be interpreted only as interpolation differences; they are the result of the full stabilized inverse solve. Thus, the figure visually confirms the table-based observation that approximation order and regularization interact in the final reconstruction.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.82\textwidth]{figures/fig_reconstruction_nonsmooth.png}
\caption{Gaussian smoothing problem with nonsmooth target $f_2(y)=|y-1/2|$ at $\delta=10^{-3}$. The test illustrates the effect of reduced solution regularity on the relative performance of the bases.}
\label{fig:reconstruction-nonsmooth}
\end{figure}

Figure~\ref{fig:reconstruction-nonsmooth} demonstrates that the nonsmooth target is a more demanding inverse test. The corner at $y=1/2$ reduces the benefit of high-order polynomial approximation and makes local support more relevant. This is why the nonsmooth experiment is included: it prevents the novelty claim from relying only on favorable smooth manufactured examples.

Figure~\ref{fig:error-vs-noise} summarizes the solution error as the data noise increases. For small noise, smooth higher-order bases can exploit the regularity of the target. As the noise level increases, however, the regularization parameter becomes more influential and the observed error reflects a combined approximation--stabilization effect. This is precisely the behavior that is missed if bases are compared only on noise-free examples.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.78\textwidth]{figures/fig_error_vs_noise.png}
\caption{Relative solution error versus noise level for the Gaussian smoothing benchmark. The comparison uses the same degrees of freedom and the same discrepancy-based regularization strategy for all bases.}
\label{fig:error-vs-noise}
\end{figure}

The trend in Figure~\ref{fig:error-vs-noise} is important because it separates approximation quality from robustness. A basis that performs well at a very small perturbation need not remain the best as $\delta$ increases. This confirms the central message of the paper: for first-kind equations, basis quality must be judged after regularization and under data perturbations.

Figure~\ref{fig:alpha-vs-noise} reports the selected value of the regularization parameter. The observed increase of $\alpha$ with the perturbation level is consistent with Morozov-type parameter selection: larger data perturbations require stronger stabilization. This diagnostic confirms that the parameter-selection rule is reacting to the noise level rather than acting as a fixed artificial damping.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.78\textwidth]{figures/fig_alpha_vs_noise.png}
\caption{Discrepancy-selected regularization parameter as a function of the relative noise level. Larger perturbations generally lead to stronger stabilization.}
\label{fig:alpha-vs-noise}
\end{figure}

Figure~\ref{fig:alpha-vs-noise} also shows that the basis affects the parameter-selection process. This is expected from the structure of the discrete inverse problem. Each basis produces a different discrete spectrum, a different distribution of small singular values, and a different representation of the noise-contaminated data. Since Morozov's discrepancy principle selects $\alpha$ by balancing the residual against the prescribed perturbation level, the selected parameter necessarily reflects both the noise level and the basis-dependent spectral response of the discretized operator. This observation is central to the paper: the basis is not only an approximation device, but also a component of the regularization mechanism.

The conditioning results in Table~\ref{tab:conditioning-results} and Figure~\ref{fig:conditioning} show the main stabilization mechanism. The unregularized matrices are severely ill-conditioned, which is typical for discretizations of compact first-kind operators. The regularized normal matrices are substantially better conditioned. This confirms that the improvement is not merely visual smoothing: the algebraic inverse problem itself becomes more stable after adding the Tikhonov term.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.78\textwidth]{figures/fig_conditioning.png}
\caption{Conditioning before and after regularization for the Gaussian smoothing benchmark at $\delta=10^{-3}$. The Tikhonov term strongly reduces the effective conditioning of the normal equations.}
\label{fig:conditioning}
\end{figure}

The conditioning plot in Figure~\ref{fig:conditioning} gives a direct algebraic explanation for the stabilizing effect. The large gap between the unregularized and regularized condition numbers shows that the Tikhonov term changes the numerical character of the linear system. Consequently, the observed reconstructions should be viewed as regularized inverse recoveries, not as ordinary collocation approximations.

Figure~\ref{fig:lcurve} gives an L-curve diagnostic for the LLMW discretization. The curve displays the trade-off between residual reduction and solution-norm growth as $\alpha$ varies. Very small values of $\alpha$ fit the data more aggressively but amplify unstable components, whereas large values over-regularize. The discrepancy-selected parameter lies in the transition region, which is the expected behavior for stable recovery from noisy data.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.70\textwidth]{figures/fig_lcurve_llmw.png}
\caption{L-curve diagnostic for the LLMW basis in the Gaussian smoothing benchmark. The curve illustrates the trade-off between the residual norm and the coefficient norm as $\alpha$ varies.}
\label{fig:lcurve}
\end{figure}

The L-curve in Figure~\ref{fig:lcurve} provides an independent consistency check on the selected parameter. The transition region represents the compromise between fitting noisy data and controlling the norm of the recovered coefficients. This supports the use of Morozov-selected regularization in the diagnostic comparison and helps justify the chosen $\alpha$ values in the tables.

The computational results in Table~\ref{tab:cost-results} add a practical basis-selection perspective. Although all bases have the same number of degrees of freedom, the assembly cost depends on support, smoothness, and evaluation cost. Consequently, the best basis is not determined solely by the smallest reconstruction error, but by the combined behavior of accuracy, conditioning, robustness, and cost.


\subsection*{Practical interpretation and synthesis}

The numerical evidence indicates that basis selection influences the regularized inverse problem through two competing mechanisms. Higher-order local bases generally improve approximation quality for smooth solutions and often produce the smallest reconstruction errors after stabilization. These gains, however, may come with increased assembly cost and a stronger dependence on the regularization parameter. In contrast, simpler bases such as Haar wavelets usually provide less accurate smooth reconstructions, but they retain attractive locality, implementation simplicity, and computational efficiency. Hybrid constructions occupy an intermediate position, combining local structure with moderate approximation capability.

Consequently, the effectiveness of a basis cannot be assessed from approximation error alone. For noisy first-kind Fredholm equations, the relevant performance measure is the combined behavior of approximation power, conditioning, discrepancy-selected regularization, residual control, noise robustness, and computational cost. This explains why a basis that is accurate in a clean-data benchmark may not be the most reliable choice once perturbations and regularization are introduced. The results therefore support the central claim of the paper: basis selection is part of the stabilization strategy for the inverse problem, not merely a representation choice.

From a practical viewpoint, the computations suggest the following guideline. When the unknown solution is expected to be smooth and high reconstruction accuracy is required, Legendre, B-spline, and Linear Legendre multiwavelet discretizations provide the most favorable balance between approximation quality and regularized recovery. When computational simplicity, low implementation cost, locality, or reduced regularity are dominant considerations, Haar and hybrid discretizations remain competitive alternatives. Thus, the preferred basis should be chosen jointly with the regularization strategy and with the expected smoothness and noise level of the data in mind.


\subsection*{Practical recommendations}
The diagnostic results lead to a more precise basis-selection rule than a standard accuracy table. The numerical evidence suggests that no single basis is uniformly optimal. For smooth solutions and moderate noise levels, Linear Legendre multiwavelets, Legendre wavelets, and B-spline wavelets generally provide the best accuracy because their higher-order local structure represents smooth variation efficiently after stabilization. For nonsmooth solutions, discontinuities, sharp localized features, or computationally inexpensive implementations, Haar and hybrid wavelets remain attractive alternatives because of their locality, simplicity, and low assembly cost. Therefore, basis selection should be guided not only by approximation order but also by the expected regularity of the solution, the noise level, the conditioning behavior of the discrete operator, and the available computational budget.

\begin{table}[htbp]
\centering
\scriptsize
\caption{Practical synthesis of the basis-dependent numerical evidence. The qualitative ratings summarize the observed behavior across accuracy, noise robustness, conditioning, and computational cost under the common regularized collocation protocol.}
\label{tab:practical-basis-summary}
\begin{tabularx}{\textwidth}{p{0.12\textwidth} p{0.13\textwidth} p{0.15\textwidth} p{0.14\textwidth} p{0.12\textwidth} X}
\toprule
Basis & Accuracy & Noise robustness & Conditioning & Cost & Recommended use \\
\midrule
LLMW & Excellent & Excellent & Very good & Moderate & Smooth inverse problems where stability and reconstruction accuracy are primary. \\
Legendre & Excellent & Very good & Very good & Moderate & Smooth recoveries requiring high-order local approximation. \\
B-spline & Excellent & Very good & Good & Moderate & Smooth reconstructions with local regularity and stable approximation. \\
Haar & Moderate & Good & Good & Very low & Nonsmooth targets, exploratory computations, and simple low-cost implementations. \\
Hybrid & Good & Good & Good & Low & Balanced accuracy--cost performance and moderately regular solutions. \\
\bottomrule
\end{tabularx}
\end{table}

Overall, the Gaussian smoothing experiment supports the main novelty claim of the paper. Recent studies have shown that wavelet bases can solve first-kind Fredholm equations and that basis choice affects approximation accuracy. The additional contribution here is to compare several compactly supported bases after the same Tikhonov--Morozov stabilization and to report the quantities that determine whether the inverse recovery is actually reliable: solution error, residual error, selected $\alpha$, condition-number reduction, L-curve behavior, noise robustness, and computational cost. In particular, two bases may give similar clean-data errors but very different conditioning, selected regularization parameters, or sensitivity to perturbations. This diagnostic layer is what distinguishes the present work from a purely approximation-based wavelet comparison and explains why basis selection for first-kind Fredholm equations must be treated as a regularized inverse-problem decision.

The observed behavior is consistent with the approximation properties of the underlying bases. Higher-order local bases such as Linear Legendre multiwavelets, Legendre wavelets, and B-spline wavelets are better able to represent smooth solution structures and therefore achieve higher reconstruction accuracy for a fixed discretization level. In contrast, Haar wavelets possess only piecewise-constant structure and consequently require additional refinement to resolve smooth features accurately. Nevertheless, their compact support, locality, and extremely low computational cost make them competitive for nonsmooth targets and large-scale exploratory computations.

An important observation is that basis selection influences not only the accuracy but also the conditioning of the discrete inverse problem. Bases producing smaller values of $\kappa(A_N)$ generally require less aggressive regularization and lead to more stable reconstructions. Conversely, highly ill-conditioned discretizations may achieve accurate noise-free approximations while becoming substantially more sensitive to data perturbations. This explains why basis rankings based solely on deterministic approximation error can differ significantly from those obtained under noisy inverse recovery.

Although all methods approximate the same continuous compact operator, their finite-dimensional representations possess different spectral distributions and conditioning properties. Since Tikhonov regularization acts directly on these discrete systems, the resulting regularized solutions may exhibit substantially different stability, discrepancy behavior, and noise-propagation characteristics. This explains why the interaction between basis construction and regularization must be evaluated at the discrete inverse-problem level, rather than inferred only from approximation order.

The numerical experiments are designed not only to compare approximation accuracy but also to investigate how basis selection affects conditioning, regularization strength, residual behavior, robustness to measurement noise, and computational efficiency within a common Tikhonov–Morozov framework.


\section{Conclusions}
\label{sec:conclusion}

This work has presented a unified regularized multiwavelet-collocation framework for ill-posed Fredholm integral equations of the first kind. Linear Legendre multiwavelets, Legendre wavelets, Haar wavelets, B-spline wavelets, and hybrid wavelets were tested under the same collocation, quadrature, Tikhonov regularization, Morozov parameter-selection, and error-measurement protocol. This common setting makes the comparison reproducible and allows the observed differences to be attributed to basis-dependent discretization and stabilization effects rather than to unrelated implementation choices.

The principal conclusion is that the approximation basis affects more than the final reconstruction error. It changes the conditioning of the discrete operator, the conditioning of the regularized normal equations, the Morozov-selected value of $\alpha$, the residual behavior, the propagation of data perturbations, and the computational cost. Thus, for first-kind Fredholm equations, basis selection should be treated as part of the regularization strategy rather than as a purely approximation-theoretic choice.

The numerical evidence gives a clear practical guideline. For smooth unknown functions and high-accuracy reconstruction, Legendre, B-spline, and Linear Legendre multiwavelet discretizations offer the most favorable balance between approximation quality and regularized recovery. When implementation simplicity, low computational cost, locality, or reduced regularity are more important, Haar and hybrid bases remain useful alternatives. The best choice is therefore problem dependent and should be made by balancing accuracy, conditioning, robustness to noise, residual control, and computational expense.

The reported computations use the stated basis sizes, noise levels, Gaussian kernel width, Morozov tolerance, and a fixed random seed for the noisy-data experiments. The manuscript package records the diagnostic code, the noise levels $10^{-4}$, $10^{-3}$, $10^{-2}$, and $5\times 10^{-2}$, the basis size $N=32$, the random seed 20250818, and the scripts used to generate the diagnostic tables and figures. This information is included to make the comparison transparent and to allow future studies to extend the same protocol to other kernels, bases, or parameter-selection rules.

The present study is restricted to compact intervals, smooth kernels, and moderate basis dimensions. Extensions to adaptive refinement strategies, higher-dimensional inverse problems, nonlinear or weakly singular integral equations, and non-Gaussian noise models constitute promising directions for future research. Overall, the results demonstrate that approximation-space selection should be regarded as an integral component of regularization design, rather than as a purely discretization-related choice, in the numerical solution of ill-posed Fredholm inverse problems.

The present study is restricted to compact intervals, smooth kernels, and moderate basis dimensions. Extensions to adaptive refinement strategies, higher-dimensional inverse problems, nonlinear integral equations, and non-Gaussian noise models constitute promising directions for future research.



% Bibliography note for citation key sathar2018numerical:
% Mohammad Hasan Abdul Sathar, Ahmad Fadly Nurullah Rasedee, Anvarjon A. Ahmedov,
% and Muhammad Asyraf Asbullah, ``Numerical Integration Based on Linear Legendre
% Multi Wavelets,'' Journal of Physics: Conference Series, 1132, 012008, 2018.
% doi:10.1088/1742-6596/1132/1/012008.
%\bibliographystyle{elsarticle-num}
%\bibliography{sample,fredholm_recent_refs_final}
%\bibliography{sample,fredholm_recent_refs_final,fredholm_reference_patch.bib}

%\bibliography{sample,fredholm_recent_refs_final,fredholm_reference_patch}
%\bibliography{sample,sample,sample.bib}
%\bibliography{sample,fredholm_recent_refs_final}
\subsection*{Declarations}

\textbf{Funding.} No funding was received for this study.

\textbf{Conflict of interest.} The author declares that they have no conflict of interest.

\textbf{Data availability.} During the current study, no experimental data were generated or analyzed. The figures were obtained from the analytical expressions reported in the manuscript.

\textbf{Code availability.} The computational scripts used to generate the figures are available from the corresponding author upon reasonable request.

\textbf{Author contributions.} Ahmed Kaffel developed the formulation, carried out the analysis, prepared the figures, interpreted the results, and wrote the manuscript.

\bibliographystyle{elsarticle-num}
\bibliography{sample,fredholm_recent_refs_final}
%\section*{Funding}  Funding not applicable
\end{document}
\end{document}

%%
%% End of file `elsarticle-template-num.tex'.
