# Convergence analysis of a numerical scheme for the porous medium equation by an energetic variational approach

The porous medium equation (PME) is a typical nonlinear degenerate parabolic equation. We have studied numerical methods for PME by an energetic variational approach in [C. Duan et al, J. Comput. Phys., 385 (2019) 13-32], where the trajectory equation can be obtained and two numerical schemes have been developed based on different dissipative energy laws. It is also proved that the nonlinear scheme, based on f log f as the total energy form of the dissipative law, is uniquely solvable on an admissible convex set and preserves the corresponding discrete dissipation law. Moreover, under certain smoothness assumption, we have also obtained the second order convergence in space and the first order convergence in time for the scheme. In this paper, we provide a rigorous proof of the error estimate by a careful higher order asymptotic expansion and two step error estimates. The latter technique contains a rough estimate to control the highly nonlinear term in a discrete W^1,∞ norm, and a refined estimate is applied to derive the optimal error order.

## Authors

• 3 publications
• 29 publications
• 76 publications
• 5 publications
06/22/2020

### A second order accurate numerical scheme for the porous medium equation by an energetic variational approach

The porous medium equation (PME) is a typical nonlinear degenerate parab...
10/31/2020

### A third order BDF energy stable linear scheme for the no-slope-selection thin film model

In this paper we propose and analyze a (temporally) third order accurate...
08/16/2020

### Structure-Preserving Numerical Methods for Nonlinear Fokker–Planck Equations with Nonlocal Interactions by an Energetic Variational Approach

In this work, we develop novel structure-preserving numerical schemes fo...
05/19/2021

### Convergence analysis of the variational operator splitting scheme for a reaction-diffusion system with detailed balance

We present a detailed convergence analysis for an operator splitting sch...
01/26/2022

### On Stability and Convergence of a Three-layer Semi-discrete Scheme for an Abstract Analogue of the Ball Integro-differential Equation

We consider the Cauchy problem for a second-order nonlinear evolution eq...
10/13/2021

### Variational and numerical analysis of a 𝐐-tensor model for smectic-A liquid crystals

We analyse an energy minimisation problem recently proposed for modellin...
02/22/2022

### Convergence Analysis of Structure-Preserving Numerical Methods Based on Slotboom Transformation for the Poisson–Nernst–Planck Equations

The analysis of structure-preserving numerical methods for the Poisson–N...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction and Background

One of the typical nonlinear degenerate parabolic equations is the porous medium equation (PME):

 ∂tf=Δx(fm), x∈Ω⊂Rd, m>1,

where is a non-negative 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 well-known 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 one-dimensional trajectory equation will be reviewed, derived by an Energetic Variational Approach [5]. We solve the following initial-boundary problem:

 (2.1) fv=−∂x(fm), x∈Ω, m>1, (2.2) f(x,0)=f0(x)≥0, x∈Ω, (2.3) ∂xf=0, x∈∂Ω, t>0, (2.4)

where is a bounded domain, is a non-negative function, is the time, is the particle position and is the velocity of particle.

The following lemma is available.

###### Lemma 2.1.

is a positive solution of (2.1)-(2.4) if and only if satisfies the corresponding energy dissipation law:

 ddt∫Ωflnfdx=−∫Ωfmfm−1|v|2dx. (2.5)

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

 ∫Ω(1+lnf)∂tfdx=−∫Ω(1+lnf)∂x(fv)dx.

Using integration by parts, in combination with (2.2), we have

 ddt∫Ωflnfdx=∫Ω∂xff(fv)dx=−∫Ωfmfm−1|v|2dx≤0. (2.6)

Subsequently, we are also able to derive (2.2) from the energy dissipation law (2.5) by EnVarA.

In addition, (2.1) is the conservation law. In the Lagrangian coordinate, its solution can be expressed by:

 f(x(X,t),t)=f0(X)∂x(X,t)∂X, (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

 Etotal:=∫Ωflnfdx. (2.8)
• Least Action Principle step.

With (2.7), the action functional in Lagrangian coordinate becomes

 A(x):=∫T∗0(−H)dt=−∫T∗0∫Ωf0(X)ln(f0(X)∂Xx)dXdt,

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

 ddε∣∣∣ε=0A(x+εy) =∫T∗0∫Ωf0(X)∂Xx⋅∂Xy dXdt =−∫T∗0∫Ω∂xf⋅˜y dxdt.

Then the conservation force turns out to be

 Fcon=δAδx=−∂xf,

in the Eulerian coordinate, and

 Fcon=−∂X(f0(X)∂Xx),

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

 Fdis:=δ12Δδv=fmfm−1v,

in the Eulerian coordinate and

 Fdis:=δ12Δδ(∂tx)=f0(X)m(f0(X)∂Xx)m−1∂tx,

in the Lagrangian coordinate.

• Force balance step.   Based on the Newton’s force balance law, we get

 f0(X)m(f0(X)∂Xx)m−1∂tx=−∂X(f0(X)∂Xx), (2.10)

in the Lagrangian coordinate, and the Darcy’s Law in the Eulerian coordinate

 fmfm−1v=−∂xf.

It is noticed that, there is an assumption that the value of initial state is positive in to make well-defined 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

 f0(X)m(f0(X)∂Xx)m−1∂tx=−∂X(f0(X)∂Xx), X∈Ω, t>0, (2.11) x|∂Ω=X|∂Ω, t>0, (2.12) x(X,0)=X, X∈Ω. (2.13)

Finally, with a substitution of (2.11) into (2.7), we obtain the solution to (2.1)-(2.4).

### 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:

 (Dhl)i−12=(li−li−1)/h, i=1,...,M, (2.14) (dhϕ)i=(ϕi+12−ϕi−12)/h, i=1,...,M−1, (2.15) (˜Dhl)i=(li+1−li−1)/2h, i=1,...,M−1, (2.16) (˜Dhl)i=(4li+1−li+2−3li)/2h, i=0, (2.17) (˜Dhl)i=(li−2−4li−1+3li)/2h, i=M, (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:

 ⟨l,g⟩:=h(12l0g0+M−1∑i=1ligi+12lMgM), (2.19) ⟨ϕ,φ⟩e:=hM−1∑i=0ϕi+12φi+12. (2.20)

The following summation by parts formula is available:

 ⟨l,dhϕ⟩=−⟨Dhl,ϕ⟩e,\ with l0=lM=0, ϕ∈CM, l∈EM. (2.21)

The inverse inequality is given by:

 ∥l∥∞≤Cm∥l∥2h1/2,  ∀l∈EM,with ∥l∥∞:=max0≤i≤M{li}%  and  ∥l∥22:=⟨l,l⟩. (2.22)

The fully discrete scheme is formulated as follows: Given the positive initial state and the particle position , find such that

 f0(Xi)m(f0(X)˜Dhxn)m−1i⋅xn+1i−xniτ=−dh[(f0(X)Dhxn+1)]i, 1≤i≤M−1, (2.23)

with and , .

It is noticed that (2.23) is still a nonlinear system which can be solved by Newton’s iteration method [5]. Then we obtain the numerical solution by

 fni=f0(X)˜Dhxni, 0≤i≤M, (2.24)

which is the discrete scheme of (2.7).

## 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 :

 W:=xe+τw(1)τ+τ2w(2)τ+h2wh, (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 point-wise condition for the exact solution, , i.e., , such that . For small , such that , and , for . As a consequence, for , we have

 DhW≥13ε0>0, (3.2)

which in turn implies that .

###### Theorem 3.2.

Assume that the initial function is positive and bounded, i.e., . Denote as the exact solution to the original PDE (2.11) (with enough regularity) and as the numerical solution to (2.23). The numerical error function is defined at a point-wise level:

 eni=xnei−xnhi, (3.3)

where , , .

Then we have

• satisfies

 ∥en∥2:=⟨en,en⟩≤C(τ+h2).
• satisfies

 ∥˜Dhen∥2≤C(τ+h2).

Moreover, the error between the numerical solution and the exact solution of equaiton (2.1)-(2.4) can be estimated by:

 ∥fnh−fne∥2≤C(τ+h2),

where is a positive constant, is the spatial step, is the time step and .

Proof: A careful Taylor expansion of the exact solution in both time and space, in terms of the numerical scheme (2.11), gives that

 f0(Xi)m(f0(Xi)˜Dhxnei)m−1xn+1ei−xneiτ=−dh(f0(X)Dhxn+1e)i+τl(1)i+τ2l(2)i+τ3l(3)i+h2g(1)i+h4g(2)i, 1≤i≤M−1, with xn+1e0=X0 , xn+1eM=XM, (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:

 f0(X)m(f0(X)∂Xxe)m−1∂tw(1)τ+m−1m(f0(X)∂Xxe)m−2∂txe⋅∂Xw(1)τ=∂X(f0(X)(∂Xxe)2∂Xw(1)τ)−l(1), w(1)τ|∂Ω=0,    w(1)τ(⋅,0)=0. (3.5)

The term is given by the following linear equation:

 f0(X)m(f0(X)∂Xxe)m−1∂tw(2)τ+m−1m(f0(X)∂Xxe)m−2∂txe⋅∂Xw(2)τ +(m−1)(m−2)2m(f0(X)∂Xxe)m−2∂Xxe(∂Xw(1)τ)2⋅∂txe+(m−1)m(f0(X)∂Xxe)m−2∂tw(1)τ⋅∂Xw(1)τ =∂X(f0(X)(∂Xxe)2∂Xw(2)τ)−∂X(f0(X)(∂Xxe)3(∂Xw(1)τ)2)−l(2), w(2)τ|∂Ω=0,    w(2)τ(⋅,0)=0. (3.6)

The term is given by the following linear equation:

 f0(X)m(f0(X)∂Xxe)m−1∂twh+(m−1)∂txem(f0(X)∂Xxe)m−2∂Xwh=∂X(f0(X)(∂Xxe)2∂Xwh)−g(1), wh|∂Ω=0,    wh(⋅,0)=0. (3.7)

Since , , are dependent only on and , we have the following estimate:

 ∥W−xe∥Hm=τ∥w(1)τ∥Hm+τ2∥w(1)τ∥Hm+h2∥wh∥Hm≤C′(τ+h2). (3.8)

With these expansion terms, the constructed approximation satisfies the numerical scheme with a higher order truncation error:

 f0(Xi)m(f0(X)˜DhWn)m−1i⋅Wn+1i−Wniτ=−dh(f0(X)DhWn+1)i+τ3l∗i+h4g∗i,  1≤i≤M−1, with\ \ Wn+10=X0,  Wn+1M=XM,  n=0,1,⋯,N−1, (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 a-priori assumption at time step :

 ∥~en∥2≤(τ114+h72). (3.10)

In turn, the following estimates become available, by making use of inverse inequalities:

 ∥˜Dh~en∥2≤C(τ74+h52), (3.11) ∥˜Dh~en∥∞≤CCm∥˜Dh~en∥2h1/2≤CCm(τ54+h2), if\ \ h=O(τ), (3.12) ∥˜Dhxnh∥∞=∥˜DhWn−˜Dh~en∥∞≤C∗+1:=C∗0, (3.13) with\ \ C∗:=∥˜DhWn∥∞,  if\ \ CCm(τ54+h2)≤1, ∥∥˜Dhxnh−˜Dhxn−1hτ∥∥∞=∥∥˜DhWn−˜DhWn−1τ−˜Dh~en−˜Dh~en−1τ∥∥∞≤~C∗t+1, (3.14) with\ \ ~C∗t:=∥∥˜DhWn−˜DhWn−1τ∥∥∞, if\ \ CCm(τ14+h)≤1. (3.15)

For , i.e., , such that , then , , if .

In turn, subtracting (3) from the numerical scheme (2.23) yields

 f0(Xi)m(f0(X)˜Dhxnh)m−1i⋅~en+1i−~eniτ+f0(Xi)m[f0(Xi)]m−1⋅Wn+1i−Wniτ⋅[(˜DhWn)m−1i−(˜Dhxnh)m−1i] =dh(f0(X)DhWn+1iDhxn+1hDh~en+1)i+τ3l∗i+h4g∗i,  1≤i≤M−1, with\ ~en+10=~en+1M=0, (3.16)

in which the form of the left term comes from the following identity:

 f0(Xi)m(f0(X)˜DhWn)m−1iWn+1i−Wniτ−f0(X)m(f0(Xi)˜Dhxnh)m−1ixn+1hi−xnhiτ =f0(Xi)τm[f0(Xi)]m−1[(˜DhWn)m−1i(Wn+1i−Wni)−(˜Dhxnh)m−1i(xn+1hi−xnhi) +(˜Dhxnh)m−1i(Wn+1i−Wni)−(˜Dhxnh)m−1i(Wn+1i−Wni)] =f0(Xi)m[f0(Xi)]m−1⋅Wn+1i−Wniτ⋅[(˜DhWn)m−1i−(˜Dhxnh)m−1i] +f0(Xi)m(f0(X)˜Dhxnh)m−1i⋅~en+1i−~eniτ.

Based on the preliminary results, taking a discrete inner product with (3.16) by gives

 2⟨αn(~en+1−~en),~en+1⟩−2τ⟨dh(f0(X)DhWn+1Dhxn+1hDh~en+1),~en+1⟩ =−2τ⟨f0(X)m[f0(X)]m−1⋅Wn+1−Wnτ⋅[(˜DhWn)m−1−(˜Dhxnh)m−1],~en+1⟩ +