¹

1 Introduction↩︎

Pseudo-parabolic equations (Sobolev-type equations) are a mathematical model for a wide class of important physical and mechanical applications, such as filtration processes in porous media [1], population aggregation dynamics [2], nonlinear long-wave propagation [3], [4], second-order unsteady fluid flows [5], fluid transport in fractured rocks [6], and non-Newtonian fluid motion [7], [8]. The inverse problems for a linear, semilinear and nonlinear pseudoparabolic equation, both by itself and with variable coefficients, have been extensively investigated; a representative selection of related works can be found in [9]–[10].

In this work, we study a backward problem for a time-fractional pseudo-parabolic equation \[\label{1461} \partial_t^\alpha u(x,t)-u_{xx}(x,t)-\mu(t)\,u_{xxt}(x,t)=f(x,t), \qquad (x,t)\in(0,l)\times(0,T),\tag{1}\] subject to the Dirichlet boundary conditions \[\label{1462} u(0,t)=u(l,t)=0,\qquad t\in[0,T],\tag{2}\] and the final-time measurement \[\label{finalT} u(x,T)=\psi(x),\qquad x\in[0,l].\tag{3}\] Here, \(u(x,t)\) denotes the temperature, \(f(x,t)\) is a given source term, \(\mu(t)\) is a pseudo-parabolic coefficient, and \(\psi\) is the measured final-time state. The operator \(\partial_t^\alpha\) denotes the Caputo fractional derivative of order \(\alpha\in(0,1)\)[11]. The backwards-in-time problem consists of reconstructing the initial state \[\label{initial} u_0(x):=u(x,0),\qquad x\in[0,l],\tag{4}\] and, consequently, the whole trajectory \(u\) on \([0,T]\) from the data \((f,\psi)\).

The equation 1 models a time-fractional pseudo-parabolic diffusion process with memory and relaxation effects, describing anomalous transport phenomena influenced by time-dependent material properties and external forcing in the space-time domain \((0,l)\times(0,T)\). In engineering applications, backwards-in-time problems are essential because they focus on reconstructing past physical states using data observed at a later time.

A pioneering contribution to the well-posedness analysis of the backward-in-time problem for the time-fractional diffusion equation is given in [12]. Further developments for backward problems in fractional diffusion-wave models include theoretical and numerical studies in bounded domains [13], as well as iterative regularisation approaches (e.g., Landweber-type methods) for identifying unknown initial data in time-space fractional settings [14]. General well-posedness frameworks for backwards-in-time problems for time-fractional diffusion and diffusion-wave equations were also established in [15], [16]. In addition, simultaneous recovery of multiple initial values for time-fractional diffusion–wave equations was investigated in [17].

For fractional pseudo-parabolic equations, backward problems and related stability issues have been analysed in [18], [19]. Regularisation methods for reconstructing the initial value in time-fractional pseudo-parabolic models were developed in [20], and robustness under random noise was further studied in [21]. Related inverse initial data problems for time-fractional pseudo-hyperbolic equations were considered in [22].

Despite recent progress, backwards-in-time analysis for time-fractional pseudo-parabolic equations remains underexplored, particularly in terms of well-posedness theory, numerical reconstruction, and noisy data treatment. To address this deficiency, we undertake a comprehensive theoretical and numerical investigation of the problem.

Compared with the constant-coefficient case, the present variable-coefficient problem has an essentially different spectral structure. If \(\mu\) is constant, the Fourier modes may be studied by using coefficient-independent damping factors or special-function representations. In contrast, when \(\mu=\mu(t)\), the \(k\)-th Fourier coefficient satisfies \[\partial_t^\alpha u_k(t)+\mu(t)\lambda_k u_k'(t)+\lambda_k u_k(t)=f_k(t),\] where the coefficient of \(u_k'\) depends on time. Therefore, the backward reconstruction cannot be reduced to a simple explicit multiplier formula. Instead, each mode must be analysed through a weakly singular Volterra equation with a time-dependent kernel. The main technical point is to obtain estimates that are uniform with respect to the Fourier mode \(k\), in particular, the positivity and lower boundedness of the reconstruction factor \(\mathcal{A}_k(T)\).

The main contributions of this paper are as follows. First, we establish classical solvability of the backward problem under explicit regularity and compatibility assumptions on the final data and source term. Second, by reducing the Fourier modes to Volterra equations with time-dependent kernels, we prove the positivity of the reconstruction factor and derive an \(L^2\)-stability estimate for the recovered initial state. Third, we construct a fully discrete finite-difference solver based on the graded-mesh L1 approximation and central spatial differences, and combine it with Tikhonov regularisation for stable initial-state reconstruction from noisy final-time data.

The remainder of the manuscript is structured as follows: Section 2 establishes the existence and uniqueness of the backward problem. Section 3 develops and analyses a stable numerical scheme for the direct model. Section 4 presents the reconstruction procedure for identifying the unknown initial state from finite-time measurements. Numerical validations, including tests with noisy data, are given in Section 5.

2 Existence and Uniqueness of the problem↩︎

In this section, we establish the existence and uniqueness of a solution to the backwards-in-time problem for the time fractional pseudo-parabolic equation 1 subject to the Dirichlet boundary conditions 2 and the final-time measurement 3 .

2.1 Definition and assumptions↩︎

Assumptions. Throughout this section, we assume:

1. \(\alpha\in(0,1)\), \(\mu\in C([0,T])\), and there exists \(\mu_0>0\) such that \(\mu(t)\ge \mu_0\) for all \(t\in[0,T]\).

2. \(\psi\in H^3(0,l)\) and \[\psi(0)=\psi(l)=0,\qquad \psi''(0)=\psi''(l)=0.\]

3. \(f\in L^2\big(0,T;H^3(0,l)\big)\) and for a.e. \(t\in(0,T)\), \[f(0,t)=f(l,t)=0,\qquad f_{xx}(0,t)=f_{xx}(l,t)=0.\]

For uniform-in-time estimates of the time-derivative series, we shall also use:

1. For every \(\delta\in(0,T)\), \(f\in L^\infty(\delta,T;H^3(0,l))\).

2.2 Existence and uniqueness↩︎

2.3 Continuous dependence on the data↩︎

In this subsection, we prove that, under assumptions [F1]–[F3], the reconstruction of the initial state \(u_0(x)=u(x,0)\) from the final-time measurement \(\psi(x)=u(x,T)\) continuously depends on the data.

3 Finite difference approximation of the direct problem↩︎

In this section, we describe the finite-difference discretisation used in the implementation for the time-fractional pseudo-parabolic direct problem \[\label{eq:direct95frac} \partial_t^\alpha u(x,t)-u_{xx}(x,t)-\mu(t)\,u_{xxt}(x,t)=f(x,t), \qquad (x,t)\in(0,l)\times(0,T),\tag{37}\] subject to the homogeneous Dirichlet boundary conditions \(u(0,t)=u(l,t)=0\). We assume that \(\mu\in C([0,T])\) and \(f\) is given. The Caputo derivative is discretised by an \(L1\)-type scheme on a graded temporal mesh, while the spatial derivatives are approximated by second-order central differences. The resulting method leads to a tridiagonal linear system at each time level[24].

3.1 Discretisations↩︎

In this subsection, we describe spatial and temporal discretisations and the full discrete scheme.

Spatial discretisation. Let \(x_i=ih\) for \(i=0,\dots,N\), where \(h=l/N\). For a grid vector \(\mathbf{v}=(v_1,\dots,v_{N-1})^\top\in\mathbb{R}^{N-1}\) we enforce the boundary conditions by setting \(v_0=v_N=0\) and define the standard second-order discrete Laplacian \(L_h:\mathbb{R}^{N-1}\to\mathbb{R}^{N-1}\) by \[\label{eq:lap95frac} (L_h\mathbf{v})_i=\frac{v_{i-1}-2v_i+v_{i+1}}{h^2},\qquad i=1,\dots,N-1.\tag{38}\] The matrix representation of \(L_h\) is symmetric tridiagonal and \((-L_h)\) is symmetric positive definite on \(\mathbb{R}^{N-1}\).

We denote the interior nodal vector of the numerical solution by \[\mathbf{u}^k:=(u_1^k,\dots,u_{N-1}^k)^\top,\qquad u_i^k\approx u(x_i,t_k),\] where \(\{t_k\}_{k=0}^M\) is the temporal grid defined below.

L1 approximation on graded meshes. Let \(0=t_0<t_1<\cdots<t_M=T\) be a time grid. In the implementation, we allow graded meshes of the form \[\label{eq:graded95mesh} t_k=T\Bigl(\frac{k}{M}\Bigr)^r,\,\;k=0,1,\dots,M,\,\, r\ge 1(\text{uniform mesh corresponds to } r=1),\tag{39}\] with \(\tau_k:=t_k-t_{k-1}\) and \(\tau_{\max}:=\max_{1\le k\le M}\tau_k\).

For each \(k\ge1\), we approximate the Caputo derivative at \(t_k\) using the \(L1\) scheme \[\label{eq:L195caputo} \partial_t^\alpha u(x_i,t_k) \approx \frac{1}{\Gamma(1-\alpha)} \sum_{j=1}^{k} \int_{t_{j-1}}^{t_j}\frac{u_t(x_i,s)}{(t_k-s)^\alpha}\,ds \;\approx\; \sum_{j=1}^{k} w_{k,j}\,\bigl(u_i^j-u_i^{j-1}\bigr),\tag{40}\] where the weights \(w_{k,j}\) depend on the mesh. For a nonuniform grid, they can be written as \[\label{eq:L195weights} w_{k,j} = \frac{1}{\Gamma(2-\alpha)}\,d_{k,j}, \quad d_{k,j} = \frac{(t_k-t_{j-1})^{1-\alpha}-(t_k-t_j)^{1-\alpha}}{t_j-t_{j-1}}, \quad 1\le j\le k.\tag{41}\]

Fully discrete scheme. Using \(u_{xxt}=(u_{xx})_t\) and the discrete Laplacian \(L_h\), we approximate at time \(t_k\): \[u_{xx}(x_i,t_k)\approx (L_h\mathbf{u}^k)_i, \qquad u_{xxt}(x_i,t_k)\approx \frac{(L_h\mathbf{u}^k)_i-(L_h\mathbf{u}^{k-1})_i}{\tau_k}.\] Let \(\boldsymbol{\mu}^k:=\mu(t_k)\) and \(\mathbf{f}^k:=\bigl(f(x_1,t_k),\dots,f(x_{N-1},t_k)\bigr)^\top.\) Combining 40 –41 with the above spatial discretisations yields, for \(k=1,\dots,M\), the linear system \[\label{eq:scheme95matrix95form} \Bigl(\frac{d_{k,k}}{\Gamma(2-\alpha)}I - \Bigl(1+\frac{\mu^k}{\tau_k}\Bigr)L_h\Bigr)\mathbf{u}^k = \mathbf{r}^k + \mathbf{f}^k - \frac{\mu^k}{\tau_k}L_h\mathbf{u}^{k-1},\tag{42}\] where \(\mathbf{r}^k\) collects the history (memory) contributions of the \(L1\) approximation, \[\label{eq:history95term} \mathbf{r}^k = \frac{d_{k,1}}{\Gamma(2-\alpha)}\mathbf{u}^0 + \frac{1}{\Gamma(2-\alpha)} \sum_{j=1}^{k-1}\bigl(d_{k,j+1}-d_{k,j}\bigr)\mathbf{u}^j .\tag{43}\] Here \(d_{k,j}\) are defined in 41 . The left-hand side matrix in 42 is tridiagonal because \(L_h\) is tridiagonal. Therefore, at each time step the system is solved efficiently by the Thomas algorithm.

Given an initial state \(u_0(x)=u(x,0)\), we set \[\boldsymbol{u_{0,h}} := \bigl(u_0(x_1),\dots,u_0(x_{N-1})\bigr)^\top\in\mathbb{R}^{N-1}, \qquad \mathbf{u}^0=\boldsymbol{u_{0,h}}.\] The direct solver 42 –43 then produces the discrete trajectory \(\{\mathbf{u}^k\}_{k=0}^M\) and, in particular, the terminal vector \(\mathbf{u}^M\) used in the inverse step.

3.2 Stability↩︎

In this subsection, we study the stability of the fully discrete method defined by 42 –43 . The Caputo derivative is approximated by the graded-mesh \(L1\) formula, with the weights \(d_{k,j}\) given in 41 . Throughout this subsection, we assume that \[\label{eq:ass95mu} \mu\in C([0,T]), \qquad 0<\mu_0\le \mu(t)\le \mu_{\max}, \qquad t\in[0,T].\tag{44}\] We use the discrete inner product and norms introduced earlier. In particular, \[(-L_h\mathbf{v},\mathbf{v})_h = \|\nabla_h\mathbf{v}\|_h^2.\]

Preliminary properties of the graded-mesh \(L1\) operator. Define the graded-mesh \(L1\) operator componentwise by \[\label{eq:disc95caputo95def} \delta_t^\alpha\mathbf{v}^k := \frac{1}{\Gamma(2-\alpha)} \sum_{j=1}^{k} d_{k,j}\bigl(\mathbf{v}^j-\mathbf{v}^{j-1}\bigr), \qquad k\ge1,\tag{45}\] where \(d_{k,j}\) are defined in 41 . For every admissible temporal mesh, and in particular for the graded mesh considered here, one has \[\label{eq:d95props} d_{k,j}>0, \qquad d_{k,1}\le d_{k,2}\le\cdots\le d_{k,k}, \qquad 1\le j\le k.\tag{46}\]

Solvability at each time level. Let \(\tau_k=t_k-t_{k-1}\) and \(\mu^k=\mu(t_k)\). The time-stepping system 42 can be written as \[\label{eq:Ak95system} \mathbf{A}^k\mathbf{u}^k=\mathbf{b}^k, \qquad \mathbf{A}^k = \frac{d_{k,k}}{\Gamma(2-\alpha)}I - \left(1+\frac{\mu^k}{\tau_k}\right)L_h,\tag{47}\] where \(\mathbf{b}^k\) contains the history term 43 , the source \(\mathbf{f}^k\), and the contribution \[-\frac{\mu^k}{\tau_k}L_h\mathbf{u}^{k-1}.\]

Let \(\mathcal{A}_h:=-L_h.\) Since \(\mathcal{A}_h\) is symmetric positive definite, there exists an orthonormal basis \(\{\boldsymbol{\varphi}_q\}_{q=1}^{N-1}\), with respect to \((\cdot,\cdot)_h\), such that \[\mathcal{A}_h\boldsymbol{\varphi}_q = \lambda_{q,h}\boldsymbol{\varphi}_q, \qquad 0<\lambda_{1,h}\le\lambda_{2,h}\le\cdots\le\lambda_{N-1,h}.\] The discrete Poincaré inequality implies \[\label{eq:lambda95lower95bound} \lambda_{1,h}^{-1/2}\le C_P,\tag{48}\] where \(C_P\) is independent of \(h\).

3.3 Convergence↩︎

In this subsection, we establish an error estimate for the fully discrete scheme 42 –43 on the graded temporal mesh introduced in Subsection 3.1.

For the graded mesh 39 , the time-step sizes are nondecreasing, and hence \[\tau_{\max}=\tau_M = T\left[ 1-\left(1-\frac{1}{M}\right)^r \right] \le \frac{rT}{M}.\] Thus, for every fixed \(r\ge1\), \[\label{eq:tau95max95graded} \tau_{\max}=O(M^{-1}).\tag{53}\]

Let \(R_h\) denote the nodal restriction operator \[R_hv:= \bigl(v(x_1),\dots,v(x_{N-1})\bigr)^\top.\] For the exact solution \(u\), define \[\mathbf{U}^k:=R_hu(\cdot,t_k), \qquad k=0,\dots,M.\] Let \(\{\mathbf{u}^k\}_{k=0}^{M}\) be the numerical solution generated by 42 –43 , and define the nodal error by \[\mathbf{e}^k:=\mathbf{U}^k-\mathbf{u}^k, \qquad k=0,\dots,M.\] Since the exact initial data are imposed in the numerical scheme, \[\label{eq:error95initial95zero} \mathbf{e}^0=\mathbf{0}.\tag{54}\]

Additional regularity assumptions. Theorem 4 establishes classical solvability under assumptions [F1]–[F4]. The convergence analysis below is conditional on additional space–time regularity of the exact solution. These stronger assumptions are not asserted to follow from [F1]–[F4]; they are imposed in order to obtain uniform consistency estimates for the central-difference approximation, the backward-difference approximation of \(u_{xxt}\), and the \(L1\) approximation of the Caputo derivative.

We assume that 44 holds and that \[\label{eq:reg95assumption95smooth} u\in C^1\bigl([0,T];H^4(0,l)\bigr) \cap C^2\bigl([0,T];H^3(0,l)\bigr), \qquad u(\cdot,t)\in H_0^1(0,l) \quad\text{for all }t\in[0,T],\tag{55}\] with \[\label{eq:reg95assumption95bound} \max_{0\le t\le T} \left( \|u(\cdot,t)\|_{H^4(0,l)} + \|u_t(\cdot,t)\|_{H^4(0,l)} + \|u_{tt}(\cdot,t)\|_{H^3(0,l)} \right) \le C_u.\tag{56}\] In particular, \[u_{xxxx t}\in C\bigl([0,T];L^2(0,l)\bigr), \qquad u_{xxtt}\in C\bigl([0,T];H^1(0,l)\bigr).\] Assumption 55 excludes the weak initial singularity that commonly occurs in time-fractional evolution problems. The treatment of such nonsmooth solutions requires time-dependent consistency estimates and a nonuniform discrete fractional Grönwall argument and is beyond the scope of the present analysis.

Consistency errors. Define the consistency errors by \[\begin{align} \boldsymbol{\xi}_\alpha^k &:= \delta_t^\alpha\mathbf{U}^k - R_h\partial_t^\alpha u(\cdot,t_k), \tag{57} \\ \boldsymbol{\xi}_x^k &:= L_h\mathbf{U}^k - R_hu_{xx}(\cdot,t_k), \tag{58} \\ \boldsymbol{\xi}_{xt}^k &:= \frac{L_h\mathbf{U}^k-L_h\mathbf{U}^{k-1}}{\tau_k} - R_hu_{xxt}(\cdot,t_k), \tag{59} \end{align}\] for \(1\le k\le M\).

For every \(v\in H^4(0,l)\), the standard central-difference estimate gives \[\label{eq:central95difference95general} \|L_hR_hv-R_hv_{xx}\|_h \le Ch^2\|v\|_{H^4(0,l)},\tag{60}\] where \(C\) is independent of \(h\). Consequently, \[\label{eq:space95consistency95bound} \|\boldsymbol{\xi}_x^k\|_h \le Ch^2, \qquad 1\le k\le M.\tag{61}\]

To estimate the consistency error associated with the pseudo-parabolic term, introduce \[\overline{u_t}^{\,k} := \frac{1}{\tau_k} \int_{t_{k-1}}^{t_k}u_t(\cdot,s)\,ds.\] Then \[\frac{\mathbf{U}^k-\mathbf{U}^{k-1}}{\tau_k} = R_h\overline{u_t}^{\,k},\] and hence \[\begin{align} \boldsymbol{\xi}_{xt}^k &= L_hR_h\overline{u_t}^{\,k} - R_hu_{xxt}(\cdot,t_k) \notag\\ &= \left( L_hR_h\overline{u_t}^{\,k} - R_h(\overline{u_t}^{\,k})_{xx} \right) + R_h\left( (\overline{u_t}^{\,k})_{xx} -u_{xxt}(\cdot,t_k) \right). \label{eq:pseudo95consistency95decomp} \end{align}\tag{62}\]

By 60 and 56 , \[\label{eq:pseudo95space95part} \left\| L_hR_h\overline{u_t}^{\,k} - R_h(\overline{u_t}^{\,k})_{xx} \right\|_h \le Ch^2.\tag{63}\]

For the second term in 62 , set \[z^k := (\overline{u_t}^{\,k})_{xx} -u_{xxt}(\cdot,t_k).\] Then \[z^k = \frac{1}{\tau_k} \int_{t_{k-1}}^{t_k} \left( u_{xxt}(\cdot,s)-u_{xxt}(\cdot,t_k) \right)\,ds.\] Using the fundamental theorem of calculus, we obtain \[\begin{align} \|z^k\|_{H^1(0,l)} &\le \frac{1}{\tau_k} \int_{t_{k-1}}^{t_k} \int_s^{t_k} \|u_{xxtt}(\cdot,\eta)\|_{H^1(0,l)} \,d\eta\,ds \\ &\le C\tau_k. \end{align}\] Moreover, the one-dimensional Sobolev embedding \(H^1(0,l)\hookrightarrow C[0,l]\) gives \[\begin{align} \|R_hz^k\|_h^2 &= h\sum_{i=1}^{N-1}|z^k(x_i)|^2 \le l\|z^k\|_{L^\infty(0,l)}^2 \le C\|z^k\|_{H^1(0,l)}^2. \label{eq:nodal95restriction95bound} \end{align}\tag{64}\] Therefore, \[\|R_hz^k\|_h\le C\tau_k.\] Combining this estimate with 63 yields \[\label{eq:pseudo95consistency95bound} \|\boldsymbol{\xi}_{xt}^k\|_h \le C\bigl(h^2+\tau_k\bigr), \qquad 1\le k\le M.\tag{65}\]

Under the smoothness assumption 55 , the standard consistency estimate for the nonuniform \(L1\) formula gives \[\label{eq:L195smooth95consistency} \max_{1\le k\le M} \|\boldsymbol{\xi}_\alpha^k\|_h \le C\tau_{\max}^{\,2-\alpha}.\tag{66}\] For consistency and convergence analyses of the \(L1\) approximation on graded temporal meshes, see, for example, [25]. The constant \(C\) in 66 is independent of \(h\) and \(M\), but may depend on \(T\), \(\alpha\), the fixed grading parameter \(r\), and the regularity constant \(C_u\).

Estimate 66 is used only under 55 . In general, it does not hold uniformly over all time levels when \[u_{tt}(t)=O(t^{\alpha-2}) \qquad\text{as }t\to0^+.\]

Error equation. Evaluating the continuous problem at the spatial grid points and subtracting the fully discrete scheme gives \[\label{eq:error95eq95graded95full} \delta_t^\alpha\mathbf{e}^k - L_h\mathbf{e}^k - \mu^k \frac{L_h\mathbf{e}^k-L_h\mathbf{e}^{k-1}}{\tau_k} = \boldsymbol{\rho}^k, \qquad 1\le k\le M,\tag{67}\] where \(\mathbf{e}^0=\mathbf{0}\) and \[\label{eq:rho95definition} \boldsymbol{\rho}^k = \boldsymbol{\xi}_\alpha^k - \boldsymbol{\xi}_x^k - \mu^k\boldsymbol{\xi}_{xt}^k.\tag{68}\] Using 44 , 61 , 65 , and 66 , we obtain \[\label{eq:rho95bound95full} \|\boldsymbol{\rho}^k\|_h \le C\left( h^2+\tau_k+\tau_{\max}^{\,2-\alpha} \right), \qquad 1\le k\le M.\tag{69}\]

4 Numerical reconstruction of the initial state↩︎

In this section, we describe the fully discrete reconstruction of the unknown initial state \(u_0=u(\cdot,0)\) from the final-time measurement \(\psi=u(\cdot,T)\) for problem 1 –3 .

The reconstruction is based on a discretise-then-regularise strategy. First, the discrete forward operator is constructed by repeated applications of the graded-mesh \(L1\) finite-difference solver developed in Section 3. The resulting finite-dimensional inverse system is then stabilised by zeroth-order Tikhonov regularisation.

4.1 Discrete forward operator and inverse relation↩︎

We use the spatial and temporal grids introduced in Subsection 3.1. For an initial vector \(\mathbf{u}_{0,h}\in\mathbb{R}^{N-1}\), the fully discrete scheme produces the trajectory \(\{\mathbf{u}^k\}_{k=0}^{M}\), with \(\mathbf{u}^0=\mathbf{u}_{0,h}\).

Let \[\boldsymbol{\psi}_h := R_h\psi = \bigl(\psi(x_1),\dots,\psi(x_{N-1})\bigr)^\top\] denote the exact nodal final-time data. When the measurements are perturbed, we write \[\boldsymbol{\psi}_h^\delta = \boldsymbol{\psi}_h+\boldsymbol{\eta}_h^\delta,\] where \(\boldsymbol{\eta}_h^\delta\) denotes the measurement error and \(\delta\ge0\) represents its level. The noise-free case corresponds to \(\delta=0\).

Let \[\mathcal{S}_{h,\tau}: \mathbb{R}^{N-1}\times\mathcal{F}_h \longrightarrow \mathbb{R}^{N-1}\] denote the discrete terminal-state map induced by the fully discrete scheme 42 –43 . Thus, for a prescribed initial vector \(\mathbf{v}\) and a discrete source \(f\), \[\mathcal{S}_{h,\tau}(\mathbf{v};f) = \mathbf{u}^M(\mathbf{v};f).\] Here and below, the dependence of the discrete operators on the temporal mesh, the fractional order, and the coefficient \(\mu\) is suppressed when no confusion can arise.

Since the fully discrete scheme is linear with respect to both the solution and the source, the terminal vector admits the decomposition \[\label{eq:affine95split} \mathbf{u}^M(\mathbf{u}_{0,h};f) = \mathbf{u}^M(\mathbf{u}_{0,h};0) + \mathbf{u}^M(\mathbf{0};f).\tag{71}\] Define the forced terminal contribution by \[\label{eq:forced95terminal} \mathbf{g}_h := \mathbf{u}^M(\mathbf{0};f),\tag{72}\] and define the homogeneous discrete forward operator \[\mathbf{F}_h\in\mathbb{R}^{(N-1)\times(N-1)}\] through \[\label{eq:homogeneous95forward} \mathbf{u}^M(\mathbf{v};0) = \mathbf{F}_h\mathbf{v}, \qquad \mathbf{v}\in\mathbb{R}^{N-1}.\tag{73}\] Consequently, \[\label{eq:discrete95affine95map} \mathbf{u}^M(\mathbf{u}_{0,h};f) = \mathbf{F}_h\mathbf{u}_{0,h}+\mathbf{g}_h.\tag{74}\]

For exact data, the discrete inverse relation is \[\label{eq:lin95system95inv95exact} \mathbf{F}_h\mathbf{u}_{0,h} = \boldsymbol{\psi}_h-\mathbf{g}_h.\tag{75}\] For perturbed measurements, we define \[\label{eq:discrete95data95vector} \mathbf{d}_h^\delta := \boldsymbol{\psi}_h^\delta-\mathbf{g}_h\tag{76}\] and seek an approximation to the solution of \[\label{eq:lin95system95inv} \mathbf{F}_h\mathbf{u}_{0,h} \approx \mathbf{d}_h^\delta.\tag{77}\]

4.1.0.1 Construction of the discrete forward operator.

Let \[\{\mathbf{e}^{(m)}\}_{m=1}^{N-1}\] be the canonical basis of \(\mathbb{R}^{N-1}\). For each \(m=1,\dots,N-1\), we solve the homogeneous fully discrete problem with \[\mathbf{u}^0=\mathbf{e}^{(m)}, \qquad f\equiv0,\] and record the corresponding terminal vector \[\mathbf{k}^{(m)} := \mathbf{u}^M(\mathbf{e}^{(m)};0).\] By linearity, \(\mathbf{k}^{(m)}\) is the \(m\)-th column of \(\mathbf{F}_h\), and therefore \[\label{eq:Fh95assembly} \mathbf{F}_h = \bigl[ \mathbf{k}^{(1)} \;\mathbf{k}^{(2)} \;\cdots \;\mathbf{k}^{(N-1)} \bigr].\tag{78}\]

This construction requires \(N-1\) homogeneous forward solves and is computationally practical for the moderate spatial dimensions used in the numerical experiments. For larger systems, explicit assembly may be replaced by a matrix-free iterative method based on products with \(\mathbf{F}_h\) and its transpose. Such an implementation generally requires the corresponding discrete adjoint operator.

4.2 Tikhonov regularisation↩︎

The continuous inverse problem is uniquely solvable and satisfies the stability estimate established in Section 2. Nevertheless, the discrete reconstruction may be affected by measurement noise, discretisation errors, and finite-precision arithmetic. We therefore use Tikhonov regularisation as a numerical stabilisation procedure.

For grid vectors, let \[(\mathbf{v},\mathbf{w})_h := h\sum_{i=1}^{N-1}v_iw_i, \qquad \|\mathbf{v}\|_h := (\mathbf{v},\mathbf{v})_h^{1/2}.\] The zeroth-order Tikhonov approximation is defined as the unique minimiser of \[\label{eq:tikhonov95fun} \mathcal{J}_\lambda(\mathbf{v}) = \frac{1}{2} \|\mathbf{F}_h\mathbf{v}-\mathbf{d}_h^\delta\|_h^2 + \frac{\lambda}{2}\|\mathbf{v}\|_h^2, \qquad \lambda>0.\tag{79}\] Since both terms contain the same factor \(h\), the first-order optimality condition is equivalent to \[\label{eq:normal95eq} \left( \mathbf{F}_h^\top\mathbf{F}_h+\lambda I \right) \mathbf{u}_{0,h}^{\lambda,\delta} = \mathbf{F}_h^\top\mathbf{d}_h^\delta.\tag{80}\] The coefficient matrix in 80 is symmetric positive definite for every \(\lambda>0\); hence, the minimiser is unique.

For the moderate dimensions considered in the numerical experiments, 80 can be solved by a dense Cholesky solver. From the viewpoint of numerical linear algebra, an equivalent augmented least-squares formulation, \[\begin{bmatrix} \mathbf{F}_h\\ \sqrt{\lambda}\,I \end{bmatrix} \mathbf{u}_{0,h}^{\lambda,\delta} \approx \begin{bmatrix} \mathbf{d}_h^\delta\\ \mathbf{0} \end{bmatrix},\] may instead be solved by QR factorisation or the singular value decomposition. This avoids the explicit squaring of the condition number associated with the formation of \(\mathbf{F}_h^\top\mathbf{F}_h\).

4.3 Choice of the regularisation parameter↩︎

The regularisation parameter \(\lambda\) determines the balance between fidelity to the measured terminal data and suppression of unstable or noise-dominated components. A value of \(\lambda\) that is too small may lead to excessive sensitivity to measurement and discretisation errors, whereas a value that is too large introduces an unnecessarily strong regularisation bias.

In the numerical experiments, \(\lambda\) is selected using the L-curve criterion. Let \[\Lambda = \{\lambda_1,\dots,\lambda_J\}\] be a logarithmically distributed set of positive candidate parameters. For each \(\lambda\in\Lambda\), we compute \(\mathbf{u}_{0,h}^{\lambda,\delta}\) from 80 and evaluate the residual and solution norms \[\begin{align} \rho(\lambda) &:= \left\| \mathbf{F}_h\mathbf{u}_{0,h}^{\lambda,\delta} -\mathbf{d}_h^\delta \right\|_h, \tag{81} \\ \eta(\lambda) &:= \left\| \mathbf{u}_{0,h}^{\lambda,\delta} \right\|_h. \tag{82} \end{align}\] The L-curve is the parametric curve \[\label{eq:lcurve} \lambda \longmapsto \bigl( \log\rho(\lambda), \log\eta(\lambda) \bigr).\tag{83}\] The selected parameter is taken near the corner of this curve, where the transition from the data-fitting regime to the regularisation-dominated regime occurs. In the discrete implementation, the corner is identified by locating a candidate parameter near the maximum curvature of the log–log L-curve.

For the numerical examples considered below, representative L-curves indicated a parameter of order \(10^{-10}\) for the noise-free tests and of order \(10^{-6}\) for the noisy-data tests. Accordingly, we use \(\lambda=10^{-10}\) in the noise-free validation experiments and \(\lambda=10^{-6}\) in the experiments with perturbed final-time data. These are empirical choices for the present test problems and are not intended as universal regularisation parameters. The values are kept fixed within each group of experiments in order to facilitate comparisons between different fractional orders and test configurations.

4.4 Reconstruction procedure↩︎

Algorithm 1 summarises the reconstruction procedure.

4.4.0.1 Verification and error measures.

To verify the reconstruction, we solve the forward problem once more using the reconstructed initial state and define the corresponding terminal state by \[\widehat{\boldsymbol{\psi}}_h^{\lambda,\delta} := \mathbf{u}^M \bigl( \mathbf{u}_{0,h}^{\lambda,\delta};f \bigr) = \mathbf{F}_h\mathbf{u}_{0,h}^{\lambda,\delta} +\mathbf{g}_h.\]

When the exact initial and terminal states are available, we set \[\mathbf{u}_{0,h}^{\mathrm{ex}} := R_hu_0, \qquad \boldsymbol{\psi}_h^{\mathrm{ex}} := R_h\psi.\] The reconstruction errors for the initial state are measured by \[\begin{align} E_{u_0,\infty} &:= \left\| \mathbf{u}_{0,h}^{\lambda,\delta} - \mathbf{u}_{0,h}^{\mathrm{ex}} \right\|_\infty = \max_{1\le j\le N-1} \left| u_{0,j}^{\lambda,\delta} -u_0(x_j) \right|, \tag{84} \\ E_{u_0,2} &:= \left\| \mathbf{u}_{0,h}^{\lambda,\delta} - \mathbf{u}_{0,h}^{\mathrm{ex}} \right\|_{2,h} \notag\\ &= \left( h\sum_{j=1}^{N-1} \left| u_{0,j}^{\lambda,\delta} -u_0(x_j) \right|^2 \right)^{1/2}. \tag{85} \end{align}\]

The corresponding errors in the reconstructed terminal state are defined by \[\begin{align} E_{\psi,\infty} &:= \left\| \widehat{\boldsymbol{\psi}}_h^{\lambda,\delta} - \boldsymbol{\psi}_h^{\mathrm{ex}} \right\|_\infty = \max_{1\le j\le N-1} \left| \widehat{\psi}_j^{\lambda,\delta} -\psi(x_j) \right|, \tag{86} \\ E_{\psi,2} &:= \left\| \widehat{\boldsymbol{\psi}}_h^{\lambda,\delta} - \boldsymbol{\psi}_h^{\mathrm{ex}} \right\|_{2,h} \notag\\ &= \left( h\sum_{j=1}^{N-1} \left| \widehat{\psi}_j^{\lambda,\delta} -\psi(x_j) \right|^2 \right)^{1/2}. \tag{87} \end{align}\] Here, \(\|\cdot\|_{2,h}\) denotes the discrete \(L^2(0,l)\)-norm.

5 Numerical validation and experiments↩︎

In this section, we present numerical experiments that demonstrate the practical behaviour of the proposed regularised reconstruction of the initial state \(u_0(x)=u(x,0)\) from the final-time measurement \(\psi(x)=u(x,T)\). All computations employ the direct solver from Section 3 and the discrete inverse formulation from Section 4.

5.1 Test setting and implementation details↩︎

Let \(0<\alpha<1\) and \((x,t)\in[0,l]\times[0,T]\). We consider the manufactured solution \[\mu(t)=1+t, \qquad u(x,t)=\bigl(1+t^{\alpha+1}\bigr)\sin\!\Bigl(\frac{\pi x}{l}\Bigr),\] together with the forcing term chosen so that \(u\) satisfies the time-fractional pseudo-parabolic model, namely \[f(x,t)=\Bigl[ \Gamma(\alpha+2)\,t +\Bigl(\frac{\pi}{l}\Bigr)^2\bigl(1+t^{\alpha+1}\bigr) +\mu(t)\Bigl(\frac{\pi}{l}\Bigr)^2(\alpha+1)t^\alpha \Bigr] \sin\!\Bigl(\frac{\pi x}{l}\Bigr).\] Accordingly, the final-time measurement and initial state are \[\psi(x)=u(x,T) =\bigl(1+T^{\alpha+1}\bigr)\sin\!\Bigl(\frac{\pi x}{l}\Bigr), \qquad u_0(x)=u(x,0) =\sin\!\Bigl(\frac{\pi x}{l}\Bigr).\]

The final-time measurement is sampled at the interior spatial nodes, as described in Section 4, and the homogeneous Dirichlet boundary conditions are imposed at all time levels.

5.2 Noise-free validation↩︎

We begin with a noise-free study to verify that the direct solver and the reconstruction procedure are mutually consistent. In this setting, the final-time measurement is sampled from the exact formula, i.e., \[\psi_i=\psi(x_i), \qquad i=1,\dots,N-1.\]

We fix \(T=1\) and refine the space–time grids according to \[(N,M)\in\{(50,50),(100,100),(200,200),(400,400)\}.\] The reconstruction is performed using the procedure described in Section 4. In accordance with the L-curve parameter-choice procedure described in Subsection 4.3, we set \(\lambda=10^{-10}\) in the noise-free experiments.

The maximum-norm and discrete \(L^2\) errors for \(u_0\) and the corresponding final-time errors for \(\psi\) are listed in Table 1. The decay of the initial-state errors under refinement confirms the convergence of the overall discrete reconstruction procedure. The final-time errors remain close to the regularisation level, demonstrating consistency between the inverse reconstruction and the direct solver.

Figure 2 shows representative reconstructions for several values of \(\alpha\) on the fixed grid \(N=M=50\). For each \(\alpha\), the left panel compares the exact initial state \(u_0(x)\) with its regularised reconstruction \(\widehat{u_0}(x)\), while the right panel compares the exact final-time measurement \(\psi(x)\) with the final-time state \(u(x,T;\widehat{u_0})\) obtained by forward propagation of the reconstructed initial state. The close agreement in both panels confirms the consistency of the inversion procedure and the direct solver.

Figure 3 compares the surface plot of the exact solution \(u(x,t)\) with the reconstructed numerical approximation obtained from the final-time measurement over the full time interval for several values of \(\alpha\). It also displays the corresponding absolute error surfaces. All simulations use the fixed grid \(N=M=200\).

5.3 Reconstruction from noisy final-time data↩︎

In practice, the final-time measurement \(\psi(x)\) is corrupted by measurement noise. To assess robustness, we perturb the discrete final-time vector by additive Gaussian noise with a prescribed relative level.

Let \(\boldsymbol{\psi}_h\in\mathbb{R}^{N-1}\) be the clean interior measurement. For \(\delta>0\), we set \[\label{eq:noise95model95alt} \boldsymbol{\psi}_h^{\delta} = \boldsymbol{\psi}_h + \sigma_\delta\,\boldsymbol{\xi}, \qquad \boldsymbol{\xi}\sim\mathcal{N}(\mathbf{0},I_{N-1}), \qquad \sigma_\delta := \delta\,\frac{\|\boldsymbol{\psi}_h\|_2}{\sqrt{N-1}}.\tag{88}\] With this choice, \[\frac{ \mathbb{E}\|\boldsymbol{\psi}_h^{\delta}-\boldsymbol{\psi}_h\|_2^2 }{ \|\boldsymbol{\psi}_h\|_2^2 } = \delta^2,\] so that \(\delta\) represents the root-mean-square relative noise level.

For each noisy data vector \(\boldsymbol{\psi}_h^\delta\), the initial state is reconstructed using the procedure in Section 4. We use the fixed grid \(N=M=100\) and regularisation parameter \(\lambda=10^{-6}.\) The errors \(E_{u_0,\infty}^{\delta}\), \(E_{u_0,2}^{\delta}\), \(E_{\psi,\infty}^{\delta}\), and \(E_{\psi,2}^{\delta}\) are computed using the definitions given in Section 4.

Representative results for \(\delta\in\{1\%,3\%,5\%\}\) are displayed in Table 2.

Table 2 quantifies the sensitivity of the reconstruction to perturbations in the final-time measurement \(\psi\). As expected for a backward problem, the reconstruction error increases as the noise level \(\delta\) grows. Nevertheless, the growth is controlled by the Tikhonov stabilisation: for \(\delta=1\%\), the recovered initial state remains accurate, with \(E_{u_0,2}^{\delta} \approx (3.7\text{--}4.1)\times10^{-2}, \quad E_{u_0,\infty}^{\delta} \approx (7.5\text{--}8.5)\times10^{-2}\) across all tested fractional orders. When the noise is increased to \(\delta=3\%\) and \(\delta=5\%\), the initial-state errors rise in a near-proportional manner. For example, \(E_{u_0,2}^{\delta}\) increases from approximately \(3.7\times10^{-2}\) to \(1.1\times10^{-1}\) and then to \(1.8\times10^{-1}\). This indicates stable dependence on the data perturbation rather than uncontrolled amplification.

A further observation is that the reported errors vary only mildly with \(\alpha\) for a fixed value of \(\delta\). This suggests that, for the present test, the regularised inversion is not overly sensitive to the fractional order and that the dominant factor affecting the accuracy is the measurement noise level. In addition, the terminal-state errors \(E_{\psi,\infty}^{\delta}\) and \(E_{\psi,2}^{\delta}\) remain comparable to the noise magnitude, showing that the reconstructed initial state produces forward solutions that remain close to the exact terminal state.

Figure 4 illustrates the influence of measurement noise on the reconstructed initial state. For the noise level \(\delta=1\%\), the reconstructed profile remains very close to the exact initial state. As the noise level increases to \(3\%\) and \(5\%\), the deviation becomes more visible, particularly in the amplitude of the reconstructed solution. Nevertheless, the main shape and spatial structure of \(u_0(x)\) are preserved in all three cases, and no significant spurious oscillations are observed. This behaviour confirms the stabilising effect of the Tikhonov regularisation and is consistent with the quantitative errors reported in Table 2

Figure 5 shows surface plots of the reconstructed numerical solution \(u(x,t)\), obtained from the final-time measurement, for several noise levels over the full time interval with \(\alpha=0.5\).

Figures 4 and 5 complement these quantitative results. Even for moderate noise, the reconstructed profiles preserve the main shape of the exact initial state, while higher noise levels primarily introduce smooth amplitude distortions rather than spurious oscillations. Overall, the experiments confirm that the proposed finite-difference-based forward solver combined with Tikhonov regularisation yields a numerically robust framework for recovering \(u_0\) from noisy final-time data.

6 Conclusion↩︎

In this work, we have provided a comprehensive analysis of the backwards-in-time problem for a time-fractional pseudo-parabolic equation with time-dependent coefficients. From a theoretical standpoint, we established the existence and uniqueness of a classical solution by means of a spectral expansion and the theory of weakly singular Volterra integral equations. We also derived stability estimates showing the continuous dependence of the initial data on the final-time measurements. Numerical treatments were developed for both the forward and inverse problems. For the direct problem, we constructed a fully discrete finite difference scheme using the \(L1\) approximation on graded meshes for the Caputo derivative. We rigorously proved that the scheme is unconditionally stable and convergent, achieving an optimal convergence rate. Building upon this efficient forward solver, we formulated the inverse reconstruction as a minimisation problem regularised by the Tikhonov method. Numerical experiments confirmed the theoretical findings, demonstrating that the proposed algorithm effectively recovers the unknown initial state even in the presence of noise in the terminal data. Future research may extend this framework to multi-dimensional domains or nonlinear source terms, where the interplay between memory effects and nonlinearity poses further challenges.

Data availability↩︎

Data will be made available on request.

Acknowledgements↩︎

This work was supported by the Ministry of Science and Higher Education of the Republic of Kazakhstan (No. AP27508473).

References↩︎