Numerical Solutions of Variable Coefficient Higher-Order Partial Differential Equations Arising in Beam Models

Abstract

In this work, an efficient and robust numerical scheme is proposed to solve the variable coefficients’ fourth-order partial differential equations (FOPDEs) that arise in Euler–Bernoulli beam models. When partial differential equations (PDEs) are of higher order and invoke variable coefficients, then the numerical solution is quite a tedious and challenging problem, which is our main concern in this paper. The current scheme is hybrid in nature in which the second-order finite difference is used for temporal discretization, while spatial derivatives and solutions are approximated via the Haar wavelet. Next, the integration and Haar matrices are used to convert partial differential equations (PDEs) to the system of linear equations, which can be handled easily. Besides this, we derive the theoretical result for stability via the Lax–Richtmyer criterion and verify it computationally. Moreover, we address the computational convergence rate, which is near order two. Several test problems are given to measure the accuracy of the suggested scheme. Computations validate that the present scheme works well for such problems. The calculated results are also compared with the earlier work and the exact solutions. The comparison shows that the outcomes are in good agreement with both the exact solutions and the available results in the literature.

1. Introduction

Small static and dynamic deflection problems can be observed properly in linear theory. For the determination of large (dynamic and static) deflection, linear theory is not beneficial and requires an accurate analysis. Linear theory admits inexact curvature in the study of beam deflections. The pioneering work in the field of thin beam theory was carried out by Bernoulli. Jacob Bernoulli then studied elastic theory and showed that the bending moment and curvature are both proportional. Later on, the Bernoulli theory was extended to loaded beams by Leonhard Euler.

The PDE of a thinand long beam is known as the Euler–Bernoulli model. The solution of this model shows the shortest distance (transverse vibration) from the beginning position in which stress and strain are linearly related. The mathematical form of the Euler–Bernoulli model [1] is described as follows:

$$\gamma \left(x\right){\partial }_{tt}u(x,t)+{\partial }_{xx}\left[\mathcal{C}\left(x\right){\partial }_{xx}u(x,t)\right]=\mathcal{E}(x,t),0<x<1,0<t\le T,$$

(1)

with appropriate initial and boundary conditions:

$$u(x,0)={\mathrm{{\rm Y}}}_1\left(x\right),{\partial }_{t}u(x,0)={\mathrm{{\rm Y}}}_2\left(x\right),0\le x\le 1,$$

(2)

$$\begin{array}{ccc}& & \left\{\begin{array}{c}u(0,t)={\epsilon }_o\left(t\right),{\partial }_{x}u(0,t)={\delta }_o\left(t\right),0<t\le T,\hfill \\ u(1,t)={\epsilon }_1\left(t\right),{\partial }_{x}u(1,t)={\delta }_1\left(t\right),0<t\le T,\hfill \end{array}\right.\hfill \end{array}$$

(3)

or

$$\begin{array}{ccc}& & \left\{\begin{array}{c}u(0,t)={\zeta }_o\left(t\right),{\partial }_{xx}u(0,t)={\mu }_o\left(t\right),0<t\le T\hfill \\ u(1,t)={\zeta }_1\left(t\right),{\partial }_{xx}u(1,t)={\mu }_1\left(t\right),0<t\le T,\hfill \end{array}\right.\hfill \end{array}$$

(4)

where the notations ${\partial }_{tt}u(x,t),{\partial }_{xx}u(x,t)$ stand for the second derivative with respect to the time and space variables, respectively. Equation (1) is the variable coefficients’ FOPDE in which $u(x,t)$ represents the beam displacement, $\gamma \left(x\right)$ is the mass per unit length, $\mathcal{C}\left(x\right)$ shows the bending stiffness of the beam, $\mathcal{E}(x,t)$ is the source term, and L is the total length. This kind of equation has widespread applications in robotics designs and large flexible space structures [2,3]. Several analytical methods have been applied to derive the closed-form solution of the governing equation. Wazwaz [4] used the Adomian decomposition method to solve the variable coefficients’ FOPDEs. Lieu [5] applied He’s variational iteration technique to explain free vibration in an Euler–Bernoulli beam. For more generic cases of initial conditions, analytical solutions are quite complicated; therefore, researchers are constantly trying to focus on the numerical solutions. A variety of numerical methods via finite difference schemes have been developed for the solution of different forms of Equation (1) such as Jain et al. [6], Evans [7], Conte [8], Richtmyer [9], and Cranial [10]. Evans et al. [11] established a stable computational method using the hopscotch algorithm. Aziz et al. [12] proposed a three-level scheme using a parametric quintic spline to solve the FOPDEs. Rashidinia [13] implemented a three-level implicit scheme coupled with a sextic spline to solve fourth-order equations. Recently, Imtiaz et al. [14] numerically studied fractional-order Korteweg–de Vries and Burger’s equations via the meshless method. Senol and his co-authors [15] investigated the Coudrey–Dodd–Gibbon–Sawada–Kotera equation with three different methods. Akinyemi [16] solved $(1+3)-$dimensional fractional reaction–diffusion trimolecular models. The same author [17] and his collaborators computed the numerical solutions of coupled nonlinear Schrödinger–Korteweg-de Vries and Maccari systems numerically. Jiwari [18] applied barycentric rational interpolation and local radial basis function algorithms for the multi-dimensional Sine–Gordon equation.

In the past few years, wavelet-based approximation techniques gained great importance for solving PDEs [19,20]. In all kinds of wavelets, the simplest family is the Haar wavelet (HW), which comprises rectangular box functions. The usage of these wavelet attracted more researchers because of its easy implementation and achievements with good results. The HW based on Haar functions introduced by Alfred Haar in 1910, which are quite simple mathematically, are discontinuous at breaking points of the interval and, hence, not differentiable. Due to this reason, the direct implementation of the HW for the solution of differential equations is not possible. For this purpose, Cattani [21,22] used interpolating splines to remove this ambiguity. Similarly, an alternative idea was given by Chen and Hasio [23]. They suggested approximating the maximum-order derivative with finite HW series. Later on, this approach was extended for the solution of different problems. Lepik [24,25] introduced a numerical method using the HW to obtain the solution of PDEs. Jiwari [26] used the HW coupled with quasi-linearization for the solution of Burgers’ equation. Mittal et al. [27] studied numerically the system of viscous Burgers’ equations with the HW. Oruç [28,29] applied the finite difference hybrid scheme combined with the HW for the solution of modified Burgers’ and KdV equations. Kumar [30] solved the system of Burgers’ equations through the finite difference HW technique. Somayeh et al. [31,32] used a semi-analytical approach for solving the Hunter–Saxton and foam drainage equations. The same author [33] implemented an HW-based scheme for the solution of the two-dimensional system of PDEs. Mittal and Pandit [34] solved unsteady squeezing nano-fluid problems via the HW. Jiwari [35] developed a hybrid numerical method consisting of the HW for the nonlinear Burger’s equation. Pandit and Mitall [36] introduced the scale-3 Haar wavelet for the numerical solutions of the fractional advection dispersion equation. The same author and his co-authors [37] defined an operational-matrix-based algorithm for computational modeling of hyperbolic-type wave equations. For more details about these wavelet, we refer to [25,38].

In the present work, a hybrid scheme consisting of the HW and finite difference is suggested to find the numerical solution of Equation (1) with the homogeneous and non-homogeneous forms with variable coefficients. The remainder of this paper is as follows. In Section 2 the main motive of the current work is given. The preliminaries of the HW and their integrals are given in Section 3. The description of the method and stability are presented in Section 4 and Section 5, respectively. For the validation of the suggested scheme, some test problems are addressed in Section 6. In Section 7 we present the detail description of the initial disturbance and noisy data. Finally, the conclusion is reported in Section 8.

2. Motivation

To develop methods for solutions of higher-order PDEs either analytical, semi-analytical, or numerical is an essential need for the physical interpretation of the problem. In the literature, several numerical techniques have been adopted for FOPDEs. According to our knowledge, a hybrid numerical method based on the HW and finite difference along with stability for FOPDEs has not been reported yet. Therefore, the main motive of this work is to determine numerical solutions of FOPDEs using finite differences and the HW.

3. Haar Wavelets and Their Integrals

The HW based on Haar functions, which were defined in 1910, belongs to a well-known class of the wavelet family known as the Daubechies wavelet. The HW is a recent mathematical tool, which became popular in the numerical study of various differential and integral equations. Initially, it was introduced for the interval $[0,1)$, but Lepik [24] extended this for any arbitrary interval $[\mathcal{A},\mathcal{B}).$

To define the HW, assume $x\in [\mathcal{A},\mathcal{B})$, then the HW family is defined for $i\ge 1$ as follows:

$$\begin{array}{ccc}& & h_1\left(x\right)=\left\{\begin{array}{cc}1,\hfill & x\in \left[\mathcal{A},\mathcal{B}\right)\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$

(5)

$$\begin{array}{ccc}& & h_i\left(x\right)=\left\{\begin{array}{cc}1,\hfill & x\in [{\xi }_1\left(i\right),{\xi }_2\left(i\right))\hfill \\ -1,\hfill & x\in [{\xi }_2\left(i\right),{\xi }_3\left(i\right)),\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.\hfill \end{array}$$

(6)

where ${\xi }_{s+1}\left(i\right)=\mathcal{A}+\left(2k+s\right)\nu \delta x,s=0,1,2\mathrm{and}\nu =\mathcal{M}/m$.

In the aforementioned equations, the interval is subdivided into $2\mathcal{M}$ intervals of equal length $\delta x=\frac{\mathcal{B}-\mathcal{A}}{2\mathcal{M}}$, where $\mathcal{M}=2^{\lambda }$ and $\lambda$ denote the maximal level of resolution. Next, two parameters $j=0,\cdots ,\lambda$ and $k=0,\cdots ,2^j-1$ are taken, which decompose the wavelet number $i=m+k+1$, where $m=2^j$. To solve $\mathrm{nth}$-order PDEs, we need the following repeated integrals:

$$\begin{array}{c}\hfill {\mathcal{R}}_{i,\beta }\left(x\right)={\int }_{\mathcal{A}}^x{\int }_{\mathcal{A}}^x\cdots {\int }_{\mathcal{A}}^x{h}_i\left(z\right)d{z}^{\beta }=\frac{1}{\left(\beta -1\right)!}{\int }_{\mathcal{A}}^x{\left(x-z\right)}^{\beta -1}{h}_i\left(z\right)dz,\end{array}$$

(7)

where

$$\beta =1,2,\cdots ,n,i=1,2,\cdots ,2\mathcal{M}.$$

Keeping in view Equations (5) and (6), the analytical expression of these integrals is given as [25]

$$\begin{array}{ccc}& & {\mathcal{R}}_{1,\beta }\left(x\right)=\frac{{\left(x-\mathcal{A}\right)}^{\beta }}{\beta !},\hfill \end{array}$$

(8)

$$\begin{array}{ccc}& & {\mathcal{R}}_{i,\beta }\left(x\right)=\left\{\begin{array}{cc}0,\hfill & x<{\xi }_1\left(i\right)\hfill \\ \frac{1}{\beta !}{\left[x-{\xi }_1\left(i\right)\right]}^{\beta }\hfill & x\in \left[{\xi }_1\left(i\right),{\xi }_2\left(i\right)\right)\hfill \\ \frac{1}{\beta !}\left[{(x-{\xi }_1\left(i\right))}^{\beta }-2({(x-{\xi }_2\left(i\right))}^{\beta }\right]\hfill & x\in [{\xi }_2\left(i\right),{\xi }_3\left(i\right))\hfill \\ \frac{1}{\beta !}\left[{(x-{\xi }_1\left(i\right))}^{\beta }-2({(x-{\xi }_2\left(i\right))}^{\beta }+{(x-{\xi }_3\left(i\right))}^{\beta }\right]\hfill & x\ge {\xi }_3\left(i\right).\hfill \end{array}\right.\hfill \end{array}$$

(9)

Function Approximation

If $w\left(x\right)\in L^2[\mathcal{A},\mathcal{B})$ is a square integrable function, then it can be approximated via HW series as:

$$w\left(x\right)=\sum _{i=1}^{2\mathcal{M}}{\alpha }_i{h}_i\left(x\right),$$

(10)

where ${\alpha }_i$ are the unknown HW coefficients. At collocation points $x\to x_l$, Equation (10) takes the following discrete form:

$$\tilde{w}=\sum _{i=1}^{2\mathcal{M}}{\alpha }_i{h}_i\left(x_l\right).$$

(11)

In matrix form, Equation (11) can be written as:

$$\Theta =\mathcal{J}\mathcal{A},$$

(12)

where $\Theta$ and $\mathcal{A}$ are $2\mathcal{M}\times 1$-dimensional matrices and $\mathcal{J}$ is a $2\mathcal{M}\times 2\mathcal{M}$-dimensional matrix. From Equation (12), one can calculate the unknown coefficients, and then, the approximation to $\mathcal{F}\left(x\right)$ can be computed using Equation (10) for different resolution levels.

4. Description of the Method

In this section of the manuscript, the proposed scheme is presented for Equation (1) with boundary conditions in the form of Equation (3). We rewrite Equation (1) as:

$$\gamma \left(x\right){\partial }_{tt}u+\mathcal{C}\left(x\right){\partial }_{xxxx}u+2{\mathcal{C}}^{\prime }\left(x\right){\partial }_{xxx}u+{\mathcal{C}}^{″}\left(x\right){\partial }_{xx}u=\mathcal{E},$$

(13)

where $C^{\prime }\left(x\right)=\frac{dC}{dx}$. Applying the finite difference to the temporal part and $(\theta$-weighted $\le \theta \le 1$ scheme, Equation (13) reduces to

$$\begin{array}{c}\gamma \left(x\right)\frac{u^{𝚥+1}-2u^{𝚥}+u^{𝚥-1}}{{\tau }^2}+\theta {\left[\mathcal{C}\left(x\right){\partial }_{xxxx}u+2{\mathcal{C}}^{\prime }{\partial }_{xxx}u+{\mathcal{C}}^{″}\left(x\right){\partial }_{xx}u\right]}^{𝚥+1}\\ +(1-\theta ){\left[\mathcal{C}\left(x\right){\partial }_{xxxx}u+2{\mathcal{C}}^{\prime }\left(x\right){\partial }_{xxx}u+{\mathcal{C}}^{″}\left(x\right){\partial }_{xx}u\right]}^{𝚥}={\mathcal{E}}^{𝚥+1},\end{array}$$

(14)

where $u^{𝚥}=u(x,t^{𝚥}),{\mathcal{E}}^{𝚥}=\mathcal{E}(x,t^{𝚥}),t^{𝚥+1}=\tau +t^{𝚥}$, and $\tau$ is the time step size. The associated boundary conditions Equation (3) is transformed to

$$\begin{array}{ccc}& & \left\{\begin{array}{c}{u}^{𝚥+1}\left(0\right)={\epsilon }_o(t^{𝚥+1}),{\partial }_x{u}^{𝚥+1}\left(0\right)={\delta }_o(t^{𝚥+1}),0<t\le T\hfill \\ u^{𝚥+1}\left(1\right)={\epsilon }_1(t^{𝚥+1}),{\partial }_x{u}^{𝚥+1}\left(1\right)={\delta }_1(t^{𝚥+1}),0<t\le T.\hfill \end{array}\right.\hfill \end{array}$$

(15)

In more simplified form, Equation (14) can be written as

$$\begin{array}{c}\gamma \left(x\right)u^{𝚥+1}+{\tau }^2\theta \mathcal{C}\left(x\right){\partial }_{xxxx}{u}^{𝚥+1}+2\theta {\tau }^2{\mathcal{C}}^{\prime }\left(x\right){\partial }_{xxx}{u}^{𝚥+1}+\theta {\tau }^2{\mathcal{C}}^{″}\left(x\right){\partial }_{xx}{u}^{𝚥+1}={\mathcal{N}}^{𝚥}+{\tau }^2{\mathcal{E}}^{𝚥+1},\end{array}$$

(16)

where

$$\begin{array}{cc}\hfill {\mathcal{N}}^{𝚥}=& \left(2u^{𝚥}-u^{𝚥-1}\right)\gamma \left(x\right)+(\theta -1){\tau }^2{\left[\mathcal{C}\left(x\right){\partial }_{xxxx}u+2{\mathcal{C}}^{\prime }\left(x\right){\partial }_{xxx}u+{\mathcal{C}}^{″}\left(x\right){\partial }_{xx}u\right]}^{𝚥}.\hfill \end{array}$$

Next, we assume the Haar wavelet approximation for the highest-order derivative as

$$\begin{array}{cc}\hfill {\partial }_{xxxx}{u}^{𝚥+1}\left(x\right)& =\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}{h}_i\left(x\right),\hfill \end{array}$$

(17)

where ${\alpha }_i^{𝚥+1}$ are unknown constants to be determined. Integration of Equation (17) four times from 0 to x leads to:

$$\begin{array}{cc}\hfill {\partial }_{xxx}{u}^{𝚥+1}\left(x\right)& =\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}{\mathcal{R}}_{i,1}\left(x\right)+{\partial }_{xxx}{u}^{𝚥+1}\left(0\right)\hfill \\ \hfill {\partial }_{xx}{u}^{𝚥+1}\left(x\right)& =\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}{\mathcal{R}}_{i,2}\left(x\right)+x{\partial }_{xxx}{u}^{𝚥+1}\left(0\right)+{\partial }_{xx}{u}^{𝚥+1}\left(0\right)\hfill \\ \hfill {\partial }_x{u}^{𝚥+1}\left(x\right)& =\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}{\mathcal{R}}_{i,3}\left(x\right)+\frac{x^2}{2!}{\partial }_{xxx}{u}^{𝚥+1}\left(0\right)+x{\partial }_{xx}{u}^{𝚥+1}\left(0\right)+{\partial }_x{u}^{𝚥+1}\left(0\right),\hfill \\ \hfill u^{𝚥+1}\left(x\right)& =\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}{\mathcal{R}}_{i,4}\left(x\right)+\frac{x^3}{3!}{\partial }_{xxx}{u}^{𝚥+1}\left(0\right)+\frac{x^2}{2!}{\partial }_{xx}{u}^{𝚥+1}\left(0\right)+x{\partial }_x{u}^{𝚥+1}\left(0\right)+u^{𝚥+1}\left(0\right).\hfill \end{array}$$

(18)

Using the boundary conditions $u^{𝚥+1}\left(1\right),{\partial }_x{u}^{𝚥+1}\left(1\right)$ in Equation (18), the unknown terms can be computed as

$$\begin{array}{cc}\hfill {\partial }_{xxx}{u}^{𝚥+1}\left(0\right)& =\varphi (t^{𝚥+1})+12\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}{\mathcal{R}}_{i,4}\left(1\right)-6\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}{\mathcal{R}}_{i,3}\left(1\right)\hfill \\ \hfill {\partial }_{xx}{u}^{𝚥+1}\left(0\right)& =\psi (t^{𝚥+1})-6\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}{\mathcal{R}}_{i,4}\left(1\right)+2\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}{\mathcal{R}}_{i,3}\left(1\right),\hfill \end{array}$$

(19)

where

$$\begin{array}{cc}\hfill \varphi \left(t^{𝚥+1}\right)& =6\left[{\delta }_1\left(t^{𝚥+1}\right)+{\delta }_o\left(t^{𝚥+1}\right)-2{\epsilon }_1\left(t^{𝚥+1}\right)+2{\epsilon }_o\left(t^{𝚥+1}\right)\right]\hfill \\ \hfill \psi \left(t^{𝚥+1}\right)& =-2\left[{\delta }_1\left(t^{𝚥+1}\right)+2{\delta }_o\left(t^{𝚥+1}\right)\right]-6\left[{\epsilon }_o\left(t^{𝚥+1}\right)-{\epsilon }_1\left(t^{𝚥+1}\right)\right].\hfill \end{array}$$

Making use of Equation (19) in Equation (18), we obtain

$$\begin{array}{cc}\hfill {\partial }_{xxx}{u}^{𝚥+1}\left(x\right)& =\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}\left[{\mathcal{R}}_{i,1}\left(x\right)-6{\mathcal{R}}_{i,3}\left(1\right)+12{\mathcal{R}}_{i,4}\left(1\right)\right]+\varphi \left(t^{𝚥+1}\right)\hfill \\ \hfill {\partial }_{xx}{u}^{𝚥+1}\left(x\right)& =\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}\left[{\mathcal{R}}_{i,2}\left(x\right)+(2-6x){\mathcal{R}}_{i,3}\left(1\right)+(12x-6){\mathcal{R}}_{i,4}\left(1\right)\right]+x\varphi \left(t^{𝚥+1}\right)+\psi \left(t^{𝚥+1}\right)\hfill \\ \hfill {\partial }_x{u}^{𝚥+1}\left(x\right)& =\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}\left[{\mathcal{R}}_{i,3}\left(x\right)+(2x-3x^2){\mathcal{R}}_{i,3}\left(1\right)+(6x^2-6x){\mathcal{R}}_{i,4}\left(1\right)\right]+\frac{x^2}{2}\varphi \left(t^{𝚥+1}\right)\hfill \\ \hfill & +x\psi \left(t^{𝚥+1}\right)+{\delta }_o\left(t^{𝚥+1}\right)\hfill \\ \hfill u^{𝚥+1}\left(x\right)& =\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}\left[{\mathcal{R}}_{i,4}\left(x\right)+(x^2-x^3){\mathcal{R}}_{i,3}\left(1\right)+(2x^3-3x^2){\mathcal{R}}_{i,4}\left(1\right)\right]+\frac{x^3}{3!}\varphi \left(t^{𝚥+1}\right)\hfill \\ \hfill & +\frac{x^2}{2}\psi \left(t^{𝚥+1}\right)+x{\delta }_o\left(t^{𝚥+1}\right)+{\epsilon }_o\left(t^{𝚥+1}\right).\hfill \end{array}$$

(20)

The proposed technique is based on the collocation approach; therefore, the collocation points are

$$x_l=\frac{l-0.5}{2\mathcal{M}},l=1,\cdots ,2\mathcal{M}.$$

Putting Equations (17) and (20) in Equation (16) and replacing x with $x_l$ leads to the following system of linear equations

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \sum _{i=1}^{2\mathcal{M}}}{\alpha }_i^{𝚥+1}[\gamma \left(x_l\right)\mathbb{A}(i,l)+{\tau }^2\theta \mathcal{C}\left(x_l\right)h_i\left(x_l\right)+2\theta {\tau }^2{\mathcal{C}}^{\prime }\left(x_l\right)\mathrm{{\rm Y}}(i,l)+{\tau }^2\theta {\mathcal{C}}^{″}\left(x_l\right)\mathrm{\Psi }(i,l)]\\ =\Lambda (i,𝚥+1)+\Omega (i,𝚥+1),\end{array}\end{array}$$

(21)

where

$$\begin{array}{cc}\hfill \mathbb{A}(i,l)=& \left[{\mathcal{R}}_{i,4}\left(x\right)+(x^2-x^3){\mathcal{R}}_{i,3}\left(1\right)+(2x^3-3x^2){\mathcal{R}}_{i,4}\left(1\right)\right]\hfill \\ \hfill \mathrm{{\rm Y}}(i,l)=& \left[{\mathcal{R}}_{i,1}\left(x\right)-6{\mathcal{R}}_{i,3}\left(1\right)+12{\mathcal{R}}_{i,4}\left(1\right)\right]\hfill \\ \hfill \mathrm{\Psi }(i,l)=& \left[{\mathcal{R}}_{i,2}\left(x\right)+(2-6x){\mathcal{R}}_{i,3}\left(1\right)+(12x-6){\mathcal{R}}_{i,4}\left(1\right)\right]\hfill \\ \hfill \Lambda (l,t^{𝚥+1})=& \mathcal{N}(l,t^{𝚥})+\mathcal{E}(l,t^{𝚥+1})\hfill \\ \hfill \Omega (i,𝚥+1)=& -\gamma \left(x_l\right)\left\{\frac{x^3}{3!}\varphi \left(t^{𝚥+1}\right)+\frac{x^2}{2}\psi \left(t^{𝚥+1}\right)+x{\delta }_o\left(t^{𝚥+1}\right)+{\epsilon }_o\left(t^{𝚥+1}\right)\right\}\hfill \\ \hfill & -2{\tau }^2\theta {\mathcal{C}}^{\prime }\left(x\right)\left\{\varphi \left(t^{𝚥+1}\right)\right\}-{\tau }^2\theta \left\{x\varphi \left(t^{𝚥+1}\right)+\psi \left(t^{𝚥+1}\right)\right\}.\hfill \end{array}$$

There are $2\mathcal{M}$ equations in Equation (21), which can be solved for $2\mathcal{M}$ unknowns iteratively. After the computation of the unknown constants, the required solution can be computed from Equation (20).

Note

To use the second kind of boundary conditions given in (Equation (4)), one may use first $u^{𝚥+1}\left(1\right),{\partial }_{xx}{u}^{𝚥+1}\left(1\right)$ in Equation (18) to obtain the following system of equations:

$$\begin{array}{cc}\hfill {\mu }_1\left(t^{𝚥+1}\right)& =\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}{\mathcal{R}}_{i,2}\left(1\right)+{\partial }_{xxx}{u}^{𝚥+1}\left(0\right)+{\mu }_0\left(t^{𝚥+1}\right)\hfill \\ \hfill {\zeta }_1\left(t^{𝚥+1}\right)& =\sum _{i=1}^{2\mathcal{M}}{\alpha }_i^{𝚥+1}{\mathcal{R}}_{i,4}\left(1\right)+\frac{1}{3!}{\partial }_{xxx}{u}^{𝚥+1}\left(0\right)+\frac{1}{2!}{\mu }_0\left(t^{𝚥+1}\right)+{\partial }_x{u}^{𝚥+1}\left(0\right)+{\zeta }_0\left(t^{𝚥+1}\right).\hfill \end{array}$$

(22)

After solving Equation (22), one can easily obtain ${\partial }_{xxx}{u}^{𝚥+1}\left(0\right)$${\partial }_x{u}^{𝚥+1}\left(0\right)$. Once these terms are calculated, then expressions for the derivatives and solutions can be extracted following the strategies given in Equations (19)–(21).

5. Stability

Here, we present the theoretical result related to the stability of the proposed technique. To derive the condition, one can deduce the following equations from Equation (20) as:

$${\partial }_{xxxx}{u}^{𝚥+1}=\mathcal{H}{\alpha }^{𝚥+1},$$

(23)

$${\partial }_{xxx}{u}^{𝚥+1}=\mathrm{{\rm Y}}{\alpha }^{𝚥+1}+{\widetilde{\mathrm{{\rm Y}}}}^{𝚥+1},$$

(24)

$${\partial }_{xx}{u}^{𝚥+1}=\mathrm{\Psi }{\alpha }^{𝚥+1}+{\tilde{\mathrm{\Psi }}}^{𝚥+1},$$

(25)

$$u^{𝚥+1}=\mathcal{F}{\alpha }^{𝚥+1}+{\tilde{\mathcal{F}}}^{𝚥+1},$$

(26)

where $\mathcal{H},\mathrm{{\rm Y}},\mathrm{\Psi }$ are differentiation matrices, $\mathcal{F}$ is the interpolation matrix, and ${\widetilde{\mathrm{{\rm Y}}}}^{𝚥+1},{\tilde{\mathrm{\Psi }}}^{𝚥+1},{\tilde{\mathcal{F}}}^{𝚥+1}$ are boundary terms given in Equation (20). Using initial condition $\frac{u^{𝚥}-u^{𝚥-1}}{\tau }={\mathrm{{\rm Y}}}_2\left(x\right)$ together with Equations (23)–(26) and $x\to x_l$ in Equation (16), we obtain

$$\begin{array}{cc}\hfill & \left[\gamma \left(x_l\right)\mathcal{F}(i,l)+{\tau }^2\theta C\left(x_l\right)\mathcal{H}(i,l)+2\theta {\tau }^2{C}^{\prime }\left(x_l\right)\mathrm{{\rm Y}}(i,l)+\theta {\tau }^2{C}^{″}\left(x_l\right)\mathrm{\Psi }(i,l)\right]{\alpha }^{𝚥+1}=[\gamma \left(x_l\right)\mathcal{F}(i,l)\hfill \\ & +(\theta -1){\tau }^2\left(C\left(x_l\right)\mathcal{H}(i,l)+2C^{\prime }\left(x_l\right)\mathrm{{\rm Y}}(i,l)+C^{″}\left(x_l\right)\mathrm{\Psi }(i,l)\right)]{\alpha }^{𝚥}+{\tau }^2{\mathcal{E}}^{𝚥+1}+\tau \gamma \left(x_l\right){\mathrm{{\rm Y}}}_2\left(x_l\right).\hfill \end{array}$$

(27)

In alternative form, Equation (27) can be written as:

$$\mathbb{M}{\alpha }^{𝚥+1}=\mathbb{N}{\alpha }^{𝚥}+{\mathbb{B}}^{𝚥+1},$$

(28)

where

$$\mathbb{M}=\gamma \left(x_l\right)\mathcal{F}(i,l)+{\tau }^2\theta C\left(x_l\right)\mathcal{H}(i,l)+2\theta {\tau }^2{C}^{\prime }\left(x_l\right)\mathrm{{\rm Y}}(i,l)+\theta {\tau }^2{C}^{″}\left(x_l\right)\mathrm{\Psi }(i,l),$$

$$\mathbb{N}=\gamma \left(x_l\right)\mathcal{F}(i,l)+(\theta -1){\tau }^2\left(C\left(x_l\right)\mathcal{H}(i,l)+2C^{\prime }\left(x_l\right)\mathrm{{\rm Y}}(i,l)+C^{″}\left(x_l\right)\mathrm{\Psi }(i,l)\right),$$

$${\mathbb{B}}^{𝚥+1}={\tau }^2{\mathcal{E}}^{𝚥+1}+\tau \gamma \left(x_l\right){\mathrm{{\rm Y}}}_2\left(x_l\right).$$

It follows from Equation (28) that

$${\alpha }^{𝚥+1}={\mathbb{M}}^{-1}\mathbb{N}{\alpha }^{𝚥}+{\mathbb{M}}^{-1}{\mathbb{B}}^{𝚥+1}.$$

(29)

Plugging Equation (29) in Equation (26), one obtains

$$\begin{array}{cc}\hfill u^{𝚥+1}=& \mathcal{F}\left({\mathbb{M}}^{-1}\mathbb{N}{\alpha }^{𝚥}+{\mathbb{M}}^{-1}{\mathbb{B}}^{𝚥+1}\right)+{\mathcal{F}}^{𝚥+1}\hfill \\ \hfill =& \mathcal{F}{\mathbb{M}}^{-1}\mathbb{N}{\alpha }^{𝚥}+\mathcal{F}{\mathbb{M}}^{-1}{\mathbb{B}}^{𝚥+1}+{\mathcal{F}}^{𝚥+1}.\hfill \end{array}$$

(30)

Using Equations (26) and (30), we have

$$\begin{array}{cc}\hfill u^{𝚥+1}& =\mathcal{F}{\mathbb{M}}^{-1}\mathbb{N}\left({\mathcal{F}}^{-1}{u}^j-{\mathcal{F}}^{-1}{\tilde{\mathcal{F}}}^{𝚥+1}\right)+\mathcal{F}{\mathbb{M}}^{-1}{\mathbb{B}}^{𝚥+1}+{\mathcal{F}}^{𝚥+1}\hfill \\ & =\mathcal{F}{\mathbb{M}}^{-1}\mathbb{N}{\mathcal{F}}^{-1}{u}^j-\mathcal{F}{\mathbb{M}}^{-1}\mathbb{N}{\mathcal{F}}^{-1}{\tilde{\mathcal{F}}}^{𝚥+1}+\mathcal{F}{\mathbb{M}}^{-1}{\mathbb{B}}^{𝚥+1}+{\mathcal{F}}^{𝚥+1}.\hfill \end{array}$$

(31)

Equation (31) shows an iterative formula between $u^{𝚥+1}$ and $u^{𝚥}$. If ${\tilde{u}}^{𝚥}$ is the approximate solution, then

$$\begin{array}{cc}\hfill {\tilde{u}}^{𝚥+1}& =\mathcal{F}{\mathbb{M}}^{-1}\mathbb{N}{\mathcal{F}}^{-1}{\tilde{u}}^j-\mathcal{F}{\mathbb{M}}^{-1}\mathbb{N}{\mathcal{F}}^{-1}{\tilde{\mathcal{F}}}^{𝚥+1}+\mathcal{F}{\mathbb{M}}^{-1}{\mathbb{B}}^{𝚥+1}+{\mathcal{F}}^{𝚥+1}.\hfill \end{array}$$

(32)

Let $e^{𝚥}=|u^{𝚥}-{\tilde{u}}^{𝚥}|$, then from Equations (31) and (32),

$$e^{𝚥+1}=\Xi e^{𝚥},$$

(33)

where $\Xi =\mathcal{F}{\mathbb{M}}^{-1}\mathbb{N}{\mathcal{F}}^{-1}$ is the amplification matrix. According to the Lax–Richtmyer criterion, the stability condition will be fulfilled if $\left|\right|\Xi \le 1\left|\right|,$ which needs the spectral radius $\rho (\Xi )\le 1$.

6. Illustrative Examples

In this section, we investigate the performance of the proposed method by solving various examples. To measure the efficiency, two error norms $L_2,L_{\infty }$ were addressed, which are defined by:

$$\begin{array}{cc}\hfill L_2& ={\left(\sum _i^{2\mathcal{M}}{\left|u_i^{ext}-{\tilde{u}}_i^{app}\right|}^2\right)}^{\frac{1}{2}}\hfill \\ \hfill L_{\infty }& =\underset{1\le i\le 2\mathcal{M}}{max}\left|u_i^{ext}-{\tilde{u}}_i^{app}\right|,\hfill \end{array}$$

where u and $\tilde{u}$ denote the exact and numerical solutions, respectively. Furthermore, we compute the convergence rate by using the following formula

$$\mathrm{Convergence} \mathrm{rate}=\frac{log(e_{\infty }^{\lambda }-e_{\infty }^{\lambda +1})}{log\left(2\right)}.$$

6.1. Problem 5.1

Consider Equation (1) with $\gamma \left(x\right)=\mathcal{C}\left(x\right)=1,$ in the following form

$${\partial }_{tt}u(x,t)+{\partial }_{xxxx}u(x,t)=({\pi }^4-1)sin\pi xcost,0\le x\le 1,0<t\le T,$$

(34)

coupled with appropriate conditions

$$u(x,0)=\sin \pi x,{\partial }_{t}u(x,0)=0,0\le x\le 1,$$

(35)

$$u(0,t)=u(1,t)=0,{\partial }_{xx}u(0,t)={\partial }_{xx}u(1,t)=0,0<t\le T.$$

(36)

The exact solution of this problem is $u(x,t)=sin\pi xcost$. This problem was solved this problem with the help of the proposed method. In

Table 1. Absolute error at different points in Problem 5.1.

Methods	Points	t	$\mathit{x}=0.1$	$\mathit{x}=0.2$	$\mathit{x}=0.3$	$\mathit{x}=0.4$	$\mathit{x}=0.5$
Present	64	0.02	6.54 × 10${}^{-7}$	1.24 × 10${}^{-6}$	1.71 × 10${}^{-6}$	2.01 × 10${}^{-6}$	2.11 × 10${}^{-6}$
	64	0.05	5.22 × 10${}^{-6}$	9.94 × 10${}^{-6}$	1.36 × 10${}^{-5}$	1.60 × 10${}^{-5}$	1.69 × 10${}^{-5}$
	64	1	7.14 × 10${}^{-4}$	1.35 × 10${}^{-3}$	1.87 × 10${}^{-3}$	2.20 × 10${}^{-3}$	2.31 × 10${}^{-3}$
	128	0.02	2.60 × 10${}^{-7}$	4.94 × 10${}^{-7}$	6.81 × 10${}^{-7}$	8.00 × 10${}^{-7}$	8.42 × 10${}^{-7}$
	128	0.05	2.59 × 10${}^{-6}$	4.94 × 10${}^{-6}$	6.80 × 10${}^{-6}$	7.99 × 10${}^{-6}$	8.40 × 10${}^{-6}$
	128	1	6.83 × 10${}^{-4}$	1.29 × 10${}^{-3}$	1.78 × 10${}^{-3}$	2.10 × 10${}^{-3}$	2.21 × 10${}^{-3}$
Mittal [39]	91	0.02	3.20 × 10${}^{-5}$	6.08 × 10${}^{-5}$	8.37 × 10${}^{-5}$	9.84 × 10${}^{-5}$	1.04 × 10${}^{-4}$
	91	0.05	3.59 × 10${}^{-5}$	6.83 × 10${}^{-5}$	9.39 × 10${}^{-5}$	1.10 × 10${}^{-4}$	1.16 × 10${}^{-4}$
	91	1	6.32 × 10${}^{-5}$	1.20 × 10${}^{-4}$	1.65 × 10${}^{-4}$	1.94 × 10${}^{-4}$	2.04 × 10${}^{-4}$
	181	0.02	3.55 × 10${}^{-6}$	6.76 × 10${}^{-6}$	9.30 × 10${}^{-6}$	1.09 × 10${}^{-5}$	1.15 × 10${}^{-5}$
	181	0.05	3.99 × 10${}^{-6}$	7.58 × 10${}^{-6}$	1.04 × 10${}^{-5}$	1.23 × 10${}^{-5}$	1.29 × 10${}^{-5}$
	181	1	7.00 × 10${}^{-6}$	1.33 × 10${}^{-5}$	1.83 × 10${}^{-5}$	2.16 × 10${}^{-5}$	2.27 × 10${}^{-5}$
Caglar [40]	121	0.02	4.80 × 10${}^{-6}$	9.70 × 10${}^{-6}$	1.40 × 10${}^{-5}$	1.90 × 10${}^{-5}$	2.40 × 10${}^{-5}$
	191	0.02	5.20 × 10${}^{-6}$	2.10 × 10${}^{-6}$	3.10 × 10${}^{-6}$	4.20 × 10${}^{-6}$	5.20 × 10${}^{-6}$
	521	0.02	4.90 × 10${}^{-7}$	9.90 × 10${}^{-7}$	1.40 × 10${}^{-6}$	1.90 × 10${}^{-6}$	2.40 × 10${}^{-6}$
Aziz et al. [12]	20	0.05	9.30 × 10${}^{-6}$	8.00 × 10${}^{-6}$	2.80 × 10${}^{-6}$	1.00 × 10${}^{-6}$	2.70 × 10${}^{-6}$
Rashidinia [13]	20	0.05	2.91 × 10${}^{-6}$	1.73 × 10${}^{-6}$	1.60 × 10${}^{-6}$	2.23 × 10${}^{-6}$	2.60 × 10${}^{-7}$
Mohammadi [41]	40	0.05	2.96 × 10${}^{-6}$	1.77 × 10${}^{-6}$	1.64 × 10${}^{-6}$	2.28 × 10${}^{-6}$	2.65 × 10${}^{-7}$

, pointwise errors are compared with existing results [12,39,40,41] for fixed $\tau =0.005$ using different time and resolution levels. From the table, it is clear that the computed results are in good agreement with those available in the literature. Moreover, we calculated the convergence rate for this problem and addressed it in

Table 2. Convergence rate of Problem 5.1 at $t=0.2$.

Problem 4.1
$\mathit{\lambda }$	$\mathit{\tau }$	${\mathit{L}}_{\infty }$	Rate
2	1/100	90,568 × 10${}^{-3}$
3	1/200	2.4428 × 10${}^{-3}$	1.8904
4	1/400	6.8223 × 10${}^{-4}$	1.8402
5	1/800	2.0398 × 10${}^{-4}$	1.7418

, which is approximately of order two. Graphical solutions in the form of 2D and 3D plots with the absolute error are shown in Figure 1 for time $t=4$. It is clear from the figure that the proposed method gives accurate solutions for a small number of collocation points and matches well the exact solution.

6.2. Problem 5.2

Consider the non-homogeneous problem of the form [41]

$${\partial }_{tt}u(x,t)+{\partial }_{xx}\left(\left(1+sin\pi x\right){\partial }_{xx}u(x,t)\right)=\mathcal{E}(x,t),0\le x\le 1,0<t\le T,$$

(37)

coupled with the initial conditions

$$u(x,0)={\partial }_{t}u(x,0)=0,0\le x\le 1,$$

(38)

and the boundary conditions

$$u(0,t)=u(1,t)=0,{\partial }_{x}u(0,t)={\partial }_{x}u(1,t)=0,0<t\le T.$$

(39)

The corresponding source term can be adjusted according to the exact solution:

$$u(x,t)=x(1-x)exp(-t)t^{2}sin4\pi x.$$

We obtained the solution of this problem in the time domain $[0,4]$. In

Table 3. Absolute error in displacement of Problem 5.2.

Methods	Points	$\mathit{\tau }$	$\mathit{t}=0.2$	$\mathit{t}=0.4$	$\mathit{t}=0.8$	$\mathit{t}=1$	$\mathit{t}=2$	$\mathit{t}=4$
Present	32	0.001	1.76 × 10${}^{-13}$	5.72 × 10${}^{-13}$	1.80 × 10${}^{-12}$	2.14 × 10${}^{-12}$	2.87 × 10${}^{-12}$	2.39 × 10${}^{-12}$
	64	0.001	1.72 × 10${}^{-13}$	5.75 × 10${}^{-13}$	1.48 × 10${}^{-12}$	1.90 × 10${}^{-12}$	2.95 × 10${}^{-12}$	6.13 × 10${}^{-13}$
[41]	100	0.01	1.78 × 10${}^{-5}$	5.85 × 10${}^{-5}$	1.57 × 10${}^{-4}$	2.00 × 10${}^{-4}$	2.95 × 10${}^{-4}$	1.60 × 10${}^{-4}$
	200	0.005	2.38 × 10${}^{-6}$	7.80 × 10${}^{-6}$	2.09 × 10${}^{-5}$	2.67 × 10${}^{-5}$	3.94 × 10${}^{-5}$	2.57 × 10${}^{-5}$

, we record the absolute error at point $x=0.5$, using different times. From table, the computations show that our results are better than those described by Mohammadi [41]. In

Table 4. Maximum error norms of Problem 5.2 at different times with $\tau =0.001$.

	$\mathit{t}=0.2$		$\mathit{t}=0.5$		$\mathit{t}=1$		$\mathit{t}=4$
$\mathit{\lambda }$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$
4	6.10 × 10${}^{-5}$	3.38 × 10${}^{-4}$	3.86 × 10${}^{-4}$	2.11 × 10${}^{-3}$	1.02 × 10${}^{-3}$	5.58 × 10${}^{-3}$	8.64 × 10${}^{-4}$	4.72 × 10${}^{-3}$
5	1.11 × 10${}^{-5}$	7.36 × 10${}^{-5}$	5.85 × 10${}^{-5}$	3.29 × 10${}^{-4}$	2.24 × 10${}^{-4}$	1.23 × 10${}^{-3}$	2.29 × 10${}^{-4}$	1.25 × 10${}^{-3}$

, different error norms are calculated, which identify that the proposed method produces good results at small resolution levels. The convergence rate of this problem was computed and is addressed in

Table 5. Convergence rate of Problem 5.2 at $t=1$.

$\mathit{\lambda }$	$\mathit{\tau }$	${\mathit{L}}_{\mathit{\infty }}$	Rate
2	1/100	1.5254 × 10${}^{-2}$
3	1/200	3.9824 × 10${}^{-3}$	1.9374
4	1/400	9.5864 × 10${}^{-4}$	2.0545
5	1/800	2.1397 × 10${}^{-4}$	2.1635

, which shows that the scheme is approximately of order two. Exact versus approximate solutions together with the absolute error are plotted at time $t=4$ in Figure 2. It is clear from the figure that the exact and approximate solutions are in good agreement.

6.3. Problem 5.3

Now, we consider the following equation:

$${\partial }_{tt}u(x,t)+{\partial }_{xx}\left(\left(1+sin\pi x\right){\partial }_{xx}u(x,t)\right)=\mathcal{E}(x,t),0\le x\le 1,0<t\le T,$$

(40)

inscribed with the initial conditions

$$u(x,0)={\partial }_{t}u(x,0)=0,0\le x\le 1,$$

(41)

and the boundary conditions

$$\begin{array}{ccc}& & \left\{\begin{array}{c}u(0,t)=0,{\partial }_{x}u(0,t)=t^2{e}^{-t},0<t\le T,\hfill \\ u(1,t)=t^2{e}^{-t},{\partial }_{x}u(1,t)=t^2{e}^{-t},0<t\le T.\hfill \end{array}\right.\hfill \end{array}$$

(42)

The exact solution of this problem is given by

$$u(x,t)=\left(x+{sin}^3\pi x\right)t^2{e}^{-t}.$$

The approximate solution of this problem was computed at different resolution levels. In

Table 6. Maximum error norms of Problem 5.3 at different times with $\tau =0.001$.

	$\mathit{t}=0.2$		$\mathit{t}=0.5$		$\mathit{t}=1$		$\mathit{t}=4$
$\mathit{\lambda }$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$
4	2.73 × 10${}^{-3}$	1.93 × 10${}^{-2}$	1.29 × 10${}^{-2}$	9.08 × 10${}^{-2}$	3.14 × 10${}^{-2}$	2.20 × 10${}^{-1}$	2.50 × 10${}^{-2}$	1.75 × 10${}^{-1}$
5	2.72 × 10${}^{-3}$	1.93 × 10${}^{-3}$	1.27 × 10${}^{-2}$	9.07 × 10${}^{-2}$	3.08 × 10${}^{-2}$	2.20 × 10${}^{-1}$	2.46 × 10${}^{-2}$	1.75 × 10${}^{-1}$

, we present the error norms for different times (t = 0.2, 0.5, 1, 4) using $\lambda =4,5$. It is observed from the table that the error norms are small, which shows the efficiency of the suggested scheme. The absolute error in displacement at $x=0.5$ matches those reported in [41] and presented in

Table 7. Absolute error in displacement at $x=0.5$ of Problem 5.3.

Methods	Points	$\mathit{\tau }$	$\mathit{t}=0.2$	$\mathit{t}=0.4$	$\mathit{t}=0.8$	$\mathit{t}=1$	$\mathit{t}=2$	$\mathit{t}=4$
Present	32	0.01	1.40 × 10${}^{-3}$	1.51 × 10${}^{-3}$	1.46 × 10${}^{-4}$	6.95 × 10${}^{-4}$	3.62 × 10${}^{-3}$	2.68 × 10${}^{-3}$
[41]	100	0.01	7.82 × 10${}^{-3}$	2.59 × 10${}^{-2}$	7.27 × 10${}^{-2}$	9.43 × 10${}^{-2}$	1.44 × 10${}^{-1}$	7.65 × 10${}^{-2}$

. It is clear from the table that the computed errors at different times are small. Moreover, the spectral radius for this problem and for previous two problems are addressed in

Table 8. Spectral radius of Problems 5.1, 5.2, and 5.3 at $t=0.2$.

$\mathit{\lambda }$	Problem  5.1	Problem 5.2	Problem 5.3
	$\rho (\Xi )$	$\rho (\Xi )$	$\rho (\Xi )$
1	0.99984	0.99937	0.99937
2	0.99984	0.99929	0.99929
3	0.99984	0.99927	0.99927
4	0.99984	0.99927	0.99927

which shows that stability condition fulfilled. Graphical solutions and the absolute error are shown in Figure 3. From the figure, one can see that the exact and approximated solutions agree mutually.

7. Initial Disturbance and Noisy Data

Here, we discuss the effect of the small perturbation in the initial condition and the noisy data. In the initial perturbation, we take the initial condition $u_0={\mathrm{{\rm Y}}}_1\left(x\right)+ϵ$ where $ϵ$ is a small number. The idea of this perturbation was taken from [42]. Simulations were performed for different values of resolution levels $\lambda$ and and for different values of $ϵ$, which are given in

Table 9. Maximum error norms of all problems for $t=1.0,$ after initial disturbances.

$\mathit{ϵ}={10}^{-2}$	Problem  5.1		Problem  5.2		Problem  5.3
$\lambda$	$L_{\infty }$	$L_2$	$L_{\infty }$	$L_2$	$L_{\infty }$	$L_2$
1	4.53 × 10${}^{-2}$	2.02 × 10${}^{-1}$	4.82 × 10${}^{-1}$	3.16 × 10${}^0$	1.30 × 10${}^{-1}$	2.27 × 10${}^{-1}$
2	2.10 × 10${}^{-2}$	9.41 × 10${}^{-2}$	1.67 × 10${}^{-2}$	8.62 × 10${}^{-2}$	4.90 × 10${}^{-2}$	3.09 × 10${}^{-1}$
3	1.46 × 10${}^{-2}$	6.55 × 10${}^{-2}$	4.14 × 10${}^{-3}$	2.18 × 10${}^{-2}$	3.40 × 10${}^{-2}$	2.26 × 10${}^{-1}$
4	1.30 × 10${}^{-2}$	5.82 × 10${}^{-2}$	1.34 × 10${}^{-3}$	5.94 × 10${}^{-3}$	3.15 × 10${}^{-2}$	2.20 × 10${}^{-1}$
$ϵ={10}^{-3}$	Problem  5.1		Problem  5.2		Problem  5.3
$\lambda$	$L_{\infty }$	$L_2$	$L_{\infty }$	$L_2$	$L_{\infty }$	$L_2$
1	3.64 × 10${}^{-2}$	1.63 × 10${}^{-1}$	4.78 × 10${}^{-1}$	3.16 × 10${}^0$	1.36 × 10${}^{-1}$	7.63 × 10${}^{-1}$
2	1.18 × 10${}^{-2}$	5.30 × 10${}^{-2}$	1.55 × 10${}^{-2}$	8.56 × 10${}^{-2}$	5.02 × 10${}^{-2}$	3.16 × 10${}^{-1}$
3	5.38 × 10${}^{-3}$	2.40 × 10${}^{-2}$	3.99 × 10${}^{-3}$	2.18 × 10${}^{-2}$	3.39 × 10${}^{-2}$	2.26 × 10${}^{-1}$
4	3.75 × 10${}^{-3}$	1.67 × 10${}^{-2}$	9.02 × 10${}^{-4}$	4.70 × 10${}^{-3}$	3.13 × 10${}^{-2}$	2.20 × 10${}^{-1}$

. From the table, we concluded that the small perturbation in the initial data produces a small change in the solution, which shows that the method is stable. Similarly, we take $u_0\left(x_i\right)={\mathrm{{\rm Y}}}_1\left(x_i\right)+{(-1)}^{\left(i\right)}\times ϵ,i=1,\cdots ,2\mathcal{M},$ for the noisy case [43]. In

Table 10. Maximum error norms of all problems for $t=1.0ϵ={10}^{-2},$ with noisy initial data.

Noise = 1%	Problem  5.1		Problem  5.2		Problem  5.3
$\mathit{\lambda }$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$	${\mathit{L}}_{\infty }$	${\mathit{L}}_{\mathbf{2}}$
1	3.55 × 10${}^{-2}$	1.58 × 10${}^{-1}$	4.77 × 10${}^{-1}$	3.16 × 10${}^0$	1.37 × 10${}^{-1}$	7.67 × 10${}^{-1}$
2	1.08 × 10${}^{-2}$	4.84 × 10${}^{-2}$	1.54 × 10${}^{-2}$	8.56 × 10${}^{-2}$	5.03 × 10${}^{-2}$	3.16 × 10${}^{-1}$
3	4.35 × 10${}^{-3}$	1.94 × 10${}^{-2}$	3.98 × 10${}^{-3}$	2.18 × 10${}^{-2}$	3.39 × 10${}^{-2}$	2.26 × 10${}^{-1}$
4	2.72 × 10${}^{-3}$	1.21 × 10${}^{-2}$	8.47 × 10${}^{-4}$	4.69 × 10${}^{-3}$	3.13 × 10${}^{-2}$	2.20 × 10${}^{-1}$

, we present the numerical results for the noisy case. From the table, it is pretty much clear that the variation in the resolution level for the noisy case still gives stable solutions.

In both cases, when we increase the resolution level, the accuracy increases, which is obvious from Table 9 and Table 10. Besides this, the absolute error for the noisy case is plotted in Figure 4, which clearly indicates the accuracy of the techniques because this agrees with the previous error when no noise exists. The simulations and graphical error indicated that the method works for both types of data.

8. Conclusions

In this work, a mixed numerical method based on the Haar wavelet coupled with the finite difference was proposed to solve the FOPDEs. First, the scheme was tested with initial data, and then, small perturbations with and with out noise were introduced. The obtained results matched with earlier work and exact solutions. Furthermore, the accuracy of the scheme was checked via computing the $L_2$ and $L_{\infty }$ error norms. It was observed that the proposed scheme works well for smooth and noisy initial data, which indicates that the method can be applied for other such problems.