An inverse source problem for a pseudoparabolic equation with memory

1.   Introduction

Statement of the problem

Let $Q_{T}: = \Omega\times [0, T]$ be a bounded rectangle, where $\Omega = (0, l)$ and $0 < l, T < \infty$. Let us consider the inverse problem of finding a pair $(y, f)$ that satisfies the following one-dimensional linear pseudoparabolic equation with a memory

subject to the initial condition

the boundary conditions

and the final overdetermination condition

where $\overline{\Omega}: = [0, l]$, $y_{0}(x), a(x), g(x, t), h(x, t)$, and $K(t)$ are given functions, and $\alpha, \beta, \gamma$ are given real numbers such that $\alpha > 0$, $\beta\geq0$, and $-\infty < \gamma < +\infty$. For simplicity, we take $|\gamma| = 1$. The studied inverse problem given by (1.1)–(1.4) involves determining, in addition to $y(x, t)$, a right-hand side coefficient $f(x)$ under the given observation (1.4). There are problems that require the determination of certain parameters, and there are obtainable solutions of direct problems when additional information is given. Such problems are usually called inverse problems. For instance, in the mathematical modeling of certain processes, certain physical parameters affecting the process may be unknown and inaccessible for direct measurement, such as when the environment is underground, or in a place with high temperatures. In these types of cases, additional information about the phenomenon can be given, allowing the recovery of the unknown parameters. Such information is usually modeled by an average value of the sought solution, or by the value at the final time. Pseudoparabolic equations with nonlocal terms such as (1.1) have various physical applications; for instance, they appear in the mathematical descriptions of heat conductions and viscous flows in materials with memory, electric signals in telegraph lines with nonlinear damping ^[1], non-Newtonian fluid dynamics ^[35], bidirectional nonlinear shallow water waves ^[32], the velocity evolution of ion-acoustic waves in collisionless plasma ^[30], population dynamics and plasma physics ^[21,33], and so on. Moreover, the physical meaning of the memory term (integral term) in the equation is important; for instance, in non-Newtonian fluid models ^[35], the elastic property of fluids has an indication. However, this term, along with its physical importance, causes mathematical difficulties in their analytical and numerical analysis, and requires additional techniques. In the absence of the memory term (i.e., $K(t) = 0$), the Eq (1.1) becomes a classical linear pseudoparabolic equation

which are used in the mathematical modeling of phenomena in fields of thermodynamics and filtration, etc.; see ^[7,8,20,28].

There are many works on the study of direct and inverse problems for these kinds of pseudoparabolic equations. Most of the papers on inverse problems for pseudoparabolic equations such as (1.5) (without a memory term) address inverse problems, which involve finding a coefficient for the right-hand side that depends on the time variable $t$; see ^{[5,6,17,23,25,31]} and others. However, the inverse problems that involve determining the right-hand side coefficients that depends on spatial variables for pseudoparabolic equations have not been studied extensively. For instance, the existence and uniqueness of classical solutions to the inverse problem for (1.5) with integral and final overdetermination conditions were established in ^[18] and ^[16], respectively. To the best of our knowledge, the inverse problems that involve determining a right-hand side coefficient that depends on the spatial variables for pseudoparabolic equations with memory, in particular, the current statement of the above inverse problem given by (1.1)–(1.4), have not yet been investigated, either analytically or numerically. Nevertheless, there are some works on the numerical analysis of the inverse problem for pseudoparabolic-type equations without a memory term; see ^{[11,12,13,14,15,26]}. For example, in ^[12,14], Huntul et al. numerically studied the inverse problems that involve reconstructing the unknown coefficients in a third-order pseudoparabolic equation from additional and nonlocal integral observations.

The rest of the paper is organized as follows. Section 2 is devoted to analytically studying the given inverse problem, i.e., the existence, uniqueness and stability of a strong generalized solution. The discretization and solution of the direct problem via an extended cubic B-spline (CBS) approach are given in Section 3, and Section 4 presents the stability analysis of the proposed scheme. The extended CBS direct solver based on the numerical method is coupled with the Tikhonov regularization method, as described in Section 5. In Section 6, the computational results for a test example are presented and discussed. Finally, Section 7 highlights the conclusions.

2.   Preliminaries

In this section, we reduce the original inverse problem given by (1.1)–(1.4) to another inverse problem with homogeneous initial conditions, whose solvability is equivalent to the solvability of an operator equation of the second kind. Throughout the entire paper, the following abbreviations are used for precisesness. The norms in the Lebesgue spaces $L^2(\Omega)$ and $L^2(Q_T)$ are denoted as follows:

and the inner product in $L^2(\Omega)$ is denoted by

For the definitions, notations of the function spaces and their properties, we refer the reader to the monographs ^[2,22]. Since the problem given by (1.1)–(1.4) is linear, we seek a solution to the original inverse problem given by (1.1)–(1.4) as follows:

where $\left \{ v(x, t), h(x, t) \right \}$ is a solution to the following direct problem:

while $\left \{ u(x, t), f(x) \right \}$ is a solution to the following inverse problem:

where $\phi(x) = a(x)-v(x, T)$, and $v(x, T)$ is a trace of a solution of the direct problem given by (2.1)–(2.4) at $T$. According to ^[28,35], the direct problem given by (2.2)–(2.4) has a unique solution

Consequently, to investigate the unique solvability of the original inverse problem given by (1.1)–(1.4), it is sufficient to study the unique solvability given by (2.5)–(2.8).

Definition 1. The pair of functions $\{{u}(x, t), {f}(x)\}$ is said to be a strong generalized solution of the inverse problem given by (2.5)–(2.8) if

and all relations given by (2.5)–(2.8) are satisfied almost everywhere in the corresponding domains.

2.1.   An equivalent operator equation

Now, we reduce the inverse problem given by (2.5)–(2.8) to an equivalent operator equation of the second kind for $f(x)$ in $L^2(\Omega)$. Here, we use the method used in the book ^[29] with few new techniques, i.e., here, we derive an operator equation by integrating the equation by $t$ in $[0, T]$ instead of taking $t = T$ in the equation as in ^[29]. The latter requires smoother solutions of the corresponding direct problem, such as $u_{xxt}\in C(0, T; L^2(\Omega))$, while $u_{xx}\in C(0, T; L^2(\Omega))$ is aready sufficient.

Let us fix an arbitrary function $f(x) \in L^2(\Omega)$. After the substitution of $f(x)$ into (2.5), the problem given by (2.5)–(2.7) can be considered to be the direct problem of determining $u(x, t)$, and it has a unique solution:

Thus, this established correspondence between $f(x)$ and $u(x, t)$ provides us with an operator

Let us consider the following operator equation of the second kind

with the operator $A$ defined by the relation

Here, $u(x, t)$ is the solution of the direct problem given by (2.5)–(2.7) with the right-hand side $F = f(x) g(x, t)$ corresponding to the fixed function $f(x) \in L^2\left (\Omega\right)$, $\eta(x)$ is a given function in $L^2\left (\Omega\right)$, and

Lemma 1. Let the following conditions be satisfied:

Then, the unique solvability of the inverse problem given by (2.5)–(2.8) is equivalent to the unique solvability of the operator (2.9), i.e., the inverse problem given by (2.5)–(2.8) is uniquely solvable, if and only if, the operator (2.9) is uniquely solvable.

Proof. (1) Let $(u, f)$ be the unique solution of the inverse problem given by (2.5)–(2.8). Integrating (2.5) by $t$ from $0$ to $T$ and using (2.6), (2.8) and (2.12), we have

It follows from (2.14) and the assumptions (2.11), (2.13), and (2.10) that (2.9) holds with $\eta$ defined by (2.13). This step completes the first part of the proof.

(2) Now, suppose that the operator equation (2.9) is uniquely solvable; let $f\in L^2 \left(\Omega \right)$ be a solution. Substituting $f\left(x\right)$ into (2.5), we consider (2.5)–(2.7) to represent a direct problem with a given right-hand side $F\left(x, t\right) = f\left(x\right)g\left(x, t\right)$, and it has a unique solution $u\left(x, t\right)\in L^2 \left(0, T; W^{2, 2} \left (\Omega \right)\cap W_0^{2, 1}\left (\Omega \right) \right)\cap W^{2, 1} \left(0, T; W^{2, 2}\left (\Omega \right) \right).$ Therefore, to prove that such a constructed pair $(u, f)$ is a strong solution of the inverse problem given by (2.5)–(2.8), it is sufficient to prove that the function ${u}$ satisfies the final overdetermination condition (2.8). Let us make an assumption for contradiction, i.e., the overdetermination condition does not hold. Suppose that

where $\phi_1(x)\neq\phi(x)$ for all $x\in\Omega$, and $\phi_1(x)\in W^{2, 2}\left(\Omega \right)$. By (2.7), we obtain

Both sides of (2.5) are integrated by the variable $t$ from $0$ to $T$, and the new overdetermination condition (2.15) is applied. By setting (2.14) and using (2.10), we obtain

Subtracting (2.9) from (2.17) and applying the condition (2.16) yields the following boundary problem for the function $E(x) = \phi(x)-\phi_1(x)$

By multiplying both sides of (2.18) by $E(x)$, and integrating from $0$ to $l$ by the variable $x$, we have

which implies that $\phi(x)\equiv \phi_1(x)$ in $\Omega$. Therefore, the function $u$ satisfies the overdetermination condition in (2.8), i.e., the pair $\left(u, f\right)$ is a unique solution of the inverse problem given by (2.5)–(2.8). □

Thus, by Lemma 1, we shall study the solvability of the operator (2.10).

Theorem 1. Assuming that in addition to (2.11)–(2.13), the following conditions are satisfied:

Then, the inverse problem given by (2.5)–(2.8) has a unique solution.

Remark 1. The condition (2.23) will be satisfied if we at least choose $T$ to be sufficiently small, i.e., the inverse problem is uniquely solvable locally in time.

For example, if $K(t) = 1$, $l = 1$, $g(x, t) = \cos (t)+\alpha \pi^2\sin(t)+\beta \pi^2 \cos (t)-\pi^2\cos (t)+ \pi^2$, and $\phi(x) = \sin (T)\sin(\pi x)$ in (2.5)–(2.8), then the exact solution of (2.5)–(2.8) corresponding to these data is the pair of functions $u(x, t) = \sin(t)\sin(\pi x)$ and $f(x) = \sin(\pi x)$. In this case, it is easy to check that the condition (2.23) holds if $\alpha = 1$, $\beta = 0.01$, and $T = 0.5$ (here $K_g = |g(x, t = 0.5)|\approx 6.898$, $g_0 = \left|\int\limits_0^{0.5}g(x, t)dt\right|\approx 11.390$, $K_0 = \sqrt{T} = 0.707$) with $\mu = 0.824$ or if $\alpha = 1$, $\beta = 2.5$, and $T = 1$ (here $K_g = |g(x, t = 0.55)|\approx28.47$, $g_0 = \left|\int\limits_0^{1}g(x, t)dt\right|\approx44.27$, $K_0 = \sqrt{T} = 1$) with $\mu = 0.9936$.

Proof. To prove this theorem, Lemma 1 is sufficient to prove the unique solvability of the operator (2.9). To this end, it is sufficient to show that the operator $A$ is defined by (2.10) and is a contraction operator in $L^2\left(\Omega\right)$. Then, by the fixed point principle, the operator (2.9) has a unique solution in $L^2(\Omega)$.

Using Holder's and Young's inequalities and the conditions (2.20)–(2.23), we have

Next, we obtain some estimates for $\left\|u_{1xx}-u_{2xx}\right\|_{L^{2}(Q_T)}^2$. By multiplying both sides of (2.5) by $u_{xx}: = u_{1xx}-u_{2xx}$, integrating over $\Omega$ and applying the formula of integrating by parts, we have

where

Now, we estimate the right hand side of (2.25)

Substituting this into (2.25), we obtain

and integrating (2.27) by $s$ from $0$ to $t$, we obtain

Omitting the first and second terms, we obtain the following from (2.28):

which, after applying Grönwall's lemma and taking $t = T$ as a result, yields

Next, omitting the first and third terms on the left-hand side of (2.28), and integrating by $s$ from $0$ to $t$, we obtain

which, by Grönwall's lemma, gives the following inequality with $t = T$

The inequalities (2.30) and (2.32) yield

Plugging the last equation with $u: = u_1-u_2$ and $f: = f_1-f_2$ into (2.24), we obtain

where $\mu: = \frac{K_gT\left(\alpha+K_0\sqrt{T}\right)}{g_0\sqrt{\alpha}} \min\left\{\sqrt{\frac{2}{\alpha}}e^{\frac{K_0^2T}{\alpha^2}}, \ \sqrt{\frac{T}{\beta}} e^{\frac{K_0^2T^2}{2\alpha\beta}}\right\}.$ Thus, by (2.23) with $\mu < 1$, the operator $A$ is a contraction, and this completes the proof of Theorem 1. □

2.2.   Stability results for the inverse problem

In this subsection, we establish the stability results for the inverse problem given by (2.5)–(2.8).

Theorem 2. If the conditions (2.10)–(2.13) hold, then the strong solution of the inverse problem given by (2.5)–(2.8) is stable, i.e., $\left(u_i, f_i\right), \, i = 1, 2$ are two solutions of the inverse problem given by (2.5)–(2.8) corresponding to the input data $\left(\phi_i, g_i\right), i = 1, 2$, then, there exists a constant $M$ such that

Remark 2. If $g_1 = g$, $g_2 = 0$, and $\phi_1 = \phi$, $\phi_2 = 0$, then we obtain the estimates for strong solutions to the inverse problem given by (2.5)–(2.8):

Remark 3. If $g_1 = g_2$ and $\phi_1 = \phi_2$, then it follows from (2.35) that the uniqueness of strong solutions to the inverse problem given by (2.5)–(2.8) holds.

Proof. Let $\left\{u_, f_i\right\}, i = 1, 2$ be two solutions of the inverse problem given by (2.5)–(2.8) corresponding to the input data $\left(\phi_i, g_i\right), i = 1, 2$. Then, we apply their differences:

We obtain the following problem:

Both sides of (2.37) are multiplied by $u_{xx}$ and integrated by $x$ over $\Omega$. Then, estimating the right-hand side, we have

As we have obtained (2.33) from (2.27), we obtain

where

Then, integrating (2.37) by $t$ from $0$ to $T$, we have

Taking the norm of both sides of the last equality and combining the result with (2.42), we obtain

It follows that

Next, combining (2.41), (2.42), and (2.45), we obtain

Now, taking the derivative of (2.37) with respect to ${u}$, we obtain

By estimating the terms on the right-hand side using the Holder and Young inequalities and combining the result with (2.45), we have

Next, multiplying (2.37) by $u_{xxt}$ and integrating over $\Omega$, we obtain

Likewise, it follows that

The inequalities (2.48), (2.46), (2.49), and (2.45) give (2.35). □

3.   The numerical scheme for the direct problem

In this section, the numerical scheme for the direct problem is described. The current approach employs the first-order backward Euler scheme to handle the time derivatives and extended CBS functions to approximate the unknown function and its spatial derivatives. Initially, the temporal domain $[0, T]$ is sliced into $N$ uniform subintervals $[t_{k}, t_{k+1}]$ such that $t_{k} = k\times \Delta t$, $k = 0, 1, 2, \cdots, N-1$, with $t_{0} = 0, t_{N} = T$ and $\Delta t = (t_{N}-t_{0})/N$. At $t = t_{k+1}$, using a theta-weighted scheme, the problem (2.5) is discretized as follows:

Moreover, the integral term in (3.1) is discretized as follows ^[9,36]:

where

Applying the above result in (3.1), we obtain

Setting $\theta = 1/2$, (3.3) takes the following shape:

By rearranging the terms, we obtain the semidiscretized form from (2.5):

Now, we subdivide the spatial domain $[0, l]$ into $M$ uniform intervals $[x_{i}, x_{i+1}]$ such that $x_{i} = i\times \Delta x$, $i = 0, 1, 2, \cdots, M$, with $x_{0} = 0, x_{M} = l$ and $\Delta x = (x_{M}-x_{0})/M$. We suppose that the extended CBS solution for the proposed problem at $t = t_{k+1}$ is given by

where $S_{p}(x)$, i.e., the extended B-spline functions, are defined as follows ^[3,34]:

Using (3.6) and (3.7), the unknown function and its first two spatial derivatives at point $(x_{i}, t_{k})$ can be approximated as follows ^[4,27]:

where $-8\leq \zeta\leq 1$ and $\varepsilon_{i}$ denotes $M+3$ unknown control points to be computed by the collocation of (2.5) and (2.7). Hence, the full discretized form of (2.5) is obtained as follows:

where

Now, applying (3.8) and (3.10) in (3.11), we have

For $i = 0, 1, 2, \cdots, M$, it follows that (3.12) gives $M+1$ linear algebraic equations involving $M+3$ unknowns. Two more equations are required to obtain a unique solution to the proposed problem, and they can be extracted from the given boundary constraints of (2.7) as follows:

Consequently, we have a system of $M+3$ equations, i.e., (3.12)–(3.14), involving $M+3$ unknown control points, $\varepsilon_{-1}(t_{k+1}), \varepsilon_{0}(t_{k+1}), \cdots, \varepsilon_{M+1}(t_{k+1})$. This system can be written in matrix form as follows:

where

The approximate solution at $t = t_{k+1}$ is obtained by solving the system of (3.15) for $\varepsilon_{i}$ and then substituting their values into (3.6). However, the control points at the initial time level must be identified before beginning any computation using (3.15). For this, the given initial condition (2.6) is transformed into the following system of equations:

where

4.   Stability analysis of the extended CBS scheme

In this section, the von Neumann approach is used to demonstrate the stability of the presented numerical method for direct computation. For simplicity, we take the full discretized form of (2.5):

where $U_{m}^{k} = U(x_{m}, t_{k})$. Applying (3.8) and (3.10) in (4.1), we obtain

where

In the context of a von Neumann stability analysis, the Fourier modes represent specific sinusoidal components resulting from the discretization process. The analysis evaluates the evolution of these Fourier modes by examining the growth factor associated with each mode during the iterative process. The growth factor in this context indicates how much an individual Fourier mode amplifies or diminishes at each iteration of the numerical algorithm. This approach helps one to determine whether the numerical scheme amplifies the errors or perturbations introduced during discretization, potentially impacting the overall stability of the solution. We assume that $\yen_{m}^k$ represents the error between the exact and approximate values of the growth factor in Fourier mode. Hence, using (4.2), the error equation is formulated as follows:

Moreover, the initial and boundary conditions are satisfied by the error equation. Now, we introduce the mesh function in Fourier form:

Hence, $\yen^{k}(x)$ in terms of the Fourier series is given by

where

By applying the norm, we obtain

Hence, by using Parseval's equality ^[19], we have

Now, we substitute $\yen_m^{k} = \varsigma_{k} e^{i \gamma m \Delta x}$ into (4.3) to obtain the following relation

where $i = \sqrt{-1}$ and $\gamma = {\frac{2\pi m }{\ell}}$. After simplification, the last equation takes the following form:

Hence,

where

Now, we demonstrate, using mathematical induction, that $|\varsigma_{k}|\leq |\varsigma_0|$ for all $k$. For $k = 0$, it follows that (4.9) takes the following form:

Assuming that the required result satisfies that $|\varsigma_{k}|\leq |\varsigma_{0}|$, we proceed as follows:

Hence,

From (4.6) and (4.10), we conclude that

Hence, the proposed scheme is numerically stable.

5.   Numerical solution of the inverse problem

In this section, our goal is to obtain simultaneously stable determinations of $f(x)$ and $u(x, t)$, satisfying (2.5)–(2.8). The inverse problem can be formulated as a nonlinear minimization of the Tikhonov regularization function given by

where $u$ satisfies the forward problem given by (2.5)–(2.7) for given $f(x)$, and $\delta$ is the nonnegative penalty parameter. The discretized form of above function is given by

The minimization of the objective function $J$ can be performed by using the MATLAB toolbox routine lsqnonlin, which does not require the user to supply the gradient of the objective function, ^[10,24]. This routine attempts to find the minimum of a sum of squares that are subject to simple bounds on the variables, starting from an initial guess.

6.   Computational results

In this section, we present a numerical test example to illustrate the accuracy and stability of the numerical methods based on the extended CB-spline technique as described in Section 3 combined with the minimization of the objective function $J$, as described in Section 5. We chose to employ the root mean square error (RMSE) in order to assess the accuracy of the numerical results, defined by

For simplicity, we take $l = 1$. The upper and lower bounds for $f(x)$ were assumed to be $10^{2}$ and $-10^2$, respectively. The inverse problem given by (2.5)–(2.8) was solved for both exact and perturbed data. The perturbed data were numerically simulated as follows:

where $\epsilon_i$ denotes the random variables generated from a Gaussian normal distribution with mean zero and standard deviation $\sigma$ given by

where $p$ represents the percentage of noise. For the perturbed data described by (6.2), $\phi(x_i)$ is replaced by $\phi^{\epsilon}(x_i)$ in (5.2). Let us investigate the inverse problem given by (2.5)–(2.8) with an unknown space-dependent coefficient $f(x)$ given by

The analytical solution $u(x, t)$ is given by

and the target output (2.8) is given by

with the following input data:

First, it can easily be checked that with this data, the condition (2.23) holds if $\alpha = 1$, $\gamma = 1$, $\beta = 24$, and $T = 1$ (here $K_g = |g(x, t = 1)|\approx 682.6$, $g_0 = \left|\int\limits_0^{1}g(x, t)dt\right|\approx 243.3$, $K_0 = \sqrt{T} = 0.657$) with $\mu = 0.99$. Hence, the uniqueness of the solution is guaranteed. The three dimensional visuals of approximate and the exact solutions are shown in Figure 1 for different values of $\Delta x$ and $\Delta t$. In this figure it can be seen that the numerical solution for $u(x, t)$ converges to the analytical solution (6.5), as the mesh size decreases.

Next, we chose to solve the inverse problem given by (2.5)–(2.8) by using the lsqnonlin minimization of the functional (5.2), with the initial guess for the vector $\underline{f}$ given by

We set $M = N = 80$ and began the investigation to determine the unknown space-dependent coefficient $f(x)$ and the solution function $u(x, t)$ in the case of exact input data, i.e., $p = 0$ in (6.3). The objective function (5.2) is depicted in Figure 2(a), where a monotonically decreasing convergence is achieved in about 20 iterations towards the low value of $O(10^{-26})$. Figure 2(b) shows the exact (i.e., (6.4)) and approximate solutions for $f(x)$ without the regularization parameter, i.e., $\delta = 0$ in (5.2). This figure shows an acceptable and stable accurate estimate for the coefficient $f(x)$, yielding RMSE$(f) =$5.6854E-04.

Next, we investigated the stability of the numerical solution with respect to various levels of $p\in\{0.01\%, 0.1\%\}$ noise in (6.3) included in the input data $\phi(x)$ in (6.2), in order to model the errors which are inherently present in any practical measurement. The exact (i.e., (6.4)) and approximate solutions for the unknown coefficient $f(x)$, with and without penalty parameter $\delta$, are shown in Figures 3 and 4. In Figures 3(a) and 4(a) it can be seen that as the noise $p$ is increased, the approximate results start to build up oscillations with RMSE$(f)\in \{0.300508, 3.058572\}$. Figures 3(b) and 4(b) illustrate the reconstructed coefficient $f(x)$ for various values of $\delta$, and one can observe that the most accurate solutions are achieved for $\delta\in\{10^{-9}, 10^{-8}\}$, yielding RMSE $(f)\in \{0.021935, 0.008473 \}$ for $p = 0.01\%$, and for $\delta\in\{10^{-8}, 10^{-7}\}$, yielding RMSE $(f)\in \{0.044552, 0.014704 \}$ for $p = 0.1\%$; see Table 1 for additional details. Finally, Figure 5 shows the absolute error between the analytical solution (6.5) and approximate solutions for $u(x, t)$ for various values of the penalty parameter $\delta$ for $p = 0.1\%$. In this figure it can be seen that the numerical solution is stable and furthermore, that its accuracy improves as $\delta > 0$.

7.   Conclusions

In this paper, the inverse source problem for a linear one-dimensional pseudoparabolic equation with a memory term has been investigated. The inverse problem that involves determining a right-hand side coefficient $f(x)$, which depends on a spatial variable and the primary dependent variable. The final overdetermined condition at the final time point has been used as additional information. The proposed inverse problem has been analyzed both theoretically and numerically. Under the appropriate assumptions on the data for the inverse problem, the existence, uniqueness and stability of a strong generalized solution have been established. For the numerical discretization, the extended CBS collocation technique has been employed as a direct solver for a specific example. The von Neumann stability analysis has also been proved. The resulting objective function was penalized by incorporating a Tikhonov regularization term to ensure stability of the numerical solutions of the inverse problem. The minimization process has been carried out iteratively through the utilization of the lsqnonlin routine in the MATLAB optimization toolbox. Finally, the generalization of the proposed CBS collocation scheme in a fourth-order pseudoparabolic equation without a memory term is an interesting topic for future research.

Use of AI tools declaration

The authors declare they have not used artificial intelligence tools in the creation of this article.

Acknowledgments

The authors extend their appreciation to the Deputyship for Research & Innovation, Ministry of Education in Saudi Arabia for funding this research work through the project number ISP-2024.

Conflict of interest

The authors declare that they have no competing interests.