Reduction to a Fredholm integral equation and numerical solution of the inverse Cauchy problem for the Schrödinger-Pauli equation

Yusif Gasimov; Abdeljalil Nachaoui; Aynura Aliyeva

doi:10.2478/ijmce-2026-0015

Introduction

The Pauli equation [1, 2] plays a fundamental role in plasma physics and in the description of atomic and subatomic particle interactions. The associated Pauli operator, introduced in quantum mechanics [2], models the dynamics of a quantum particle with spin subjected to a magnetic field. Depending on the electromagnetic and external fields, the resulting differential equation is linear and may be time-dependent. In unbounded domains or infinite strips, the equation can be solved numerically with relative ease. A wide variety of numerical approaches have been proposed, including finite difference methods [3–7], polynomial-based schemes [8,9], pseudo-spectral methods [10–13], and neural-network-based techniques [14].

More recently, the Pauli equation posed with boundary conditions has been investigated in [15], where the authors proposed a trigonometric series representation of the solution. To the best of our knowledge, however, the Cauchy problem associated with this equation has not yet been studied. The present work aims to fill this gap by formulating an algorithm for solving the corresponding inverse problem and by numerically validating its convergence.

We focus on the simplest nontrivial physical setting, namely an electron interacting solely with an external homogeneous magnetic field $\vec{B} \in ℝ^{2}$ \vec B \in {^2}. The system is governed by the Pauli equation 1 $P Ψ = f,$ P\Psi = f, where $P = [{(- i \nabla - a)}^{2} + V] \cdot J + σ B .$ P = \,\,\left[ {{{( - i\nabla - a)}^2} + V} \right]\,\,\cdot\,J + \sigma B.

A convenient way to “complexify” (1) is to restrict the spatial variables to a bounded domain Ω ⊂ ℝ² and to impose complex-valued Robin-type boundary conditions, 2 $α \frac{\partial Ψ}{\partial n} + A Ψ = 0 on \partial Ω,$ \alpha {{\partial \Psi } \over {\partial n}} + A\Psi = 0\quad {\rm{ on }}\quad \partial \Omega , where n denotes the outward unit normal to ∂Ω, α ∈ ℝ, and A is a 2 × 2 complex-valued matrix. The solution Ψ of (1) is therefore required to satisfy the boundary constraint (2).

When B = 0, equation (1) reduces to the classical Schrödinger equation, which has been extensively studied over the past century in the context of atomic structure and particle interactions (see, e.g., [1,2,16]). In contrast, much less attention has been paid to the Pauli equation posed on bounded domains. Although the problem remains linear, the imposed boundary conditions introduce significant analytical and numerical challenges that limit the applicability of standard methods.

For Dirichlet boundary conditions (α = 0 and A = I), the authors of [15] constructed a trigonometric series solution to the Pauli equation and established both existence and uniqueness results. In the present work, we address a more challenging situation in which the boundary matrix A is known only on an accessible portion Γ₁ of the boundary ∂Ω and must be reconstructed on an inaccessible part Γ₀ using measurements of Ψ on Γ₁. This leads to an inverse Cauchy problem for the Pauli equation.

Inverse problems of this type arise in numerous practical applications, where unknown initial or boundary data, or inaccessible parameters, must be identified from indirect measurements [17–28]. Such problems are typically ill-posed in the sense of Hadamard [29], meaning that small perturbations in the data may produce large errors in the solution. This motivates the development of stable and robust numerical techniques [20,23,30–48].

A large body of research has focused on parameter identification problems. For instance, [49] studied the recovery of time-independent coefficients in a fractional diffusion model from final-time data, while [50] and [22] addressed hydraulic parameter estimation in porous media and Richards equation, respectively. In heat transfer, the identification of time-dependent thermal conductivity was investigated in [51].

Nonlocal and inverse source problems have also been widely examined, with applications ranging from groundwater dynamics [52, 53] to viscoelasticity and food processing [54]. Inverse Cauchy problems constitute another important class, particularly for elliptic equations. Iterative and relaxed methods introduced by Nachaoui [30] have been extended to nonlinear elliptic problems [33], elasticity [55], convection–diffusion equations [56], biharmonic equations [57], and hyperbolic models such as the telegraph equation [58]. Boundary element and meshless approaches have also been successfully employed [27, 34, 47, 59–65].

The remainder of this paper is organized as follows. Section 2 presents preliminary material related to the Pauli problem. The formulation of the inverse Cauchy problem is developed in Section 3. The main contribution, namely the reduction of this problem to a Fredholm equation and its approximation by Haar wavelets, is found in Section 4. Section 5 is devoted to numerical experiments that illustrate the performance and convergence of the proposed method. Finally, Section 6 concludes the paper.

Structure of the two-dimensional Pauli operator

The Pauli equation extends the classical Schrödinger equation by incorporating the spin of a quantum particle. In particular, it describes the dynamics of a spin- $\frac{1}{2}$ {1 \over 2} particle subjected to an external magnetic field, which constitutes its main distinction from the standard Schrödinger operator (see, e.g., [1, 2,16]).

In two spatial dimensions, the Schrödinger-Pauli operator is defined as $P = P (a, V) \cdot J + σ B,$ P = P(a,V)\cdotJ + \sigma B, where $J = (\begin{array}{l} 1 & 0 \\ 0 & 1 \end{array}), σ = (\begin{matrix} 1 & 0 \\ 0 & - 1 \end{matrix})$ J = \left( {\matrix{ 1 \hfill & 0 \hfill \cr 0 \hfill & 1 \hfill \cr } } \right),\quad \sigma = \left( {\matrix{ 1 & 0 \cr 0 & { - 1} \cr } } \right) are diagonal matrices, and $P (a, V) = {(- i \nabla - a)}^{2} + V, \nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial y}) .$ P(a,V) = {( - i\nabla - a)^2} + V,\quad \nabla = \left( {{\partial \over {\partial x}},{\partial \over {\partial y}}} \right).

Here, i denotes the imaginary unit, a = (a₁,a₂) is the magnetic vector potential, and V is the electric potential. The magnetic field B, orthogonal to the plane, is generated by the vector potential a and is given by $B = B_{3} = a_{2} \frac{\partial}{\partial x} - a_{1} \frac{\partial}{\partial y} .$ B = {B_3} = {a_2}{\partial \over {\partial x}} - {a_1}{\partial \over {\partial y}}.

A direct computation shows that the Pauli operator admits a diagonal representation, 3 $P = (\begin{matrix} {(- i \nabla - a)}^{2} + a_{2} \frac{\partial}{\partial x} - a_{1} \frac{\partial}{\partial y} + V & 0 \\ 0 & {(- i \nabla - a)}^{2} - a_{2} \frac{\partial}{\partial x} + a_{1} \frac{\partial}{\partial y} + V \end{matrix}) = (\begin{matrix} P_{1} & 0 \\ 0 & P_{2} \end{matrix}),$ P = \left( {\matrix{ {{{( - i\nabla - a)}^2} + {a_2}{\partial \over {\partial x}} - {a_1}{\partial \over {\partial y}} + V} & 0 \cr 0 & {{{( - i\nabla - a)}^2} - {a_2}{\partial \over {\partial x}} + {a_1}{\partial \over {\partial y}} + V} \cr } } \right) = \left( {\matrix{ {{P_1}} & 0 \cr 0 & {{P_2}} \cr } } \right), where 4 $P_{1} = - Δ + (2 i a_{1} + a_{2}) \frac{\partial}{\partial x} + (2 i a_{2} - a_{1}) \frac{\partial}{\partial y} + a^{2} + V,$ {P_1} = - \Delta + \left( {2i{a_1} + {a_2}} \right){\partial \over {\partial x}} + \left( {2i{a_2} - {a_1}} \right){\partial \over {\partial y}} + {a^2} + V, 5 $P_{2} = - Δ + (2 i a_{1} - a_{2}) \frac{\partial}{\partial x} + (2 i a_{2} + a_{1}) \frac{\partial}{\partial y} + a^{2} + V .$ {P_2} = - \Delta + \left( {2i{a_1} - {a_2}} \right){\partial \over {\partial x}} + \left( {2i{a_2} + {a_1}} \right){\partial \over {\partial y}} + {a^2} + V.

In what follows, we assume without loss of generality that V = 0 and L = 1. We consider equation (1) posed in the domain Ω and seek the unknown spinor Ψ = (Ψ₁, Ψ₂) corresponding to a given source term f = (f₁, f₂). Using (3)–(5), the Schrödinger-Pauli equation decouples into the system 6 ${\begin{array}{l} - {ΔΨ}_{1} + (2 i a_{1} + a_{2}) \frac{\partial Ψ_{1}}{\partial x} + (2 i a_{2} - a_{1}) \frac{\partial Ψ_{1}}{\partial y} + a^{2} Ψ_{1} = f_{1}, (x, y) \in Ω, \\ - {ΔΨ}_{2} + (2 i a_{1} - a_{2}) \frac{\partial Ψ_{2}}{\partial x} + (2 i a_{2} + a_{1}) \frac{\partial Ψ_{2}}{\partial y} + a^{2} Ψ_{2} = f_{2}, (x, y) \in Ω . \end{array}$ \left\{ {\matrix{ { - \Delta {\Psi _1} + \left( {2i{a_1} + {a_2}} \right){{\partial {\Psi _1}} \over {\partial x}} + \left( {2i{a_2} - {a_1}} \right){{\partial {\Psi _1}} \over {\partial y}} + {a^2}{\Psi _1} = {f_1},\,\,\,\,(x,y) \in \Omega \,,} \hfill \cr { - \Delta {\Psi _2} + \left( {2i{a_1} - {a_2}} \right){{\partial {\Psi _2}} \over {\partial x}} + \left( {2i{a_2} + {a_1}} \right){{\partial {\Psi _2}} \over {\partial y}} + {a^2}{\Psi _2} = {f_2},(x,y) \in \Omega .} \hfill \cr } } \right.

Separating the real and imaginary parts yields the equivalent system 7 ${\begin{array}{l} - Δ Ψ_{1} + a_{2} \frac{\partial Ψ_{1}}{\partial x} - a_{1} \frac{\partial Ψ_{1}}{\partial y} + a^{2} Ψ_{1} = f_{1}, \\ a_{1} \frac{\partial Ψ_{1}}{\partial x} + a_{2} \frac{\partial Ψ_{1}}{\partial y} = 0, \\ - Δ Ψ_{2} - a_{2} \frac{\partial Ψ_{2}}{\partial x} + a_{1} \frac{\partial Ψ_{2}}{\partial y} + a^{2} Ψ_{2} = f_{2}, \\ a_{1} \frac{\partial Ψ_{2}}{\partial x} + a_{2} \frac{\partial Ψ_{2}}{\partial y} = 0, \end{array}$ \left\{ {\matrix{ { - \Delta {\Psi _1} + {a_2}{{\partial {\Psi _1}} \over {\partial x}} - {a_1}{{\partial {\Psi _1}} \over {\partial y}} + {a^2}{\Psi _1} = {f_1},} \hfill \cr {{a_1}{{\partial {\Psi _1}} \over {\partial x}} + {a_2}{{\partial {\Psi _1}} \over {\partial y}} = 0,} \hfill \cr { - \Delta {\Psi _2} - {a_2}{{\partial {\Psi _2}} \over {\partial x}} + {a_1}{{\partial {\Psi _2}} \over {\partial y}} + {a^2}{\Psi _2} = {f_2},} \hfill \cr {{a_1}{{\partial {\Psi _2}} \over {\partial x}} + {a_2}{{\partial {\Psi _2}} \over {\partial y}} = 0,} \hfill \cr } } \right. to be solved together with the boundary conditions (2).

Assuming a₁ ≠ 0 and a₂ ≠ 0, the constraint equations in (7) imply 8 ${\begin{array}{l} \frac{\partial Ψ_{1}}{\partial y} = - \frac{a_{1}}{a_{2}} \frac{\partial Ψ_{1}}{\partial x}, \\ \frac{\partial Ψ_{2}}{\partial x} = - \frac{a_{2}}{a_{1}} \frac{\partial Ψ_{2}}{\partial y}, \end{array}$ \left\{ {\matrix{ {{{\partial {\Psi _1}} \over {\partial y}} = - {{{a_1}} \over {{a_2}}}{{\partial {\Psi _1}} \over {\partial x}},} \hfill \cr {{{\partial {\Psi _2}} \over {\partial x}} = - {{{a_2}} \over {{a_1}}}{{\partial {\Psi _2}} \over {\partial y}},} \hfill \cr } } \right. which allow further simplification of the governing equations. Substituting (8) into (7)_1,3 leads to 9 ${\begin{array}{l} - Δ Ψ_{1} + \frac{a^{2}}{a_{2}} \frac{\partial Ψ_{1}}{\partial x} + a^{2} Ψ_{1} = f_{1}, \\ - Δ Ψ_{2} + \frac{a^{2}}{a_{1}} \frac{\partial Ψ_{2}}{\partial y} + a^{2} Ψ_{2} = f_{2}, \end{array}$ \left\{ {\matrix{ { - \Delta {\Psi _1} + {{{a^2}} \over {{a_2}}}{{\partial {\Psi _1}} \over {\partial x}} + {a^2}{\Psi _1} = {f_1},} \hfill \cr { - \Delta {\Psi _2} + {{{a^2}} \over {{a_1}}}{{\partial {\Psi _2}} \over {\partial y}} + {a^2}{\Psi _2} = {f_2},} \hfill \cr } } \right. which provides a simplified formulation suitable for the analysis and numerical treatment developed in the subsequent sections.

Formulation of the inverse Cauchy problem

Let Ω be the unit square $Ω = (0, 1) \times (0, 1) \subset ℝ^{2} .$ \Omega = (0,1)\,\, \times \,\,(0,1) \subset {^2}.

In accordance with the boundary condition (2), we impose homogeneous Dirichlet conditions on three sides of the boundary, 10 ${\begin{array}{l} Ψ_{1} (x, 0) = Ψ_{1} (0, y) = Ψ_{1} (1, y) = 0, \\ Ψ_{2} (x, 0) = Ψ_{2} (0, y) = Ψ_{2} (1, y) = 0, \end{array}$ \left\{ {\matrix{ {{\Psi _1}(x,0) = {\Psi _1}(0,y) = {\Psi _1}(1,y) = 0,} \hfill \cr {{\Psi _2}(x,0) = {\Psi _2}(0,y) = {\Psi _2}(1,y) = 0,} \hfill \cr } } \right. together with the overdetermination conditions 11 $- \frac{\partial Ψ_{1}}{\partial y} (x, 0) = g_{1} (x),$ - {{\partial {\Psi _1}} \over {\partial y}}(x,0) = {g_1}(x), 12 $- \frac{\partial Ψ_{2}}{\partial y} (x, 0) = g_{2} (x) .$ - {{\partial {\Psi _2}} \over {\partial y}}(x,0) = {g_2}(x).

Due to the symmetry of the boundary conditions and the equivalence of the governing equations in (9), it is sufficient to restrict the analysis to a single scalar component. We therefore consider the following inverse Cauchy problem: determine Ψ₁ and the regularization parameter α such that 13 ${\begin{array}{l} - Δ Ψ_{1} + \frac{a^{2}}{a_{2}} \frac{\partial Ψ_{1}}{\partial x} + a^{2} Ψ_{1} = f_{1}, & (x, y) \in Ω, \\ Ψ_{1} (x, 0) = Ψ_{1} (0, y) = Ψ_{1} (1, y) = 0, \\ - \frac{\partial Ψ_{1}}{\partial y} (x, 0) = g_{1} (x), \\ \frac{\partial Ψ_{1}}{\partial y} (x, 1) + α Ψ_{1} (x, 1) = 0. \end{array}$ \left\{ {\matrix{ { - \Delta {\Psi _1} + {{{a^2}} \over {{a_2}}}{{\partial {\Psi _1}} \over {\partial x}} + {a^2}{\Psi _1} = {f_1},} \hfill & {(x,y) \in \Omega ,} \hfill \cr {{\Psi _1}(x,0) = {\Psi _1}(0,y) = {\Psi _1}(1,y) = 0,} \hfill & {} \hfill \cr { - {{\partial {\Psi _1}} \over {\partial y}}(x,0) = {g_1}(x),} \hfill & {} \hfill \cr {{{\partial {\Psi _1}} \over {\partial y}}(x,1) + \alpha {\Psi _1}(x,1) = 0.} \hfill & {} \hfill \cr } } \right.

The main idea is to express Ψ₁(x, y) analytically in terms of the unknown trace $φ (x) = Ψ_{1} (x, 1),$ \varphi (x) = {\Psi _1}(x,1), by means of a separation-of-variables (Fourier) expansion. Differentiation with respect to y and evaluation at y = 0, then, yields a Fredholm integral equation of the first kind for φ, with the measured data g₁(x) appearing on the right-hand side. Once φ has been reconstructed, it can be substituted back into the analytical representation to recover Ψ₁(x, y) throughout the domain.

Derivation of the Fredholm integral equation

We first consider the boundary value problem 14 ${\begin{array}{l} - Δ Ψ_{1} + \frac{a^{2}}{a_{2}} \frac{\partial Ψ_{1}}{\partial x} + a^{2} Ψ_{1} = f_{1}, & (x, y) \in {(0, 1)}^{2}, \\ Ψ_{1} (x, 0) = Ψ_{1} (0, y) = Ψ_{1} (1, y) = 0, \\ Ψ_{1} (x, 1) = φ (x), \end{array}$ \left\{ {\matrix{ { - \Delta {\Psi _1} + {{{a^2}} \over {{a_2}}}{{\partial {\Psi _1}} \over {\partial x}} + {a^2}{\Psi _1} = {f_1},} \hfill & {(x,y) \in {{(0,1)}^2},} \hfill \cr {{\Psi _1}(x,0) = {\Psi _1}(0,y) = {\Psi _1}(1,y) = 0,} \hfill & {} \hfill \cr {{\Psi _1}(x,1) = \varphi (x),} \hfill & {} \hfill \cr } } \right. supplemented with the additional condition $- \frac{\partial Ψ_{1}}{\partial y} (x, 0) = g_{1} (x) .$ - {{\partial {\Psi _1}} \over {\partial y}}(x,0) = {g_1}(x).

We begin with the homogeneous case f₁ = 0 and introduce the notation $c = \frac{a^{2}}{a_{2}} .$ c = {{{a^2}} \over {{a_2}}}.

Seeking separable solutions of the form Ψ₁(x, y) = X(x)Y(y) leads to $\frac{X^{''} - c X^{'}}{X} = - λ, \frac{Y^{''} - a^{2} Y}{Y} = λ,$ {{{X^{\prime \prime }} - c{X^\prime }} \over X} = - \lambda ,\quad {{{Y^{\prime \prime }} - {a^2}Y} \over Y} = \lambda , and hence to the system 15 ${\begin{matrix} X^{''} - c X^{'} + λ X & = & 0, \\ Y^{''} - (λ + a^{2}) Y & = & 0 . \end{matrix}$ \left\{ {\matrix{ {{X^{\prime \prime }} - c{X^\prime } + \lambda X = 0,} \hfill \cr {{Y^{\prime \prime }} - \left( {\lambda + {a^2}} \right)Y = 0.} \hfill \cr } } \right.

Now, let’s consider eigenvalue problem in the x-direction.

The equation $X^{''} - c X^{'} + λ X = 0,$ {X^{\prime \prime }} - c{X^\prime } + \lambda X = 0, subject to the Dirichlet conditions X(0) = X(1) = 0, can be rewritten in self-adjoint form by introducing the transformation $X (x) = e^{c x / 2} ϕ (x) .$ X(x) = {e^{cx/2}}\phi (x).

A direct computation yields $ϕ^{''} + (λ - \frac{c^{2}}{4}) ϕ = 0, ϕ (0) = ϕ (1) = 0.$ {\phi ^{\prime \prime }} + \left( {\lambda - {{{c^2}} \over 4}} \right)\phi = 0,\quad \phi (0) = \phi (1) = 0.

This classical Sturm–Liouville problem admits the eigenpairs $λ_{n} = (^{n π) 2} + \frac{c^{2}}{4}, X_{n} (x) = e^{c x / 2} \sin (n π x), n \geq 1.$ {\lambda _n} = {(n\pi )^2} + {{{c^2}} \over 4},\quad {X_n}(x) = {e^{cx/2}}\sin (n\pi x),\quad n \ge 1.

Now let’s consider the solution in the y-direction. For this purpose for each n ≥ 1, define 16 $γ_{n}^{2} = (^{n π) 2} + \frac{c^{2}}{4} + a^{2} > 0.$ \gamma _n^2 = {(n\pi )^2} + {{{c^2}} \over 4} + {a^2} > 0.

The corresponding y-dependent equation $Y^{''} - γ_{n}^{2} Y = 0$ {Y^{\prime \prime }} - \gamma _n^2Y = 0 with the condition Y(0) = 0 admits the solution $Y_{n} (y) = \sinh (γ_{n} y) .$ {Y_n}(y) = \sinh \left( {{\gamma _n}y} \right).

By superposition, the solution of (14) can be written as $Ψ_{1} (x, y) = \sum_{n = 1}^{\infty} A_{n} e^{c x / 2} \sin (n π x) \sinh (γ_{n} y) .$ {\Psi _1}(x,y) = \sum\limits_{n = 1}^\infty {{A_n}} {e^{cx/2}}\sin (n\pi x)\sinh \left( {{\gamma _n}y} \right).

Imposing the boundary condition Ψ₁(x, 1) = φ(x) yields $φ (x) = \sum_{n = 1}^{\infty} A_{n} e^{c x / 2} \sin (n π x) \sinh (γ_{n}) .$ \varphi (x) = \sum\limits_{n = 1}^\infty {{A_n}} {e^{cx/2}}\sin (n\pi x)\sinh \left( {{\gamma _n}} \right).

Multiply by e^−cx/2 sin(πx), integrate on [0, 1] and using orthogonality of the sine functions, we obtain $A_{n} = \frac{2}{\sinh (γ_{n})} \int_{0}^{1} φ (s) e^{- c s / 2} \sin (n π s) d s .$ {A_n} = {2 \over {\sinh \left( {{\gamma _n}} \right)}}\int_0^1 \varphi (s){e^{ - cs/2}}\sin (n\pi s)ds.

Hence, 17 $Ψ_{1} (x, y) = \sum_{n = 1}^{\infty} \frac{2 \sinh (γ_{n} y)}{\sinh (γ_{n})} (\int_{0}^{1} φ (s) e^{- c s / 2} \sin (n π s) d s) e^{c x / 2} \sin (n π x) .$ {\Psi _1}(x,y) = \sum\limits_{n = 1}^\infty {{{2\sinh \left( {{\gamma _n}y} \right)} \over {\sinh \left( {{\gamma _n}} \right)}}} \left( {\int_0^1 \varphi (s){e^{ - cs/2}}\sin (n\pi s)ds} \right){e^{cx/2}}\sin (n\pi x).

Now we can construct the Fredholm integral equation. Differentiating (17) with respect to y and evaluating at y = 0, the overdetermination condition yields $g_{1} (x) = - \sum_{n = 1}^{\infty} \frac{2 γ_{n}}{\sinh (γ_{n})} e^{c x / 2} \sin (n π x) \int_{0}^{1} φ (s) e^{- c s / 2} \sin (n π s) d s .$ {g_1}(x) = - \sum\limits_{n = 1}^\infty {{{2{\gamma _n}} \over {\sinh \left( {{\gamma _n}} \right)}}} {e^{cx/2}}\sin (n\pi x)\int_0^1 \varphi (s){e^{ - cs/2}}\sin (n\pi s)ds.

Interchanging summation and integration leads to the Fredholm integral equation of the first kind 18 $g_{1} (x) = \int_{0}^{1} K (x, s) φ (s) d s,$ {g_1}(x) = \int_0^1 K (x,s)\varphi (s)ds, with kernel $K (x, s) = - 2 e^{\frac{c}{2} (x - s)} \sum_{n = 1}^{\infty} \frac{γ_{n}}{\sinh (γ_{n})} \sin (n π x) \sin (n π s) .$ K(x,s) = - 2{e^{{c \over 2}(x - s)}}\sum\limits_{n = 1}^\infty {{{{\gamma _n}} \over {\sinh \left( {{\gamma _n}} \right)}}} \sin (n\pi x)\sin (n\pi s).

This operator is compact and severely ill-posed.

To stabilize the inversion, we adopt Lavrentiev regularization and consider the second-kind Fredholm equation 19 $α φ_{α} (x) + \int_{0}^{1} K (x, t) φ_{α} (t) d t = g_{1} (x), x \in [0, 1],$ \alpha {\varphi _\alpha }(x) + \int_0^1 K (x,t){\varphi _\alpha }(t)dt = {g_1}(x),\quad x \in [0,1], where α > 0 is the regularization parameter.

4.1

Truncation as regularization

The kernel of the Fredholm integral equation is given by $K (x, s) = - 2 e^{\frac{c}{2} (x - s)} \sum_{n = 1}^{\infty} β_{n} \sin (n π x) \sin (n π s), β_{n} = \frac{γ_{n}}{\sinh (γ_{n})} .$ K(x,s) = - 2{e^{{c \over 2}(x - s)}}\sum\limits_{n = 1}^\infty {{\beta _n}} \sin (n\pi x)\sin (n\pi s),\quad {\beta _n} = {{{\gamma _n}} \over {\sinh \left( {{\gamma _n}} \right)}}.

We will use another type of regularization which truncates this series; this is justified by the rapid decrease of the β_n coefficients. Indeed, for large n, γ_n ~ nπ and $\sinh (γ_{n}) ~ \frac{1}{2} e^{γ_{n}}$ \sinh \left( {{\gamma _n}} \right)\~{1 \over 2}{e^{{\gamma _n}}}, yielding $β_{n} ~ 2 n π e^{- n π},$ {\beta _n}\~2n\pi {e^{ - n\pi }}, demonstrating exponential decay as e^−nπ. We show in the numerical results section that β_n decreases towards zero after a few terms.

Therefore, we will use the truncated kernel after N terms: $K_{N} (x, t) = - 2 e^{\frac{c}{2} (x - t)} \sum_{n = 1}^{N} β_{n} \sin (n π x) \sin (n π t) .$ {K_N}(x,t) = - 2{e^{{c \over 2}(x - t)}}\sum\limits_{n = 1}^N {{\beta _n}} \sin (n\pi x)\sin (n\pi t).

Note that Truncating at N

Introduces negligible approximation error due to exponential decay.
Acts as spectral cutoff regularization, filtering high-frequency noise.
Replaces the ill-posed first-kind equation with a finite-rank approximation that can be solved stably after applying Lavrentiev regularization.

Thus, the truncated kernel K_N(x, s) is both accurate and computationally efficient, with the exponential decay of β_n providing rigorous justification for the truncation.

The regularized equation (19) is then solved numerically by replacing K with K_N and using a Haar wavelet approximation.

4.2

The Haar wavelets approximation

Haar wavelets, introduced by Haar [66], form a simple orthogonal basis of piecewise constant functions. Owing to their compact support and computational efficiency, they have been widely employed in the numerical solutions of ordinary and partial differential equations [67–74]. In this work, the unknown boundary function is expanded in a Haar series and the resulting integral equation is discretized by a collocation strategy.

In the Haar wavelet approximation, the dimension of the basis is determined by the maximum resolution level (the finest scale of detail that can be captured with the basis). If the maximum resolution level is denoted J, the total number of basis functions (the dimension) is 2M with M = 2^J.

On an interval [a, b], the Haar basis consists of the scaling function $h_{1} (x) = {\begin{array}{l} 1, & x \in [a, b), \\ 0, & otherwise, \end{array}$ {h_1}(x) = \left\{ {\matrix{ {1,} \hfill & {x \in [a,b),} \hfill \cr {0,} \hfill & {{\rm{ otherwise,}}} \hfill \cr } } \right. the mother wavelet $h_{2} (x) = {\begin{array}{l} 1, & x \in [a, \frac{a + b}{2}), \\ - 1, & x \in [\frac{a + b}{2}, b), \\ 0, & otherwise, \end{array}$ {h_2}(x) = \left\{ {\matrix{ {1,} \hfill & {x \in \left[ {a,{{a + b} \over 2}} \right),} \hfill \cr { - 1,} \hfill & {x \in \left[ {{{a + b} \over 2},b} \right),} \hfill \cr {0,} \hfill & {{\rm{ otherwise,}}} \hfill \cr } } \right. and the wavelets $h_{i} (x) = {\begin{array}{l} 1, & x \in [ξ_{1} (i), ξ_{2} (i)), \\ - 1, & x \in [ξ_{2} (i), ξ_{3} (i)), & i = m + k + 1, \\ 0, & otherwise, \end{array}$ {h_i}(x) = \left\{ {\matrix{ {1,} \hfill & {x \in \left[ {{\xi _1}(i),{\xi _2}(i)} \right),} \hfill & {} \hfill \cr { - 1,} \hfill & {x \in \left[ {{\xi _2}(i),{\xi _3}(i)} \right),} \hfill & {i = m + k + 1,} \hfill \cr {0,} \hfill & {{\rm{ otherwise,}}} \hfill & {} \hfill \cr } } \right. where j = 0, ⋯, J, m = 2^j, k = 0, ⋯, m – 1, and $ξ_{1} (i) = a + \frac{k (b - a)}{m}, ξ_{2} (i) = a + \frac{(k + 0.5) (b - a)}{m}, ξ_{3} (i) = a + \frac{(k + 1) (b - a)}{m} .$ {\xi _1}(i) = a + {{k(b - a)} \over m},\quad {\xi _2}(i) = a + {{(k + 0.5)(b - a)} \over m},\quad {\xi _3}(i) = a + {{(k + 1)(b - a)} \over m}.

The interval [a, b] is partitioned into 2M equal subintervals of length Δx = (b – a)/(2M).

Any function u ∈ L²([a, b]) can be approximated by the truncated Haar series 20 $u (x) \approx \sum_{i = 1}^{2 M} u_{i} h_{i} (x), u_{i} = \int_{a}^{b} u (x) h_{i} (x) d x .$ u(x) \approx \sum\limits_{i = 1}^{2M} {{u_i}} {h_i}(x),\quad {u_i} = \int_a^b u (x){h_i}(x)dx.

Setting a = 0 and b = 1, we approximate the solution of the regularized Fredholm equation (19) by $φ_{α} (t) \approx \sum_{i = 1}^{2 M} c_{i} h_{i} (t),$ {\varphi _\alpha }(t) \approx \sum\limits_{i = 1}^{2M} {{c_i}} {h_i}(t), where c_i are unknown coefficients. The collocation points are chosen as the midpoints $x_{j} = (j - \frac{1}{2}) Δ, Δ = \frac{1}{2 M}, j = 1, \dots, 2 M .$ {x_j} = \left( {j - {1 \over 2}} \right)\Delta ,\quad \Delta = {1 \over {2M}},\quad j = 1, \cdots ,2M.

Substituting the Haar expansion into the regularized equation (19) and enforcing the equality at each collocation point yields the linear system $\sum_{i = 1}^{2 M} (α h_{i} (x_{j}) + K_{j, i}) c_{i} = g_{1} (x_{j}), j = 1, \dots, 2 M,$ \sum\limits_{i = 1}^{2M} {\left( {\alpha {h_i}\left( {{x_j}} \right) + {K_{j,i}}} \right)} {c_i} = {g_1}\left( {{x_j}} \right),\quad j = 1, \cdots ,2M, where 21 $K_{j, i} = \int_{0}^{1} K_{N} (x_{j}, t) h_{i} (t) d t .$ {K_{j,i}} = \int_0^1 {{K_N}} \left( {{x_j},t} \right){h_i}(t)dt.

This system can be written in matrix form as $A c = g,$ A{\bf{c}} = {\bf{g}}, where A ∈ ℝ^2M×2M, c, g ∈ ℝ^2M defined by $A_{j, i} = α h_{i} (x_{j}) + K_{j, i}, c = (c_{i}), g = (g_{1} (x_{j})) .$ {A_{j,i}} = \alpha {h_i}\left( {{x_j}} \right) + {K_{j,i}},\quad {\bf{c}} = \left( {{c_i}} \right),\quad {\bf{g}} = \left( {{g_1}\left( {{x_j}} \right)} \right).

Since each Haar function is piecewise constant with compact support, the integrals in (21) can be evaluated exactly, avoiding numerical quadrature.

Now we consider efficient evaluation of the matrix entries. Using the truncated kernel $K_{N} (x, t) = - 2 e^{- A x} \sum_{n = 1}^{N} β_{n} \sin (w_{n} x) e^{A t} \sin (w_{n} t),$ {K_N}(x,t) = - 2{e^{ - Ax}}\sum\limits_{n = 1}^N {{\beta _n}} \sin \left( {{w_n}x} \right){e^{At}}\sin \left( {{w_n}t} \right), where $A = - \frac{c}{2}, w_{n} = n π, λ_{n} = \frac{c^{2}}{4} + n^{2} π^{2}, β_{n} = \frac{γ_{n}}{\sinh (γ_{n})},$ A = - {c \over 2},\quad {w_n} = n\pi ,\quad {\lambda _n} = {{{c^2}} \over 4} + {n^2}{\pi ^2},\quad {\beta _n} = {{{\gamma _n}} \over {\sinh \left( {{\gamma _n}} \right)}}, we obtain, for each collocation point x_j, $K_{j, i} = - 2 e^{- A x_{j}} \sum_{n = 1}^{N} β_{n} \sin (w_{n} x_{j}) I_{n, i},$ {K_{j,i}} = - 2{e^{ - A{x_j}}}\sum\limits_{n = 1}^N {{\beta _n}} \sin \left( {{w_n}{x_j}} \right){I_{n,i}}, with $I_{n, i} = \int_{supp (h_{i})} e^{A t} \sin (w_{n} t) h_{i} (t) d t .$ {I_{n,i}} = \int_{{\mathop{\rm supp}\nolimits} \left( {{h_i}} \right)} {{e^{At}}} \sin \left( {{w_n}t} \right){h_i}(t)dt.

Since h_i is supported on [ξ₁(i), ξ₃(i)] and changes sign at ξ₂(i), we have $I_{n, i} = \int_{ξ_{1} (i)}^{ξ_{2} (i)} e^{A t} \sin (w_{n} t) d t - \int_{ξ_{2} (i)}^{ξ_{3} (i)} e^{A t} \sin (w_{n} t) d t .$ {I_{n,i}} = \int_{{\xi _1}(i)}^{{\xi _2}(i)} {{e^{At}}} \sin \left( {{w_n}t} \right)dt - \int_{{\xi _2}(i)}^{{\xi _3}(i)} {{e^{At}}} \sin \left( {{w_n}t} \right)dt.

Using the antiderivative $\int e^{A t} \sin (w_{n} t) d t = \frac{e^{A t}}{λ_{n}} (A \sin (w_{n} t) - w_{n} \cos (w_{n} t)),$ \int {{e^{At}}} \sin \left( {{w_n}t} \right)dt = {{{e^{At}}} \over {{\lambda _n}}}\left( {A\sin \left( {{w_n}t} \right) - {w_n}\cos \left( {{w_n}t} \right)} \right), explicit expressions for I_n,i are obtained.

Case i = 1 (scaling function). For h 1 supported on [0, 1], $I_{n, 1} = \frac{w_{n}}{λ_{n}} (1 - {(- 1)}^{n} e^{A}) .$ {I_{n,1}} = {{{w_n}} \over {{\lambda _n}}}\left( {1 - {{( - 1)}^n}{e^A}} \right).

Case i = 2 (mother wavelet). For h₂ supported on [0, 1] with sign change at 1/2, $I_{n, 2} = \frac{2 e^{A / 2}}{λ_{n}} (A \sin (\frac{w_{n}}{2}) - w_{n} \cos (\frac{w_{n}}{2})) + \frac{w_{n}}{λ_{n}} (1 + {(- 1)}^{n} e^{A}) .$ {I_{n,2}} = {{2{e^{A/2}}} \over {{\lambda _n}}}\left( {A\sin \left( {{{{w_n}} \over 2}} \right) - {w_n}\cos \left( {{{{w_n}} \over 2}} \right)} \right) + {{{w_n}} \over {{\lambda _n}}}\left( {1 + {{( - 1)}^n}{e^A}} \right).

Case i ≥ 3 (wavelets). For i = m + k + 1, $ξ_{1} = \frac{k}{m}, ξ_{2} = \frac{k + 0.5}{m}, ξ_{3} = \frac{k + 1}{m},$ {\xi _1} = {k \over m},\quad {\xi _2} = {{k + 0.5} \over m},\quad {\xi _3} = {{k + 1} \over m}, we obtain $\begin{matrix} I_{n, i} & = & \frac{2 e^{A ξ_{2}}}{λ_{n}} (A \sin (w_{n} ξ_{2}) - w_{n} \cos (w_{n} ξ_{2})) \\ - \frac{e^{A ξ_{1}}}{λ_{n}} (A \sin (w_{n} ξ_{1}) - w_{n} \cos (w_{n} ξ_{1})) \\ - \frac{e^{A ξ_{3}}}{λ_{n}} (A \sin (w_{n} ξ_{3}) - w_{n} \cos (w_{n} ξ_{3})) . \end{matrix}$ \matrix{ {{I_{n,i}}} & = & {{{2{e^{A{\xi _2}}}} \over {{\lambda _n}}}\left( {A\sin \left( {{w_n}{\xi _2}} \right) - {w_n}\cos \left( {{w_n}{\xi _2}} \right)} \right)} \cr {} & {} & { - {{{e^{A{\xi _1}}}} \over {{\lambda _n}}}\left( {A\sin \left( {{w_n}{\xi _1}} \right) - {w_n}\cos \left( {{w_n}{\xi _1}} \right)} \right)} \cr {} & {} & { - {{{e^{A{\xi _3}}}} \over {{\lambda _n}}}\left( {A\sin \left( {{w_n}{\xi _3}} \right) - {w_n}\cos \left( {{w_n}{\xi _3}} \right)} \right).} \cr }

Finally, the entries of the collocation matrix are given explicitly by 22 $A_{j, i} = α h_{i} (x_{j}) - 2 e^{- A x_{j}} \sum_{n = 1}^{N} β_{n} \sin (n π x_{j}) I_{n, i},$ {A_{j,i}} = \alpha {h_i}\left( {{x_j}} \right) - 2{e^{ - A{x_j}}}\sum\limits_{n = 1}^N {{\beta _n}} \sin \left( {n\pi {x_j}} \right){I_{n,i}}, which completes the construction of the discrete system.

Numerical results

In this section, we present numerical results for the inverse Cauchy problem described in Section 3. The Haar wavelet method described in Section 4 is implemented with careful attention to the efficient evaluation of the matrix entries using the explicit formulas (22).

First, we fix the physical parameters as follows: $\begin{matrix} a_{1} = 1, a_{2} = 1, a^{2} = a_{1}^{2} + a_{2}^{2} = 2, \\ c = \frac{a^{2}}{a_{2}} = 2, A = - \frac{c}{2} = - 1. \end{matrix}$ \matrix{ {{a_1} = 1,\quad {a_2} = 1,\quad {a^2} = a_1^2 + a_2^2 = 2,} \cr {c = {{{a^2}} \over {{a_2}}} = 2,\quad A = - {c \over 2} = - 1.} \cr }

For each n ≥ 1, we compute: $γ_{n} = \sqrt{{(n π)}^{2} + \frac{c^{2}}{4} + a^{2}} = \sqrt{{(n π)}^{2} + 3}, β_{n} = \frac{γ_{n}}{\sinh (γ_{n})} .$ {\gamma _n} = \sqrt {{{(n\pi )}^2} + {{{c^2}} \over 4} + {a^2}} = \sqrt {{{(n\pi )}^2} + 3} ,\quad {\beta _n} = {{{\gamma _n}} \over {\sinh \left( {{\gamma _n}} \right)}}.

Table 1 and Fig. 1 show the computed values of γ_n and β_n for n = 1,2, ⋯, 20. The rapid decay of β_n indicates that the kernel series converges quickly. For n ≥ 4, β_n < 10⁻⁶, suggesting that truncating at N = 5 captures essentially all significant contributions. Therefore, we set N = 5 for all subsequent calculations.

Table 1

Values of γ_n and β_n for selected n.

n	γ_n	β_n	n	γ_n	β_n
1	4.3973	0.0251	7	22.4306	3.48 × 10⁻¹¹
2	7.0963	0.0009	8	25.7158	1.09 × 10⁻¹²
3	9.9713	3.16 × 10⁻⁵	9	29.0538	3.42 × 10⁻¹⁴
4	12.9614	1.06 × 10⁻⁶	10	32.4397	1.07 × 10⁻¹⁵
5	16.0442	3.44 × 10⁻⁸	15	49.9741	3.10 × 10⁻²³
6	19.2042	1.10 × 10⁻⁹	20	68.3345	8.70 × 10⁻³¹

We first consider solutions containing a single mode.

5.1

Case m = 1 (Fundamental mode)

The exact solution is: $Ψ_{1}^{exact} (x, y) = \frac{\sinh (γ_{1} y)}{\sinh (γ_{1})} e^{x} \sin (π x) .$ \Psi _1^{{\rm{exact }}}(x,y) = {{\sinh \left( {{\gamma _1}y} \right)} \over {\sinh \left( {{\gamma _1}} \right)}}{e^x}\sin (\pi x).

No Noise Case (δ = 0) (a) L-curve analysis: Figure 2 shows the L-curve for J = 3 (M = 8). The corner of the L-curve corresponds to the optimal regularization parameter α_opt = 10⁻³.

(b) Convergence study:

The convergence study presented in Table 2 examines the performance of the Haar wavelet method as the resolution level J increases. We observe that

The relative error ℰ decreases consistently with increasing resolution, from 6.21 × 10⁻² at J = 1 to 1.98 × 10⁻³ at J = 6, representing an overall error reduction of approximately 31 ×. This demonstrates the method’s numerical convergence.
The method achieves satisfactory accuracy (ℰ < 0.2%) with moderate resolution (J = 6), making it computationally feasible for practical applications.

Table 2

Convergence study for m = 1, δ = 0.

J	M = 2^J	Relative Error ℰ	Convergence Rate	Error Reduction
1	2	6.21e-02	–	–
2	4	3.15e-02	0.98	1.97
3	8	1.58e-02	1.00	1.99
4	16	7.92e-03	1.00	2.00
5	32	3.96e-03	1.00	2.00
6	64	1.98e-03	1.00	2.00

(c) Boundary trace reconstruction: Figure 3 compares the exact boundary trace φ^exact(x) = e^x sin(πx) with the reconstructed trace. The reconstruction is visually indistinguishable from the exact solution.

(d) Full solution reconstruction: Figure 4 shows the exact and reconstructed solutions over Ω. The maximum pointwise error is approximately 2.5 × 10⁻³, occurring near the inaccessible boundary y = 1.

(e) Error profile: Figure 5 shows the absolute error along the vertical line x = 0.5. The error grows from the accessible boundary (y = 0) toward the inaccessible boundary (y = 1), as expected for this ill-posed inverse problem.

Noisy Data Cases (a) Relative errors: Table 3 shows the reconstruction errors for different noise levels with J = 3. The optimal regularization parameter α increases with noise level to provide appropriate stabilization.

Table 3

Performance comparison and error increase factors for different noise levels m = 1 J = 3

Noise Level δ	Relative Error ℰ(δ)	Optimal α	Error Increase Factor	Noise-to-Error Ratio	α Increase Factor
0.000	0.016	0.001	1.00	–	1.00
0.001	0.018	0.001	1.13	0.056	1.00
0.005	0.026	0.003	1.63	0.200	3.00
0.010	0.035	0.005	2.19	0.286	5.00
0.050	0.089	0.018	5.56	0.562	18.00
0.100	0.142	0.032	8.88	0.704	32.00

The performance comparison shows several important characteristics of the proposed method’s robustness to noise. The error increase factor is calculated as: $Error Increase Factor = \frac{ℰ (δ)}{ℰ (0)},$ {\rm{Error Increase Factor }} = {{{\cal E}({\bf{\delta }})} \over {{\cal E}(0)}}, where ℰ(0) = 0.016 is the noise-free reference error.

From this table, we observe also that:

The reconstruction error increases nonlinearly with the noise level:

For 0.1% noise: Error increases by factor 1.13× compared to noise-free For 0.5% noise: Error increases by factor 1.63 × compared to noise-free For 1% noise: Error increases by factor 2.19 × compared to noise-free For 5% noise: Error increases by factor 5.56 × compared to noise-free For 10% noise: Error increases by factor 8.88 × compared to noise-free

This controlled error growth demonstrates the method’s strong stability properties even under extreme noise conditions.
The Noise-to-Error Ratio decreases from 0.056 (for δ = 0.001) to 0.704 (for δ = 0.100), suggesting diminishing returns in noise handling as noise levels increase.
The optimal regularization parameter α increases approximately linearly with noise level: $α \approx 0.32 δ (empirical relationship from the data)$ \alpha \approx 0.32\delta {\rm{ (empirical relationship from the data)}}

This linear relationship indicates consistent regularization behavior across different noise levels.
For noise up to (δ = 0.010): The error increases by a factor of at most 2 compared to the noise-free case.

Note that the nonlinear relationship between noise level and reconstruction error is characteristic of ill-posed inverse problems, where small perturbations in input data can lead to disproportionately large errors in the solution.

The error remains consistently subordinate to the input noise level up to 5% noise. For noise levels ≤ 1%, the method achieves good precision with errors below 3.5%, rapidly approaching the theoretical noise-free performance.

(b) Boundary trace reconstruction with noise: Figure 6 compares the exact boundary trace φ^exact(x) and the reconstructed trace φ_a(x) for J = 3 and various noise levels. This figure provides comprehensive visual evidence of the method’s noise handling capabilities:

For 10% noise: Despite substantial noise contamination, the reconstructed solution maintains correct phase and amplitude characteristics with smooth recovery of the underlying harmonic pattern.

For 5% noise: the figures show good noise suppression with reconstructed solutions closely tracking the exact profile. Minor residual oscillations are well within acceptable limits for practical applications.

For 1%, 0.5% and 0.1% noise: the figures demonstrate virtually perfect reconstruction where computed and exact solutions are visually identical, confirming superb performance under high-quality measurement conditions.

This controlled error growth demonstrates the method’s strong stability properties even under extreme noise conditions.

5.2

Case multi-mode

We now consider a more challenging case with a solution containing multiple frequencies: 23 $Ψ_{1}^{exact} (x, y) = \frac{\sinh (γ_{1} y)}{\sinh (γ_{1})} e^{x} \sin (π x) + \frac{\sinh (γ_{4} y)}{\sinh (γ_{4})} e^{x} \sin (4 π x) .$ \Psi _1^{{\rm{exact }}}(x,y) = {{\sinh \left( {{\gamma _1}y} \right)} \over {\sinh \left( {{\gamma _1}} \right)}}{e^x}\sin (\pi x) + {{\sinh \left( {{\gamma _4}y} \right)} \over {\sinh \left( {{\gamma _4}} \right)}}{e^x}\sin (4\pi x).

This solution combines low-frequency (sin(πx)) and high-frequency (sin(4πx)) components, making it a more realistic test case.

No Noise Case (δ = 0) For J = 3, the optimal regularization parameter selected via the L-curve is α_opt = 5 × 10⁻³. The relative reconstruction error is 0.047, which is larger than for either single mode alone due to the interaction between modes.

Figure 7 shows the boundary trace reconstruction. The method successfully captures both the overall shape and the high-frequency oscillations, though some smoothing is evident in the peaks of the high-frequency component.

Here below we consider performance analysis for multi-mode solution.

5.3

Error analysis with noise

Table 4 presents the performance comparison including error increase factors for the multi-mode solution at different noise levels. Table 5 gives the comparison of error increase factors single mode (m = 1) vs multi-mode.

Table 4

Performance comparison and error increase factors for multi-mode solution with different noise levels, J = 3.

Noise Level δ	Relative Error ℰ(δ)	Optimal α	Error Increase Factor	Noise-to-Error Ratio	α Increase Factor
0.000	0.047	0.005	1.00	–	1.00
0.001	0.048	0.005	1.02	0.021	1.00
0.005	0.051	0.007	1.09	0.098	1.40
0.010	0.058	0.010	1.23	0.172	2.00
0.050	0.112	0.032	2.38	0.446	6.40
0.100	0.183	0.056	3.89	0.546	11.20

Table 5

Comparison of error increase factors: Single mode (m = 1) vs Multi-mode.

Noise Level	Error Increase Factor	Error Increase Factor	Ratio
δ	Single mode (m = 1)	Multi-mode	(Multi/Single)
0.001	1.13	1.02	0.90
0.005	1.63	1.09	0.67
0.010	2.19	1.23	0.56
0.050	5.56	2.38	0.43
0.100	8.88	3.89	0.44

5.4

Comparison with single mode results

To understand the differences in noise sensitivity, we compare the multi-mode results with the single mode (m = 1) case from Table 3:

We observe that

The multi-mode solution exhibits significantly lower sensitivity to noise compared to the single mode solution. For example, with 10% noise (δ = 0.100), the error increases by a factor of 3.89 for the multimode solution versus 8.88 for the single mode solution.
While error still increases nonlinearly with noise level, the growth is less dramatic for the multi-mode case:
- 0.1% noise → 2% error increase (factor: 1.02)
- 0.5% noise → 9% error increase (factor: 1.09)
- 1% noise → 23% error increase (factor: 1.23)
- 5% noise → 138% error increase (factor: 2.38)
- 10% noise → 289% error increase (factor: 3.89)
The noise-free reconstruction error for the multi-mode solution (ℰ(0) = 0.047) is approximately three times higher than for the single mode solution (ℰ(0) = 0.016). This reflects the inherent difficulty of reconstructing solutions with multiple frequency components.
The optimal regularization parameter α for the multi-mode solution is consistently higher than for the single mode solution at equivalent noise levels. This suggests that stronger regularization is necessary to handle the high-frequency components effectively.

The observed behavior can be explained by several factors:

Regularization Effects: The stronger regularization required for multi-mode solutions (higher a values) provides enhanced noise suppression, albeit at the potential cost of smoothing high-frequency details. This trade-off appears beneficial for noise robustness.

Error Composition: The total reconstruction error for the multi-mode solution comprises both approximation errors from high-frequency components and noise-induced errors. The relative contribution of additive noise may be smaller when substantial approximation errors are already present.

Frequency Interactions: Errors in reconstructing different frequency components may interact in ways that partially compensate or average out random noise contributions, reducing overall sensitivity to measurement noise.

Figure 8 shows reconstructions for various noise levels. With 10% noise, the high-frequency oscillations become difficult to distinguish, though the low-frequency component is still reasonably captured.

5.5

Advantages of the proposed method

The proposed approach offers several distinct advantages over conventional numerical methods such as finite elements, finite differences, or functional optimization techniques. Unlike these methods, which typically require iterative procedures, mesh generation, and the solution of multiple large-scale systems, our method is direct and mesh-free. The solution is obtained by solving a single linear system with a matrix of very small size (e.g., 2M × 2M with M = 2^J and typically J ≤ 6). This avoids the computational overhead associated with mesh refinement, remeshing, or iterative updates.

For instance, solving the regularized Fredholm integral equation of the second kind using a finite difference approach necessitates an iterative method starting from an interpolation polynomial obtained by solving a system of equations [75]. Such iterative schemes can be computationally expensive and sensitive to the choice of initial guess. In contrast, our method computes the solution directly via matrix inversion or a fast linear solver, making it exceptionally fast and suitable for real-time or large-scale inverse problems where computational efficiency is critical.

Furthermore, the mesh-free nature of the Haar wavelet discretization circumvents mesh-related issues such as distortion, adaptivity, or singularities, which are common challenges in finite element or finite difference formulations. This simplicity, combined with the explicit formulas derived for the matrix entries, ensures that the method is both accurate and easy to implement.

Conclusion

This work introduces an original and comprehensive numerical framework for reconstructing inaccessible boundary conditions in the inverse Cauchy problem associated with the two-dimensional Pauli equation. By exploiting the diagonal structure of the Pauli operator, the problem is reduced to a scalar boundary condition task. Using separation of variables, we reformulate the obtained problem as a Fredholm integral equation of the first kind, which we regularize using the Lavrentiev method to mitigate its ill-posedness. This work introduces an original and comprehensive numerical framework for reconstructing inaccessible boundary conditions in the inverse Cauchy problem associated with the two-dimensional Pauli equation. By exploiting the diagonal structure of the Pauli operator, the problem is reduced to a scalar boundary condition task. Using separation of variables, we reformulate the obtained problem as a Fredholm integral equation of the first kind, which we regularize using the Lavrentiev method to mitigate its ill-posedness.

The numerical scheme, based on Haar wavelet discretization combined with an explicit evaluation of the kernel integrals, proved both efficient and accurate. In the noise-free setting, the method exhibited first-order convergence with respect to the resolution level, achieving a relative error below 0.2% at moderate discretization. The L-curve criterion provided a reliable mechanism for selecting the optimal regularization parameter.

The explicit formulas derived for the matrix entries avoid numerical quadrature, ensuring computational efficiency and making the method suitable for practical implementation. The results confirm that the proposed Haar wavelet approach, combined with Lavrentiev regularization, constitutes a stable and convergent numerical tool for this class of ill-posed inverse problems.

Future work could explore the extension of this methodology to time-dependent Pauli systems, three-dimensional configurations, or problems involving more complex electromagnetic potentials. Additionally, adaptive wavelet strategies or alternative regularization techniques might further improve the reconstruction of high-frequency features in highly noisy environments.

Declarations

Reduction to a Fredholm integral equation and numerical solution of the inverse Cauchy problem for the Schrödinger-Pauli equation

Full Article

Paradigm

My account