1 Introduction and Background
One of the typical nonlinear degenerate parabolic equations is the porous medium equation (PME):
where is a nonnegative scalar function of space () and the time , and is a constant larger than 1. It has been applied in many physical and biological models, such as an isentropic gas flow through a porous medium, the viscous gravity currents, nonlinear heat transfer and image processing [18], etc.
It is well known that the PME is degenerate at points where . In turn, the PME has many special features: the finite speed of propagation, the free boundary, a possible waiting time phenomenon [5, 18]. Various numerical methods have been studied for the PME, such as finite difference approach [8], tracking algorithm method [3], a local discontinuous Galerkin finite element method [24], Variational Particle Scheme (VPS) [23] and an adaptive moving mesh finite element method [13]. Many theoretical analyses have been derived in the existing literature [1, 12, 14, 16, 17, 18], etc.
Relevant detailed descriptions can be found in a recent paper [5], in which the numerical methods for the PME were constructed by an Energetic Variational Approach (EnVarA) to naturally keep the physical laws, such as the conservation of mass, energy dissipation and force balance. Meanwhile, based on different dissipative energy laws, two different numerical schemes have been studied. In more details, based on the total energy form and , a fully discrete nonlinear scheme and a linear numerical scheme could be appropriately designed for the trajectory equation, respectively. It has also been proved that the former one is uniquely solvable on an admissible convex set, and both schemes preserve the corresponding discrete dissipation law. Numerical experiments have demonstrated that both schemes have yielded a good approximation for the solution without oscillation and the free boundary. The notable advantage is that the waiting time problem could be naturally treated, which has been a wellknown difficult issue for all the existing methods. In addition, under certain smoothness assumption, the second order convergence in space and the first order convergence in time have been reported for both schemes in [5]. The aim of the paper is to provide a rigorous proof of the optimal rate convergence analysis for the nonlinear scheme. On the other hand, the highly nonlinear nature of the trajectory equation makes the convergence analysis every challenging. To overcome these difficulties, we use a higher order expansion technique to ensure a higher order consistency estimate, which is needed to obtain a discrete bound of the numerical solution. Similar ideas have been reported in earlier literature for incompressible fluid equations [6, 7, 15, 21], while the analysis presented in this work turns out to be more complicated, due to the lack of a linear diffusion term in the trajectory equation of the PME. In addition, we have to carry out two step estimates to recover the nonlinear analysis:

Step 1 A rough estimate for the discrete derivative of numerical solution, namely () at time , to control the nonlinear term;

Step 2 A refined estimate for the numerical error function to obtain an optimal convergence order.
Different from a standard error estimate, the rough estimate controls the nonlinear term, which is an effective approach to handle the highly nonlinear term.
This paper is organized as follows. The trajectory equation of the PME and the numerical scheme are outlined in Sec. 2.1 and Sec. 2.2, respectively. Subsequently, the proof of optimal rate convergence analysis is provided in Sec. 3. Finally we present a simple numerical example to demonstrate the convergence rate of the numerical scheme in Sec. 4.
2 Trajectory equation and the numerical scheme
In this section, we review the trajectory equation and the corresponding numerical scheme.
2.1 Trajectory equation of the PME
In this part, the onedimensional trajectory equation will be reviewed, derived by an Energetic Variational Approach [5]. We solve the following initialboundary problem:
(2.1)  
(2.2)  
(2.3)  
(2.4) 
where is a bounded domain, is a nonnegative function, is the time, is the particle position and is the velocity of particle.
The following lemma is available.
Lemma 2.1.
Proof: We first prove the energy dissipation law (2.5) if is the solution of (2.1)(2.4). Multiplying by and integrating on both sides of (2.1), we get
Using integration by parts, in combination with (2.2), we have
(2.6) 
In addition, (2.1) is the conservation law. In the Lagrangian coordinate, its solution can be expressed by:
(2.7) 
where is the positive initial data and is the deformation gradient in one dimension.
Based on an Energetic Variational Approach, we can obtain the trajectory equation.

Energy Dissipation Law.
The total energy of the PME is
(2.8) 
Least Action Principle step.
With (2.7), the action functional in Lagrangian coordinate becomes
where is a given terminal time and is the free energy depending on . Thus for any test function and , taking the variational of with respect to , we have
Then the conservation force turns out to be
in the Eulerian coordinate, and
in the Lagrangian coordinate.

Maximal Dissipation Principle step.
Define the entropy production . Taking the variational of with respect to the velocity and , we obtain the dissipation force
in the Eulerian coordinate and
in the Lagrangian coordinate.

Force balance step. Based on the Newton’s force balance law, we get
(2.10) 
in the Lagrangian coordinate, and the Darcy’s Law in the Eulerian coordinate
It is noticed that, there is an assumption that the value of initial state is positive in to make welldefined in (2.5). More details can be found in [5].
Since then, we should first settle the initial and boundary conditions for (2.10). From (2.1) and (2.4), we have , for . This means that the particles lying on boundary will stay there forever, so a Dirichlet boundary condition should be subject to as , for . As a result, the trajectory problem becomes
2.2 Numerical scheme of the trajectory equation
Let , where , is the final time, and the grid points are given by , . Let be the left point of and be the spatial step, . Denote by , where takes on integer and half integer values. Let and be the spaces of functions whose domains are and , respectively. In component form, these functions are identified via , , for , and , , for .
The difference operator , , and can be defined as:
(2.14)  
(2.15)  
(2.16)  
(2.17)  
(2.18) 
respectively.
Let with its boundary set . Then is a closed convex set. Its physical meaning indicates that particles are arranged in the order without twisting or exchanging in .
A few more notations have to be introduced. Let , and , . We define the inner product on space and respectively as:
(2.19)  
(2.20) 
The following summation by parts formula is available:
(2.21) 
The inverse inequality is given by:
(2.22) 
The fully discrete scheme is formulated as follows: Given the positive initial state and the particle position , find such that
(2.23) 
with and , .
3 Convergence analysis
In this section, the second order spatial convergence and the first order temporal convergence will be theoretically justified for the numerical scheme (2.23). We first introduce a higher order approximate expansion of the exact solution, since a consistency estimate (second order in space and first order in time) is not able to control the discrete norm of the numerical solution. Also see the related works in the earlier literature [2, 6, 7, 9, 10, 15, 19, 20, 21, 22], etc.
Lemma 3.1.
Assume a higher order approximate solution of the exact solution :
(3.1) 
where , , . Then there exists a small , such that , , i.e., , where and are the time step and the spatial mesh sizes, respectively.
Proof: Because of a pointwise condition for the exact solution, , i.e., , such that . For small , such that , and , for . As a consequence, for , we have
(3.2) 
which in turn implies that .
Theorem 3.2.
Proof: A careful Taylor expansion of the exact solution in both time and space, in terms of the numerical scheme (2.11), gives that
(3.4) 
where , , , , , with only dependent on the exact solution.
To perform a higher order consistency analysis for an approximate solution of the exact solution, we have to construct the approximation as in (3.1).
The term is given by the following linear equation:
(3.5) 
The term is given by the following linear equation:
(3.6) 
The term is given by the following linear equation:
(3.7) 
Since , , are dependent only on and , we have the following estimate:
(3.8) 
With these expansion terms, the constructed approximation satisfies the numerical scheme with a higher order truncation error:
(3.9) 
where , are dependent only on , , , , and the derivatives of , , .
Then we define , , . In other words, instead of a direct comparison between the numerical solution and exact PDE solution, we evaluate the numerical error between the numerical solution and the constructed solution . The higher order truncation error enables us to obtain a required of the numerical solution, which is necessary in the nonlinear convergence analysis.
Note that the discrete norm at time step . We make the following apriori assumption at time step :
(3.10) 
In turn, the following estimates become available, by making use of inverse inequalities:
(3.11)  
(3.12)  
(3.13)  
(3.14)  
(3.15) 
For , i.e., , such that , then , , if .
In turn, subtracting (3) from the numerical scheme (2.23) yields
(3.16) 
in which the form of the left term comes from the following identity:
Based on the preliminary results, taking a discrete inner product with (3.16) by gives
Comments
There are no comments yet.