Numerical solutions for the fourth order extended fisher-kolmogorov equation with high accuracy by d
* Author to whom correspondence should be addressed.
Sigma Journal of Engineering and Natural Sciences 2018, Vol. 9, Suppl. 3, pp. 273-284; doi.org/10.62051/ytu.sigma-journal-of-engineering-and-natural-sciences-numerical-solutions-for-the-fourth-order-extended-fisher-kolmogorov-equation-wit
Abstract
Keywords: Partial differential equations; differential quadrature method; EFK equation; modified cubic B- Splines.
1. Introduction
In this study, we have investigated numerical solutions of fourth order extended FisherKolmogorov (EFK) equation by using modified cubic B-spline Differential Quadrature Method (MCB-DQM). The EFK equation is given in the following form 𝑢𝑡 + 𝛾𝑢𝑥𝑥𝑥𝑥 − 𝑢𝑥𝑥 + 𝜓(𝑢) = 0, 𝑥 ∈ [𝑎, 𝑏], 𝑡 ∈ [0, 𝑇] where 𝑢 ∈ [𝑎, 𝑏] × [0, 𝑇], 𝜓(𝑢) = 𝑢3 − 𝑢 and 𝛾 > 0.
Coullet et al.[1] and van Saarlos [2-4] and Dee and van Saarlos [5] introduced the EFK equation given in the Eq. (1) by adding a stabilizing fourth-order derivative for the value of the 𝛾 ≠ 0. For the value of the 𝛾 = 0 the Eq. (1) turns into the standard Fisher-Kolmogorov (FK) equation.
Corresponding Author: e-mail: alibashan@gmail.com, tel: (372) 291 17 88 / 1788
A. Başhan, Y. Uçar, N.M. Yağmurlu, A. Esen / Sigma J Eng & Nat Sci 9 (3), 273-284, 2018
EFK equ ation has many applications in several scientific field such as spatiotemporal chaos [1] pattern formation in bi-stable systems [5], propagation of domain walls in liquid crystals [6], phase transition near a Lifshitz point [7,8]. EFK equation (1) has been solved with various methods by many researchers. Among others, Danumjaya and Pani [9] used orthogonal cubic spline collocation method, Mittal and Arora [10] used quantic B-spline collocation, Mittal and Dahiya [11] used quintic B-spline differential quadrature method. In this paper, we investigated numerical solutions of EFK equation via modified cubic B-spline differential quadrature method since it has high accurate solutions. In recent years, DQM which is first introduced by Bellman et al. [12] 1972, has had wide application areas because of its using considerably less number of grid points. DQM is a method in which partial derivative of a function with respect to a coordinate direction is expressed as a linear weighted sum of all the functional values at all mesh points along that direction[13]. Many researchers have developed various types of DQM by using different base functions such as Legendre polynomials and spline functions [12,14-17], Lagrange interpolation polynomials [1820], Lagrange interpolated trigonometric polynomials [21], Hermite polynomials [22], radial basis functions [23], harmonic functions [24], Sinc functions [25,26], B-spline functions [27-33]. In the present study, MCB-DQM is going to be applied to obtain approximate solutions of the EFK equation. MCB-DQM has been preferred for solving fourth-order EFK equation since it has high accurate solutions, low storage allocation and simple applicability to equation.
2. Modified Cubic B-Spline DQM
We are going to consider the Eq. (1) with the boundary conditions taken from: 𝑈(𝑎, 𝑡) = 𝑔1 (𝑡),
with the initial condition 𝑈(𝑥, 0) = 𝑢0 (𝑥), 𝑎 ≤ 𝑥 ≤ 𝑏.
Let us take the uniform grid distribution 𝑎 = 𝑥1 < 𝑥2 < ⋯ < 𝑥𝑁 = 𝑏 of a finite interval [𝑎, 𝑏] into consideration. Provided that any given function 𝑈(𝑥) is smooth enough over the solution domain, its derivatives with respect to 𝑥 at a grid point 𝑥𝑖 can be approximated by a linear summation of all the functional values in the solution domain, namely, 𝑑 (𝑟) 𝑈
where 𝑟 represents the order of derivative, 𝑤𝑖𝑗 ’s denote the weighting coefficients of the 𝑟 th order derivative approximation, and 𝑁 represents the number of nodal points in the given (𝑟) solution domain. Here, the index 𝑗 indicates the fact that 𝑤𝑖𝑗 is the corresponding weighting coefficient of the value of the function 𝑈(𝑥𝑗 ). In this work, we are going to need the second and fourth order derivatives of the function 𝑈(𝑥). So, we are going to begin to find the value of the equation (4) for the 𝑟 = 1. The main idea behind DQM approximation is to find out the corresponding weighting (𝑟) coefficients 𝑤𝑖𝑗 by means of a set of base functions spanning the problem domain. While determining the corresponding weighting coefficients, different basis may be used. In the present study, we will try to compute weighting coefficients with modified cubic B-spline basis. Let ∁𝑚 (𝑥) be the cubic B-splines with knots at the points 𝑥𝑖 where the uniformly distributed N grid points are taken as 𝑎 = 𝑥1 < 𝑥2 < ⋯ < 𝑥𝑁 = 𝑏 on the ordinary real axis. Then, the cubic B-splines {∁0 , ∁1 , … , ∁𝑁+1 } form a basis for the functions defined over [𝑎, 𝑏]. The cubic B-splines ∁𝑚 (𝑥) are defined by the relationships:
Numerical Solutions for the Fourth Order Extended … / Sigma J Eng & Nat Sci 9 (3), 273-284, 2018
(𝑥 − 𝑥𝑚−2 )3 , [𝑥𝑚−2 , 𝑥𝑚−1 ] (𝑥 − 𝑥𝑚−2 )3 − 4(𝑥 − 𝑥𝑚−1 )3 , [𝑥𝑚−1 , 𝑥𝑚 ] 1 ∁𝑚 (𝑥) = ℎ3 (𝑥𝑚+2 − 𝑥)3 − 4(𝑥𝑚+1 − 𝑥)3 , [𝑥𝑚 , 𝑥𝑚+1 ] 3 (𝑥𝑚+2 − 𝑥) , [𝑥𝑚+1 , 𝑥𝑚+2 ] {0 , 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 where ℎ = 𝑥𝑚 − 𝑥𝑚−1 for all 𝑚 [34]. Using the modified cubic B-splines results in a diagonally dominant matrix system of equations. This structur e has great importance for the stability analysis. Modification of cubic Bsplines can be carried out differently. Among others, Mittal and Jain [35] have introduced modified cubic B-splines at the grid points as follows ∅1 (𝑥) = ∁1 (𝑥) + 2∁0 (𝑥) ∅2 (𝑥) = ∁2 (𝑥) − ∁0 (𝑥) ∅𝑠 (𝑥) = ∁𝑠 (𝑥), for 𝑠 = 3,4, … , 𝑁 − 2 ∅𝑁−1 (𝑥) = ∁𝑁−1 (𝑥) − ∁𝑁+1 (𝑥) ∅𝑁 (𝑥) = ∁𝑁 (𝑥) + 2∁𝑁+1 (𝑥)
where ∅𝑘 , (𝑘 = 1, 2, … , 𝑁) forms a basis functions over the [𝑎, 𝑏] domain.
2.1. Weighting Coefficients of the First Order Derivative
From Eq. (4) with value of 𝑟 = 1, we have obtained the following equation (1)
∅′𝑘 (𝑥𝑖 ) = ∑𝑁 𝑗=1 𝑤𝑖𝑗 ∅𝑘 (𝑥𝑗 ) for 𝑖 = 1, 2, … , 𝑁; 𝑘 = 1, 2, … , 𝑁
For the first grid point 𝑥1 , from Eq. (6) we get an equation in the following form (1)
∅′𝑘 (𝑥1 ) = ∑𝑁 𝑗=1 𝑤1𝑗 ∅𝑘 (𝑥𝑗 ) for 𝑘 = 1, 2, … , 𝑁
and by using the value of modified cubic basis functions 6 0
0 𝑤1,𝑁−1 0 [ 6 ] 𝑤 (1) [ 1,𝑁 ] [ 0 ] the above equation system is obtained. Similarly, by using the value of modified cubic basis functions at the grid points 𝑥𝑖 , (2 ≤ 𝑖 ≤ 𝑁 − 1), respectively, (1)
the above equation system is obtained. For the last grid point 𝑥𝑁
A. Başhan, Y. Uçar, N.M. Yağmurlu, A. Esen / Sigma J Eng & Nat Sci 9 (3), 273-284, 2018
the above equation system is obtained. So, weighting coefficients 𝑤𝑖,𝑗 which are related to the 𝑥𝑖 (𝑖 = 1, 2, … , 𝑁) are found quite easily by solving system of Eqs. (8)-(10) with Thomas algorithm.
2.2. Weighting Coefficients of the Second Order Derivative
This method depends on the first order weighting coefficients when obtaining weighting coefficients of the second order derivatives. By the Shu's recurrence formula, the second order weighting coefficients are determined for 𝑖 = 1, 2, … , 𝑁 and 𝑗 = 1, 2, … , 𝑁 as below [13]: (2)
2.3. Weighting Coefficients of the Fourth Order Derivative
This method depends on the second order weighting coefficients to obtain the weighting coefficients of the fourth order derivatives. By the matrix multiplication approach, the fourth order weighting coefficients are determined as below [13]: [𝐴(𝑚) ] = [𝐴(1) ][𝐴(𝑚−1) ] = [𝐴(𝑚−1) ][𝐴(1) ] (𝑚−1)
where [𝐴 ] and [𝐴 ] are the weighting coefficients matrices of the (𝑚 − 1) − 𝑡ℎ and 𝑚 − 𝑡ℎ order derivatives, respectively. Although Eq. (13) looks simple, it involves more arithmetic operations as compared to Eqs. (11) and (12). It is noted that the calculation of each weighting coefficient by Eq. (13) involves N multiplications and (N-1) additions, i.e. a total of (2N-1) arithmetic operations. On the other hand, Shu's recurrence relationship (11) involves two multiplications, one division and one subtraction, i.e. a total of four arithmetic operations [13].
3. Numerical Discretization
The EFK equation of the form 𝑈𝑡 + 𝛾𝑈𝑥𝑥𝑥𝑥 − 𝑈𝑥𝑥 + 𝑈 3 − 𝑈 = 0
with the boundary conditions (2) and the initial condition (3) is rewritten as 𝑈𝑡 = −𝛾𝑈𝑥𝑥𝑥𝑥 + 𝑈𝑥𝑥 − 𝑈 3 + 𝑈.
Then, the differential quadrature derivative approximations of the second and the fourth orders have been used in Eq. (15) 𝑑𝑈(𝑥𝑖 ) 𝑑𝑡
and ordinary differential equation (16) is obtained. Then, the ordinary differential equation given by Eq. (16) is integrated with respect to time. Here, we have preferred strong stability-
Numerical Solutions for the Fourth Order Extended … / Sigma J Eng & Nat Sci 9 (3), 273-284, 2018
preserving low storage Runge-Kutta43 method[36] due to its advantages such as accuracy, stability and memory allocation properties.
4. Numerical Examples And Stability
Here, we have obtained the numerical solutions of the EFK by the MCB-DQM. The accuracy of the numerical method is checked using the error norms 𝐿2 and 𝐿∞ , respectively: 𝑁
𝐿∞ = ‖𝑈 − 𝑢‖∞ ≅ max|𝑈𝑗 − 𝑢𝑗 |, 𝑗 = 1,2, … 𝑁 − 1. 𝑗
Since the analytical solution of EFK equation does not exist, newly obtained numerical solution is compared with those solutions obtained when grid number is taken N=160 instead of exact solution. Stability analysis of a numerical method for a nonlinear differential equation requires the determination of eigenvalues of coefficient matrices. With the numerical discretization of partial differential equation EFK, it turns into an ordinary differential equation. The stability of a time-dependent problem: 𝜕𝑈 𝜕𝑡
with the proper initial and boundary conditions, where 𝑙 is a spatial differential operator. After discretization with DQM, Eq. (17) is reduced into a set of ordinary differential equations in time as follows 𝑑{𝑈} 𝑑𝑡
where {𝑢} is an unknown vector of the functional values at the grid points except the left and right boundary points, {𝑏} is a vector containing the non-homogenous part and the boundary conditions and 𝐴 is the coefficient matrix. The stability of a numerical scheme for numerical integration of Eq. (18) depends on the stability of the ordinary differential Eq. (18). If the ordinary differential Eq. (18) is not stable, numerical methods may not generate converged solutions. The stability of Eq. (18) is related to the eigenvalues of the matrix 𝐴, since its exact solution is directly determined by the eigenvalues of the matrix 𝐴. When all 𝑅𝑒(𝜆𝑖 ) ≤ 0 for all it is enough to show the stability of the exact solution of {𝑢} as 𝑡 → ∞ where 𝑅𝑒(𝜆𝑖 ) denotes the real part of the (4) (2) eigenvalues 𝜆𝑖 of the matrix 𝐴. The matrix 𝐴 in Eq. (18) is determined as 𝐴𝑖𝑗 = −𝛾𝑤𝑖𝑗 + 𝑤𝑖𝑗 − 3 𝛼𝑖 + 𝛼𝑖 where 𝛼𝑖 = 𝑈(𝑥𝑖 , 𝑡) [13]. The eigenvalues of matrix 𝐴 should be in the stability region as shown in Figure 1 [37].
A. Başhan, Y. Uçar, N.M. Yağmurlu, A. Esen / Sigma J Eng & Nat Sci 9 (3), 273-284, 2018
Figure 1. Stability regions of fourth order SSPRK eigenvalues.
4.1. Test Problem 1
The first test problem has initial condition as follows: 𝑈(𝑥, 0) = −𝑠𝑖𝑛(𝜋𝑥)
with boundary conditions 𝑈(𝑥0 , 𝑡) = 𝑈(𝑥𝑁 , 𝑡) = 0, at the domain −4 ≤ 𝑥 ≤ 4.
We fix the number of grid points N=81 and time increment ∆𝑡 = 0.0001 when 𝛾 = 0, 𝛾 = 0.0001 and 𝛾 = 0.1, respectively. Numerical simulations are given in Figures 2-4.
Numerical Solutions for the Fourth Order Extended … / Sigma J Eng & Nat Sci 9 (3), 273-284, 2018
Figure 4. Simulations for 𝛾 = 0.1. Table 1. 𝐿2 and 𝐿∞ error norms at 𝑡 = 0.2. N 20 40 80
Present (MCB-DQM) 𝐿2 𝐿∞ 0.01881 0.01839 0.00230 0.00222 0.00024 0.00023
Quin. Coll. [10] 𝐿2 𝐿∞ 0.01158 0.00551 0.00282 0.00134 0.00057 0.00028
Quin. DQM[11] 𝐿2 𝐿∞ 0.02135 0.00116 0.00222 0.00122 0.00031 0.00015
As it seen straightforwardly from Figure 2 and Figure 3 that the behaviour of solutions for 𝛾 = 0 and 𝛾 = 0.0001 are similar to each other. Except for 𝛾 = 0.1 simulations given in Figure 4 that solutions decline to 0 very rapidly because of stabilizing behaviour of EFK. The calculated and compared values of the error norms 𝐿2 and 𝐿∞ are given in Table 1. As it is seen straightforward from comparison with earlier works [10, 11] given in Table 1 that the present error norms 𝐿2 and 𝐿∞ are acceptably good.
A. Başhan, Y. Uçar, N.M. Yağmurlu, A. Esen / Sigma J Eng & Nat Sci 9 (3), 273-284, 2018
4.3. Test Problem 2
The second test problem has initial condition as follows: 𝑈(𝑥, 0) = 10−3 𝑒𝑥𝑝(−𝑥 2 )
with boundary conditions 𝑈(𝑥0 , 𝑡) = 𝑈(𝑥𝑁 , 𝑡) = 1, at the domain −4 ≤ 𝑥 ≤ 4.
The behaviours of solutions for time running up from 𝑡 = 0.25 to 𝑡 = 4.5 are given in Figure 5. As it seen from Figure 5 that approximate solution of 𝑈 declines as time increases and eventually it comes close to the value 1. Numerical results are agreeable to those given in [9].
4.3. Test Problem 3
Our last test problem has initial condition 𝑈(𝑥, 0) = −10−3 𝑒𝑥𝑝(−𝑥 2 )
with boundary conditions 𝑈(𝑥0 , 𝑡) = 𝑈(𝑥𝑁 , 𝑡) = −1, at the domain −4 ≤ 𝑥 ≤ 4.
The behaviours of solutions for time running up from 𝑡 = 0.25 to 𝑡 = 4.5 are given in Figure 6. As it seen from Figure 6 that the approximate solution of 𝑈 decline by time increases and eventually it come close to the value -1. Numerical results are agreeable to the as given in [9].
Numerical Solutions for the Fourth Order Extended … / Sigma J Eng & Nat Sci 9 (3), 273-284, 2018
Figure 6. Simulations for 𝛾 = 0.0001. A matrix stability analysis is also carried out for the MCB-DQM. We have used matlab r2013b program to obtain the eigenvalues of the coefficient matrix for the first test problem. Eigenvalues of the suggested method for N=21, N=41, N=81 and N=161 are presented in Figure 7. None of them has imaginary parts. All the eigenvalues are in convenience with stability criteria [37].
Figure 7. Eigenvalues for various number of grid points. Also, maximum absolute value of eigenvalues for various value of 𝛾 = 0, 𝛾 = 0.0001 and 𝛾 = 0.1 for different number of grid points N=21 , N=41, N=81 and N=161 are tabulated in Table 2. All results are pure real eigenvalues when 𝛾 = 0, 𝛾 = 0.0001 and 𝛾 = 0.1. As it seen from
A. Başhan, Y. Uçar, N.M. Yağmurlu, A. Esen / Sigma J Eng & Nat Sci 9 (3), 273-284, 2018
Table 2 when the number of the grid points are increased, absolute value of eigenvalue grows, so time step should be decrease to obtain the stable solution. Table 2. Maximum absolute value of eigenvalues at various number of grid points. Grid Number 𝑀𝑎𝑥|𝑅𝑒(𝜆)| 𝑀𝑎𝑥|𝐼𝑚(𝜆)| Grid Number 𝑀𝑎𝑥|𝑅𝑒(𝜆)| 𝑀𝑎𝑥|𝐼𝑚(𝜆)| Grid Number 𝑀𝑎𝑥|𝑅𝑒(𝜆)| Max|Im(λ)|
5. Conclusion
In this study, we have implemented MCB-DQM for numerical solution of EFK equation. To obtain the second order weighting coefficients, we used Shu's recurrence formulae and to obtain the fourth order weighting coefficients we used the matrix multiplication approach. The performance and accuracy of the method have been shown by calculating the error norms 𝐿2 and 𝐿∞ . As can be observed by the comparison between the obtained values of the error norms of the present method and earlier works, MCB-DQM results are acceptable good. The obtained results show that MCB-DQM can be used to produce reasonable accurate numerical solutions of the EFK equation. As the results of investigation of eigenvalues, the present method is stable. So, MCBDQM is a reliable one for getting the numerical solutions of some physically important nonlinear problems.
Share and Cite
BAŞHAN, A.; UÇAR, Y.; YAĞMURLU, N.M.; ESEN, A. Numerical solutions for the fourth order extended fisher-kolmogorov equation with high accuracy by d. Sigma Journal of Engineering and Natural Sciences 2018, Vol. 9, pp. 273-284. https://doi.org/10.62051/ytu.sigma-journal-of-engineering-and-natural-sciences-numerical-solutions-for-the-fourth-order-extended-fisher-kolmogorov-equation-wit