A discontinuous Galerkin method for nonlinear parabolic equations and gradient flow problems with interaction potentials

We consider a class of time dependent second order partial differential equations governed by a decaying entropy. The solution usually corresponds to a density distribution, hence positivity (non-negativity) is expected. This class of problems covers important cases such as Fokker-Planck type equations and aggregation models, which have been studied intensively in the past decades. In this paper, we design a high order discontinuous Galerkin method for such problems. If the interaction potential is not involved, or the interaction is defined by a smooth kernel, our semi-discrete scheme admits an entropy inequality on the discrete level. Furthermore, by applying the positivity-preserving limiter, our fully discretized scheme produces non-negative solutions for all cases under a time step constraint. Our method also applies to two dimensional problems on Cartesian meshes. Numerical examples are given to confirm the high order accuracy for smooth test cases and to demonstrate the effectiveness for preserving long time asymptotics.


Introduction
In this paper, we propose a discontinuous Galerkin method for solving the initial value problem, ρ(x, 0) = ρ 0 (x). (1.1) Here ρ = ρ(x, t) ≥ 0 is a time-dependent density. H is strictly increasing and W (x) = W (−x) is symmetric. We assume f (ρ) ≥ 0 and it attains zero if and only if ρ = 0. f either increases strictly with respect to ρ, or it satisfies the property f (ρ) ρ ≤ C for some fixed constant C.

1.2)
Many classic problems can be included in this setting, such as the heat equation, the porous media equation, the Fokker-Planck equation and so on.
In the second case, several authors take the interaction term W * ρ into account, while imposing f (ρ) = ρ. Hence the equation takes the form of a continuity equation, and the velocity field is determined by the gradient of a potential function.
Here, H is a density of internal energy, V is a confinement potential, and W is an interaction potential. It is used to model, for example, interacting gases [16,43], granular flows [4,5] and aggregation behaviors in biology [40,42]. This equation is also related with the gradient flow for the Wasserstein metric on the space of probability measures [2]. Here u is defined as that in (1.3) and I is referred to as the entropy dissipation.
Indeed, (1.5) has provided much insight into the problem and has helped people to study the dynamics of (1.2) and (1.3), see, for example [15,16,43]. Hence it is desirable to develop numerical schemes mimicking a similar entropy-entropy dissipation relationship in the discrete sense.
Another challenge for developing the numerical schemes is to ensure the non-negativity of the numerical density without violating the mass conservation. It is not only for the preservation of physical meanings, but also for the well-posedness of the initial value problem.
For example, in (1.3), the entropy may not necessarily decay if ρ admits negative values.
Numerical schemes addressing both of these concerns have been studied intensively very recently. In [6], the authors designed a second order finite volume scheme for (1.2). Later on, a direct discontinuous Galerkin method has been proposed by Liu and Wang in [37]. Their scheme achieves high order accuracy but the preservation of non-negativity only holds for limited cases. As for (1.3), a variety of numerical methods have been developed, including a mixed finite element method [8], a finite volume method [12], a particle method [14], a method of evolving diffeomorphisms [17] and a blob method [22] (for H = 0, V = 0).
In this paper, we design a high order discontinuous Galerkin (DG) method for (1.1), which covers both (1.2) and (1.3). The DG method is a class of finite element methods using spaces of discontinuous piecewise polynomials and is especially suitable for solving hyperbolic conservation laws. Coupled with strong stability preserving Runge-Kutta (SSP-RK) time discretization and suitable limiters (the so-called RKDG method, developed by Cockburn et al. in [26,25,24,23,27]), the method captures shocks effectively and achieves high order accuracy in smooth regions [29]. The method has also been generalized for problems involving diffusion and higher order derivatives, for example, the local DG method [3,28], the ultra-weak DG method [21] and the direct DG method [38].
Our idea is to formally treat (1.1) as a classical conservation law and apply the techniques there to overcome the challenges. The main ingredients for our schemes are: 1. Legendre-Gauss-Lobatto quadrature rule for numerical integration, 2. positivity-preserving limiter with SSP-RK time discretization.
The quadrature is used to stabilize the semi-discrete scheme. Such approach has already been studied in different contexts, such as in spectral collocation methods [34] and in nodal discontinuous Galerkin methods [35]. Recently, based on the methodology established in [11,31,32], Chen and Shu proposed a unified framework for designing entropy stable DG scheme for hyperbolic conservation laws using suitable quadrature rules [20]. Their approach is also related with the summation-by-part technique in finite difference methods. The positivitypreserving limiter with provable high order accuracy is firstly designed by Zhang and Shu in [45] to numerically ensure the maximum principle of scalar hyperbolic conservation laws.
Then the methodology has been generalized for developing the bound-preserving schemes for various systems. We refer to [46] and the references therein for more details. Our approach of implementing the positivity-preserving limiter for parabolic problems is mainly inspired by the recent work of Zhang on compressible Navier-Stokes equation [44], in which the author considers the conservative form of the problem and introduces the diffusion flux to handle the second order derivatives.
Based on these techniques, we propose a discontinuous Galerkin scheme for (1.1) such that 1. the semi-discrete scheme satisfies an entropy inequality for smooth W , 2. the fully discretized scheme is positivity-preserving.
Our method also has other desired properties. It achieves high order accuracy, conserves the total mass and preserves numerical steady states. Special care is needed for the case of non-smooth interaction potential, due to the fact that one should adopt exact integration to calculate the convolution, see [12], as the Gauss-Lobatto quadrature is no longer accurate.
The remaining part of the paper is organized as follows. In Section 2, we design the numerical method for one dimensional problems. We firstly introduce the notations and briefly discuss our motivation on deriving the scheme. Then it follows with the semi-discrete scheme and the discrete entropy inequality. The next part is on the time discretization and the positivity-preserving property of the fully discretized scheme. Finally we outline the matrix formulation and the algorithm flowchart. Section 3 is organized similarly for two dimensional problems on Cartesian meshes, while the implementation details are omitted.
Then in Section 4 and Section 5, we present numerical examples for one dimensional and two dimensional problems respectively. Finally conclusions are drawn in Section 6.
2 Numerical method: one dimensional case

Notations and motivations
For now, we focus on the one dimensional case of (1.1) In general, the problem can be either defined on a connected compact domain with proper boundary conditions, or it can involve the whole real line with solutions vanishing at the infinity. In our numerical scheme, we will always choose Ω to be a connected interval. For simplicity, the periodic or compactly supported boundary conditions are applied. But we remark that our approach can be extended to more general types of boundary conditions, for example, with zero-flux boundaries.
) and I = ∪ N i=1 I i be a regular partition of the domain Ω. Denote and h = max i h i . We will seek a numerical solution in the discontinuous piecewise polynomial space, Here P k (I i ) is the space of k-th order polynomials on I i . Note that the functions in V h can be double-valued at the cell interfaces. Hence the notations v + h and v − h are introduced for the right limit and the left limit of v h . For a function s = s(ρ) or s = s(ρ, x), we denote by s h = s(ρ h ) or s h = s(ρ h , x) respectively. Furthermore, the notation s + h stands for s(ρ + h ) or s(ρ + h , x), and s − h stands for s(ρ − h ) or s(ρ − h , x). To define the DG method, we split the original problem into the following system, By applying the DG approximation, we obtain the following scheme as a rudiment. We seek dxdy, while f u and ξ are numerical fluxes.
By setting ϕ h (x) = 1 I i (x), one will get the following evolution equation for the cell This is very similar to that in hyperbolic conservation laws. In particular, for k = 0 and u ≡ 1, by using the so called monotone flux, (2.3) will become a monotone scheme, which satisfies many desirable properties.
In order to achieve an entropy-entropy dissipation relationship as that for the exact solution, one may hope to set ϕ h = ξ h in (2.2a) and ψ h = u h f h in (2.2b). Unluckily, neither of them falls into the test function space. A natural attempt is to use the usual L 2 projection to enforce this property, as that in [37]. However, the projection will change the values at the cell interfaces, and the desired form in (2.3) will be violated. This will cause trouble when one seeks to preserve the non-negativity of the solution. Inspired by the recent work of Chen and Shu in [20], we introduce a suitable quadrature to overcome this difficulty. Moreover, the Gauss-Lobatto quadrature is used to preserve the values at the cell ends.
Let us denote by {x r i } k+1 r=1 the k+1 Gauss-Lobatto quadrature points on I i and {w r } k+1 r=1 the k +1 Gauss-Lobatto quadrature weights on [−1, 1]. In particular On each cell I i , the operator I returns the k-th order polynomial interpolating at {x r i } k+1 r=1 .
We will use the notation for the Gauss-Lobatto quadrature. As a convention, ∼ Ω stands for i ∼ I i . We will finally come to the fully discretized scheme, for which we denote by τ the time step length and λ i = τ h i .

Semi-discrete scheme and entropy inequality
Our scheme is to replace the integrals in (2.2) by the Gauss-Lobatto quadrature. In other words, we seek ρ h , u h ∈ V h , such that for any test functions Here, when the interaction potential W is smooth, we set While for non-smooth W , the quadrature may not achieve sufficient accuracy of the convolution. Hence the exact integration is applied The numerical fluxes are chosen in the following way, ξ is the central flux. When f (ρ) = ρ, one can choose g(ρ) = ρ and f u is the Lax-Friedrich flux.
Remark 2.1. Although we formally require that ϕ h and ψ h are taken from the finite element space, our scheme (2.4) actually can not distinguish a function from its interpolation polynomial at the Gauss-Lobatto points. Hence, (2.4) will also hold for "test functions" outside V h . This facilitates our proof of the discrete entropy inequality. This fact can also be justified through the matrix formulation, which is presented in Section 2.4.
This semi-discrete scheme satisfies the following entropy inequality.
Theorem 2.1. For smooth interaction kernel W , assume that the semi-discrete scheme (2.4) has a solution, then it satisfies a similar entropy-entropy dissipation relationship, as that in is the discrete entropy andĨ is the associated discrete entropy dissipation. Proof. Using the symmetry of W , we have d dt Hence, . Note is a polynomial of degree 2k − 1 and the Gauss-Lobatto quadrature with k + 1 points is exact. Hence we can replace the quadrature with the exact integral and integrate by parts. Then we obtain By using the same trick, we change the exact integral back to the Gauss-Lobatto quadrature, and apply the scheme (2.4b) to obtain ≥ 0, which completes the proof of the first claim.
It is easy to see that the entropy inequality will hold as long as the coefficients α i+ 1 2 are non-negative. Let us assume now that α i+ 1 2 > 0 for all i and ρ h is a non-negative stationary state of the semi-discrete scheme (2.4), namely Then from (2.6) we deduce that both terms in the right-hand side must vanish, that is Therefore each term in the summations will be zero. On the one hand, if α i+ 1 2 > 0, then Here ξ ± h = ξ is guaranteed by the continuity of ξ h (implied by that of ρ h ). Therefore, Iξ h is constant on each i ∈ Λ. Due to the fact that ξ h is continuous globally, all these constants must be the same and the piecewise polynomial interpolation of ξ h is constant on J.

Time discretization and preservation of positivity
The semi-discrete scheme itself does not guarantee the positivity of the numerical solution. If no special treatment is applied, one may produce nonsense density with negative values and the problem can become illposed. Hence we adopt the methodology developed by Zhang and Shu in [45], which enforces the positivity of the solution without violating the mass conservation. Their idea is to incorporate a positivity-preserving limiter into the strong stability preserving Runge-Kutta (SSP-RK) time discretization. Under certain time step constraints, each Euler forward step preserves the positivity of the cell average (referred to as the weak positivity in the literature). Then one can scale the solution, without affecting spatial accuracy, to ensure the point-wise non-negativity. SSP-RK time discretization will preserve the non-negativity of the solution in the Euler forward steps.

First order Euler forward in time
Let us firstly consider the Euler forward time stepping. We use the superscript "pre" for the solution obtained by Euler forward method before applying the positivity-preserving limiter.
The time discretization of (2.4a) becomes Here, in the constraint of λ i , we consider 0 0 := +∞ by convention.
Proof. We drop all subscripts h in this proof. Take ϕ = 1 in (2.7), we havē Note that ρ n is a polynomial of degree k, the Gauss-Lobatto quadrature is exact for evaluating the cell averageρ n i . More specifically, we havē The superscripts n will also be omitted for simplicity in the rest. Hencē The first term is automatically non-negative, since the weights w r ≥ 0 and the nodal values The positivity of the last two terms is guaranteed by our choice of α and g. One } to ensure the second and the third term to be being 0, the corresponding term is also 0 and there is nothing to impose. Hence we introduce the notation 0 0 := +∞.) Thereforeρ n+1,pre i ≥ 0 under the prescribed time step constraint.
1. According to the definition of α and g, ( will always be nonnegative. 2. Although the original equation can be parabolic, we have incorporated the second order derivative into u, such that one can formally treat it as a hyperbolic problem. This technique is introduced by Zhang for the compressible Navier-Stokes equation [44]. In non-negative at the previous time step (at the Gauss-Lobatto quadrature points), as long as the time step is smaller than a threshold, the cell average at next time step will remain non-negative. In order to close the loop, one would need to ensure the nodal values at the quadrature points of the next time step are also non-negative. This indeed can be achieved by applying a scaling limiter, which luckily does not affect the spatial accuracy. We refer to [44] for more details.
where ρ(x, t n+1 ) is the exact solution at time t n+1 and C k is a constant depending only on the polynomial degree k.
Remark 2.3. Our scheme only uses the nodal values at the Gauss-Lobatto quadrature points, hence we only need to ensure the non-negativity at these nodes. One can also squash the solution polynomials so that the solution is non-negative everywhere on the domain. The proof will still go through.
Theorem 2.2. With the scaling limiter in Lemma 2.2, the Euler forward time discretization of the semi-discrete scheme is positivity-preserving, provided the time step restriction specified in Lemma 2.1 is satisfied.

High order time discretization
The SSP-RK method will be used for time discretization. We refer readers to [33] for more details. Since the time step scales like τ = Ch 2 , the Euler forward method will be sufficient for piecewise linear elements to achieve overall second order accuracy. For k = 2, 3, we will use the second order SSP-RK scheme For k = 4, 5, the third order SSP-RK scheme is used The positivity-preserving limiter should be applied immediately after each Euler forward stage. As one can see, the SSP-RK schemes (2.8) and (2.9) can be rewritten as convex combinations of the Euler forward steps. Since each Euler forward step preserves the positivity, the numerical density at the next time level will remain non-negative. We also mention several other properties of the fully discretized scheme, whose proofs are omitted. Such properties also hold for two dimensional cases.

Matrix formulation and implementation
At the end of this section, we would like to introduce the matrix formulation of our numerical scheme and outline the flowchart of the algorithm.

Matrix formulation
The derivation of the matrix formulation is similar to that in Section 3.1 of [20]. We refer to that paper for more details.
On each cell, the unknown function can be represented as Here ζ i is the mapping from I i to [−1, 1]. To determine ρ, it suffices to identify the coefficients The matrix formulation can be written as follows.
] T and f u * i is defined similarly for f u. We remind the readers that one should replace I i by

Algorithm flowchart
For simplicity, we only consider the Euler forward time stepping. The algorithm with SSP-RK time discretization can be implemented by repeating the following flowchart in each stage.
• If {ρ i } is a set of non-negative numbers. Apply the positivity preserving limiter to obtain {( ρ i ) n+1 } and enter the next time level.
• Otherwise halve the time step τ and restart from 1.  N ). The idea is that, for each fixed i, the convolution can be evaluated by, If the convolution kernel is periodic, then (2.11) can be formulated as the multiplication of an (N × (k + 1)) × (N × (k + 1)) block circulant matrix and a vector. The FFT acceleration for such problems is standard.
Although W is not periodic most of the time, ρ is usually a (numerically) compactly supported function. One only needs to evaluate ξ precisely on the same interval. Hence we can simply extend the problem to a larger domain to adopt the previous procedure. 3. In our numerical tests, both for one dimensional and two dimensional cases, we will use a sufficiently small time step to avoid the cell average attaining negative values.
Also, g(ρ) = ρ will be used to define the numerical flux, unless otherwise stated.

Numerical method: two dimensional case
In this section, we apply our method to solve two dimensional problems on Cartesian meshes.
Here Ω = I × J is a rectangular domain and the periodic boundary conditions are applied.
The mesh size is denoted by . The finite element spaces are defined as Here Q k (I i × J j ) is the tensor product space of P k (I i ) and P k (J j ).
The semi-discrete DG scheme is formulated as follows. One needs to find ρ h ∈ V h and , y)dy Here, when the interaction potential W is smooth, we set While for non-smooth W , the quadrature may not achieve sufficient accuracy. Hence the exact integration is applied The numerical fluxes are chosen in the following way, For smooth W , one can obtain an entropy inequality as we have done for one dimensional problems.
Theorem 3.1. For smooth interaction kernel W , assume that the semi-discrete scheme defined by (3.1) and (3.2) has a solution, then it satisfies the following entropy inequality.
is the discrete entropy andĨ Proof. We will focus on the entropy-entropy dissipation relationship and the proof of the second part of the theorem is omitted. Using the symmetry of W , we have d dt Hence, , y)dy is a polynomial of degree 2k − 1 with respect to x. Hence the Gauss-Lobatto quadrature with k + 1 nodes is exact. We replace the quadrature with the exact integral, integrate by parts and then change back to the quadrature. The same argument also applies to the second integral. One can then obtain, , y)dy Use the scheme (3.2) one can get , y)dy By our choices of g, the strict monotonicity of H and the fact that V and W are singlevalued, the last term is non-positive, which gives (3.3).

Time discretization and preservation of positivity
It suffices to ensure the positivity-preserving property of the Euler forward scheme. The high order case is automatically taken care of by SSP-RK time discretization.
The first step is to show that, provided the solution at the current time level is nonnegative, the cell average at next time level will also be non-negative, if a specific time step restriction is satisfied.
Here, in the constraint of λ x i and λ x j , we formally denote by 0 0 := +∞.
Proof. As before, we drop all the subscripts h in this proof. The superscript n will also be omitted for simplicity. Take ϕ = 1 in (3.1), we havē Let λ x i = τ h x i and λ y j = τ h y j , then we havē ) .
The first term is automatically non-negative, since the weights w r , w s ≥ 0 and the nodal values ρ(x r i , y s j ) ≥ 0. The positivity of the last four terms is guaranteed by our choice of α and g. One only needs (3.4) to ensure the second and the third term to be non-negative.
And as before, one can check the convention 0 0 = +∞ does make sense. Henceρ n+1 i,j ≥ 0 under the prescribed time step constraint.
Then, as we have done in the one dimensional case, a scaling limiter is applied to sure the numerical polynomial solution takes non-negative values at the quadrature points. Hence the assumption in Lemma 3.1 is met and the fully discretized scheme is positivity-preserving.

One dimensional numerical tests 4.1 Accuracy tests
In this part, we examine the accuracy of the numerical schemes with P 1 , P 2 , P 3 and P 4 elements. The error is measured in the discrete norms.
The problem has an exact solution u(x, t) = 1 + sin(x + t). In this test, f (ρ) = ρ, H (ρ) = 0, V (x) = x and W (x) = 0. To be consistent at the boundaries, one needs to manually impose ξ(π) = π and ξ(−π) = −π. (This will gives u ≡ 1 and the scheme is equivalent to the usual upwinding DG method with a mass lumping treatment.) We compute up to t = 2 and the time step is τ = 0.02h 2 .
Due to our choice of the initial condition, the solution has point vacuum and the numerical solution may become negative in its neighborhood. We perform numerical tests without and with the positivity-preserving limiter and the results are listed in Table 4.1 and Table 4.2 respectively. As one can see, without the limiter, the convergence rate is optimal. The rate degenerates a little bit for the P 4 scheme when one applies the limiter.
ρ(x, 0) = 2 + sin(x), with periodic boundary conditions. The exact solution to the problem is ρ(x, t) = 2 + e −t sin(x). The decomposition of the equation into the desired form is not unique. Let us consider two test cases, Note that for both of the cases, the schemes are nonlinear, although the original problem is linear. We use the time step τ = 0.01h 2 to compute to t = 2. Error tables are given in Table 4.3 and Table 4.4 respectively. According to our numerical results, we see that different choices of decomposition lead to negligible difference. For both of the tests, P 2 and P 4 schemes are of the optimal rate of convergence, but the order for P 1 and P 3 schemes seems to be reduced.
for problems with interaction potentials. The order of accuracy seems to be optimal for odd k, while the order degenerates for even k. See Table 4.5 and Table 4

Fokker-Planck type equations
which is used to model the flow of a gas through a porous interface. The equation fits our model with H(ρ) = 1 m−1 ρ m and V (x) = x 2 2 . The property of the equation is studied by Carrillo and Toscani in [19] using an entropy approach. They have proved that the equation converges to a unique steady state given by a Barenblatt-Pattle type formula, Here the constant C is determined by ensuring the mass conservation. Furthermore, the relative entropy E(t|∞) = E(ρ(t)) − E(ρ ∞ ) decays exponentially, E(t|∞) ≤ E(0|∞)e −2t and the rate −2 is sharp. We particularly choose m = 2 in our numerical test.
ρ(x, 0) = max{1 − |x|, 0}, with periodic boundary conditions. The stationary solution is We  zero. We then evaluate the relative entropy using that of the numerical steady state as a reference. But we can only plot up to a certain time, before the relative entropy becomes slightly negative, since the entropy decay of the semi-discrete scheme may not be preserved after applying the time discretization and the limiter. We stop the plotting after negative relative entropy appears, not only because it can not be depicted in the logarithm scale, but also because the unresolved tail should be regarded as a discretization error. If we choose N = 320 while keeping τ = 0.005h 2 , the decay will continue further.
For the symmetric initial condition, our numerical tests indicate the decay rates is around e −6.7t . Indeed, symmetric initial data converge faster to equilibrium than the sharp rate since they preserve the invariance of the center mass, see [13] for more details. We then test the problem under the same settings, except for the initial condition shifted to the right ρ(x, 0) = max{1 − |x − 1 2 |, 0} and the final time set to t = 10. The corresponding plot of E(t|∞) is given in Figure 4.3 with the exponential decay rate −2, which coincides with the    result in [19]. Similar numerical test can be found in [6].
Here κ = 1 corresponds to boson gases and κ = −1 relates to fermion gases. The long time asymptotics of the one dimensional model has been studied in [18]   For both of the test cases, we use numerical steady states as references to calculate the relative entropy. The result for the boson case is given in Figure 4.4a. As one can see, the decay rate is around −2.6. While for the fermion case, which is exhibited in Figure 4.4b, the relative entropy decays at a slower rate of −1.44.  mass will evolve a singularity at the origin, which has been confirmed numerically in [6] and [37]. In this test, we repeat the numerical experiment in [6] and [37], setting The initial datum is chosen as We test with both the subcritical case with M = 1 and the supercritical case with M = 10.
The P 4 elements are used in our numerical scheme and we compute on the domain [− 6,6] with N = 120. For M = 1, the time step is chosen as τ = 0.003h 2 and for M = 10, it is τ = 0.0005h 2 . And we use g(ρ) = f (ρ) when defining the numerical flux. According to the numerical results in Figure 4.5, our scheme does capture the asymptotics of the equation.   Here 0 ≤ ρ ≤ 1. ν > 0 and m > 1 are parameters to be specified. The convolution kernel W in this example is nonlocal and smooth. Under this setting, the attraction effect is weak and the solution would either end up with a steady state or spread out in the whole domain with bounded initial data. The compactly supported steady state is of special interest, due to its application for modeling the biological aggregation, such as flocks and swarms. Indeed, such stationary solution can be reached for m > 2 with arbitrary ν. While for 1 < m ≤ 2, the long time behavior of the solution can be sophisticated. We refer to [9] and [10] for details.

Aggregation models
For our numerical test, we focus on the specific setting, We apply periodic boundary conditions and use a mesh with N = 120 for computation. The     When they get closer, the attraction becomes strong. A sudden decay of the relative entropy occurs when the two bumps merge. After that, the contribution of the interaction potential to the total energy becomes small. The equation is again dominated by the diffusion term, and the relative entropy decays exponentially as we have seen before.
We omit the plots for m = 2. And for m = 3, the diffusion is relatively weak and the exponential decay after the steep drop is hard to observe. Hence we only plot the entropy in the normal scale. But still we can see the sharp drop when the bumps merge in Figure 4.8.
The initial stage featured with the weak long-range-interaction is referred as metastability.
To convince the readers that the disconnected profile in Figure 4.10b is indeed the stationary solution, we plot ρ and ξ in Figure 4.11b. As one can observe, ρ∂ x ξ ≈ 0. Hence ∂ t ρ ≈ 0 and ρ will be trapped in this steady state. Therefore, the observation in Figure   4.10 confirms our previous claim, that different initial density distributions may end up with steady states with distinct connectivity. Let us remark that such phenomenon has also been explored numerically in [12].
Example 5.2 (dumbbell model for polymers). The dumbbell model is widely used to describe the rheological behavior of dilute polymer solutions. In this model, the polymer molecular is treated as a dumbbell made of two beads jointed by a spring. We will consider the simplest case, in which the flow is homogeneous and the scaling constant is set to 1.
Then the configuration probability density is governed by the Fokker-Plank equation, ∂ t ρ(x, t) = ∇ · (ρ∇(U − 1 2 xKx)) + ∆ρ.   Here x = (x, y) corresponds to the direction vector of the molecule, while U is the spring potential and the 2 × 2 matrix K is the velocity gradient of the background flow. For the incompressible flow, Tr(K) = 0. In our numerical test, we consider the finitely extensible nonlinear elastic (FENE) model. The potential U is given by It is close to the Hookean potential when |x| r, while the distance between the two beads are restricted within r. Rigorously, one should consider the equation on the ball {x : |x| ≤ r}, and the singularity near the boundary will cause challenges both analytically and numerically, see [30,36,41,39] and the references therein. While in our numerical test, we only consider a simpler case, that the solution seems to be supported within the ball and it hardly reaches the boundaries. More specifically, we solve (5.   Patlak-Keller-Segel system is a mathematical model to describe the motion of the organism.
The equation can be rewritten in a compact way ∂ t ρ = ∇ · (ρ∇ (log(ρ) + W * ρ)) , W (x, y) = 1 2π log( x 2 + y 2 ) (x, y) ∈ R 2 , t > 0. (5.5) Such system has been studied intensively in the past decades. It has been shown that the behavior of the equation (5.5) is determined by its initial mass (see [7], for example). If the initial value M is smaller than a critical value M c = 8π, then the solution will exist globally.
Otherwise, if M lies beyond M c , the solution will blow up in a finite time, which is referred as chemotactic collapse.

Concluding remarks
In this paper, we develop a high order DG method for solving a class of parabolic equations and gradient flow problems with interaction potentials. Such equations are governed by an entropy-entropy dissipation relationship and are featured with non-negative solutions.
By applying the Gauss-Lobatto quadrature rule, our numerical scheme admits an entropy inequality for problems with smooth interaction kernels. Furthermore, with the SSP-RK time discretization and the positivity-preserving limiter, the fully discretized scheme preserves the non-negativity of the numerical density. It also conserves mass, and preserves numerical steady states for certain problems. We also apply the method to two dimensional problems on Cartesian meshes. Numerical examples are given to demonstrate the performance of the scheme.