We consider the scalar conservation law in one space dimension with a genuinely nonlinear flux. We assume that an appropriate velocity function depending on the entropy solution of the conservation law is given for the comprising particles, and study their corresponding trajectories under the flow. The differential equation that each of these trajectories satisfies depends on the entropy solution of the conservation law which is typically discontinuous in both time and space variables. The existence and uniqueness of these trajectories are guaranteed by the Filippov theory of differential equations. We show that such a Filippov solution is compatible with the front tracking and vanishing viscosity approximations in the sense that the approximate trajectories given by either of these methods converge uniformly to the trajectories corresponding to the entropy solution of the scalar conservation law. For certain classes of flux functions, illustrated by traffic flow, we prove the H\"older continuity of the particle trajectories with respect to the initial field or the flux function. We then consider the inverse problem of recovering the initial field or the flux function of the scalar conservation law from discrete pointwise measurements of the particle trajectories. We show that the above continuity properties translate to the stability of the Bayesian regularised solutions of these inverse problems with respect to appropriate approximations of the forward map. We also discuss the limitations of the situation where the same inverse problems are considered with pointwise observations made from the entropy solution itself.
In PDE-constrained optimization, one aims to find design parameters that minimize some objective, subject to the satisfaction of a partial differential equation. A major challenges is computing gradients of the objective to the design parameters, as applying the chain rule requires computing the Jacobian of the design parameters to the PDE's state. The adjoint method avoids this Jacobian by computing partial derivatives of a Lagrangian. Evaluating these derivatives requires the solution of a second PDE with the adjoint differential operator to the constraint, resulting in a backwards-in-time simulation. Particle-based Monte Carlo solvers are often used to compute the solution to high-dimensional PDEs. However, such solvers have the drawback of introducing noise to the computed results, thus requiring stochastic optimization methods. To guarantee convergence in this setting, both the constraint and adjoint Monte Carlo simulations should simulate the same particle trajectories. For large simulations, storing full paths from the constraint equation for re-use in the adjoint equation becomes infeasible due to memory limitations. In this paper, we provide a reversible extension to the family of permuted congruential pseudorandom number generators (PCG). We then use such a generator to recompute these time-reversed paths for the heat equation, avoiding these memory issues.
The problem of optimal recovering high-order mixed derivatives of bivariate functions with finite smoothness is studied. On the basis of the truncation method, an algorithm for numerical differentiation is constructed, which is order-optimal both in the sense of accuracy and in terms of the amount of involved Galerkin information.
Anomalous diffusion in the presence or absence of an external force field is often modelled in terms of the fractional evolution equations, which can involve the hyper-singular source term. For this case, conventional time stepping methods may exhibit a severe order reduction. Although a second-order numerical algorithm is provided for the subdiffusion model with a simple hyper-singular source term $t^{\mu}$, $-2<\mu<-1$ in [arXiv:2207.08447], the convergence analysis remain to be proved. To fill in these gaps, we present a simple and robust smoothing method for the hyper-singular source term, where the Hadamard finite-part integral is introduced. This method is based on the smoothing/ID$m$-BDF$k$ method proposed by the authors [Shi and Chen, SIAM J. Numer. Anal., to appear] for subdiffusion equation with a weakly singular source term. We prove that the $k$th-order convergence rate can be restored for the diffusion-wave case $\gamma \in (1,2)$ and sketch the proof for the subdiffusion case $\gamma \in (0,1)$, even if the source term is hyper-singular and the initial data is not compatible. Numerical experiments are provided to confirm the theoretical results.
We study an optimal control problem governed by elliptic PDEs with interface, which the control acts on the interface. Due to the jump of the coefficient across the interface and the control acting on the interface, the regularity of solution of the control problem is limited on the whole domain, but smoother on subdomains. The control function with pointwise inequality constraints is served as the flux jump condition which we called Neumann interface control. We use a simple uniform mesh that is independent of the interface. The standard linear finite element method can not achieve optimal convergence when the uniform mesh is used. Therefore the state and adjoint state equations are discretized by piecewise linear immersed finite element method (IFEM). While the accuracy of the piecewise constant approximation of the optimal control on the interface is improved by a postprocessing step which possesses superconvergence properties; as well as the variational discretization concept for the optimal control is used to improve the error estimates. Optimal error estimates for the control, suboptimal error estimates for state and adjoint state are derived. Numerical examples with and without constraints are provided to illustrate the effectiveness of the proposed scheme and correctness of the theoretical analysis.
We present a general theory to quantify the uncertainty from imposing structural assumptions on the second-order structure of nonstationary Hilbert space-valued processes, which can be measured via functionals of time-dependent spectral density operators. The second-order dynamics are well-known to be elements of the space of trace-class operators, the latter is a Banach space of type 1 and of cotype 2, which makes the development of statistical inference tools more challenging. A part of our contribution is to obtain a weak invariance principle as well as concentration inequalities for (functionals of) the sequential time-varying spectral density operator. In addition, we introduce deviation measures in the nonstationary context, and derive estimators that are asymptotically pivotal. We then apply this framework and propose statistical methodology to investigate the validity of structural assumptions for nonstationary response surface data, such as low-rank assumptions in the context of time-varying dynamic fPCA and principle separable component analysis, deviations from stationarity with respect to the square root distance, and deviations from zero functional canonical coherency.
We consider the weighted least squares spline approximation of a noisy dataset. By interpreting the weights as a probability distribution, we maximize the associated entropy subject to the constraint that the mean squared error is prescribed to a desired (small) value. Acting on this error yields a robust regression method that automatically detects and removes outliers from the data during the fitting procedure, by assigning them a very small weight. We discuss the use of both spline functions and spline curves. A number of numerical illustrations have been included to disclose the potentialities of the maximal-entropy approach in different application fields.
Solving multiphysics-based inverse problems for geological carbon storage monitoring can be challenging when multimodal time-lapse data are expensive to collect and costly to simulate numerically. We overcome these challenges by combining computationally cheap learned surrogates with learned constraints. Not only does this combination lead to vastly improved inversions for the important fluid-flow property, permeability, it also provides a natural platform for inverting multimodal data including well measurements and active-source time-lapse seismic data. By adding a learned constraint, we arrive at a computationally feasible inversion approach that remains accurate. This is accomplished by including a trained deep neural network, known as a normalizing flow, which forces the model iterates to remain in-distribution, thereby safeguarding the accuracy of trained Fourier neural operators that act as surrogates for the computationally expensive multiphase flow simulations involving partial differential equation solves. By means of carefully selected experiments, centered around the problem of geological carbon storage, we demonstrate the efficacy of the proposed constrained optimization method on two different data modalities, namely time-lapse well and time-lapse seismic data. While permeability inversions from both these two modalities have their pluses and minuses, their joint inversion benefits from either, yielding valuable superior permeability inversions and CO2 plume predictions near, and far away, from the monitoring wells.
Transition amplitudes and transition probabilities are relevant to many areas of physics simulation, including the calculation of response properties and correlation functions. These quantities can also be related to solving linear systems of equations. Here we present three related algorithms for calculating transition probabilities. First, we extend a previously published short-depth algorithm, allowing for the two input states to be non-orthogonal. Building on this first procedure, we then derive a higher-depth algorithm based on Trotterization and Richardson extrapolation that requires fewer circuit evaluations. Third, we introduce a tunable algorithm that allows for trading off circuit depth and measurement complexity, yielding an algorithm that can be tailored to specific hardware characteristics. Finally, we implement proof-of-principle numerics for models in physics and chemistry and for a subroutine in variational quantum linear solving (VQLS). The primary benefits of our approaches are that (a) arbitrary non-orthogonal states may now be used with small increases in quantum resources, (b) we (like another recently proposed method) entirely avoid subroutines such as the Hadamard test that may require three-qubit gates to be decomposed, and (c) in some cases fewer quantum circuit evaluations are required as compared to the previous state-of-the-art in NISQ algorithms for transition probabilities.
Iterative refinement (IR) is a popular scheme for solving a linear system of equations based on gradually improving the accuracy of an initial approximation. Originally developed to improve upon the accuracy of Gaussian elimination, interest in IR has been revived because of its suitability for execution on fast low-precision hardware such as analog devices and graphics processing units. IR generally converges when the error associated with the solution method is small, but is known to diverge when this error is large. We propose and analyze a novel enhancement to the IR algorithm by adding a line search optimization step that guarantees the algorithm will not diverge. Numerical experiments verify our theoretical results and illustrate the effectiveness of our proposed scheme.
The solutions of scalar ordinary differential equations become more complex as their coefficients increase in magnitude. As a consequence, when a standard solver is applied to such an equation, its running time grows with the magnitudes of the equation's coefficients. It is well known, however, that scalar ordinary differential equations with slowly-varying coefficients admit slowly-varying phase functions whose cost to represent via standard techniques is largely independent of the magnitude of the equation's coefficients. This observation is the basis of most methods for the asymptotic approximation of the solutions of ordinary differential equations, including the WKB method. Here, we introduce two numerical algorithms for constructing phase functions for scalar ordinary differential equations inspired by the classical Levin method for the calculation of oscillatory integrals. In the case of a large class of scalar ordinary differential equations with slowly-varying coefficients, their running times are independent of the magnitude of the equation's coefficients. The results of extensive numerical experiments demonstrating the properties of our algorithms are presented.