Method of approximations for the convection-dominated anomalous diffusion equation in a rectangular plate using high-resolution compact discretization.

The present method describes the high-resolution compact discretization method for the numerical solution of the nonlinear fractal convection-diffusion model on a rectangular plate by employing the Hausdorff distance metric. Estimation of anomalous diffusion is formulated by averaging forward and backward mesh stencils. The higher-order fractional derivatives are appropriately approximated on a minimum mesh stencil and subsequently considered for designing a numerical method that falls in the scope of expanded accuracy. Compact discretization is an efficient technique for partial differential equations; however, studies that apply high-resolution scheme for fractional-order systems are still uninvestigated. A second and fourth-order numerical method for the fractional-order convection-dominated anomalous diffusion equation in two dimensions is constructed for practical applications. Convergence of high-order method is obtained for the nonlinear partial differential equations employing Hausdorff fractal distance metric. The numerical simulations with fractal Graetz-Nusselt equation, fractal Poisson equation, fractal Schrödinger equation, and anomalous diffusion equations with variable and constant coefficients are considered to illustrate the utility of the numerical method in the context of local fractional partial differential equations.•The paper demonstrates a computational method for the fractal convection-diffusion model on a rectangular plate.•Two numerical methods of order two and four for the mildly nonlinear fractional-order convection-dominated anomalous diffusion equations are proposed.•The high-resolution scheme is computationally efficient and makes use of minimal data storage.Method name: High-order method for 2D convection-dominated anomalous diffusion equation, Graetz-Nusselt equation, Poisson equation, and Schrödinger equation in fractal media.

Specifications table

Method basics and direct submitted research

The differential system identifying normal diffusion is inadequate in characterizing many complex circumstances, for example, the diffusion process in heterogeneous or inhomogeneous medium (porous media). Fractal derivative by the Hausdorff distance metric is an alternative approach to describing the complex interplay between convection and anomalous diffusion. Anomalous diffusion describes many physical scenarios, such as crowded systems, for instance, diffusion through porous media or diffusion of protein within cells; sub-diffusion emerged as an estimation of macromolecular crowding in the cytoplasm. Anomalous diffusion is observed in many systems, including telomeres in the nucleus of cells, ultra-cold atoms, plasma membrane ion channels, mixing in the interstellar medium, colloidal particles in the cytoplasm [1], worm-like micellar solutions [2], moisture transport in cement-based substances [3]. The fractal diffusion equations' introduction analyzes anomalous diffusion's characterization and qualitative behaviour. The presence of convective phenomenon extends it to the fractal convection-diffusion equation. In recent years, local fractional or fractal partial differential equations (PDEs) in fractal spacetime have emerged as a hot topic for engineers, mathematicians, and numerical analysts. The non-uniqueness in the definitions of fractional derivatives, with unknown practical applications, restricts engineers from selecting the proper fractional derivative formula [4]. Among many available fractional derivatives, the formula described in [5], [6], [7], [8] resembles mathematical correctness, practical relevance, and physical foundation. The scientific explanation of magnetic resonance imaging [9], water conservancy [10], economic problems [11], and viscoelastic materials [12] is analyzed by differential equations on non-Euclidean Hausdorff fractal distance. He et al. [13] presented the geometrical explanation and understanding of practical applications of fractional calculus. The fractal derivative discussed by Chen [14] is a local differential operator and computationally efficient compared with the global fractional derivative that characterizes the Mittag–Leffler function and Levy stable distribution. The fractal diffusion equations describe the underlying long familiar stretched exponential decay of Kohlrausch–Williams–Watts and stretched Gaussian statistics [15]. The equivalence of fractal and Hausdorff derivatives is proposed by Li and Ostoja-Starzewski [16]. Most of the work presents a qualitative analysis of fractional and fractal differential, but their quantitative analysis is still immature.

This paper considers a two-dimensional nonlinear local fractional convection-diffusion equation on a rectangular region derived by the Hausdorff fractal distances metric. The fractal Poisson equation described in [17] is extended here for the mildly nonlinear model in a bounded domain defined byand

The main challenge in the numerical treatment of (1)-(2) lies in the solutions' instability. A sharp change in the solution behavior is observed with a slight change in the fractal parameter . The discrete replacement of partial derivatives on the Hausdorff fractal distances metric yields a simple but computationally inefficient method and often fails to resemble the theoretical convergence rate. This happens due to the presence of several singular coefficients. The asymptotic expansion of singular terms makes the computation feasible up to a reasonable convergence order of the numerical method. Although, the most accurate solution in a short computing time requires the scheme to fall in the range of enhanced accuracy and implements a comparatively smaller number of mesh points. Compact discretization with high accuracy is one way to overcome these difficulties. The construction of high-resolution compact discretization restricts the number of mesh points used in the computational stencil to a maximum of three in each spatial direction, irrespective of the order of computational schemes. Moreover, it brings a block-tridiagonal system that can be efficiently solved numerically. Applications of compact formulation in fractional differential equations have been discussed earlier in [18]. Solution techniques with compact discretization combined with B-spline and preconditioned methods for Fisher's equation, advection-diffusion equations, and diffusion-convection equation have been described in [19], [20], [21].

Preliminaries and local fractional difference operators

This section presents definitions, difference operators, and theorems to be used later. The fractal partial derivative is local in nature and does not involve integral convolution like other fractional derivatives such as the Riemann-Liouville, the Caputo-Fabrizio, the Grünwald-Letnikov, the Hilfer fractional, the Atangana-Baleanu, and the Hadamard fractional derivative [4]. It is noticed that, in the fractal space-time dimension, the trajectories associated with quantum particles are continuous yet not differentiable [22,23]. Therefore, fractal space-time demands redefining the concepts of velocity and acceleration [24]. The fractal partial derivatives and define the gradient of on the fractal media, which might be continuous but not differentiable anywhere. A detailed geometrical description of fractal derivatives in one dimension is explained by He [25]. Numerical simulation of bed‐load transport model in a heterogeneous sand bed by the fractal derivative is recently conducted by Nie et al. [26]. The local fractional or fractal partial derivative for modeling and analyzing the physical phenomenon through the Hausdorff derivative and the properties of their difference operators are described here.

The fractal partial derivatives of a function concerning the fractal measure and are defined as

Let , and . Then, one obtains

In the limiting case, the fractal partial derivatives along -direction is obtained as

On a similar line, one can quickly find

If is a continuous function, then there exists constant and , such thatand since , and , for .

Let the spatial domain is partitioned with uniformly spaced mesh points , , , having step-size along -axis and along -axis. We shall denote the exact solution value and point-wise partial derivative values as

Similar to the classical forward and backward difference operators, Eqs. (3)–(5) yields the following discretization formula

The compact operators (10)-(13) derived on Hausdorff fractal distance is a generalization of first and second central-difference operators on Euclidean distance.

The operators and respectively approximate the fractal partial derivatives and with -accuracy.

Proof. The Taylor series expansion with the assumptions , yields

In a similar manner, the operators and approximate the fractal partial derivatives and with -accuracy, since

The -order fractal partial derivatives can be expressed in terms of first-order partial derivatives, since

The -order fractal partial derivatives can be expressed as a linear combination of first and second-order partial derivatives. That is,

The operators , and , involve evaluations at four adjacent meshes and one central mesh point. Therefore, it is easy to estimate the nonlinear elliptic PDEs in fractal media numerically, whose exact solution value is not known.

Two-dimensional fractal elliptic PDEs

The general form of mildly nonlinear fractal PDEs for the function is given by subject to the boundary data , , , and . Upon using Theorem 1 and Eqs. (14)–(16), a second-order discretization to (21) leads to the system of nonlinear equations

The discrete replacement of Eq. (21) by the method (22) provides an -accurate solution values.

Proof. The proof follows from the application of series expansion to (14) - (16) and utilizing the Eq. (21).

. The approximated solution values of the fractal PDEs (21) can be obtained using the scheme (22) after associating the boundary data (2). But their computational cost is very high for getting precise solution values. It requires significantly small mesh steps and more iterations for convergence to the accurate solution values. Therefore, an increased convergence rate scheme will help determine a more precise solution with optimal computing time.

. The exact solution to the nonlinear fractal PDEs (21) is not known for arbitrary choice of , even if . Therefore, we need to formulate a high-resolution computational scheme for linear and nonlinear fractal PDEs separately.

Operator compact high-resolution computational method

We shall initially describe a high-resolution scheme for the nonlinear elliptic PDEsand later, the scheme will be generalized for to solve the fractal PDEs (21). The high-order discretization of the nonlinear PDEs (23) is obtained with the help of three-point compact operators, defined by

These operators and their composite help approximate partial derivatives of various orders. We will use them to design the algorithm that numerically estimates the fractal PDEs (21) with high accuracy. The construction of a high-resolution method requires the continuity condition . Now, let us define

Then, the high-resolution scheme for the numerical approximation of PDEs (23) is given bywhereand

The scheme (35) has a local truncation error of , which yields -accuracy to the solution values. One needs to assign the boundary data (2) for the algorithmic implementation, and the iterative method helps determine approximate solution values. The computational scheme for regular elliptic PDEs in two dimensions has been obtained in the past [27,28].

The application of the method (35) to the fractal PDEs (21), after replacing the fractional derivatives by the formula (6)-(7) fails, if the domain contains zero (singular point). The partial derivatives replacement of fractal derivatives leads to the coefficients that exhibit singularity over the domain . Thus, the high-resolution scheme for fractal elliptic PDEs requires further mathematical analysis. Hence, the method for linear and nonlinear equations will be designed separately to resolve singularities.

Computational method for linear fractal elliptic PDEs

The nonhomogeneous local fractional elliptic PDEs employing Hausdorff fractal distance over the domain takes the general form

The application of fractal partial derivatives (6) and (7) transforms the Eq. (36) into the singular PDEswhere

The elliptic character of PDEs (37) remains intact since , and the terms of appears with power having even multiple, as a result . Upon substituting the following linear formin the scheme (35), we obtainwhereand the values of coefficients are given in Appendix A. The local truncation error associated with the scheme (39) is and renders fourth-order solutions to the linear fractal PDEs (36). The application of the computational method (35) to the PDEs (37) appears with the terms like , in case , , the evaluation of leads to zero-divisors at . To overcome this difficulty, we use a third-order series expansionand

The implementation of series expansion to all such terms yields a modified scheme (39) that is free from the singular terms and estimates the approximate solution of fractal PDEs (36) after removing the local truncation errors. A detailed Maple code generation for linear fractal PDEs is reported in Appendix B, illustrating compact operator replacement and removal of singularities.

Computational method for nonlinear fractal elliptic PDEs

The computational scheme for nonlinear PDEs employing Hausdorff fractal distance involves many higher-order terms, and therefore, their computational scheme requires different treatment. To illustrate the method, we consider the following nonlinear convection-dominated anomalous diffusion equation

Using the formula (6) and (7), the nonlinear fractal PDEs (40) results in a variable coefficient singular elliptic PDEwhere

Upon substituting , in the formula (35), and performing similar treatment of third-order series expansion to each coefficient , etc., we obtain

The values of coefficients are given in Appendix A. The scheme (42) is consistent with the fractal PDEs (40) and yields fourth-order numerical solution values. The nonlinear scheme's numerical solution and computational convergence rate in (42) are executed upon removing the local truncation error terms and incorporating the boundary data. In practice, the convergence rate is computed by opting -spatial mesh step-size as a constant multiple of -spatial mesh step or vice versa. Maple code generation for nonlinear fractal PDEs can be obtained in a manner similar to linear fractal PDEs.

Convergence analysis

The approximation in the solution values depends on the fractal parameter, type of discretization, size of mesh steps, method order, stability, and nature of convergence. This section will describe the solution behaviour with diminishing mesh step-size and their dependency over the fractal parameter. Since the fractal PDE (21) is elliptic, therefore , as a result, both and are of the same sign. In the following analysis, we shall assume them to be positive. The nonlinear fractal PDEs (21) can be written aswhere

Let and denotes the exact and approximate solution and be the point-wise error. The high-order relation (30) corresponding to the exact solution can be expressed aswhile the expression corresponding to approximate solution is given by

The application of the Mean-value theorem yields

In a similar manner, we can obtainwhereand expressions of , , , and , are obtained from (26)-(29), (32) and (33) respectively, upon the symbolic replacement of to . The linear dependence of spatial mesh steps and , being real numbers, suggest to choose . Now, the error equation corresponding to the high-order scheme (35) is given by

An term involved in the error Eq. (49) is neglected for error analysis since approaches zero faster than the term . The error Eq. (49) in the matrix-vector form yieldswhere is the matrix of point-wise errors, is the vector of local truncation errors, and is the block-tri-diagonal matrix. Non-zero elements of the matrix F are obtained as

For the sufficiently small mesh steps and , the block of sub-matrices in the matrix has non-zero values. Further, is an irreducible matrix since the graph of matrix is strongly connected. The connected character can be seen by defining a directed path corresponding to non-zero entries of the matrix . It is observed that there exists a directed path connecting arbitrary couple of nodes and corresponding to non-zero entries of the matrix , [29,30]. Fig. 1 illustrate the strongly connected graph on a mesh network. Let be the row element sum of the matrix , then, in the limiting case of mesh step-size , we find

Fig. 1

Strongly connected graph on a mesh network.

Since, and has the same sign (assumed positive), it follows that the matrix is monotone [31]. Consequently, the matrix is monotone and irreducible, thus invertible and . Let be the ()th element of . Now, define the following matrix and vector normsand

Let denote the vector with one as each element; then, by using the algebraic relationship , we can achieve

Using the relation (51), and applying Taylor's series expansion, one can determine the upper bounds on each non-zero element of the matrix . It is computed as

Upon employing the above inequalities on the elements of the matrix , one can obtain the following bonds on the vector of point-wise error from the Eq. (50) in the following manner

That is,

The relation (52) corroborates the (fourth-order) convergence to the computational method for the approximate numerical solution of fractal PDEs (21). The convergence analysis in case both and are simultaneously negative can be carried on a similar line.

The parameter corresponds to the coefficient in the Eq. (48), and its value depends upon , and . The quantity remains non-negative since and existence of unique solution ensure .

Applications and numerical simulations

In this section, we shall present the efficiency of the computational method of some fractal PDEs that appears in mathematical physics. The -accurate scheme (22) provides a convergence rate of second-order while the -accurate scheme (39) or (42) yields a fourth-order convergence rate. The new high-resolution scheme is designed with the influence of only eight neighboring mesh points and one central mesh point. Such a formulation consumes minimum computational resources and is memory space efficient. The maximum-absolute-error and numerical convergence rate help analyze deviation in the approximate and exact solution value. These measures are defined as

The Newton-Raphson method with zero vector as the initial solution guess and taking the error tolerance solve the nonlinear fractal PDEs, while the Gauss-Seidel iterative method has been applied to solve the linear difference equations [32]. In each case, the boundary data is obtained from the exact solution as a test procedure with . Maple 2017 and Python 3.0 are used for symbolic and numerical computations on MacBook Pro's Catalina operating system with 16GB of memory and a 2.6GHz 6-Core Intel i7 processor.

helps determine the electric potential for a given charge distribution, . The forcing function is the ratio of total volume charge density and permittivity of the fractal medium. The electric potential is computed using the proposed computational schemes and compared with the exact solution . The -errors and convergence rate is presented in Table 1 for and 1. In each case, the convergence order by the scheme is close to four, while the scheme yields an order less than two. The high-order scheme exhibits smaller -errors, and numerical simulations with other values of shows similar accuracies in the electric potential. Fig. 2 presents the changing nature of electric potential for and Fig. 3 depicts the effect of fractal parameter on implicit electric potential , for .

Table 1

Maximum-absolute-errors and convergence order for .

M	ℓ_∞	Θ_∞	ℓ_∞	Θ_∞
α=1	O(h+k)^2-scheme		O(h+k)^4-scheme
4	7.48e-02	—	2.29e-02	—
8	2.66e-02	1.5	2.11e-03	3.4
16	8.48e-03	1.7	1.79e-04	3.6
α=0.8
4	4.66e-02	—	4.67e-03	—
8	1.62e-02	1.5	4.25e-04	3.5
16	4.72e-03	1.8	4.04e-05	3.4
α=0.6
4	1.84e-02	—	3.74e-02	—
8	6.35e-03	1.5	1.83e-03	4.4
16	1.84e-03	1.8	8.11e-04	4.5

Fig. 2

Changing nature of electric potential for and 1.

Fig. 3

Effect of fractal parameter on implicit solution plot for .

The fractal Poisson's equation

Fractal equation for an electron confined on a rectangular plate is given by

The approximate values of the wavefunction are obtained by and scheme and compared with the exact wavefunction

The -errors and convergence rate for are presented in Table 2 for various mesh points. Fig. 4 presents the changing nature of wavefunction for different fractal parameter and 1.

Table 2

Maximum-absolute-errors and convergence order with -scheme.

M	ℓ_∞	Θ_∞	ℓ_∞	Θ_∞
	α=1.0		α=0.9
4	2.31e-00	–	2.03e-00	—
8	6.06e-01	1.9	4.41e-01	2.2
16	7.83e-02	3.0	5.23e-02	3.1
32	5.76e-03	3.8	3.84e-03	3.8
	α=0.8		α=0.7
4	1.67e-00	—	1.27e-00	–
8	3.00e-01	2.5	2.08e-01	2.6
16	3.30e-02	3.2	1.96e-02	3.4
32	2.44e-03	3.8	2.73e-03	2.8

Fig. 4

Changing nature of wavefunction for and 1.

The fractal convection-diffusion equation describes the joint effect of convection and anomalous diffusion

The source function is obtained from the mass concentration over the domain . Table 3 enumerates the -errors and convergence rate of exact and approximate concentration for different mesh arrangements. Fig. 5 presents the changing nature of mass concentration for , and 1. Fig. 6 depicts the effect of fractal parameter on implicit electric potential , for .

Table 3

Maximum-absolute-errors and convergence order for .

M	ℓ_∞	Θ_∞	ℓ_∞	Θ_∞
α=1	O(h+k)^2-scheme		O(h+k)^4-scheme
4	6.92e-03	—	3.25e-03	—
8	2.07e-03	1.7	3.55e-04	3.2
16	5.96e-04	1.8	5.08e-05	2.8
32	1.59e-04	1.9	4.41e-06	3.5
α=0.8
4	6.92e-03	—	3.99e-03	—
8	2.07e-03	1.7	3.94e-04	3.3
16	7.59e-04	1.9	9.76e-05	2.0
32	2.02e-04	1.9	1.55e-05	2.7

Fig. 5

Changing nature in mass concentration for and 1.

Fig. 6

Effect of fractal parameter on implicit solution plot for .

The description of anomalous heat flux exchange is governed by the fractal Graetz-Nusselt equation

Here the Peclet number is the ratio of advective transport rate to anomalous diffusive transport rate [33]. The errors and convergence rate in approximate and exact solutionis reported in Table 4 for , and different mesh concentrations using the method.

Table 4

Maximum-absolute-errors and convergence order with -scheme.

M	ℓ_∞	Θ_∞	ℓ_∞	Θ_∞	ℓ_∞	Θ_∞
	α=1.0		α=0.8		α=0.6
4	1.04e-03	–	1.15e-03	—	5.23e-04	–
8	9.24e-05	3.5	8.61e-05	3.7	1.04e-04	2.3
16	1.02e-05	3.2	9.16e-06	3.2	1.68e-05	2.6
32	8.40e-07	3.6	1.50e-06	2.6	1.61e-06	3.4

The complex interaction of convective phenomenon and anomalous diffusion often takes the fractal convection-diffusion equation in nonlinear form. We consider the simulation of two models, one with constant and another with variable coefficients.

Case 1: Constant coefficients anomalous diffusionwith the concentration of mass transfer as .

Case 2: Variable coefficients anomalous diffusionwith the mass transfer concentration , and . The numerical simulations in each case are conducted over the rectangular domain and the forcing term is obtained from the concentration of mass transfer (case 1). Table 5 presents the maximum absolute errors in case 1 for using the scheme. The errors in case 2 are reported in Table 6 for and . The numerical rate corroborates the fourth-order convergence of the method for . Fig. 7 presents the changing nature of mass concentration for and Fig. 8 depicts the effect of fractal parameter on implicit mass concentration , for in case 2.

Table 5

Maximum-absolute-errors and convergence order with -scheme at .

M	l_∞	Θ_∞
4	1.86e-02	–
8	1.49e-03	3.6
16	1.80e-04	3.1
32	2.90e-05	2.6

Table 6

Maximum-absolute-errors and convergence order with -scheme.

M	ℓ_∞	Θ_∞	ℓ_∞	Θ_∞
	α=0.8		α=0.6
4	2.04e-02	–	1.09e-01	—
8	2.20e-03	3.2	1.39e-02	3.0
16	3.64e-04	2.6	1.75e-03	3.0
32	4.58e-05	3.0	1.49e-04	3.5

Fig. 7

Changing nature in mass concentration for .

Fig. 8

Effect of fractal parameter on implicit solution plot for .

Conclusion

The fractional-order elliptic equations in two dimensions with nonlinear anomalous diffusion is considered. The fractal PDEs inherently appear with singular coefficients, and their computational scheme often yields inaccuracies and instability in approximate solution values. Therefore, designing a computationally efficient technique that results in a precise solution value with the involvement of a minimum number of mesh stencils is essential. The minimum mesh stencil saves computational resources and yields an algorithm that makes use of space-time tradeoffs. We proposed an approximation method based on nine-point meshes and fractional compact operators, leading to a generalized scheme applied to fractional order elliptic equations. The computational method for convection-dominated anomalous diffusion shows a fourth-order convergence of the approximate solution values. If the mesh width is further diminished, see Table 1, Table 2, Table 3, Table 4, Table 5, Table 6, the -errors further improve with a change in fractional values of . The new results will provide valuable insight for the numerical analysis of nonlinear parabolic PDEs in fractal media.

Declaration of Competing Interest

The authors of this article declare that they have no conflict of interest.

Subject Area;	Mathematical and Computational Sciences
More specific subject area;	Computational method in anomalous diffusion
Method name;	High-order method for the convection-dominated anomalous diffusion equation
Name and reference of original method;	He, J., 2018. Fractal calculus and its geometrical explanation, Results in Physics, 10: 272–276.https://doi.org/10.1016/j.rinp.2018.06.011
Resource availability;	N.A.

c_10	[{k^2Rl,myy+h^2(Rl,mxx+2Vl,mx)+12R_l,m}Q_l,m−k^2Rl,my(2Ql,my+S_l,m)]P_l,m−h^2Q_l,m(Rl,mx+V_l,m)(2Pl,mx+R_l,m)
c_01	[{h^2Sl,mxx+k^2(Sl,myy+2Vl,my)+12S_l,m}P_l,m−h^2Sl,mx(2Pl,mx+R_l,m)]Q_l,m−k^2P_l,m(Sl,my+V_l,m)(2Ql,my+S_l,m)
c_20	12Q_l,mPl,m2+[Q_l,m{Pl,mxx−2Rl,mx−V_l,m}h^2+{Pl,myyQ_l,m−Pl,my(2Ql,my+S_l,m)}k^2]P_l,m−h^2Q_l,m(2Pl,mx+R_l,m)(Pl,mx−R_l,m)
c_02	12P_l,mQl,m2+[P_l,m{Ql,myy−2Sl,my−V_l,m}k^2+{Ql,mxxP_l,m−Ql,mx(2Pl,mx+R_l,m)}h^2]Q_l,m−k^2P_l,m(2Ql,my+S_l,m)(Ql,my−S_l,m)
c_11	h^2Q_l,m[2P_l,mSl,mx−S_l,m(2Pl,mx+R_l,m)]+k^2P_l,m[2Q_l,mRl,my−R_l,m(2Ql,my+S_l,m)]
c_12	P_l,m(2h^2Ql,mx−k^2R_l,m)−h^2Q_l,m(2Pl,mx+R_l,m)
c_21	Q_l,m(2k^2Pl,my−h^2S_l,m)−k^2P_l,m(2Ql,my+S_l,m)
c_22	P_l,mQ_l,m(k^2P_l,m+h^2Q_l,m)
t_1	[{k^2Rl,myy+h^2(Rl,mxx+2Xl,mx)+12R_l,m}Q_l,m−Rl,myk^2(2Ql,my+S_l,m)]P_l,m−h^2Q_l,m[2V_l,m(Xl,mxu_00+Yl,mx)+(2Pl,mx+R_l,m)(Rl,mx+X_l,m)]
t_2	[{Sl,mxxh^2+(Sl,myy+2Xl,my)k^2+12S_l,m}P_l,m−Sl,mxh^2(2Pl,mx+R_l,m)]Q_l,m−P_l,mk^2[2W_l,m(Xl,myu_00+Yl,my)+(2Ql,my+S_l,m)(Sl,my+X_l,m)]
t_3	h^2Q_l,m{P_l,mVl,mxx−2(Rl,mx+X_l,m)V_l,m−(2Pl,mx+R_l,m)Vl,mx}+k^2P_l,m{Q_l,mVl,myy−Vl,my(2Ql,my+S_l,m)}+12P_l,mQ_l,mV_l,m
t_4	k^2P_l,m{Q_l,mWl,myy−2(Sl,my+X_l,m)W_l,m−(2Ql,my+S_l,m)Wl,my}+Q_l,mh^2{P_l,mWl,mxx−Wl,mx(2Pl,mx+R_l,m)}+12P_l,mQ_l,mW_l,m,
t_5	h^2Q_l,m{(Pl,mxx−2Rl,mx−X_l,m)P_l,m−(2Pl,mx+R_l,m)(Pl,mx−R_l,m)}−2hQ_l,mδx(1){2P_l,mVl,mx−(Pl,mx+2R_l,m)V_l,m}−k^2P_l,m{(2Ql,my+S_l,m)Pl,my−Pl,myyQ_l,m}−2kP_l,mPl,myW_l,mδy(1)+4Q_l,m(Vl,m2(δx(1))^2+3Pl,m2)
t_6	k^2P_l,m[(Ql,myy−2Sl,my−X_l,m)Q_l,m−(2Ql,my+S_l,m)(Ql,my−S_l,m)]−2kP_l,mδy(1){2Q_l,mWl,my−(Ql,my+2S_l,m)W_l,m}−h^2Q_l,m{(2Pl,mx+R_l,m)Ql,mx−P_l,mQl,mxx}−2hQ_l,mQl,mxV_l,mδx(1)+4P_l,m(Wl,m2(δy(1))^2+3Ql,m2)
t_7	P_l,mQ_l,m(Yl,myyk^2+Yl,mxxh^2+12Y_l,m)−P_l,mYl,myk^2(2Ql,my+S_l,m)−Q_l,mYl,mxh^2(2Pl,mx+R_l,m)
t_8	P_l,mQ_l,m(Xl,myyk^2+Xl,mxxh^2+12X_l,m)−P_l,mXl,myk^2(2Ql,my+S_l,m)−Q_l,mXl,mxh^2(2Pl,mx+R_l,m)
t_9	k^2P_l,m{2Q_l,mVl,my−(2Ql,my+S_l,m)V_l,m}−h^2Q_l,mS_l,mV_l,m
t_10	h^2Q_l,m{2P_l,mWl,mx−(2Pl,mx+R_l,m)W_l,m}−k^2P_l,mR_l,mW_l,m
t_11	h^2Q_l,m{2P_l,mSl,mx−(2Pl,mx+R_l,m)S_l,m}+k^2P_l,m{2Q_l,mRl,my−(2Ql,my+S_l,m)R_l,m}
t_12	P_l,mk[{2Pl,myQ_l,m−P_l,m(2Ql,my+S_l,m)}k^2−Q_l,mS_l,mh^2]−2P_l,mW_l,m(P_l,mk^2+Q_l,mh^2)δy(1)
t_13	hQ_l,m[h^2{2P_l,mQl,mx−Q_l,m(2Pl,mx+R_l,m)}−P_l,mR_l,mk^2]−2Q_l,mV_l,m(P_l,mk^2+Q_l,mh^2)δx(1)