We present a space-time ultra-weak discontinuous Galerkin discretization of the linear Schr\"odinger equation with variable potential. The proposed method is well-posed and quasi-optimal in mesh-dependent norms for very general discrete spaces. Optimal $h$-convergence error estimates are derived for the method when test and trial spaces are chosen either as piecewise polynomials, or as a novel quasi-Trefftz polynomial space. The latter allows for a substantial reduction of the number of degrees of freedom and admits piecewise-smooth potentials. Several numerical experiments validate the accuracy and advantages of the proposed method.
The hazard function represents one of the main quantities of interest in the analysis of survival data. We propose a general approach for modelling the dynamics of the hazard function using systems of autonomous ordinary differential equations (ODEs). This modelling approach can be used to provide qualitative and quantitative analyses of the evolution of the hazard function over time. Our proposal capitalises on the extensive literature of ODEs which, in particular, allow for establishing basic rules or laws on the dynamics of the hazard function via the use of autonomous ODEs. We show how to implement the proposed modelling framework in cases where there is an analytic solution to the system of ODEs or where an ODE solver is required to obtain a numerical solution. We focus on the use of a Bayesian modelling approach, but the proposed methodology can also be coupled with maximum likelihood estimation. A simulation study is presented to illustrate the performance of these models and the interplay of sample size and censoring. Two case studies using real data are presented to illustrate the use of the proposed approach and to highlight the interpretability of the corresponding models. We conclude with a discussion on potential extensions of our work and strategies to include covariates into our framework.
We describe a `discretize-then-relax' strategy to globally minimize integral functionals over functions $u$ in a Sobolev space satisfying prescribed Dirichlet boundary conditions. The strategy applies whenever the integral functional depends polynomially on $u$ and its derivatives, even if it is nonconvex. The `discretize' step uses a bounded finite-element scheme to approximate the integral minimization problem with a convergent hierarchy of polynomial optimization problems over a compact feasible set, indexed by the decreasing size $h$ of the finite-element mesh. The `relax' step employs sparse moment-SOS relaxations to approximate each polynomial optimization problem with a hierarchy of convex semidefinite programs, indexed by an increasing relaxation order $\omega$. We prove that, as $\omega\to\infty$ and $h\to 0$, solutions of such semidefinite programs provide approximate minimizers that converge in $L^p$ to the global minimizer of the original integral functional if this is unique. We also report computational experiments that show our numerical strategy works well even when technical conditions required by our theoretical analysis are not satisfied.
We propose a matrix-free parallel two-level-deflation preconditioner combined with the Complex Shifted Laplacian preconditioner(CSLP) for the two-dimensional Helmholtz problems. The Helmholtz equation is widely studied in seismic exploration, antennas, and medical imaging. It is one of the hardest problems to solve both in terms of accuracy and convergence, due to scalability issues of the numerical solvers. Motivated by the observation that for large wavenumbers, the eigenvalues of the CSLP-preconditioned system shift towards zero, deflation with multigrid vectors, and further high-order vectors were incorporated to obtain wave-number-independent convergence. For large-scale applications, high-performance parallel scalable methods are also indispensable. In our method, we consider the preconditioned Krylov subspace methods for solving the linear system obtained from finite-difference discretization. The CSLP preconditioner is approximated by one parallel geometric multigrid V-cycle. For the two-level deflation, the matrix-free Galerkin coarsening as well as high-order re-discretization approaches on the coarse grid are studied. The results of matrix-vector multiplications in Krylov subspace methods and the interpolation/restriction operators are implemented based on the finite-difference grids without constructing any coefficient matrix. These adjustments lead to direct improvements in terms of memory consumption. Numerical experiments of model problems show that wavenumber independence has been obtained for medium wavenumbers. The matrix-free parallel framework shows satisfactory weak and strong parallel scalability.
The Helmholtz equation is related to seismic exploration, sonar, antennas, and medical imaging applications. It is one of the most challenging problems to solve in terms of accuracy and convergence due to the scalability issues of the numerical solvers. For 3D large-scale applications, high-performance parallel solvers are also needed. In this paper, a matrix-free parallel iterative solver is presented for the three-dimensional (3D) heterogeneous Helmholtz equation. We consider the preconditioned Krylov subspace methods for solving the linear system obtained from finite-difference discretization. The Complex Shifted Laplace Preconditioner (CSLP) is employed since it results in a linear increase in the number of iterations as a function of the wavenumber. The preconditioner is approximately inverted using one parallel 3D multigrid cycle. For parallel computing, the global domain is partitioned blockwise. The matrix-vector multiplication and preconditioning operator are implemented in a matrix-free way instead of constructing large, memory-consuming coefficient matrices. Numerical experiments of 3D model problems demonstrate the robustness and outstanding strong scaling of our matrix-free parallel solution method. Moreover, the weak parallel scalability indicates our approach is suitable for realistic 3D heterogeneous Helmholtz problems with minimized pollution error.
For terminal value problems of fractional differential equations of order $\alpha \in (0,1)$ that use Caputo derivatives, shooting methods are a well developed and investigated approach. Based on recently established analytic properties of such problems, we develop a new technique to select the required initial values that solves such shooting problems quickly and accurately. Numerical experiments indicate that this new proportional secting technique converges very quickly and accurately to the solution. Run time measurements indicate a speedup factor of between 4 and 10 when compared to the standard bisection method.
Invariant finite-difference schemes for the one-dimensional shallow water equations in the presence of a magnetic field for various bottom topographies are constructed. Based on the results of the group classification recently carried out by the authors, finite-difference analogues of the conservation laws of the original differential model are obtained. Some typical problems are considered numerically, for which a comparison is made between the cases of a magnetic field presence and when it is absent (the standard shallow water model). The invariance of difference schemes in Lagrangian coordinates and the energy preservation on the obtained numerical solutions are also discussed.
Partial differential equations (PDEs) with uncertain or random inputs have been considered in many studies of uncertainty quantification. In forward uncertainty quantification, one is interested in analyzing the stochastic response of the PDE subject to input uncertainty, which usually involves solving high-dimensional integrals of the PDE output over a sequence of stochastic variables. In practical computations, one typically needs to discretize the problem in several ways: approximating an infinite-dimensional input random field with a finite-dimensional random field, spatial discretization of the PDE using, e.g., finite elements, and approximating high-dimensional integrals using cubatures such as quasi-Monte Carlo methods. In this paper, we focus on the error resulting from dimension truncation of an input random field. We show how Taylor series can be used to derive theoretical dimension truncation rates for a wide class of problems and we provide a simple checklist of conditions that a parametric mathematical model needs to satisfy in order for our dimension truncation error bound to hold. Some of the novel features of our approach include that our results are applicable to non-affine parametric operator equations, dimensionally-truncated conforming finite element discretized solutions of parametric PDEs, and even compositions of PDE solutions with smooth nonlinear quantities of interest. As a specific application of our method, we derive an improved dimension truncation error bound for elliptic PDEs with lognormally parameterized diffusion coefficients. Numerical examples support our theoretical findings.
The discovery of equations with knowledge of the process origin is a tempting prospect. However, most equation discovery tools rely on gradient methods, which offer limited control over parameters. An alternative approach is the evolutionary equation discovery, which allows modification of almost every optimization stage. In this paper, we examine the modifications that can be introduced into the evolutionary operators of the equation discovery algorithm, taking inspiration from directed evolution techniques employed in fields such as chemistry and biology. The resulting approach, dubbed directed equation discovery, demonstrates a greater ability to converge towards accurate solutions than the conventional method. To support our findings, we present experiments based on Burgers', wave, and Korteweg--de Vries equations.
We revisit the relation between the gradient-flow equations and Hamilton's equations in information geometry. By regarding the gradient-flow equations as Huygens' equations in geometric optics, we have related the gradient flows to the geodesic flows induced by the geodesic Hamiltonian in an appropriate Riemannian geometry. The original evolution parameter $\textit{t}$ in the gradient-flow equations is related to the arc-length parameter in the associated Riemannian manifold by Jacobi-Maupertuis transformation. As a by-product, it is found the relation between the gradient-flow equation and replicator equations.
Ghost, or fictitious points allow to capture boundary conditions that are not located on the finite difference grid discretization. We explore in this paper the impact of ghost points on the stability of the explicit Euler finite difference scheme in the context of the diffusion equation. In particular, we consider the case of a one-touch option under the Black-Scholes model. The observations and results are however valid for a much wider range of financial contracts and models.