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

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.


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: γ(x)∂ tt u(x, t) + ∂ xx [C(x)∂ xx u(x, t)] = E (x, t), 0 < x < 1, 0 < t ≤ T, with appropriate initial and boundary conditions: or where the notations ∂ tt u(x, t), ∂ 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, γ(x) is the mass per unit length, C(x) shows the bending stiffness of the beam, 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 closedform 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 semianalytical 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 nonhomogeneous 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 Sections 4 and 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.

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.

Haar Wavelets and Their Integrals
The HW based on Haar functions, which were defined in 1910, belongs to a wellknown 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 [A, B).
Keeping in view Equations (5) and (6), the analytical expression of these integrals is given as [25] Function Approximation is a square integrable function, then it can be approximated via HW series as: where α i are the unknown HW coefficients. At collocation points x → x l , Equation (10) takes the following discrete form:w In matrix form, Equation (11) can be written as: where Θ and A are 2M × 1-dimensional matrices and J is a 2M × 2M-dimensional matrix. From Equation (12), one can calculate the unknown coefficients, and then, the approximation to F (x) can be computed using Equation (10) for different resolution levels.

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: where C (x) = dC dx . Applying the finite difference to the temporal part and (θ-weighted ≤ θ ≤ 1 scheme, Equation (13) reduces to where and τ is the time step size. The associated boundary conditions Equation (3) is transformed to In more simplified form, Equation (14) can be written as where Next, we assume the Haar wavelet approximation for the highest-order derivative as where α +1 i are unknown constants to be determined. Integration of Equation (17) four times from 0 to x leads to: Using the boundary conditions u +1 (1), ∂ x u +1 (1) in Equation (18), the unknown terms can be computed as where Making use of Equation (19) in Equation (18), we obtain The proposed technique is based on the collocation approach; therefore, the collocation points are Putting Equations (17) and (20) in Equation (16) and replacing x with x l leads to the following system of linear equations where There are 2M equations in Equation (21), which can be solved for 2M 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 (1), ∂ xx u +1 (1) in Equation (18) to obtain the following system of equations: After solving Equation (22), one can easily obtain ∂ xxx u +1 (0) ∂ x u +1 (0). Once these terms are calculated, then expressions for the derivatives and solutions can be extracted following the strategies given in Equations (19)-(21).

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: where H, Υ, Ψ are differentiation matrices, F is the interpolation matrix, andΥ +1 ,Ψ +1 ,F +1 are boundary terms given in Equation (20). Using initial condition u  −u −1 τ = Υ 2 (x) together with Equations (23)-(26) and x → x l in Equation (16), we obtain In alternative form, Equation (27) can be written as: where It follows from Equation (28) that Plugging Equation (29) in Equation (26), one obtains Using Equations (26) and (30), we have Equation (31) shows an iterative formula between u +1 and u  . Ifũ  is the approximate solution, theñ Let e  = |u  −ũ  |, then from Equations (31) and (32), where Ξ = F M −1 NF −1 is the amplification matrix. According to the Lax-Richtmyer criterion, the stability condition will be fulfilled if ||Ξ ≤ 1||, which needs the spectral radius ρ(Ξ) ≤ 1.

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 ∞ were addressed, which are defined by: where u andũ denote the exact and numerical solutions, respectively. Furthermore, we compute the convergence rate by using the following formula Convergence rate = log(e λ ∞ − e λ+1 ∞ ) log(2) .

Problem 5.1
Consider Equation (1) with γ(x) = C(x) = 1, in the following form coupled with appropriate conditions The exact solution of this problem is u(x, t) = sin πx cos t. This problem was solved this problem with the help of the proposed method. In Table 1 Table 2, 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.

Problem 5.2
Consider the non-homogeneous problem of the form [41] ∂ tt u(x, t) coupled with the initial conditions and the boundary conditions The corresponding source term can be adjusted according to the exact solution: u(x, t) = x(1 − x) exp(−t)t 2 sin 4πx.
We obtained the solution of this problem in the time domain [0, 4]. In Table 3, 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, 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, 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.     Now, we consider the following equation:

Methods Points
inscribed with the initial conditions and the boundary conditions The exact solution of this problem is given by The approximate solution of this problem was computed at different resolution levels. In Table 6, we present the error norms for different times (t = 0.2, 0.5, 1, 4) using λ = 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. 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 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.

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 = Υ 1 (x) + where is a small number. The idea of this perturbation was taken from [42]. Simulations were performed for different values of resolution levels λ and and for different values of , which are given in Table 9. 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 (x i ) = Υ 1 (x i ) + (−1) (i) × , i = 1, · · · , 2M, for the noisy case [43]. In Table 10, 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 Tables 9 and 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.

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 ∞ 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.