Numerical Simulation that Provides Non-oscillatory Solutions for Porous Medium Equation and Error Approximation of Boussinesq’s Equation

So far, so many works have been done in a different way to find the exact track that can grasp the true solution for the one-dimensional porous media equation (PME). For instance, Monika studied using relaxation [1], Q. Zhang and Z. Wu. have already done the similar work by local degenerate Galarkin (LDG) method [12] and so on. Still, this is a challenge to find an appropriate scheme that can track the true solution when adiabatic exponent increases monotonically. In this paper, we have studied numerical result for where we have used Explicit-Implicit Finite Difference Method (EIFDM). Since so far PME is a degenerate parabolic equation and analytically the existence and uniqueness occur weakly only in the Sobolev sense, it is very hard to track the true solution numerically. Our main objective is to study numerically the PME with mixed boundary conditions and shown the result is helpful to track the Barenblatt’s self similar solution and its interface when adiabatic exponent larger than 3 that provides much less error. This paper will show a possibility that Finite Difference Method (FDM) is also helpful rather the Finite Elements Method to track the interface in the simulation with an appropriate initial guess. Also checked L1 ,L2 and L∞-error for Boussinesq’s equation which is a fundamental equation of groundwater flow, hopefully, the simulated results can help when this equation is useful in the practical world. Finally, all studied results are given to show the advantage of the θ-scheme method in the simulation of the PME and its capability to capture accurately sharp interfaces without oscillation.


Introduction
The aim of this paper is to study numerical approximation of the nonlinear diffusion equation.
and usually it is called the porous medium equation (PME), with due attention paid to its closest relatives. The default settings are:ρ = ρ(x, t) is a non-negative scalar function of space x ∈ R d and time t ∈ R, the space dimension is d ≥ 1, and γ is a constant larger than 1 and ∆ represents the Laplacian form. The equation can be posed for all x ∈ R d and 0 < t < ∞, and then initial conditions are needed to determine the solutions; but it is quite often posed, especially in practical problems, in a bounded sub-domain Ω ∈ R d for 0 < t < T , and then determination of a unique solution asks for boundary conditions as well as initial conditions. [4,6] This equation is one of the simplest examples of a nonlinear evolution equation of parabolic type. It appears in the description of different natural phenomena, and its theory and properties depart strongly from the heat equation, u t = ∆u, its most famous relative. Hence the interest of its analytical and numerical study is an important task for both for the pure mathematicians and the applied statisticians.
Many of those mathematical models are described by a (system of) partial differential equations (PDEs). Well established laws of nature and their mathematical counterparts have led to most of the development of suitable PDEs, e.g., heat conduction, fluid dynamics or deformation of solids [7,4,11,9]. However, sometimes the development and understanding of a particular mathematical model can only be tackled by a trial and error approach, which is a two-step procedure. In a first step, the simulated results are compared to experimental data and in a second step, these comparisons are used to modify the mathematical model, i.e., most of the case it is the underlying PDE. Hence, numerical simulations of PDEs are of tremendous importance when trying to understand real world processes. On the one hand, in the field of numerical ODEs highly valuable methods and results exist which are of practical use for solving time-dependent PDEs, something which is often not fully exploited by numerical PDE researchers [9].
Although many problems can be solved by Euler's method or the Crank-Nicolson method, better alternatives are often available which can significantly reduce the computational 26 Numerical Simulation that Provides Non-oscillatory Solutions for Porous Medium Equation and Error Approximation of Boussinesq's Equation effort needed to solve practical problems [7]. On the other hand, many numerical ODE researchers are unaware of the vast amount of highly interesting results on discretization methods for PDEs [9]. Moreover, when solving PDEs, discretizations in space and time have to be matched, and different spatial discretizations may require different temporal discretizations. The phenomena modeled by partial differential equations become increasingly complicated, and so do the partial differential equations themselves. Often, one wishes a model to capture different aspects of a situation, for instance both convective transport and dispersive oscillations on a small scale. These different aspects of the model are then reflected in a partial differential equation, which may contain terms (operators) that are mathematically very different, making these models hard to analyze, both theoretically and numerically. (See [1,9,10,11,15,18,19,20,21,22,24]).
In order to study numerically any kind of problem we need a stable and consistent numerical scheme [7,11,9]. Besides, a good scheme has to reproduce all of the important features of the original model, which arise from the physical background of the problem. First of all the scheme has to preserve the non negativity of solutions as we deal with densities and concentrations. Then, if we consider bounded domains with no-flux boundary conditions, the numerical approximation has to conserve the total mass. Conservation laws with reaction terms are characterized by a special balance between the fluxes and the sources, which can lead to non constant stationary solutions. Their preservation is essential and in the case of geometric sources, containing for example space derivatives, it is impossible by using standard schemes treating the flux and the source at different time steps. Moreover, the presence of the vacuum states brings another difficulty as many schemes produce oscillations at the interface between the regions, where the density is strictly positive and where it vanishes. In order to study the behaviour of the boundary, the approximate solution not only has to be free from oscillations, but also it has to deal with tracking the movement of exact front. Using an appropriate numerical scheme for parabolic model of porous medium type that can satisfy all these properties is not an easy task (see [3,12]).
In this paper, we focus our attention on degenerate parabolic equations such as the porous medium equation (PME). It is well known that the degeneracy implies a finite speed of propagation [4,6]. Equivalently it means that when the initial density has compact support, it will remain compact for all times [14]. In this case, the classical smooth solution may not always exist in general, even if the initial solution is smooth. It is necessary to consider the weak energy solution, whose behavior causes many difficulties for a good numerical simulation. For example, the weak solution may lose its classical derivative at some (interface) points, and the sharp interface of support may propagate with finite speed if the initial data have compact support. Enamored of these interesting facts, there have been many works on the simu-lation for the non smooth solution of the PME, for example, the finite different method by Graveleau and Jamet [25], the interface tracking algorithm by DiBenedetto and Hoff [26], and the relaxation scheme referred to in [1].
In numerical point of view, implicit θ-scheme possesses several properties to make it very attractive for practical computations, such as parallelization, adaptivity, and simple treatment of boundary conditions. The most important properties of this method is its strong stability and high-order accuracy; as a result, it is very good at capturing discontinuous jumps and sharp transient layers. Explicit relaxation scheme has a good L 2 and L ∞ -stability for the standard nonlinear diffusion equation, including the PME for short time simulation where as for long time time simulation explicit-implicit Crank-Nicolson scheme is better in the sense of L 2 and L ∞stability [1]. We compared between fully implicit Backward Euler scheme and explicit-implicit Crank-Nicolson scheme and measured the accuracy by the L 1 , L 2 and L ∞ -norm.
There are two main components in this paper. The first is the simulation by the finite difference method (FDM) for smooth solutions of the PME (see [6]). For this, we design a non-negativity-preserving the physical relevancy of the numerical solutions (Q. Zhang and Z. Wu. has been already done by local degenerate Galarkin (LDG) method, see [12], however, his works was only for initial value problem, where we have considered boundary value problem with homogenous Neuman boundary condition). The given numerical results verify the above advantage of the FDM method that it has the ability to capture sharp interfaces accurately without numerical oscillation (see the Figure 7). As a comparison, time integration method based of non-negativity is also considered to show the importance to avoid the numerical difficulties. Here we want to mention shortly that, before us, researchers mostly considered either Dirichlet or Neuman boundary conditions on both sides. So,our work clearly much more effective when people wants to use PME to study cell migrations and for the treatment of cornea where mathematician usually use PME to formulate a mathematical model for the blood flow in the vessels of humanł or animalł eyes and it is required to use mixed boundary conditions over there.
The second component is an analysis for the nonnegativity preservation principle of the considered θ-method, i.e., in each cell of numerical discretization of the Barenblatt solution [8] for the PME remains non-negative for different adiabatic exponents (see the Figure 7). At the same time, this paper will point out, anyone who works on PME, need to take care of the oscillation for the the PME, for any choice of the adiabatic parameter ensures the stability not only in the L 2 sense but also in the sense of L ∞ , then it would be guaranteed that the approximation to the desired solution response will lie within the illustrated bounds.
To learn more about numerical methods for parabolic equations including PME, we invite the reader to take a closer look at one or several of the references [7, 9, (1), a solution of self-similar form can be recovered in a similar way as for the linear diffusion equation. We therefore look for a solution of the form in one dimension, to the one-dimensional PME where, P (ρ) = ρ γ , γ > 1.
First of all, the conservation of the total mass easily implies and therefore we have, β = α.
In the new variablesρ and ζ = xt −α , we obtain, which requires the choice in order to get rid of the time variable. The above can be written as In the domain whereρ > 0, we can imposê , where, C is any arbitrary constant. Clearly, the above gives problems in case,ρ becomes negative. Since,ρ = 0 on an interval solves the equation forρ, we can introduce the following solution where, ρ + = max(ρ, 0) and equation (5) is called Barenblatt solution. Such a solution has compact support on the interval which grows as t grows. Moreover, it can be easily proven that B (.; 0) is a multiple of the Dirac delta distribution, exactly as it is the case for the Gaussian solution of the linear diffusion equation. Analyzing the differences with the linear diffusion case, we clearly can see that the Barenblatt solution is not smooth, opposite to the Gaussian solution (G) to the linear diffusion equation, which is C ∞ . More precisely, the Barenblatt solution has possible lack of smoothness at the boundary of its support, where it eventually features discontinuities in a space derivative of a certain order depending on the adiabatic exponent γ. At any other points, B is C ∞ . A key difference between G and B is that the support of G becomes unbounded immediately after t = 0 (in finite speed of propagation) whereas the support of B is compact for all times, it travels with finite speed.

Numerical Schemes
In this section we are going to present numerical schemes based on finite difference method. In particular, we consider two different kind of the implicit method (pure implicit scheme and Crank-Nicolson implicit-explicit scheme) on an uniform grid with the mesh size ∆x = h and time step ∆t = k. In one dimension on an interval [0, L]×[0, T ] we define nodes x i , i = 1; . . . ; s such that x 1 = 0 and x s = L. Denote by ρ n i an approximation of a function ρ(x, t) at node x i at time t n , where, t ∈ [0, T ] and for any function f (ρ) denote f (ρ n i ) = f n i for simplicity. This notation is going to be used throughout the paper. Since we know, in θ− scheme, where, θ = 1, scheme is called backward Euler scheme which is fully implicit scheme and θ = 1/2, scheme is called Crank-Nicolson scheme which is implicit-explicit scheme. we want to compare numerical error between these two schemes and for this reason, first we are going to present θ− scheme for PME.
Numerical Simulation that Provides Non-oscillatory Solutions for Porous Medium Equation and Error Approximation of Boussinesq's Equation When the parameter γ increases, the Barenblatt solution tends to vary more slowly than Gaussian Solution inside its support, and it tends to be steeper near the interface of the support. We have to notice that in this case the presence of nonlinear diffusion complicated the scheme. The system of linear equations in the previous case is replaced by nonlinear one. It has to be solved at each time iteration, which increases computational cost. We are going to use Newton method. It requires calculation of a Jacobian for a s − 2-dimensional vector not only at each time iteration but also at each iteration inside Newton method until reaching the desired accuracy. To construct the scheme for (6) we need at first space discretization of the diffusion terms. We are going to use conservative, second order, centered discretizations. The approximation for diffusion is given at node x i by Using a forward difference in time and a central difference in space for the equation (6),

θ− scheme
Introducing vectors W ∈ R s such that θ − scheme takes the form, The above results in a nonlinear system of simultaneous equations, As mentioned before the presence of non linear stiff terms requires at each time step the solution of a system of nonlinear equations. Let us define for every vector w ∈ R s−2 a function where, f (ρ n ) := ρ n + (1 − θ)kW n .

Simulation
In the numerical point of view, solving non-linear PDE is based on iterative process. There are several methods can apply to solve non-linear PDE, such as, bisection method, line search technique, gradient search technique, Newton-Raphshon technique etc. For the simulation we choose Newton-Raphshon method. In Newton method initial guessing is an important task to have faster convergence. To understand this, let's study this case numerically for simple linear and non-linear functions that why it is important to consider for the non-linear PDE because a time dependent PDE model basically If we know time depended PDE's exact solution, it is easy to guess the initial condition should be same at t = 0 and the function ρ(x, t) becomes ρ(x) if it provides a finite result and it will be conversed very fast. Now, to understand steady-sate PDE function, first consider a linear function: The root of this function is 1/10. Figure 2 represents the importance of initial guess for the Newton-Rapshon method. For this function, we see, if we choose initial guess except −4 < x 0 < 4 with respect to tolerance = 1e − 15, it takes one iteration. On the other hand if we choose initial guess in between −4 to 4 i.e, −4 < x 0 < 4, say, x 0 = −3, then it will take 0 iteration.
The roots of this function are x ≈ 0 and x ≈ 1.40441. table 1 represents the 6 numerical results and the number of iteration for the convergence to the root for the function, f (x) = x 2 − 2sin(x). Figure 3 represents the plot of the table 1, also we see, for this non-linear function, number of iteration iterations are  But, this is not helpful as always, specially when smoothness comes forward to deal with it. When f (x) is non-linear, an appropriate initial guess in "Newton's"-methods if Jacobian provides a stable eigenvalues of the f (x), theoretically it is obvious that, in the time depended solution will also be stable. So, from this knowledge, our concern is to find appropriate initial guess for PME (1) and thus we took initial guess the value which is at time t 0 = 0.01 for the true Barenblatt solution.
After solving the PME using Newton method, we are interested now to check the error in each grid points for the exponent γ > 1, to get better accuracy first suppose, we know the true Barenblatt solution (B n ) at time t = t n . Let E n denote the error in the calculation with time and space, that is in each grid point, as computed using the exact solution.
Let's suppose that E n is a scalar as same as chapter 2 that we did for linear diffusion, typically some norm of the error over the grid points, i.e., E n = ρ numerical − ρ exact where ρ numerical is the numerical solution vector and ρ exact is the true solution evaluated at each grid point by using time steps and space steps. So, At each grid point there will be some error which is important to check. So, first we have checked error on the basis of L 1 , L 2 and L ∞ norm for the exponent and usually it is called the Boussinesq's equation and observed the qualitative behavior of error for the exponent γ = 2 on each grid point with grid refinement to have better result. Now, presenting some numerical results to show the effectiveness of the θ − method. To do that, we begin our simulation for the Barenblatt solution of the PME (6), where the initial condition is taken as the Barenblatt solution at t = 1, and the boundary condition is u x (L, t) = 0, x ∈ [0, 25], for, t > 1. For this, we have divide the computation domain into N = 251 uniform cells, and plot in Figure 4, the density profile for different exponents, γ = 2, 3, 4, 5, 6 using space step ∆x = 0.1 for both schemes. All results are compared with the Barenblatt solution. we see that both implicit schemes are able to capture the location of the front correctly without noticeable oscillations near the interface, but need deep insight to check the accuracy for the smooth part of those solutions and decide which one is better for the approximation.
For this region, we choose the parameter of PME for γ = 2, because, since it's not locally integrable but if we can measure error by the L ∞ sense then it can also possible to have better approximation for the larger adiabatic exponent. Started from big space and then refined mesh size ∆x, from 0.25 and end up refining at ∆x ≈ 0.03. The accuracy table 2 and 3 are given for the error between numerical and Barenblatt solution of PME (1) with exponent parameter γ = 2 for the Backward Euler scheme and Crank-Nicolson scheme, where considered space, x ∈ [0, 15] and time, t ∈ [0, 5]. It shows that the implicit-explicit Crank scheme is better than Backward Eular scheme. And comparison results are plotted in Figure 4. Table 2. Estimated error for the PME, numerical simulation has been done by using Backward Euler scheme from initial time t 0 = 0.2 to the final time T = 5, step sizes ∆t = ∆x, C = 1 for the Barenblatt   Using, L 2 − norm now, we are going to measure error to know at final time response and see how our simulation is differ from exact solution at each final time. But, we also want the best fit to a specified numerical response in an L ∞ sense rather than the L 2 sense. This minimizes the maximum error of the approximation. If the L ∞ − norm of the error is small, then we are guaranteed that the approximation to our desired solution response will lie within the illustrated bounds.
As mention before, If the L ∞ − norm of the error is small, then we are guaranteed that the approximation to our desired solution response will lie within the illustrated bounds. Figure 6 represents the individual error for both schemes, where we can see, both schemes satisfy this condition, however, Figure 5 represents that Crank-Nicolson scheme even more accurate than Backward Euler implicit scheme for PME. Now pay more attention to the movement of the numerical interface, to check whether the Crank-Nicolson schemes has the ability to capture the true interface accurately. The position of the numerical interface is detected as follows: we scan the numerical solution in the space x from 12.5 to 14, where x ∈ [0, 15], and find the right endpoint of this cell as the numerical interface. We plot in Figure 7