To address fourth-order non-linear equations, which involve a differential equation alongside four constraints, numerical techniques leverage these constraints to discretize the equation effectively. Specifically, the Quintic Hermite techniques we re employed. Quintic Hermite splines inherently include the discretized form of the differential equation at collocation points along with two boundary constraints for the interpolating polynomial. To obtain such a scheme, the fourth-order differential equation was split into a second-order form by introducing a twice continuously differentiable function, thereby reducing the constraints to two for each interpolating polynomial. To execute this, the following theorem is considered a necessary part of implementing the proposed technique: Theorem 1:

The Hermite interpolation problem of order n ∈ N is applicable to two-point boundary value problems. ¹⁹

Consider the non-linear EFK type Equation (1). To implement the proposed technique, we introduce a twice continuously differentiable function v ( ζ , t ) such that:- v ( ζ , t ) = − u ζζ ( ζ , t ) .

As investigated in, ⁵⁵ the induced coupled system of differential equations of our governing problem is given as v + u ζζ = 0 , (3) u t − β v ζζ + v + f ( u ) = 0 , (4) where f ( u ) = u − u m for m ∈ ℝ , subject to the following initial and boundary conditions: u ( ζ , 0 ) = u 0 ( ζ ) , (5) u ( a , t ) = g 1 ( t ) , u ( b , t ) = g 2 ( t ) . (6) v ( a , t ) = 0 , v ( b , t ) = 0 . (7)

These equations are the transformed forms of EFK into a system comprising one linear equation and one semi-linear equation that involves two functions, u and v.

The hybrid technique employs semi-discretization methods to discretize both the time and spatial directions. This two-step approach was utilized for a second-order coupled system of singularly perturbed partial differential equations (PDEs):

a. Time discretization

For temporal discretization, we utilized the Crank-Nicolson scheme, a specific instance of the weighted finite difference method, as investigated in. ⁵⁶ This scheme takes advantage of the nearest neighboring points to maintain continuity in behavior, avoiding any time discontinuities (see ⁵⁷ ). To apply this technique to the reduced coupled system of the EFK equation, we define the following terms.

Let us consider that both v n + 1 and u n + 1 are represent the approximating polynomial function at the current time, t n + 1 . Also, let ∆ t = T N , N > 0 , for t n = n ∆ t for n = 0 , 1 , 2 , … , N , be the cell of the time step, and N be the maximum of the sub-interval in the temporal direction. Since, for any function F, we have F n + 1 2 = F n + 1 + F n 2 (8)

Hence, using the above condition, the semi-discretized forms of our coupled systems are given as v n + 1 2 = − u ζζ n + 1 2 u n + 1 − u n ∆ t = β v ζζ n + 1 2 − v n + 1 2 − ( u n − u ) n + 1 2 (9) with the corresponding boundary conditions given in Eqs. (6) and (7), respectively.

b. Spatial discretization

Let a = ζ 0 < ζ 1 < ζ 2 < … < ζ M = b , where M is the maximum number of subintervals of I in the spatial direction. This indicates that we have a non-uniform grid spacing ζ i ∈ I = [ a , b ] for which h i = ζ i + 1 − ζ i . Hence, our subintervals are given by I i = [ ζ i , ζ i + 1 ] . For each sub-interval, I i , six collocation points were selected, including two boundary points. These collocation points are the roots of a shifted Legendre polynomial of order six, which is not equidistant from the mean. Details of these collocation points can be found in. ⁵⁸ In the spatial direction, we utilized the quintic Hermite interpolating polynomials p i ( ζ ) , p i ¯ ( ζ ) and p i ̿ ( ζ ) that introduced in ¹⁹ to approximate the trial function within each partition I i of the solution domain. The approximate function is defined as follows: u ( ζ ) = u ( ζ i ) ∑ i = 1 2 P i ( ζ ) + u ′ ( ζ i ) ∑ i = 1 2 p i ¯ ( ζ ) + u ′ ′ ( ζ i ) ∑ i = 1 2 p i ̿ ( ζ ) (10a) where p i ( ζ ) is the interpolated function u ( ζ ) at each nodal point of sub-interval I i and p i ¯ ( ζ ) and p i ̿ ( ζ ) are the interpolated slope and curvature at the nodal points, respectively. As in, ¹⁹ they are, respectively, given by p i ( ζ ) = { 6 ( ζ i + 1 − ζ ) 5 ( ζ i + 1 − ζ i ) 5 − 15 ( ζ i + 1 − ζ ) 4 ( ζ i + 1 − ζ i ) 4 + 10 ( ζ i + 1 − ζ ) 3 ( ζ i + 1 − ζ i ) 3 , ζ i ≤ ζ ≤ ζ i + 1 6 ( ζ − ζ i − 1 ) 5 ( ζ i − ζ i − 1 ) 5 − 15 ( ζ − ζ i − 1 ) 4 ( ζ i − ζ i − 1 ) 4 + 10 ( ζ − ζ i − 1 ) 3 ( ζ i − ζ i − 1 ) 3 , ζ i − 1 ≤ ζ ≤ ζ i 0 , otherwise ; p i ¯ ( ζ ) = { 3 ( ζ i + 1 − ζ ) 5 ( ζ i + 1 − ζ i ) 4 − 7 ( ζ i + 1 − ζ ) 4 ( ζ i + 1 − ζ i ) 3 + 4 ( ζ i + 1 − ζ ) 3 ( ζ i + 1 − ζ i ) 2 , ζ i ≤ ζ ≤ ζ i + 1 − 3 ( ζ − ζ i − 1 ) 5 ( ζ i − ζ i − 1 ) 4 + 7 ( ζ − ζ i − 1 ) 4 ( ζ i − ζ i − 1 ) 3 − 4 ( ζ − ζ i − 1 ) 3 ( ζ i − ζ i − 1 ) 2 , ζ i − 1 ≤ ζ ≤ ζ i 0 , otherwise ; p i ̿ ( ζ ) = { 2 ( ζ i + 1 − ζ ) 5 ( ζ i + 1 − ζ i ) 3 − ( ζ i + 1 − ζ ) 4 ( ζ i + 1 − ζ i ) 2 + 2 ( ζ i + 1 − ζ ) 3 ( ζ i + 1 − ζ i ) , ζ i ≤ ζ ≤ ζ i + 1 2 ( ζ − ζ i − 1 ) 5 ( ζ i − ζ i − 1 ) 3 + ( ζ − ζ i − 1 ) 4 ( ζ i − ζ i − 1 ) 3 + 2 ( ζ − ζ i − 1 ) 3 ( ζ i − ζ i − 1 ) 2 , ζ i − 1 ≤ ζ ≤ ζ i 0 , otherwise . (10b)

To apply orthogonal collocation within each subinterval I i , a new variable ω is introduced, such that ω = ζ − ζ i h i where ω varies from 0 to 1 as ζ varies from ζ i to ζ i + 1 . This transformation reduces the polynomial structure in Eq. (10) into the tabular form shown in Table 1. The Quintic Hermit polynomials are presented in the interval [ ζ i , ζ i + 1 ] in Table 1, and the values of the Hermit splines and their derivatives are presented in Table 2 at ω = { 0 , 1 } . Hence, using this idea, the approximating function at each time t n + 1 in each subinterval is taken as u n + 1 ( ω ) = ∑ i = 1 6 c i n + 1 H i ( ω ) v n + 1 ( ω ) = ∑ i = 1 6 d i n + 1 H i ( ω ) (11) where c i n + 1 and d i n + 1 are constants of approximations.

The process of selecting collocation points is crucial in the Hermit collocation methods. To obtain these colocation points, we use the zeros of the orthogonal polynomials called Jacobi polynomials. The Legendre and Chebyshev polynomials are special cases of Jacobi polynomials. Hence, in this study, we utilized the zeros of the shifted Legendre polynomials as collocation points.

However, at each of these collocation points, the semi-discretization of the aforementioned equation by substituting Eq. (8) into Eq. (9) in the coupled system, we obtain v n + 1 + v n 2 = − [ u ζζ n + 1 + u ζζ n 2 ] u n + 1 − u n = β ∆ t [ v ζζ n + 1 + v ζζ n 2 ] − ∆ t [ v n + 1 + v n 2 ] − ∆ t 2 [ f n + 1 ( u ) + f n ( u ) ] . (12)

As investigated in, ⁵⁹ the quasi-linearization of f ( u ) = u m − u where m ∈ ℝ has been obtained using the Taylor series approximation, and it is given as u m ( t + ∆ t ) = u m ( t ) + ∆ t ( u m ) t + O ( ∆ t ) 2 = u m ( t ) + ∆ tm u m − 1 u t + O ( ∆ t ) 2 (13)

By using finite difference approximation of u ( t + ∆ t ) with order two, we have: u ( t + ∆ t ) = u ( t ) + ∆ t u t + O ( ∆ t ) 2

Then, neglecting ( ∆ t ) 2 , this equation gives us: u t = u ( t + ∆ t ) − u ( t ) ∆ t (14)

Substituting Eq. (14) into Eq. (13), at t = t n , we obtain u m ( t n + ∆ t ) = u m ( t n ) + m u m − 1 ( u ( t n + ∆ t ) − u ( t n ) ) + O ( ∆ t ) 2 = u m ( t n ) + m u m − 1 u ( t n + ∆ t ) + m u m ( t n ) + O ( ∆ t ) 2 .

By approximating the value of u ( t n + 1 ) from the above difference equation, we obtain: u ( t n + 1 ) = u m ( t n ) + m u m − 1 u ( t n + ∆ t ) + m u m ( t n ) + O ( ∆ t ) 2 . (15)

Using Eq. (15) in Eq. (12) after neglecting the second-order truncation term, we obtain v n + 1 + u ζζ n + 1 2 = − ( v n + u ζζ n ) 2 u n + 1 − u n = β ∆ t ( v ζζ n + 1 + v ζζ n 2 ) − ∆ t ( v n + 1 + v n 2 ) + ∆ t ( u n + 1 + u n 2 ) − ∆ t 2 ( mu n + 1 ( u m − 1 ) n − ( m − 2 ) u mn ) .

By rearranging the terms, we obtain the coupled system of equations given by: v n + 1 + u ζζ n + 1 = − v n − u ζζ n u n + 1 + ∆ t 2 ( v n + 1 + mu n + 1 ( u m − 1 ) n − u n + 1 − β v ζζ n + 1 ) = u n + ∆ t 2 ( β v ζζ n + ( m − 2 ) u mn + u n − v n ) (16)

Since, from Eq. (11), the discretization of both u n + 1 and v n + 1 in the spatial direction for 0 ≤ ω ≤ 1 is given by u n + 1 = ∑ i = 1 6 c i n H i ( ω ) ; v n + 1 = ∑ i = 1 6 d i n H i ( ω ) . (17)

Since we know that The approximate solution lies in each subinterval I i , where H i ( ω ) are quintic Hermite splines defined over the interval [0, 1]. Therefore, using this concept and Eq. (17), the system of collocation equations obtained after discretization is explained as follows: ∑ i = 1 6 { c 4 ( r − 1 ) + i n + 1 K ji + d 4 ( r − 1 ) + i n + 1 η ji } = ∑ i = 1 6 { c 4 ( r − 1 ) + i n ν ji + d 4 ( r − 1 ) + i n α ji } (18) ∑ i = 1 6 { c 4 ( r − 1 ) + i n + 1 ρ ji + d 4 ( r − 1 ) + i n + 1 γ ji } = ∑ i = 1 6 { c 4 ( r − 1 ) + i n χ ji + d 4 ( r − 1 ) + i n λ ji } + ( m − 2 ) Δ t 2 [ ∑ i = 1 6 c 4 ( r − 1 ) + i n H i ( ω ) ] m (19) where r = 1 ( 1 ) M , j = 1 ( 1 ) 4 , c i , & d i are the unknown coefficients to be determined. The constant parameters are defined as in, ¹⁹ and are given by ρ ji = H i ( ω j ) − Δt 2 H i ( ω j ) + mΔt 2 H i ( ω j ) [ ∑ i = 1 6 c 4 ( r − 1 ) + i n H i ( ω j ) ] m ; K ji = 1 2 h 2 H ′ ′ i ( ω j ) , η ji = 1 2 H i ( ω j ) , α ji = − 1 2 H i ( ω j ) ; γ ji = − βΔt 2 h 2 H ′ ′ i ( ω j ) + Δt 2 H i ( ω j ) ; λ ji = βΔt 2 h 2 H ′ ′ i ( ω j ) − Δt 2 H i ( ω j ) ; χ ji = H i ( ω j ) + Δt 2 H i ( ω j ) , ν ji = − 1 2 h 2 H ′ ′ i ( ω j ) . (20)

The total number of variables was 8M + 4, and the collocation equations were 8M. The additional four variables are computed from the boundary conditions of the system of equations. Therefore, by introducing a vector of unknown variables, the compact form of the system of difference equations is [ A ] Y n + 1 = [ B ] Y n + f n (21)

where A and B are 8 M × 8 M diagonally dominated matrices with each block A r and B r for r = 2 , … , M − 1 of order 8 × 12 and they are given as: Y n = [ d n ( i ) c n ( i ) ] 8 M × 1 , d n = [ d 1 n d 2 n ⋮ d 4 M n ] 4 M × 1 , c n = [ c 1 n c 2 n ⋮ c 4 M n ] 4 M × 1 , A = [ A 1 A 2 ⋱ A r ⋱ A M ] 8 M × 8 M , B = [ B 1 B 2 ⋱ B r ⋱ B N ] 8 M × 8 M ; and the source term f n contains both the nonlinear term ( m − 2 ) Δ t 2 [ ∑ i = 1 6 c 4 ( r − 1 ) + i n H i ( ω ) ] m and the non-homogeneous boundary data points of the functional values. It is given by: fn = [ 0 ( m − 2 ) Δ t 2 [ ∑ i = 1 6 c i n H i ( ω j ) ] m ⋮ 0 ( m − 2 ) Δ t 2 [ ∑ i = 1 6 c 4 ( r − 1 ) + i n H i ( ω j ) ] m ⋮ 0 ( m − 2 ) Δ t 2 [ ∑ i = 1 6 c 4 ( M − 1 ) + i n H i ( ω j ) ] m ]

Each block matrix, A r and B r also contained a sub-block matrix of order 2 × 2 . The block matrices A r and B r having the following structures: A r = [ [ η ji K ji γ ji ρ ji ] ij ] , B r = [ [ α ji ν ji λ ji χ ji ] ij ] where i = 1 , 2 , 3 , 4 and j = 1 , 2 , … , 6 except for the case of r = 1 and M due to known boundary data. This means that for known boundary constraints, all block matrixes , A 1 , A M , B 1 , and B r have an order 8 × 10 . Hence, from Eq. (21), Y n + 1 is 8 M × 1 column vectors given by. Y n + 1 = [ A − 1 B ] Y n + [ A − 1 f n ] (22)

After determining the unknown parameters c i , & d i from the coupled system of difference equations in Eq. (22), we can determine the discretized functional values u n + 1 and v n + 1 in the spatial direction by applying Eq. (17).

However, in the numerical scheme above, the presence of a non-linear term may lead to error propagation during the iteration process to solve the problem. Therefore, we mitigated this issue by linearizing our coupled system of difference equations.

Let u be considered a locally constant value, denoted as μ, for quasi-linearization of the non-linear term. The implemented technique involves operators Ψ and Φ, which operate on u and v, given as Ψ u , v and Φ u , v respectively. This operator approximates two sets of values for each iteration at the selected collocation points for the two approximate solutions, ŭ and v ˘ . This means that the iterative scheme at each time step approximates a function comprising two component functional values U = { u , v } . To obtain this approximate function using the scheme, we employed the component Hermite interpolating polynomial, defined as follows: H = { H 1 ˘ , H 2 ˘ } on C 6 [ a , b ] . Let u be the bounded functional value. The operators can then be linearized by considering: μ = max m , n [ u ( t n ) ] m − 1 .

We assume that U ˘ = { u ˘ , v ˘ } is an approximation of U = { u , v } . As investigated in, ¹⁹ we applied the Crank-Nicolson discretization scheme to our operators, which are given as follows: Ψ u , v ( ω i ) = P 1 ( u j ( ω i ) , v j ( ω i ) ) ; Φ u , v ( ω i ) = P 2 ( u j ( ω i ) , v j ( ω i ) ) ; (23a) where Ψ ( U ( ω i ) ) ≈ Ψ u , v ( ω i ) = v j + 1 ( ω i ) + u ζζ j + 1 ( ω i ) Φ ( U ( ω i ) ≈ Φ u , v ( ω i ) ) = u j + 1 ( ω i ) − ∆ t 2 ( β v ζζ j + 1 ( ω i ) − v j + 1 ( ω i ) + ( 1 − mμ ) u j + 1 ( ω i ) ) ; P 1 ( u j ( ω i ) , v j ( ω i ) ) ≅ − v j ( ω i ) − u ζζ j ( ω i ) ; P 2 ( u j ( ω i ) , v j ( ω i ) ) ≅ u j ( ω i ) + ∆ t 2 ( β v ζζ j ( ω i ) − v j ( ω i ) + ( 1 − mμ ) u j ( ω i ) ) j = 0 , 1 , 2 , … , N & i = 0 , 1 , 2 , … , M .

The matrix notation for the parameters and collocation equations can be expressed as follows: A Y j + 1 = [ Ψ u , v Φ u , v ] , B Y j = [ P 1 j P 2 j ] , Y j = [ d i j c i j ] (23b) where Y j contains the coefficients of Hermite polynomials. In this approximation technique, both A and B are symmetrical and diagonally dominant. Consequently, every symmetric and diagonally dominant matrix is positive definite. Thus, by employing the Gauss-Seidel iterative method, we can determine the unknown coefficients of Hermite polynomials.

Because, to discretize Eq. (1), the Crank-Nicolson scheme is applied to discretize the temporal derivative, which has a second-order accurate numerical technique, as indicated in Eq. (13). Subsequently, the quintic hermite spline collocation technique was employed in the spatial direction. The degree of the interpolated quintic hermite splines function on the collocation points is six. Thus, the numerical method resulting from this interpolation has 6th-order-accurate as shown in Eq. (10) in a specific direction. Therefore, by combining these cases, the method proposed in Eq. (16) exhibits a six-order accuracy in the spatial direction and a second-order accuracy in the temporal direction. This implies that the order of accuracy of the proposed method is O ( h 6 + ∆ t 2 ) .

One advantage of the RE technique is that it enhances the accuracy of the approximate solutions for given differential equation (see ^{52, 60} ). It is also used to accelerate the convergence of the underlying method in solving the provided equation (see ⁵³ ).

However, the goal of this extrapolation is to apply the difference scheme to two consecutive grids and then combine the resulting solutions to achieve higher-order approximations (see ^{54, 61} ). To do this, let us consider a set containing consecutive grid points, which are called Ω h i ∆ t coarser and Ω h i ∆ t / 2 finer points that are defined on the grid points of the temporal direction. Hence, at each and every one of these grid points, we can obtain a numerical solution of Eq. (1) using the following difference equations: u h i ∆ t = u m ( t n ) + m u m − 1 u ( t n + ∆ t ) + m u m ( t n ) + C ∆ t 2 + R ( h i , ∆ t ) , ( x i , t n ) ∈ Ω h i ∆ t . (24a) u h i ∆ t / 2 = u m ( t n ) + m u m − 1 u ( t n + ∆ t ) + m u m ( t n ) + C 2 ∆ t 2 + R ( h i , ∆ t / 2 ) , ( x i , t n ) ∈ Ω h i ∆ t / 2 (24b) where R ( h i , ∆ t ) and R ( h i , ∆ t / 2 ) are the remainders of order O ( h 6 + ∆ t 4 ) and C is a constant independent of h and ∆ t . To eliminate C ∆ t 2 , Eq. (24a) by − 1 / 15 and Eq. (24b) by 16 / 15 and combining the results, we obtain a more accurate approximate solution using the formulas given by u ( t n ) = 16 u h i ∆ t / 2 − u h i ∆ t 15 + O ( h 6 + ∆ t 4 ) , ( x i , t n ) ∈ Ω h i ∆ t (25a)

This can be extended to both v and u . Based on Eq. (16), we define the values of V and U as: V n = v n + u ζζ n U n = u n + ∆ t 2 ( β v ζζ n + ( m − 2 ) u mn + u n − v n )

Now using Richardson Extrapolation Formula, we have: V ( t n ) = 2 p V h i n / 2 − V h i n 2 p − 1 + O ( h 6 + ∆ t p ) U ( t n ) = 2 p U h i n / 2 − U h i n 2 p − 1 + O ( h 6 + ∆ t p ) where p is the order of convergence, ∆ t / 2 ≈ n / 2 , V ( t n ) and U ( t n ) represent the extrapolated values of v and u at time t n , V h i n / 2 and U h i n / 2 are their approximations values using a smaller time step ∆ t / 2 , V h i n and U h i n are the approximations using the original time step ∆ t respectively. In our case, the order of convergence for the time variable was = 4 . Using this and cutting off the truncation error, the extrapolation for V is given by V ( t n ) = 16 V h i n / 2 − V h i n 15 (25b)

Again, by eliminating the truncation error, the extrapolation of U is given by U ( t n ) = 16 U h i n / 2 − U h i n 15 (25c) where V ( t n ) and U ( t n ) represent the improved approximations for velocity and displacement, respectively, based on the values computed at the original and halved time steps.

These schemes demonstrate that the order of accuracy of the proposed scheme increases by two in the temporal direction after applying the Richardson extrapolation process. This means that the accuracy of the present scheme improved to the sixth order in the spatial direction and fourth order in the temporal direction following the Richardson extrapolation. The order of accuracy can be expressed as O ( h 6 + ∆ t 4 ) . This demonstrates that applying Richardson extrapolation doubles the convergence rate of the proposed scheme compared to that of the original scheme.