Time-Splitting Procedures for the Numerical Solution of the 2D Advection-Diffusion Equation

Abstract

We perform a spectral analysis of the dispersive and dissipative properties of two time-splitting procedures, namely, locally one-dimensional (LOD) Lax-Wendroff and LOD (1, 5) [9] for the numerical solution of the 2D advection-diffusion equation. We solve a 2D numerical experiment described by an advection-diffusion partial differential equation with specified initial and boundary conditions for which the exact solution is known. Some errors are computed, namely, the error rate with respect to the 𝐿₁ norm, dispersion and dissipation errors. Lastly, an optimization technique is implemented to find the optimal value of temporal step size that minimizes the dispersion error for both schemes when the spatial step is chosen as 0.025, and this is validated by numerical experiments.

1. Introduction

The advection-diffusion equation is a parabolic partial differential equation combining the diffusion and advection (convection) equations, which describes physical phenomena where particles, energy, or other physical quantities are transferred inside a physical system due to two processes: diffusion and advection [1]. The numerical solution of advection-diffusion equation plays an important role in many fields of science and engineering. These include the transport of air and groundwater pollutants, oil reservoir flow [2], heat transfer in draining film, flow through porous media, the dispersion of pollutants in rivers and streams, long range transport of pollutants in the atmosphere, thermal pollution in river systems, and dispersion of dissolved salts in groundwater [3].

The 3D advection-diffusion equation is given by 𝜕𝑢𝜕𝑡+𝛽_𝑥𝜕𝑢𝜕𝑥+𝛽_𝑦𝜕𝑢𝜕𝑦+𝛽_𝑧𝜕𝑢𝜕𝑧=𝛼_𝑥𝜕²𝑢𝜕𝑥²+𝛼_𝑦𝜕²𝑢𝜕𝑦²+𝛼_𝑧𝜕²𝑢𝜕𝑧², (1) where 𝛽_𝑥, 𝛽_𝑦, and 𝛽_𝑧 are the velocity components of advection in the directions of 𝑥, 𝑦, and 𝑧, respectively, and 𝛼_𝑥, 𝛼_𝑦, and 𝛼_𝑧 are the coefficients of diffusivity in the 𝑥-, 𝑦-, and 𝑧-directions, respectively.

This study deals with the 2D advection-diffusion equation, 𝜕𝑢𝜕𝑡+𝛽_𝑥𝜕𝑢𝜕𝑥+𝛽_𝑦𝜕𝑢𝜕𝑦=𝛼_𝑥𝜕²𝑢𝜕𝑥²+𝛼_𝑦𝜕²𝑢𝜕𝑦², (2) where 0<𝑡≤𝑇, in the domain 0≤𝑥≤1, 0≤𝑦≤1, with specified initial and boundary conditions.

Dehghan [3] proposed two time-splitting procedures for the solution of the two-dimensional transport equation. The time-splitting procedure used is the locally-one dimensional (LOD) in proceeding from one time step to the next step. LOD replaces the complicated multidimensional partial differential equations by a sequence of solutions of simpler one-dimensional partial differential equations. The originality in this work is that we perform a spectral analysis of the dispersion and dissipation properties of the two schemes at some values of the temporal and spatial step sizes. 3D plots of the relative phase error per unit time step (RPE) and the modulus of the amplification factor (AFM) versus phase angles in 𝑥- and 𝑦-directions are obtained. We then use optimization strategies to compute the optimal values of the temporal step size for the two schemes when the spatial step size is 0.025. We then validate this result by performing a 2D numerical experiment with specified initial and boundary conditions. Lastly, we have the conclusion and references.

2. Numerical Dispersion and Dissipation

Dissipation reduces the amplitude of sinusoids in a Fourier series. This is caused by the presence of derivatives like 𝑢_𝑥𝑥 and −𝑢_{𝑥𝑥𝑥𝑥} in the modified equation [4]. On the other hand, the amplitude of sinusoids in a Fourier series is increased by antidissipation. Derivatives like −𝑢_𝑥𝑥 and 𝑢_{𝑥𝑥𝑥𝑥} in the modified equation are generally antidissipative. Dispersion affects the speed of sinusoids in a Fourier series causing phase lag or phase lead and is caused due to the presence of odd-order derivatives in the modified equation [4].

The modulus of the amplification factor (AFM) is a measure of the stability of a scheme and it is also used to measure the dissipative characteristics of the scheme. If the modulus of the amplification factor is equal to one, a disturbance neither grows nor damps [5]. If the modulus of the amplification factor is greater than one, then the scheme is unstable [6]; if it is less than one, damping occurs [7]. The partial differential equation given by (2) is dissipative in nature due to the terms 𝑢_𝑥𝑥 and 𝑢_𝑦𝑦.

The relative phase error (RPE) is a measure of the dispersive characteristics of a scheme. The relative phase error of a scheme approximating the 1D advection-diffusion equation is given by 1RPE=−ℑ(𝜉)𝑐𝜔arctanℜ(𝜉), (3) where 𝑐 is the Courant number, 𝜔 is phase angle, 𝜉 is the amplification factor of the numerical scheme approximating the 1D advection-diffusion equation, and ℜ(𝜉) and ℑ(𝜉) are the real and imaginary parts of 𝜉, respectively [8].

We extend the work on the relative phase error in [8] for the case of the 2D advection-diffusion equation. The relative phase error for a numerical scheme approximating (2) is obtained on substituting 𝑢 by exp(𝐼(𝜔₁𝑡−𝜃₁𝑥−𝜃₂𝑦)) [9], where 𝐼=√−1, 𝑡 is the time, 𝜃₁ and 𝜃₂ are wavenumbers in the directions of 𝑥 and 𝑦, respectively, and 𝜔₁ is the dispersion relation. Thus, we get 𝜕𝑢𝜕𝑡=𝐼𝜔₁exp𝐼𝜔₁𝑡−𝜃₁𝑥−𝜃₂𝑦,𝜕𝑢𝜕𝑥=−𝐼𝜃₁exp𝐼𝜔₁𝑡−𝜃₁𝑥−𝜃₂𝑦,𝜕𝑢𝜕𝑦=−𝐼𝜃₂exp𝐼𝜔₁𝑡−𝜃₁𝑥−𝜃₂𝑦,𝜕²𝑢𝜕𝑥²=−𝜃₁²exp𝐼𝜔₁𝑡−𝜃₁𝑥−𝜃₂𝑦,𝜕²𝑢𝜕𝑦²=−𝜃₂²exp𝐼𝜔₁𝑡−𝜃₁𝑥−𝜃₂𝑦. (4) On substituting (4) into (2), we obtain 𝜔₁=𝛽_𝑥𝜃₁+𝛽_𝑦𝜃₂+𝛼_𝑥𝜃₁²𝐼+𝛼_𝑦𝜃₂²𝐼. (5)

The exact phase velocity is ℜ(𝜔₁)/wavenumber and we obtain 𝛽exactphasevelocity=^𝑥𝜃₁+𝛽_𝑦𝜃₂.wavenumber (6)

The amplification factor can be written as 𝜉=𝜉₁+𝐼𝜉₂, where 𝜉₁ and 𝜉₂ are the real and imaginary parts of 𝜉, respectively. Also, we can express 𝜉 as 𝜉=exp(−𝑏𝑘), where 𝑘 is time step and 𝑏 is exponential growth rate [9]. Thus, we obtain 1𝑏=𝑘𝜉log¹−𝐼𝜉₂𝜉²₁+𝜉²₂. (7)

The numerical phase velocity is calculated as ℑ(𝑏)/wavenumber and we get −()𝜉numericalphasevelocity=1/𝑘arctan₂/𝜉₁.wavenumber (8) The relative phase error is the ratio of the numerical phase velocity to the exact phase velocity [10]. It is calculated as 1RPE=−𝑘𝛽_𝑥𝜃₁+𝛽_𝑦𝜃₂𝜉arctan²𝜉₁. (9) The phase angles in the directions of 𝑥 and 𝑦 are given by 𝜔_𝑥=Δ𝑥𝜃₁ and 𝜔_𝑦=Δ𝑦𝜃₂, respectively. Hence, the relative phase error is 1RPE=−𝑐_𝑥𝜔_𝑥+𝑐_𝑦𝜔_𝑦𝜉arctan²𝜉₁, (10) where 𝑐_𝑥=𝛽_𝑥𝑘/Δ𝑥 and 𝑐_𝑦=𝛽_𝑦𝑘/Δ𝑦 are the Courant numbers in the directions of 𝑥 and 𝑦, respectively.

3. Time-splitting Procedures and Numerical Experiments

The domain we consider is 𝑥,𝑦∈[0,1]. We divide the spatial interval [0,1] along 𝑥- and 𝑦-directions into 𝑁 and 𝑀 nodes, respectively, such that (𝑁−1)Δ𝑥=1 and (𝑀−1)Δ𝑦=1 and also divide the time interval [0,𝑇] into 𝑃 grid points such that (𝑃−1)Δ𝑡=𝑇. Then the grid points (𝑥_𝑖,𝑦_𝑗,𝑡_𝑛) are defined by 𝑥_𝑖𝑦=𝑖Δ𝑥,𝑖=1,2,…,𝑁,_𝑗𝑡=𝑗Δ𝑦,𝑗=1,2,…,𝑀,_𝑛=𝑛Δ𝑡,𝑛=1,2,…,𝑃. (11) For simplification we take Δ𝑥=Δ𝑦=ℎ. Let Δ𝑡=𝑘; then the parameters ℎ and 𝑘 represent the spatial and temporal grid spacing, respectively. We denote the approximated value of 𝑢 at the grid point (𝑖,𝑗,𝑛) by 𝑢^𝑛_𝑖,𝑗.

3.1. Time-Splitting Procedures

Since one-dimensional schemes are easier to use than two-dimensional schemes, (2) is split into the following two one-dimensional equations: 12𝜕𝑢𝜕𝑡+𝛽_𝑥𝜕𝑢𝜕𝑥=𝛼_𝑥𝜕²𝑢𝜕𝑥²,12𝜕𝑢𝜕𝑡+𝛽_𝑦𝜕𝑢𝜕𝑦=𝛼_𝑦𝜕²𝑢𝜕𝑦². (12)(13) Each of (12) and (13) can be solved over half of a time step to be used for the complete 2D advection-diffusion equation, using the procedures developed for the 1D advection-diffusion equation.

Some work on time-splitting procedures can be found in [3, 11]. In this paper, we refer to [3] on how a 2D advection-diffusion equation is converted into two 1D advection-diffusion equations using the locally one-dimensional (LOD) time-splitting procedure. Solving (12) and (13) in each half time step is equivalent to solving the following equations over a full-time step: 𝜕𝑢𝜕𝑡+𝛽_𝑥𝜕𝑢𝜕𝑥=𝛼_𝑥𝜕²𝑢𝜕𝑥²,𝜕𝑢𝜕𝑡+𝛽_𝑦𝜕𝑢𝜕𝑦=𝛼_𝑦𝜕²𝑢𝜕𝑦². (14)(15) Then we can use the schemes used for solving the 1D advection-diffusion equation to solve (14) and (15).

3.2. Numerical Experiments

In [12], different explicit and implicit finite differences schemes are used to solve the 1D advection-diffusion equation 𝜕𝑢+𝜕𝑡𝜕𝑢𝜕𝜕𝑥=0.01²𝑢𝜕𝑥². (16) Three values of the cell Reynolds number, namely, 𝑅_Δ=2,4, and 8, are used for the numerical experiments. In [3], the time-splitting procedure is used to solve the two-dimensional transport equation in the region 0≤𝑥≤1, 0≤𝑦≤1, with 𝛼_𝑥=𝛼_𝑦=𝛼=0.01,𝛽_𝑥=𝛽_𝑦=𝛽=0.8 (17) and with the following boundary and initial conditions: 𝑔₀()=1𝑥,𝑡−()4𝑡+1exp𝑥−0.8𝑡−0.5²()−()0.014𝑡+10.8𝑡+0.5²(),𝑔0.014𝑡+1₁()=1𝑥,𝑡−()4𝑡+1exp𝑥−0.8𝑡−0.5²()−()0.014𝑡+10.5−0.8𝑡²(),ℎ0.014𝑡+1₀=1𝑦,𝑡−()4𝑡+1exp0.8𝑡+0.5²()−0.014𝑡+1𝑦−0.8𝑡−0.5²(),ℎ0.014𝑡+1₁=1𝑦,𝑡−()4𝑡+1exp0.5−0.8𝑡²()−0.014𝑡+1𝑦−0.8𝑡−0.5²(),𝑓−()0.014𝑡+1𝑥,𝑦=exp𝑥−0.5²−0.01𝑦−0.5²,0.01 (18)(19)(20)(21)(22) for which the exact solution is 𝑢=1𝑥,𝑦,𝑡−()4𝑡+1exp𝑥−0.8𝑡−0.5²()−0.014𝑡+1𝑦−0.8𝑡−0.5²().0.014𝑡+1 (23) We consider two time-splitting procedures LOD Lax-Wendroff and LOD (1,5) to solve (2) with 𝛼_𝑥=𝛼_𝑦=𝛼=0.01 and 𝛽_𝑥=𝛽_𝑦=𝛽=0.8, subject to the boundary conditions (18)–(22) at time 𝑇=0.3.

We consider two values of 𝑅_Δ; 𝑅_Δ=2 and 𝑅_Δ=4. Since 𝑐=0.8𝑘/ℎ and 𝑠=0.01𝑘/ℎ², we have 𝑅_Δ=80ℎ. Thus for 𝑅_Δ=2 and 4, we have ℎ=0.025 and ℎ=0.05, respectively. When ℎ=0.025, we have 𝑐=32𝑘 and 𝑠=16𝑘 and for ℎ=0.05 we have 𝑐=16𝑘 and 𝑠=4𝑘.

3.3. Quantification of Errors from Numerical Results

In this subsection, we describe how errors from numerical results can be quantified into dispersion and dissipation by a technique devised by Takacs [13].

The total mean square error in the 1D case [13] is calculated as 𝐸_TMS=1𝑁^𝑁_𝑖=1𝑢_𝑖−𝑣_𝑖², (24) where 𝑢_𝑖 is the exact solution, 𝑣_𝑖 is the numerical solution at a grid point 𝑖, and 𝑁 is the number of grid points. The total mean square error can be expressed as 1𝑁^𝑁_𝑖=1𝑢_𝑖−𝑣_𝑖²=1𝑁^𝑁_𝑖=1𝑢_𝑖−𝑢²+1𝑁^𝑁_𝑖=1𝑣_𝑖−𝑣²+2𝑁^𝑁_𝑖=1𝑢_𝑖2𝑢+𝑁^𝑁_𝑖=1𝑣_𝑖1𝑣−𝑁^𝑁_𝑖=1(𝑢)²−1𝑁^𝑁_𝑖=1(𝑣)²−2𝑁^𝑁_𝑖=1𝑢_𝑖𝑣_𝑖. (25)

The right hand side of (25) can be rewritten as 𝜎²(𝑢)+𝜎²(𝑣)+2(𝑢)²+2(𝑣)²−(𝑢)²−(𝑣)²−2𝑁^𝑁_𝑖=1𝑢_𝑖𝑣_𝑖, (26) where 𝜎²(𝑢) and 𝜎²(𝑣) denote the variance of 𝑢 and 𝑣, respectively, and 𝑢 and 𝑣 denote the mean values of 𝑢 and 𝑣, respectively. Then we have 𝐸_TMS=𝜎²(𝑢)+𝜎²(𝑣)+(𝑢)²−2𝑢𝑣+(𝑣)²+2𝑢2𝑣−𝑁^𝑁_𝑖=1𝑢_𝑖𝑣_𝑖=𝜎²(𝑢)+𝜎²(𝑣)+(𝑢−𝑣)²1−2𝑁^𝑁_𝑖=1𝑢_𝑖𝑣_𝑖−𝑢𝑣. (27) Therefore, 1𝑁^𝑁_𝑖=1𝑢_𝑖−𝑣_𝑖²=𝜎²(𝑢)+𝜎²(𝑣)+(𝑢−𝑣)²−2Cov(𝑢,𝑣), (28) where Cov(𝑢,𝑣)=(1/𝑁)(∑^𝑁_𝑖=1𝑢_𝑖𝑣_𝑖−𝑢𝑣).

The total mean square error can be expressed as 1𝑁^𝑁_𝑖=1𝑢_𝑖−𝑣_𝑖²=(𝜎(𝑢)−𝜎(𝑣))²+(𝑢−𝑣)²+21−𝜌𝜎(𝑢)𝜎(𝑣), (29) where 𝜌=Cov(𝑢,𝑣)/(𝜎(𝑢)𝜎(𝑣)) is the coefficient of correlation.

The term 2(1−𝜌)𝜎(𝑢)𝜎(𝑣) measures the dispersion error and the term (𝜎(𝑢)−𝜎(𝑣))²+(𝑢−𝑣)² measures the dissipation error.

We extend the work on quantification of errors in [14, 15] for the 2D case. The total mean square error for the 2D case is calculated as 𝐸_TMS=1𝑁𝑀^𝑀𝑁_𝑗=1𝑖=1𝑢_𝑖,𝑗−𝑣_𝑖,𝑗², (30) where 𝑢_𝑖,𝑗 and 𝑣_𝑖,𝑗 are the exact and numerical solutions at a grid point (𝑖,𝑗), respectively, and 1𝑁𝑀^𝑀𝑁_𝑗=1𝑖=1𝑢_𝑖,𝑗−𝑣_𝑖,𝑗²=1𝑁𝑀^𝑀𝑁_𝑗=1𝑖=1𝑢²_𝑖,𝑗−2𝑢_𝑖,𝑗𝑣_𝑖,𝑗+𝑣²_𝑖,𝑗=1𝑁𝑀^𝑀𝑁_𝑗=1𝑖=1𝑢²_𝑖,𝑗+^𝑀𝑁_𝑗=1𝑖=1𝑣²_𝑖,𝑗−2^𝑀𝑁_𝑗=1𝑖=1𝑢_𝑖,𝑗𝑣_𝑖,𝑗. (31) Since ^𝑀𝑁_𝑗=1𝑖=1𝑢_𝑖,𝑗−𝑢²=^𝑀𝑁_𝑗=1𝑖=1𝑢²_𝑖,𝑗−2𝑢_𝑖,𝑗𝑢+𝑢²,_{𝑀𝑁𝑗=1𝑖=1}𝑣_𝑖,𝑗−𝑣²=^𝑀𝑁_𝑗=1𝑖=1𝑣²_𝑖,𝑗−2𝑣_𝑖,𝑗𝑣+𝑣², (32) we have 𝐸_TMS=1𝑁𝑀^𝑀𝑁_𝑗=1𝑖=1𝑢_𝑖,𝑗−𝑢²+^𝑀𝑁_𝑗=1𝑖=1𝑣_𝑖,𝑗−𝑣²+2^𝑀𝑁_𝑗=1𝑖=1𝑢_𝑖,𝑗𝑢+2^𝑀𝑁_𝑗=1𝑖=1𝑣_𝑖,𝑗𝑣−2^𝑀𝑁_𝑗=1𝑖=1𝑢_𝑖,𝑗𝑣_𝑖,𝑗−^𝑀𝑁_𝑗=1𝑖=1𝑢²−^𝑀𝑁_𝑗=1𝑖=1𝑣²=𝜎²(𝑢)+𝜎²(𝑣)+2𝑢²+2𝑣²−2𝑁𝑀^𝑀𝑁_𝑗=1𝑖=1𝑢_𝑖,𝑗𝑣_𝑖,𝑗−𝑢²−𝑣²=𝜎²(𝑢)+𝜎²(𝑣)+𝑢²−2𝑢𝑣+𝑣²+2𝑢2𝑣−𝑁𝑀^𝑀𝑁_𝑗=1𝑖=1𝑢_𝑖,𝑗𝑣_𝑖,𝑗=𝜎²(𝑢)+𝜎²(𝑣)+(𝑢−𝑣)²1−2𝑁𝑀^𝑀𝑁_𝑗=1𝑖=1𝑢_𝑖,𝑗𝑣_𝑖,𝑗−𝑢𝑣=𝜎²(𝑢)+𝜎²(𝑣)+(𝑢−𝑣)²−2Cov(𝑢,𝑣). (33)

Hence, 𝐸_TMS=(𝜎(𝑢)−𝜎(𝑣))²+(𝑢−𝑣)²+21−𝜌𝜎(𝑢)𝜎(𝑣), (34) where 𝜌=Cov(𝑢,𝑣)/𝜎(𝑢)𝜎(𝑣). The dissipation error is (𝜎(𝑢)−𝜎(𝑣))²+(𝑢−𝑣)² and the dispersion error is 2(1−𝜌)𝜎(𝑢)𝜎(𝑣).

The error rate with respect to 𝐿₁ norm for 𝑥,𝑦∈[0,1] is calculated as 𝐸_num=1(𝑁−1)(𝑀−1)^𝑀𝑁_𝑗=1𝑖=1𝑢_𝑖,𝑗−𝑣_𝑖,𝑗. (35)

4. Construction of the LOD Lax-Wendroff Procedure

We use the following approximations in the first time step of the LOD procedure [3]: 𝜕𝑢𝜕𝑡^𝑛_𝑖,𝑗≃𝑢^𝑛+(1/2)_𝑖,𝑗−𝑢^𝑛_𝑖,𝑗,Δ𝑡𝜕𝑢𝜕𝑥^𝑛_𝑖,𝑗≃1−𝛾_𝑥𝑢^𝑛_𝑖,𝑗−𝑢^𝑛_{𝑖−1,𝑗}+𝛾_𝑥𝑢^𝑛_𝑖+1,𝑗−𝑢^𝑛_𝑖,𝑗,𝜕Δ𝑥₂𝑢𝜕𝑥^2𝑛_𝑖,𝑗≃𝑢^𝑛_𝑖+1,𝑗−2𝑢^𝑛_𝑖,𝑗−𝑢^𝑛_{𝑖−1,𝑗}(Δ𝑥)², (36) where 𝛾_𝑥 is the spatial weighting factor in the direction of 𝑥. The following approximations are used for the second step of the LOD procedure: 𝜕𝑢𝜕𝑡^𝑛+(1/2)_𝑖,𝑗≃𝑢^𝑛+1_𝑖,𝑗−𝑢^𝑛+(1/2)_𝑖,𝑗,Δ𝑡𝜕𝑢𝜕𝑦^𝑛+(1/2)_𝑖,𝑗≃1−𝛾_𝑦𝑢^𝑛+(1/2)_𝑖,𝑗−𝑢^𝑛+(1/2)_{𝑖,𝑗−1}+𝛾_𝑦𝑢^𝑛+(1/2)_𝑖,𝑗+1−𝑢^𝑛+(1/2)_𝑖,𝑗,𝜕Δ𝑦₂𝑢𝜕𝑦^2𝑛+(1/2)_𝑖,𝑗≃𝑢^𝑛+(1/2)_𝑖,𝑗+1−2𝑢^𝑛+(1/2)_𝑖,𝑗−𝑢^𝑛+(1/2)_{𝑖,𝑗−1}Δ𝑦², (37) where 𝛾_𝑦 is the spatial weighting factor in the direction of 𝑦. By using the following relationships 𝑐_𝑥=𝛽_𝑥𝑘ℎ,𝑠_𝑥=𝛼_𝑥𝑘ℎ²,𝛾_𝑥=1−𝑐_𝑥2,𝑐_𝑦=𝛽_𝑦𝑘ℎ,𝑠_𝑦=𝛼_𝑦𝑘ℎ²,𝛾_𝑦=1−𝑐_𝑦2, (38) the finite difference formula for the first half time step of the LOD Lax-Wendroff procedure is given by 𝑢^𝑛+(1/2)_𝑖,𝑗=122𝑠_𝑥+𝑐_𝑥+𝑐²_𝑥𝑢^𝑛_{𝑖−1,𝑗}+1−2𝑠_𝑥−𝑐²_𝑥𝑢^𝑛_𝑖,𝑗+122𝑠_𝑥−𝑐_𝑥+𝑐²_𝑥𝑢^𝑛_𝑖+1,𝑗. (39) At the second half time step the finite difference is given by 𝑢^𝑛+1_𝑖,𝑗=122𝑠_𝑦+𝑐_𝑦+𝑐²_𝑦𝑢^𝑛+(1/2)_{𝑖,𝑗−1}+1−2𝑠_𝑦−𝑐²_𝑦𝑢^𝑛+(1/2)_𝑖,𝑗+122𝑠_𝑦−𝑐_𝑦+𝑐²_𝑦𝑢^𝑛+(1/2)_𝑖,𝑗+1, (40) where 𝑢^𝑛+(1/2)_{𝑖,𝑗−1}=122𝑠_𝑥+𝑐_𝑥+𝑐²_𝑥𝑢^𝑛_{𝑖−1,𝑗−1}+1−2𝑠_𝑥−𝑐²_𝑥𝑢^𝑛_{𝑖,𝑗−1}+122𝑠_𝑥−𝑐_𝑥+𝑐²_𝑥𝑢^𝑛_{𝑖+1,𝑗−1},𝑢^𝑛+(1/2)_𝑖,𝑗+1=122𝑠_𝑥+𝑐_𝑥+𝑐²_𝑥𝑢^𝑛_{𝑖−1,𝑗+1}+1−2𝑠_𝑥−𝑐²_𝑥𝑢^𝑛_𝑖,𝑗+1+122𝑠_𝑥−𝑐_𝑥+𝑐²_𝑥𝑢^𝑛_{𝑖+1,𝑗+1}. (41) We obtain a single expression for 𝑢^𝑛+1_𝑖,𝑗 in a complete time step as follows: 𝑢^𝑛+1_𝑖,𝑗=142𝑠_𝑦+𝑐_𝑦+𝑐²_𝑦2𝑠_𝑥+𝑐_𝑥+𝑐²_𝑥𝑢^𝑛_{𝑖−1,𝑗−1}+122𝑠_𝑦+𝑐_𝑦+𝑐²_𝑦1−2𝑠_𝑥−𝑐²_𝑥𝑢^𝑛_{𝑖,𝑗−1}+142𝑠_𝑦+𝑐_𝑦+𝑐²_𝑦2𝑠_𝑥−𝑐_𝑥+𝑐²_𝑥𝑢^𝑛_{𝑖+1,𝑗−1}+121−2𝑠_𝑦−𝑐²_𝑦2𝑠_𝑥+𝑐_𝑥+𝑐²_𝑥𝑢^𝑛_{𝑖−1,𝑗}+1−2𝑠_𝑦−𝑐²_𝑦1−2𝑠_𝑥−𝑐²_𝑥𝑢^𝑛_𝑖,𝑗+121−2𝑠_𝑦−𝑐²_𝑦2𝑠_𝑥−𝑐_𝑥+𝑐²_𝑥𝑢^𝑛_𝑖+1,𝑗+142𝑠_𝑦−𝑐_𝑦+𝑐²_𝑦2𝑠_𝑥+𝑐_𝑥+𝑐²_𝑥𝑢^𝑛_{𝑖−1,𝑗+1}+122𝑠_𝑦−𝑐_𝑦+𝑐²_𝑦1−2𝑠_𝑥−𝑐²_𝑥𝑢^𝑛_𝑖,𝑗+1+142𝑠_𝑦−𝑐_𝑦+𝑐²_𝑦2𝑠_𝑥−𝑐_𝑥+𝑐²_𝑥𝑢^𝑛_{𝑖+1,𝑗+1}. (42) To find the modified equation of the scheme, we first find Taylor’s expansion of each term in (42) about 𝑢^𝑛_𝑖,𝑗. The Taylor series expansion of 𝑢^𝑛+1_𝑖,𝑗 is given by 𝑢^𝑛+1_𝑖,𝑗=𝑢+𝑘𝑢_𝑡+𝑘²2𝑢_𝑡𝑡+𝑘³6𝑢_𝑡𝑡𝑡+𝑘⁴𝑢4!^{𝑡𝑡𝑡𝑡}+⋯. (43)

The Taylor series expansions for some grid points about 𝑢^𝑛_𝑖,𝑗 are given as follows: 𝑢^𝑛_{𝑖−1,𝑗−1}=𝑢−ℎ𝑢_𝑥−ℎ𝑢_𝑦+ℎ²2𝑢_𝑥𝑥+2𝑢_𝑥𝑦+𝑢_𝑦𝑦−ℎ³6𝑢_𝑥𝑥𝑥+3𝑢_𝑥𝑥𝑦+3𝑢_𝑥𝑦𝑦+𝑢_𝑦𝑦𝑦+ℎ⁴𝑢4!^{𝑥𝑥𝑥𝑥}+4𝑢_{𝑥𝑥𝑥𝑦}+6𝑢_{𝑥𝑥𝑦𝑦}+4𝑢_{𝑥𝑦𝑦𝑦}+𝑢_{𝑦𝑦𝑦𝑦}𝑢+⋯,_{𝑛𝑖−1,𝑗}=𝑢−ℎ𝑢_𝑥+ℎ²2𝑢_𝑥𝑥−ℎ³6𝑢_𝑥𝑥𝑥+ℎ⁴𝑢4!^{𝑥𝑥𝑥𝑥}𝑢+⋯,_{𝑛𝑖,𝑗−1}=𝑢−ℎ𝑢_𝑦+ℎ²2𝑢_𝑦𝑦−ℎ³6𝑢_𝑦𝑦𝑦+ℎ⁴𝑢4!^{𝑦𝑦𝑦𝑦}𝑢+⋯,_{𝑛𝑖−1,𝑗+1}=𝑢−ℎ𝑢_𝑥+ℎ𝑢_𝑦+ℎ²2𝑢_𝑥𝑥−2𝑢_𝑥𝑦+𝑢_𝑦𝑦−ℎ³6𝑢_𝑥𝑥𝑥−3𝑢_𝑥𝑥𝑦+3𝑢_𝑥𝑦𝑦−𝑢_𝑦𝑦𝑦+ℎ⁴𝑢4!^{𝑥𝑥𝑥𝑥}−4𝑢_{𝑥𝑥𝑥𝑦}+6𝑢_{𝑥𝑥𝑦𝑦}−4𝑢_{𝑥𝑦𝑦𝑦}+𝑢_{𝑦𝑦𝑦𝑦}𝑢+⋯,_{𝑛𝑖+1,𝑗−1}=𝑢+ℎ𝑢_𝑥−ℎ𝑢_𝑦+ℎ²2𝑢_𝑥𝑥−2𝑢_𝑥𝑦+𝑢_𝑦𝑦−ℎ³6−𝑢_𝑥𝑥𝑥+3𝑢_𝑥𝑥𝑦−3𝑢_𝑥𝑦𝑦+𝑢_𝑦𝑦𝑦+ℎ⁴𝑢4!^{𝑥𝑥𝑥𝑥}−4𝑢_{𝑥𝑥𝑥𝑦}+6𝑢_{𝑥𝑥𝑦𝑦}−4𝑢_{𝑥𝑦𝑦𝑦}+𝑢_{𝑦𝑦𝑦𝑦}𝑢+⋯,_{𝑛𝑖+1,𝑗}=𝑢+ℎ𝑢_𝑥+ℎ²2𝑢_𝑥𝑥+ℎ³6𝑢_𝑥𝑥𝑥+ℎ⁴𝑢4!^{𝑥𝑥𝑥𝑥}𝑢+⋯,_{𝑛𝑖,𝑗+1}=𝑢+ℎ𝑢_𝑦+ℎ²2𝑢_𝑦𝑦+ℎ³6𝑢_𝑦𝑦𝑦+ℎ⁴𝑢4!^{𝑦𝑦𝑦𝑦}𝑢+⋯,_{𝑛𝑖+1,𝑗+1}=𝑢+ℎ𝑢_𝑥+ℎ𝑢_𝑦+ℎ²2𝑢_𝑥𝑥+2𝑢_𝑥𝑦+𝑢_𝑦𝑦+ℎ³6𝑢_𝑥𝑥𝑥+3𝑢_𝑥𝑥𝑦+3𝑢_𝑥𝑦𝑦+𝑢_𝑦𝑦𝑦+ℎ⁴𝑢4!^{𝑥𝑥𝑥𝑥}+4𝑢_{𝑥𝑥𝑥𝑦}+6𝑢_{𝑥𝑥𝑦𝑦}+4𝑢_{𝑥𝑦𝑦𝑦}+𝑢_{𝑦𝑦𝑦𝑦}+⋯. (44) Using (2), we convert the temporal derivatives 𝑢_𝑡𝑡, 𝑢_𝑡𝑡𝑡, and 𝑢_{𝑡𝑡𝑡𝑡} into spatial derivatives. We thus have 𝑢_𝑡𝑡=𝛼²_𝑥𝑢_{𝑥𝑥𝑥𝑥}+𝛼²_𝑦𝑢_{𝑦𝑦𝑦𝑦}+𝛼_𝑥𝛼_𝑦𝑢_{𝑥𝑥𝑦𝑦}+𝑢_{𝑦𝑦𝑥𝑥}−2𝛼_𝑥𝛽_𝑥𝑢_𝑥𝑥𝑥−2𝛼_𝑦𝛽_𝑦𝑢_𝑦𝑦𝑦−𝛼_𝑥𝛽_𝑦𝑢_𝑥𝑥𝑦+𝑢_𝑦𝑥𝑥−𝛼_𝑦𝛽_𝑥𝑢_𝑥𝑦𝑦+𝑢_𝑦𝑦𝑥+𝛽²_𝑥𝑢_𝑥𝑥+𝛽²_𝑦𝑢_𝑦𝑦+𝛽_𝑥𝛽_𝑦𝑢_𝑥𝑦+𝑢_𝑦𝑥𝑢+⋯,_𝑡𝑡𝑡=3𝛼_𝑥𝛽²_𝑥𝑢_{𝑥𝑥𝑥𝑥}+3𝛼_𝑦𝛽²_𝑦𝑢_{𝑦𝑦𝑦𝑦}+𝛼_𝑥𝛽²_𝑦𝑢_{𝑦𝑦𝑥𝑥}+𝑢_{𝑦𝑥𝑥𝑦}+𝑢_{𝑥𝑥𝑦𝑦}+𝛼_𝑦𝛽²_𝑥𝑢_{𝑥𝑥𝑦𝑦}+𝑢_{𝑦𝑦𝑥𝑥}+𝑢_{𝑥𝑦𝑦𝑥}+𝛼_𝑦𝛽_𝑥𝛽_𝑦2𝑢_{𝑦𝑦𝑦𝑥}+𝑢_{𝑦𝑦𝑥𝑦}+2𝑢_{𝑥𝑦𝑦𝑦}+2𝑢_{𝑦𝑥𝑦𝑦}+𝛼_𝑥𝛽_𝑥𝛽_𝑦2𝑢_{𝑦𝑥𝑥𝑥}+2𝑢_{𝑥𝑥𝑥𝑦}+𝑢_{𝑥𝑥𝑦𝑥}+2𝑢_{𝑥𝑦𝑥𝑥}−𝛽³_𝑥𝑢_𝑥𝑥𝑥−𝛽_𝑥𝛽²_𝑦𝑢_𝑦𝑦𝑥+𝑢_𝑦𝑥𝑦+𝑢_𝑥𝑦𝑦−𝛽³_𝑦𝑢_𝑦𝑦𝑦−𝛽²_𝑥𝛽_𝑦𝑢_𝑥𝑥𝑦+𝑢_𝑥𝑦𝑥+𝑢_𝑦𝑥𝑥𝑢+⋯,_{𝑡𝑡𝑡𝑡}=𝛽_𝑥𝛽³_𝑥𝑢_{𝑥𝑥𝑥𝑥}+𝛽²_𝑥𝛽_𝑦𝑢_{𝑥𝑥𝑦𝑥}+𝑢_{𝑥𝑦𝑥𝑥}+𝑢_{𝑦𝑥𝑥𝑥}+𝛽_𝑥𝛽²_𝑦𝑢_{𝑦𝑦𝑥𝑥}+𝑢_{𝑦𝑥𝑦𝑥}+_{𝑥𝑦𝑦𝑥}+𝛽³_𝑦𝑢_{𝑦𝑦𝑦𝑥}+𝛽_𝑦𝛽³_𝑦𝑢_{𝑦𝑦𝑦𝑦}+𝛽²_𝑦𝛽_𝑥𝑢_{𝑦𝑦𝑥𝑦}+𝑢_{𝑦𝑥𝑦𝑦}+𝑢_{𝑥𝑦𝑦𝑦}+𝛽_𝑦𝛽²_𝑥𝑢_{𝑥𝑥𝑦𝑦}+𝑢_{𝑥𝑦𝑥𝑦}+𝑢_{𝑦𝑥𝑥𝑦}+𝛽³_𝑥𝑢_{𝑥𝑥𝑥𝑦}+⋯. (45) Then on substituting the Taylor series expansions of each term of the difference scheme we get the following modified equation for the LOD Lax-Wendroff 𝑢_𝑡+𝛽_𝑥𝑢_𝑥+𝛽_𝑦𝑢_𝑦=𝛼_𝑥𝑢_𝑥𝑥+𝛼_𝑦𝑢_𝑦𝑦−16𝛽_𝑥ℎ²1−𝑐²_𝑥−6𝑠_𝑥𝑢_𝑥𝑥𝑥−16𝛽_𝑦ℎ²1−𝑐²_𝑦−6𝑠_𝑦𝑢_𝑦𝑦𝑦−ℎ³24−2𝑠_𝑥−𝑐²_𝑥+12𝑠_𝑥𝑐²_𝑥+12𝑠²_𝑥+𝑐⁴_𝑥𝑢_{𝑥𝑥𝑥𝑥}−ℎ³24−2𝑠_𝑦−𝑐²_𝑦+12𝑠_𝑦𝑐²_𝑦+12𝑠²_𝑦+𝑐⁴_𝑦𝑢_{𝑦𝑦𝑦𝑦}+⋯. (46) The scheme is second order accurate in space and the leading error terms are dispersive in nature (presence of odd-order derivatives 𝑢_𝑥𝑥𝑥 and 𝑢_𝑦𝑦𝑦). As the time and spatial increments go to zero, the modified equation (46) reduces to its original equation, that is, (2). Hence LOD Lax-Wendroff is consistent.

We now study the spectral analysis of the dispersive and dissipative properties of the scheme for the case 𝑐_𝑥=𝑐_𝑦=𝑐 and 𝑠_𝑥=𝑠_𝑦=𝑠. To obtain the amplification factor we use the Von Neumann stability analysis by substituting 𝑢^𝑛_𝑖,𝑗 by 𝜉^𝑛exp(𝐼(𝑖𝜔_𝑥+𝑗𝜔_𝑦)) in (42), where 𝐼=√−1. We thus have 1𝜉=42𝑠_𝑦+𝑐_𝑦+𝑐²_𝑦2𝑠_𝑥+𝑐_𝑥+𝑐²_𝑥exp−𝐼𝜔_𝑥+𝜔_𝑦+122𝑠_𝑦+𝑐_𝑦+𝑐²_𝑦1−2𝑠_𝑥−𝑐²_𝑥exp−𝐼𝜔_𝑦+142𝑠_𝑦+𝑐_𝑦+𝑐²_𝑦2𝑠_𝑥−𝑐_𝑥+𝑐²_𝑥exp𝐼𝜔_𝑥−𝜔_𝑦+121−2𝑠_𝑦−𝑐²_𝑦2𝑠_𝑥+𝑐_𝑥+𝑐²_𝑥exp−𝐼𝜔_𝑥+1−2𝑠_𝑦−𝑐²_𝑦1−2𝑠_𝑥−𝑐²_𝑥+121−2𝑠_𝑦−𝑐²_𝑦2𝑠_𝑥−𝑐_𝑥+𝑐²_𝑥exp𝐼𝜔_𝑥+142𝑠_𝑦−𝑐_𝑦+𝑐²_𝑦2𝑠_𝑥+𝑐_𝑥+𝑐²_𝑥exp−𝐼𝜔_𝑥−𝜔_𝑦+122𝑠_𝑦−𝑐_𝑦+𝑐²_𝑦1−2𝑠_𝑥−𝑐²_𝑥exp𝐼𝜔_𝑦+142𝑠_𝑦−𝑐_𝑦+𝑐²_𝑦2𝑠_𝑥−𝑐_𝑥+𝑐²_𝑥exp𝐼𝜔_𝑥+𝜔_𝑦. (47) The real and imaginary parts of the amplification factor are given by Real(𝜉)=1+4𝑠𝑐²cos𝜔_𝑥cos𝜔_𝑦+𝑐⁴cos𝜔_𝑥cos𝜔_𝑦+4𝑠²cos𝜔_𝑥cos𝜔_𝑦−4𝑠𝑐²cos𝜔_𝑦−𝑐²sin𝜔_𝑥sin𝜔_𝑦+4𝑠²+4𝑠𝑐²−4𝑠𝑐²cos𝜔_𝑥+𝑐²cos𝜔_𝑥−𝑐⁴cos𝜔_𝑥+4𝑠cos𝜔_𝑥+2𝑠cos𝜔_𝑦+2𝑠cos𝜔_𝑥−4𝑠²cos𝜔_𝑦−4𝑠²cos𝜔_𝑥−2𝑐²+𝑐²cos𝜔_𝑦−𝑐⁴cos𝜔_𝑦+𝑐⁴(𝜉)−4𝑠,Imag=−𝑐³cos𝜔_𝑥sin𝜔_𝑦+2𝑠𝑐sin𝜔_𝑥+2𝑠𝑐sin𝜔_𝑦−2𝑠𝑐cos𝜔_𝑥sin𝜔_𝑦−2𝑠𝑐sin𝜔_𝑥cos𝜔_𝑦−𝑐sin𝜔_𝑥+𝑐³sin𝜔_𝑥−𝑐sin𝜔_𝑦+𝑐³sin𝜔_𝑦−𝑐³sin𝜔_𝑥cos𝜔_𝑦, (48) respectively. The modulus of the amplification factor, AFM, is obtained as AFM=(ℜ(𝜉))²+(ℑ(𝜉))². (49) We find the region of stability using the approach of Hindmarsh et al. [16] and Sousa [17]. We consider the case when 𝜔_𝑥=𝜔_𝑦=𝜋 and 𝜔_𝑥→0 and 𝜔_𝑦→0. When 𝜔_𝑥=𝜔_𝑦=𝜋, (49) gives 𝜉=√1+16𝑠𝑐²−4𝑐²+4𝑐⁴+16𝑠²−8𝑠. (50) From the Von Neumann stability analysis, the scheme is stable if and only if |𝜉|≤1. Thus, we get 1+16𝑠𝑐²−4𝑐²+4𝑐⁴+16𝑠²−8𝑠≤1. (51)

On simplification, we get 𝑠+2𝑐²−14≤14, (52) which reduces to 𝑠≤1−𝑐²2. (53)

When 𝜔_𝑥→0 and 𝜔_𝑦→0, we use the following approximations: cos𝜔_𝑥𝜔≃1−^2𝑥2,cos𝜔_𝑦𝜔≃1−^2𝑦2,sin𝜔_𝑥≃𝜔_𝑥,sin𝜔_𝑦≃𝜔_𝑦. (54) We consider (49) and use the approximations in (54) to obtain 𝜉²≃1−2𝑠𝜔²_𝑥+𝜔²_𝑦. (55)

Thus LOD Lax-Wendroff is stable if −2𝑠≤0. Therefore, we have 𝑠≥0. (56)

On combining (53) and (56), we obtain the region of stability for the LOD Lax-Wendroff procedure as 0≤𝑠≤1−𝑐²2. (57) We choose 𝛼=0.01 and 𝛽=0.8 [3]. For ℎ=0.025 from (57), we have 0≤0.01𝑘(0.025)²≤1−(0.8𝑘/0.025)²2. (58) On solving for 𝑘, we get 0≤𝑘≤0.0193. (59)

Therefore, for ℎ=0.025, the stability region for the LOD Lax-Wendroff procedure is 0≤𝑘≤0.0193. We choose 𝑘 such that 0.3/𝑘 is an integer as for our numerical experiments, 𝑇=0.3. Then 0.3/𝑘 gives the number of time intervals. We choose 𝑘=0.005,0.01,3/160, and 1/60.

For ℎ=0.05, we have 𝑐=16𝑘 and 𝑠=4𝑘. When 𝜔_𝑥=𝜋 and 𝜔_𝑦=𝜋, using (51) we get 0≤1−32𝑘+1638𝑘³−768𝑘²+262144𝑘⁴≤1, (60) which gives 0≤𝑘≤0.0487985. (61) Therefore, for ℎ=0.05, LOD Lax-Wendroff is stable if 0≤𝑘≤0.0487985. We then choose some values of 𝑘∈(0,0.048) for the numerical experiments. for our numerical experiments.

3D plots of the modulus of the amplification factor versus phase angles in 𝑥- and 𝑦-directions for some different values of ℎ and 𝑘 are depicted in Figures 1 and 2. 2D plots of the modulus of the amplification factor versus 𝜔_𝑥, when 𝜔_𝑦=0, are shown in Figure 3. The scheme is in general less dissipative at ℎ=0.05 as compared to ℎ=0.025. Out of the eight combinations of values of ℎ and 𝑘, the scheme is least dissipative when ℎ=0.05 and 𝑘=0.005 and it is most dissipative when ℎ=0.025 and 𝑘=0.01.

(a) ℎ=0.025,𝑘=0.005

(b) ℎ=0.025,𝑘=0.01

(c) ℎ=0.025,𝑘=1/60

(d) ℎ=0.025,𝑘=3/160

Figure 1 

Plots of modulus of amplification factor versus phase angle in 𝑥-direction, 𝜔𝑥, versus phase angle in 𝑦-direction, 𝜔𝑦, at ℎ=0.025 with some different values of 𝑘 for the LOD Lax-Wendroff scheme.

(a) ℎ=0.05,𝑘=0.005

(b) ℎ=0.05,𝑘=0.01

(c) ℎ=0.05,𝑘=0.02

(d) ℎ=0.05,𝑘=0.03

Figure 2 

Plots of modulus of amplification factor versus phase angle in 𝑥-direction, 𝜔𝑥, versus phase angle in 𝑦-direction, 𝜔𝑦, at ℎ=0.05 with some different values of 𝑘 for the LOD Lax-Wendroff scheme.

(a) ℎ=0.025

(b) ℎ=0.05

Figure 3 

Plots of modulus of amplification factor versus 𝜔𝑥 when 𝜔𝑦=0 at ℎ=0.025 and ℎ=0.05 with some different values of 𝑘 for the LOD Lax-Wendroff scheme.

Figures 4 and 5 show the 3D plots of the relative phase error versus 𝜔_𝑥 versus 𝜔_𝑦 for some different values of ℎ and 𝑘. Figure 6 shows the 2D plots of the relative phase error versus 𝜔_𝑥, for the case 𝜔_𝑦=0 at ℎ=0.025 and ℎ=0.05. For ℎ=0.025, we observe phase lag behaviour at 𝑘=0.01 and 𝑘=0.005 and phase lead behaviour at 𝑘=1/60 and 𝑘=3/160. For ℎ=0.05, we have phase lag behaviour at 𝑘=0.005,0.01,0.02, and 0.03. The scheme is least dispersive when ℎ=0.05; 𝑘=0.03 and ℎ=0.025; 𝑘=0.01. In the following section, we consider the LOD (1,5) scheme.

(a) ℎ=0.025,𝑘=0.005

(b) ℎ=0.025,𝑘=0.01

(c) ℎ=0.025,𝑘=1/60

(d) ℎ=0.025,𝑘=3/160

Figure 4 

Plots of relative phase error versus phase angle in 𝑥-direction (𝜔𝑥) versus phase angle in 𝑦-direction (𝜔−𝑦), at ℎ=0.025 with some different values of 𝑘 for the LOD Lax-Wendroff scheme.

(a) ℎ=0.05,𝑘=0.005

(b) ℎ=0.05,𝑘=0.01

(c) ℎ=0.05,𝑘=0.02

(d) ℎ=0.05,𝑘=0.03

Figure 5 

Plots of relative phase error versus 𝜔𝑥 versus 𝜔𝑦, at ℎ=0.05 with some different values of 𝑘 for the LOD Lax-Wendroff scheme.

(a) ℎ=0.025

(b) ℎ=0.05

Figure 6 

Plots of relative phase error versus 𝜔𝑥 when 𝜔𝑦=0 at ℎ=0.025 and ℎ=0.05 with some different values of 𝑘 for the LOD Lax-Wendroff scheme.

5. LOD (1,5) Explicit Procedure

In this procedure the following approximations are used in the first half time step [3]: 𝜕𝑢𝜕𝑡^𝑛_𝑖,𝑗≃𝑢^𝑛+(1/2)_𝑖,𝑗−𝑢^𝑛_𝑖,𝑗,Δ𝑡𝜕𝑢𝜕𝑥^𝑛_𝑖,𝑗≃12𝑠_𝑥+2𝑐²_𝑥−3𝑐_𝑥−2𝑢12^{𝑛𝑖+2,𝑗}−𝑢^𝑛_𝑖,𝑗+2Δ𝑥12𝑠_𝑥+2𝑐²_𝑥+3𝑐_𝑥−2𝑢12^{𝑛𝑖,𝑗}−𝑢^𝑛_{𝑖−2,𝑗}−𝑐2Δ𝑥_2𝑥+6𝑠_𝑥−4𝑢12^{𝑛𝑖+1,𝑗}−𝑢^𝑛_{𝑖−1,𝑗},𝜕2Δ𝑥₂𝑢𝜕𝑥^2𝑛_𝑖,𝑗≃−𝑐⁴_𝑥+4𝑐²_𝑥−12𝑠²_𝑥−12𝑠_𝑥𝑐²_𝑥+8𝑠_𝑥6𝑠_𝑥