For the solution of the cubic nonlinear Schr\"odinger equation in one space dimension, we propose and analyse a fully discrete low-regularity integrator. The scheme is explicit and can easily be implemented using the fast Fourier transform with a complexity of $\mathcal{O}(N\log N)$ operations per time step, where $N$ denotes the degrees of freedom in the spatial discretisation. We prove that the new scheme provides an $\mathcal{O}(\tau^{\frac32\gamma-\frac12-\varepsilon}+N^{-\gamma})$ error bound in $L^2$ for any initial data belonging to $H^\gamma$, $\frac12<\gamma\leq 1$, where $\tau$ denotes the temporal step size. Numerical examples illustrate this convergence behavior.
We consider a Johnson-N\'ed\'elec FEM-BEM coupling, which is a direct and non-symmetric coupling of finite and boundary element methods, in order to solve interface problems for the magnetostatic Maxwell's equations with the magnetic vector potential ansatz. In the FEM-domain, equations may be non-linear, whereas they are exclusively linear in the BEM-part to guarantee the existence of a fundamental solution. First, the weak problem is formulated in quotient spaces to avoid resolving to a saddle point problem. Second, we establish in this setting well-posedness of the arising problem using the framework of Lipschitz and strongly monotone operators as well as a stability result for a special type of non-linearity, which is typically considered in magnetostatic applications. Then, the discretization is performed in the isogeometric context, i.e., the same type of basis functions that are used for geometry design are considered as ansatz functions for the discrete setting. In particular, NURBS are employed for geometry considerations, and B-Splines, which can be understood as a special type of NURBS, for analysis purposes. In this context, we derive a priori estimates w.r.t. h-refinement, and point out to an interesting behavior of BEM, which consists in an amelioration of the convergence rates, when a functional of the solution is evaluated in the exterior BEM-domain. This improvement may lead to a doubling of the convergence rate under certain assumptions. Finally, we end the paper with a numerical example to illustrate the theoretical results, along with a conclusion and an outlook.
In this technical note, we consider a dynamic linear, cantilevered rectangular plate. The evolutionary PDE model is given by the fourth order plate dynamics (via the spatial biharmonic operator) with clamped-free-free-free boundary conditions. We additionally consider damping/dissipation terms, as well as non-conservative lower order terms arising in various applications. Dynamical numerical simulations are achieved by way of a finite difference spatial approximation with a MATLAB time integrator. The rectangular geometry allows the use of standard 2D spatial finite differences, while the high spatial order of the problem and mixed clamped-free type boundary conditions present challenges. Dynamic energies are also computed. The relevant code is presented, with discussion of the model and context.
In molecular dynamics, penalized overdamped Langevin dynamics are used to model the motion of a set of particles that follow constraints up to a parameter $\varepsilon$. The most used schemes for simulating these dynamics are the Euler integrator in $\mathbb{R}^d$ and the constrained Euler integrator. Both have weak order one of accuracy, but work properly only in specific regimes depending on the size of the parameter $\varepsilon$. We propose in this paper a new consistent method with an accuracy independent of $\varepsilon$ for solving penalized dynamics on a manifold of any dimension. Moreover, this method converges to the constrained Euler scheme when $\varepsilon$ goes to zero. The numerical experiments confirm the theoretical findings, in the context of weak convergence and for the invariant measure, on a torus and on the orthogonal group in high dimension and high codimension.
We consider the numerical taxonomy problem of fitting a positive distance function ${D:{S\choose 2}\rightarrow \mathbb R_{>0}}$ by a tree metric. We want a tree $T$ with positive edge weights and including $S$ among the vertices so that their distances in $T$ match those in $D$. A nice application is in evolutionary biology where the tree $T$ aims to approximate the branching process leading to the observed distances in $D$ [Cavalli-Sforza and Edwards 1967]. We consider the total error, that is the sum of distance errors over all pairs of points. We present a deterministic polynomial time algorithm minimizing the total error within a constant factor. We can do this both for general trees, and for the special case of ultrametrics with a root having the same distance to all vertices in $S$. The problems are APX-hard, so a constant factor is the best we can hope for in polynomial time. The best previous approximation factor was $O((\log n)(\log \log n))$ by Ailon and Charikar [2005] who wrote "Determining whether an $O(1)$ approximation can be obtained is a fascinating question".
We study the problem of the nonparametric estimation for the density $\pi$ of the stationary distribution of a $d$-dimensional stochastic differential equation $(X_t)_{t \in [0, T]}$. From the continuous observation of the sampling path on $[0, T]$, we study the rate of estimation of $\pi(x)$ as $T$ goes to infinity. One finding is that, for $d \ge 3$, the rate of estimation depends on the smoothness $\beta = (\beta_1, ... , \beta_d)$ of $\pi$. In particular, having ordered the smoothness such that $\beta_1 \le ... \le \beta_d$, it depends on the fact that $\beta_2 < \beta_3$ or $\beta_2 = \beta_3$. We show that kernel density estimators achieve the rate $(\frac{\log T}{T})^\gamma$ in the first case and $(\frac{1}{T})^\gamma$ in the second, for an explicit exponent $\gamma$ depending on the dimension and on $\bar{\beta}_3$, the harmonic mean of the smoothness over the $d$ directions after having removed $\beta_1$ and $\beta_2$, the smallest ones. Moreover, we obtain a minimax lower bound on the $\mathbf{L}^2$-risk for the pointwise estimation with the same rates $(\frac{\log T}{T})^\gamma$ or $(\frac{1}{T})^\gamma$, depending on the value of $\beta_2$ and $\beta_3$.
We provide a control-theoretic perspective on optimal tensor algorithms for minimizing a convex function in a finite-dimensional Euclidean space. Given a function $\Phi: \mathbb{R}^d \rightarrow \mathbb{R}$ that is convex and twice continuously differentiable, we study a closed-loop control system that is governed by the operators $\nabla \Phi$ and $\nabla^2 \Phi$ together with a feedback control law $\lambda(\cdot)$ satisfying the algebraic equation $(\lambda(t))^p\|\nabla\Phi(x(t))\|^{p-1} = \theta$ for some $\theta \in (0, 1)$. Our first contribution is to prove the existence and uniqueness of a local solution to this system via the Banach fixed-point theorem. We present a simple yet nontrivial Lyapunov function that allows us to establish the existence and uniqueness of a global solution under certain regularity conditions and analyze the convergence properties of trajectories. The rate of convergence is $O(1/t^{(3p+1)/2})$ in terms of objective function gap and $O(1/t^{3p})$ in terms of squared gradient norm. Our second contribution is to provide two algorithmic frameworks obtained from discretization of our continuous-time system, one of which generalizes the large-step A-HPE framework and the other of which leads to a new optimal $p$-th order tensor algorithm. While our discrete-time analysis can be seen as a simplification and generalization of~\citet{Monteiro-2013-Accelerated}, it is largely motivated by the aforementioned continuous-time analysis, demonstrating the fundamental role that the feedback control plays in optimal acceleration and the clear advantage that the continuous-time perspective brings to algorithmic design. A highlight of our analysis is that we show that all of the $p$-th order optimal tensor algorithms that we discuss minimize the squared gradient norm at a rate of $O(k^{-3p})$, which complements the recent analysis.
We consider the numerical approximation of the inertial Landau-Lifshitz-Gilbert equation (iLLG), which describes the dynamics of the magnetization in ferromagnetic materials at subpicosecond time scales. We propose and analyze two fully discrete numerical schemes: The first method is based on a reformulation of the problem as a linear constrained variational formulation for the linear velocity. The second method exploits a reformulation of the problem as a first order system in time for the magnetization and the angular momentum. Both schemes are implicit, based on first-order finite elements, and generate approximations satisfying the unit-length constraint of iLLG at the vertices of the underlying mesh. For both methods, we prove convergence of the approximations towards a weak solution of the problem. Numerical experiments validate the theoretical results and show the applicability of the methods for the simulation of ultrafast magnetic processes.
Fully implicit Runge-Kutta (IRK) methods have many desirable properties as time integration schemes in terms of accuracy and stability, but high-order IRK methods are not commonly used in practice with numerical PDEs due to the difficulty of solving the stage equations. This paper introduces a theoretical and algorithmic preconditioning framework for solving the systems of equations that arise from IRK methods applied to linear numerical PDEs (without algebraic constraints). This framework also naturally applies to discontinuous Galerkin discretizations in time. Under quite general assumptions on the spatial discretization that yield stable time integration, the preconditioned operator is proven to have condition number bounded by a small, order-one constant, independent of the spatial mesh and time-step size, and with only weak dependence on number of stages/polynomial order; for example, the preconditioned operator for 10th-order Gauss IRK has condition number less than two, independent of the spatial discretization and time step. The new method can be used with arbitrary existing preconditioners for backward Euler-type time stepping schemes, and is amenable to the use of three-term recursion Krylov methods when the underlying spatial discretization is symmetric. The new method is demonstrated to be effective on various high-order finite-difference and finite-element discretizations of linear parabolic and hyperbolic problems, demonstrating fast, scalable solution of up to 10th order accuracy. The new method consistently outperforms existing block preconditioning approaches, and in several cases, the new method can achieve 4th-order accuracy using Gauss integration with roughly half the number of preconditioner applications and wallclock time as required using standard diagonally implicit RK methods.
Fully implicit Runge-Kutta (IRK) methods have many desirable accuracy and stability properties as time integration schemes, but high-order IRK methods are not commonly used in practice with large-scale numerical PDEs because of the difficulty of solving the stage equations. This paper introduces a theoretical and algorithmic framework for solving the nonlinear equations that arise from IRK methods (and discontinuous Galerkin discretizations in time) applied to nonlinear numerical PDEs, including PDEs with algebraic constraints. Several new linearizations of the nonlinear IRK equations are developed, offering faster and more robust convergence than the often-considered simplified Newton, as well as an effective preconditioner for the true Jacobian if exact Newton iterations are desired. Inverting these linearizations requires solving a set of block 2x2 systems. Under quite general assumptions, it is proven that the preconditioned 2x2 operator's condition number is bounded by a small constant close to one, independent of the spatial discretization, spatial mesh, and time step, and with only weak dependence on the number of stages or integration accuracy. Moreover, the new method is built using the same preconditioners needed for backward Euler-type time stepping schemes, so can be readily added to existing codes. The new methods are applied to several challenging fluid flow problems, including the compressible Euler and Navier Stokes equations, and the vorticity-streamfunction formulation of the incompressible Euler and Navier Stokes equations. Up to 10th-order accuracy is demonstrated using Gauss IRK, while in all cases 4th-order Gauss IRK requires roughly half the number of preconditioner applications as required by standard SDIRK methods.
In order to avoid the curse of dimensionality, frequently encountered in Big Data analysis, there was a vast development in the field of linear and nonlinear dimension reduction techniques in recent years. These techniques (sometimes referred to as manifold learning) assume that the scattered input data is lying on a lower dimensional manifold, thus the high dimensionality problem can be overcome by learning the lower dimensionality behavior. However, in real life applications, data is often very noisy. In this work, we propose a method to approximate $\mathcal{M}$ a $d$-dimensional $C^{m+1}$ smooth submanifold of $\mathbb{R}^n$ ($d \ll n$) based upon noisy scattered data points (i.e., a data cloud). We assume that the data points are located "near" the lower dimensional manifold and suggest a non-linear moving least-squares projection on an approximating $d$-dimensional manifold. Under some mild assumptions, the resulting approximant is shown to be infinitely smooth and of high approximation order (i.e., $O(h^{m+1})$, where $h$ is the fill distance and $m$ is the degree of the local polynomial approximation). The method presented here assumes no analytic knowledge of the approximated manifold and the approximation algorithm is linear in the large dimension $n$. Furthermore, the approximating manifold can serve as a framework to perform operations directly on the high dimensional data in a computationally efficient manner. This way, the preparatory step of dimension reduction, which induces distortions to the data, can be avoided altogether.