Jon Wilkening

C. H. Arthur Cheng, Rafael Granero-Belinchon, Steve Shkoller et al.

We develop a rigorous asymptotic derivation for two mathematical models of water waves that capture the full nonlinearity of the Euler equations up to quadratic and cubic interactions, respectively. Specifically, letting epsilon denote an asymptotic parameter denoting the steepness of the water wave, we use a Stokes expansion in epsilon to derive a set of linear recursion relations for the tangential component of velocity, the stream function, and the water wave parameterization. The solution of the water waves system is obtained as an infinite sum of solutions to linear problems at each epsilon^k level, and truncation of this series leads to our two asymptotic models, that we call the quadratic and cubic h-models. Using the growth rate of the Catalan numbers (from number theory), we prove well-posedness of the h-models in spaces of analytic functions, and prove error bounds for solutions of the h-models compared against solutions of the water waves system. We also show that the Craig-Sulem models of water waves can be obtained from our asymptotic procedure and that their WW2 model is well-posed in our functional framework. We then develop a novel numerical algorithm to solve the quadratic and cubic h-models as well as the full water waves system. For three very different examples, we show that the agreement between the model equations and the water waves solution is excellent, even when the wave steepness is quite large. We also present a numerical example of corner formation for water waves.

Jon Wilkening, Antoine Cerfon, Matt Landreman

We examine the merits of using a family of polynomials that are orthogonal with respect to a non-classical weight function to discretize the speed variable in continuum kinetic calculations. We consider a model one-dimensional partial differential equation describing energy diffusion in velocity space due to Fokker-Planck collisions. This relatively simple case allows us to compare the results of the projected dynamics with an expensive but highly accurate spectral transform approach. It also allows us to integrate in time exactly, and to focus entirely on the effectiveness of the discretization of the speed variable. We show that for a fixed number of modes or grid points, the non-classical polynomials can be many orders of magnitude more accurate than classical Hermite polynomials or finite-difference solvers for kinetic equations in plasma physics. We provide a detailed analysis of the difference in behavior and accuracy of the two families of polynomials. For the non-classical polynomials, if the initial condition is not smooth at the origin when interpreted as a three-dimensional radial function, the exact solution leaves the polynomial subspace for a time, but returns (up to roundoff accuracy) to the same point evolved to by the projected dynamics in that time. By contrast, using classical polynomials, the exact solution differs significantly from the projected dynamics solution when it returns to the subspace. We also explore the connection between eigenfunctions of the projected evolution operator and (non-normalizable) eigenfunctions of the full evolution operator, as well as the effect of truncating the computational domain.

Michael Westdickenberg, Jon Wilkening

Both the porous medium equation and the system of isentropic Euler equations can be considered as steepest descents on suitable manifolds of probability measures in the framework of optimal transport theory. By discretizing these variational characterizations instead of the partial differential equations themselves, we obtain new schemes with remarkable stability properties. We show that they capture successfully the nonlinear features of the flows, such as shocks and rarefaction waves for the isentropic Euler equations. We also show how to design higher order methods for these problems in the optimal transport setting using backward differentiation formula (BDF) multi-step methods or diagonally implicit Runge-Kutta methods.

Charles L. Epstein, Jon Wilkening

We study the classical Kimura diffusion operator defined on the n-simplex, $$L^{Kim}=\sum_{1\leq i,j\leq n+1}x_ix_j\partial_{x_i}\partial_{x_j}$$ We give novel constructions for the basis of eigenpolynomials, and the solution to the inhomogeneous Dirichlet problem, which are well adapted to numerical applications. Our solution of the Dirichlet problem is quite explicit and provides a precise description of the singularities that arise along the boundary.

Eugenia Kim, Jon Wilkening

The dynamics of the magnetic distribution in a ferromagnetic material is governed by the Landau-Lifshitz equation, which is a nonlinear geometric dispersive equation with a nonconvex constraint that requires the magnetization to remain of unit length throughout the domain. In this article, we present a mass-lumped finite element method for the Landau-Lifshitz equation. This method preserves the nonconvex constraint at each node of the finite element mesh, and is energy nonincreasing. We show that the numerical solution of our method for the Landau-Lifshitz equation converges to a weak solution of the Landau-Lifshitz-Gilbert equation using a simple proof technique that cancels out the product of weakly convergent sequences. Numerical tests for both explicit and implicit versions of the method on a unit square with periodic boundary conditions are provided for structured and unstructured meshes.

Saad Qadeer, Jon Wilkening

The computation of the Dirichlet-Neumann operator for the Laplace equation is the primary challenge for the numerical simulation of the ideal fluid equations. The techniques used commonly for 2D fluids, such as conformal mapping and boundary integral methods, fail to generalize suitably to 3D. In this study, we address this problem by developing a Transformed Field Expansion method for computing the Dirichlet-Neumann operator in a cylindrical geometry with a variable upper boundary. This technique reduces the problem to a sequence of Poisson equations on a flat geometry. We design a fast and accurate solver for these sub-problems, a key ingredient being the use of Zernike polynomials for the circular cross-section instead of the traditional Bessel functions. This lends spectral accuracy to the method as well as allowing significant computational speed-up. We rigorously analyze the algorithm and prove its applicability to a wide class of problems before demonstrating its effectiveness numerically.

Matthias Morzfeld, Xuemin Tu, Jon Wilkening et al.

Implicit sampling is a weighted sampling method that is used in data assimilation, where one sequentially updates estimates of the state of a stochastic model based on a stream of noisy or incomplete data. Here we describe how to use implicit sampling in parameter estimation problems, where the goal is to find parameters of a numerical model, e.g.~a partial differential equation (PDE), such that the output of the numerical model is compatible with (noisy) data. We use the Bayesian approach to parameter estimation, in which a posterior probability density describes the probability of the parameter conditioned on data and compute an empirical estimate of this posterior with implicit sampling. Our approach generates independent samples, so that some of the practical difficulties one encounters with Markov Chain Monte Carlo methods, e.g.~burn-in time or correlations among dependent samples, are avoided. We describe a new implementation of implicit sampling for parameter estimation problems that makes use of multiple grids (coarse to fine) and BFGS optimization coupled to adjoint equations for the required gradient calculations. The implementation is "dimension independent", in the sense that a well-defined finite dimensional subspace is sampled as the mesh used for discretization of the PDE is refined. We illustrate the algorithm with an example where we estimate a diffusion coefficient in an elliptic equation from sparse and noisy pressure measurements. In the example, dimension\slash mesh-independence is achieved via Karhunen-Loève expansions.

Jeremiah Birrell, Jon Wilkening, Johann Rafelski

We present a novel method to solve the spatially homogeneous and isotropic relativistic Boltzmann equation. We employ a basis set of orthogonal polynomials dynamically adapted to allow for emergence of chemical non-equilibrium. Two time dependent parameters characterize the set of orthogonal polynomials, the effective temperature $T(t)$ and phase space occupation factor $Υ(t)$. In this first paper we address (effectively) massless fermions and derive dynamical equations for $T(t)$ and $Υ(t)$ such that the zeroth order term of the basis alone captures the particle number density and energy density of each particle distribution. We validate our method and illustrate the reduced computational cost and the ability to easily represent final state chemical non-equilibrium by studying a model problem that is motivated by the physics of the neutrino freeze-out processes in the early Universe, where the essential physical characteristics include reheating from another disappearing particle component ($e^\pm$-annihilation).