Search

We develop a new finite difference scheme for the Maxwell-Stefan diffusion system. The scheme is conservative, energy stable and positivity-preserving. These nice properties stem from a variational structure and are proved by reformulating the finite difference scheme into an equivalent optimization problem. The solution to the scheme emerges as the minimizer of the optimization problem, and as a consequence energy stability and positivity-preserving properties are obtained.

Many problems in solid state physics and quantum chemistry require the solution of the Schrodinger equation in the presence of an electromagnetic field. In this paper, we construct, analyze, and nu...

In this work, we design, analyze, and numerically test an invariant preserving discontinuous Galerkin method for solving the nonlinear Camassa--Holm equation. This model is integrable and admits peakon solitons. The proposed numerical method is high order accurate, and preserves two invariants, momentum and energy, of this nonlinear equation. The $L^2$-stability of the scheme for general solutions is a consequence of the energy preserving property. The numerical simulation results for different types of solutions of the Camassa--Holm equation are provided to illustrate the accuracy and capability of the proposed method.

In this paper, we present entropy satisfying schemes for solving an integro-differential equation that describes the evolution of a population structured with respect to a continuous trait. In [P.-E. Jabin and G. Raoul, J. Math. Biol., 63 (2011), pp. 493--517] solutions are shown to converge toward the so-called evolutionary stable distribution (ESD) as time becomes large, using the relative entropy. At the discrete level, the ESD is shown to be the solution to a quadratic programming problem and can be computed by any well-established nonlinear programing algorithm. The schemes are then shown to satisfy the entropy dissipation inequality on the set where initial data are positive and the numerical solutions tend toward the discrete ESD in time. An alternative algorithm is presented to capture the global ESD for nonnegative initial data, which is made possible due to the mutation mechanism built into the modified scheme. A series of numerical tests are given to confirm both accuracy and the entropy satisf...

We design and analyze up to third order accurate discontinuous Galerkin (DG) methods satisfying a strict maximum principle for Fokker--Planck equations. A procedure is established to identify an effective test set in each computational cell to ensure the desired bounds of numerical averages during time evolution. This is achievable by taking advantage of the two parameters in the numerical flux and a novel decomposition of weighted cell averages. Based on this result, a scaling limiter for the DG method with first order Euler forward time discretization is proposed to solve the one-dimensional Fokker--Planck equations. Strong stability preserving high order time discretizations will keep the maximum principle. It is straightforward to extend the method to two and higher dimensions on rectangular meshes. We also show that a modified limiter can preserve the strict maximum principle for DG schemes solving Fokker--Planck equations. As a consequence, the present schemes preserve steady states. Numerical tests...

In this work, we propose a high-resolution alternating evolution (AE) scheme to solve Hamilton--Jacobi equations. The construction of the AE scheme is based on an alternating evolution system of the Hamilton--Jacobi equation, following the idea previously developed for hyperbolic conservation laws. A semidiscrete scheme derives directly from a sampling of this system on alternating grids. Higher order accuracy is achieved by a combination of high order nonoscillatory polynomial reconstruction from the obtained grid values and a time discretization with matching accuracy. Local AE schemes are made possible by choosing the scale parameter $\epsilon$ to reflect the local distribution of waves. The AE schemes have the advantage of easy formulation and implementation and efficient computation of the solution. For the first local AE scheme and the second order local AE scheme with a limiter, we prove the numerical stability in the sense of satisfying the maximum principle. Numerical experiments for a set of Ham...

In this paper, we propose an entropy satisfying conservative method to solve the Fokker–Planck equation of the finitely extensible nonlinear elastic dumbbell model for polymers, subject to homogeneous fluids. Both semidiscrete and fully discrete schemes satisfy all three desired properties—(i) mass conservation, (ii) positivity preserving, and (iii) entropy satisfying—in the sense that these schemes satisfy discrete entropy inequalities for both the physical entropy and the quadratic entropy. These ensure that the computed solution is a probability density and the schemes are entropy stable and preserve the equilibrium solutions. We also prove convergence of the numerical solution to the equilibrium solution as time becomes large. Zero flux at boundary is naturally incorporated, and boundary behavior is resolved sharply. Both one- and two-dimensional numerical results are provided to demonstrate the good qualities of the scheme and the effects of some canonical homogeneous flows.

The finitely extensible nonlinear elastic (FENE) dumbbell model consists of the incompressible Navier--Stokes equation for the solvent and the Fokker--Planck equation for the polymer distribution. In such a model, the polymer elongation cannot exceed a limit $m_0$ which yields all interesting features of solutions near this limit. This work is concerned with the sharpness of boundary conditions in terms of the dimensionless parameter $b=\frac{Hm_0^2}{k_BT}$. Through a careful analysis of the Fokker--Planck operator coupled with the Navier--Stokes equation, we establish a local well-posedness for the full coupled FENE dumbbell model under a class of Dirichlet-type boundary conditions dictated by the parameter $b$. For each $b>0$ we identify a sharp boundary requirement for the underlying density distribution, while the sharpness follows from the existence result for each specification of the boundary behavior. It is shown that the probability density governed by the Fokker--Planck equation approaches zero ...

The alternating evolution (AE) system of Liu [J. Hyperbolic Differ. Equ., 5 (2008), pp. 421-447], $\partial_t u+ \partial_x f(v) =\frac{1}{\epsilon}(v-u)$, $\partial_t v + \partial_x f(u) =\frac{1}{\epsilon}(u-v)$, serves as a refined description of systems of hyperbolic conservation laws $\partial_t \phi+ \partial_x f(\phi)=0$, $\phi(x, 0)=\phi_0(x)$. The solution of conservation laws is precisely captured when two components take the same initial value as $\phi_0$, or is approached by two components exponentially fast when $\epsilon \downarrow 0$ if two initial states are sufficiently close. This nice property enables us to construct novel shock capturing schemes by sampling the AE system on alternating grids. In this paper we develop a class of local AE schemes by taking advantage of the AE system. Our approach is based on an average of the AE system over a hypercube centered at $x$ with vertices at $x\pm\Delta x$. The numerical scheme is then constructed by sampling the averaged system over alternating grids. Higher order accuracy is achieved by a combination of high order nonoscillatory polynomial reconstruction from the obtained averages and a matching Runge-Kutta solver in time discretization. Local AE schemes are made possible by letting the scale parameter $\epsilon$ reflect the local distribution of nonlinear waves. The AE schemes have the advantage of easier formulation and implementation, and efficient computation of the solution. For the first and second order local AE schemes applied to one-dimensional scalar conservation laws, we prove the numerical stability in the sense of satisfying the maximum principle and total variation diminishing property. The formulation procedure of AE schemes in multiple dimensions is given, followed by both the first and second order AE schemes for two-dimensional conservation laws. Numerical experiments for both scalar conservation laws and compressible Euler equations are presented to demonstrate the high order accuracy and capacity of these AE schemes.

Computation of high frequency solutions to wave equations is important in many applications, and notoriously difficult in resolving wave oscillations. Gaussian beams are asymptotically valid high frequency solutions concentrated on a single curve through the physical domain, and superposition of Gaussian beams provides a powerful tool for generating more general high frequency solutions to PDEs. An alternative way to compute Gaussian beam components such as phase, amplitude, and Hessian of the phase is to capture them in phase space by solving Liouville-type equations on uniform grids. In this work we review and extend recent constructions of asymptotic high frequency wave fields from computations in phase space. We give a new level set method of computing the Hessian and higher derivatives of the phase. Moreover, we prove that the kth order phase space–based Gaussian beam superposition converges to the original wave field in $L^2$ at the rate of $\epsilon^{\frac{k}{2}-\frac{n}{4}}$ in dimension n.

A new discontinuous Galerkin finite element method for solving diffusion problems is introduced. Unlike the traditional local discontinuous Galerkin method, the scheme called the direct discontinuous Galerkin (DDG) method is based on the direct weak formulation for solutions of parabolic equations in each computational cell and lets cells communicate via the numerical flux $\widehat{u_x}$ only. We propose a general numerical flux formula for the solution derivative, which is consistent and conservative; and we then introduce a concept of admissibility to identify a class of numerical fluxes so that the nonlinear stability for both one-dimensional (1D) and multidimensional problems are ensured. Furthermore, when applying the DDG scheme with admissible numerical flux to the 1D linear case, $k$th order accuracy in an energy norm is proven when using $k$th degree polynomials. The DDG method has the advantage of easier formulation and implementation and efficient computation of the solution. A series of numerical examples are presented to demonstrate the high order accuracy of the method. In particular, we study the numerical performance of the scheme with different admissible numerical fluxes.

We consider the microscopic equation of finite extensible nonlinear elasticity (FENE) models for polymeric fluids under a steady flow field. It is shown that for the underlying Fokker–Planck type of equations, any preassigned distribution on the boundary will become redundant once the nondimensional number $\text{{\it Li\/}} := \frac{Hb}{k_BT} \geq 2$, where H is the elasticity constant, $\sqrt{b}$ is the maximum dumbbell extension, T is the temperature, and $k_B$ is the usual Boltzmann constant. Moreover, if the probability density function is regular enough for its trace to be defined on the sphere $|m| = \sqrt{b}$, then the trace is necessarily zero when $\text{{\it Li\/}} > 2$. These results are consistent with our numerical simulations as well as some recent well-posedness results by preassuming a zero boundary distribution.

We provide a complete description of the critical threshold phenomenon for the two-dimensional localized Euler-Poisson equations, introduced by the authors in [Comm. Math. Phys., 228 (2002), pp. 435-466]. Here, the questions of global regularity vs. finite-time breakdown for the two-dimensional (2D) restricted Euler-Poisson solutions are classified in terms of precise explicit formulae, describing a remarkable variety of critical threshold surfaces of initial configurations. In particular, it is shown that the 2D critical thresholds depend on the relative sizes of three quantities: the initial density, the initial divergence, and the initial spectral gap, that is, the difference between the two eigenvalues of the 2 × 2 initial velocity gradient.

The phenomena of concentration and cavitation and the formation of $\delta$-shocks and vacuum states in solutions to the Euler equations for isentropic fluids are identified and analyzed as the pressure vanishes. It is shown that, as the pressure vanishes, any two-shock Riemann solution to the Euler equations for isentropic fluids tends to a $\delta$-shock solution to the Euler equations for pressureless fluids, and the intermediate density between the two shocks tends to a weighted $\delta$-measure that forms the $\delta$-shock. By contrast, any two-rarefaction-wave Riemann solution of the Euler equations for isentropic fluids is shown to tend to a two-contact-discontinuity solution to the Euler equations for pressureless fluids, whose intermediate state between the two contact discontinuities is a vacuum state, even when the initial data stays away from the vacuum. Some numerical results exhibiting the formation process of $\delta$-shocks are also presented.

In this work we consider a convolution model for nonlinear conservation laws. Due to the delicate balance between the nonlinear convection and the nonlocal forcing, this model allows for narrower shock layers than those in the viscous Burgers' equation and yet exhibits the conditional finite time breakdown as in the damped Burgers' equation. We show the critical threshold phenomenon by presenting a lower threshold for the breakdown of the solutions and an upper threshold for the global existence of the smooth solution. The threshold condition depends only on the relative size of the minimum slope of the initial velocity and its maximal variation. We show the exact blow-up rate when the slope of the initial profile is below the lower threshold. We further prove the L1 stability of the smooth shock profile, provided the slope of the initial profile is above the critical threshold.

In this paper, we prove a global error estimate for a relaxation scheme approximating scalar conservation laws. To this end, we decompose the error into a relaxation error and a discretization error. Including an initial error $\omega(\ep)$ we obtain the rate of convergence of $\sqrt{\ep}$ in L1 for the relaxation step. The estimate here is independent of the type of nonlinearity. In the discretization step a convergence rate of $\sqrt{\Del x} $ in L1 is obtained. These rates are independent of the choice of initial error $\omega(\ep)$. Thereby, we obtain the order 1/2 for the total error.

We show the discrete lip+-stability for a relaxation scheme proposed by Jin and Xin [ Comm. Pure Appl. Math., 48 (1995), pp. 235--277] to approximate convex conservation laws. Equipped with the lip+-stability we obtain global error estimates in the spaces W s,p for $-1 \leq s \leq 1/p$, $1 \leq p \leq \infty$ and pointwise error estimates for the approximate solution obtained by the relaxation scheme. The proof uses the framework introduced by Nessyahu and Tadmor [SIAM J. Numer. Anal., 29 (1992), pp. 1505--1519]. We also show a maximum principle for the relaxation scheme when the initial data are in an equilibrium state.