In this paper, we study numerical approximations for stochastic differential equations (SDEs) that use adaptive step sizes. In particular, we consider a general setting where decisions to reduce step sizes are allowed to depend on the future trajectory of the underlying Brownian motion. Since these adaptive step sizes may not be previsible, the standard mean squared error analysis cannot be directly applied to show that the numerical method converges to the solution of the SDE. Building upon the pioneering work of Gaines and Lyons, we shall instead use rough path theory to establish convergence for a wide class of adaptive numerical methods on general Stratonovich SDEs (with sufficiently smooth vector fields). To the author's knowledge, this is the first error analysis applicable to standard solvers, such as the Milstein and Heun methods, with non-previsible step sizes. In our analysis, we require the sequence of adaptive step sizes to be nested and the SDE solver to have unbiased "L\'evy area" terms in its Taylor expansion. We conjecture that for adaptive SDE solvers more generally, convergence is still possible provided the method does not introduce "L\'evy area bias". We present a simple example where the step size control can skip over previously considered times, resulting in the numerical method converging to an incorrect limit (i.e. not the Stratonovich SDE). Finally, we conclude with a numerical experiment demonstrating a newly introduced adaptive scheme and showing the potential improvements in accuracy when step sizes are allowed to be non-previsible.
The starting point of this paper is a collection of properties of an algorithm that have been distilled from the informal descriptions of what an algorithm is that are given in standard works from the mathematical and computer science literature. Based on that, the notion of a proto-algorithm is introduced. The thought is that algorithms are equivalence classes of proto-algorithms under some equivalence relation. Three equivalence relations are defined. Two of them give bounds between which an appropriate equivalence relation must lie. The third lies in between these two and is likely an appropriate equivalence relation. A sound method is presented to prove, using an imperative process algebra based on ACP, that this equivalence relation holds between two proto-algorithms.
In this second part of our two-part paper, we extend to multiple spatial dimensions the one-dimensional, fully conservative, positivity-preserving, and entropy-bounded discontinuous Galerkin scheme developed in the first part for the chemically reacting Euler equations. Our primary objective is to enable robust and accurate solutions to complex reacting-flow problems using the high-order discontinuous Galerkin method without requiring extremely high resolution. Variable thermodynamics and detailed chemistry are considered. Our multidimensional framework can be regarded as a further generalization of similar positivity-preserving and/or entropy-bounded discontinuous Galerkin schemes in the literature. In particular, the proposed formulation is compatible with curved elements of arbitrary shape, a variety of numerical flux functions, general quadrature rules with positive weights, and mixtures of thermally perfect gases. Preservation of pressure equilibrium between adjacent elements, especially crucial in simulations of multicomponent flows, is discussed. Complex detonation waves in two and three dimensions are accurately computed using high-order polynomials. Enforcement of an entropy bound, as opposed to solely the positivity property, is found to significantly improve stability. Mass, total energy, and atomic elements are shown to be discretely conserved.
In this paper we introduce a multilevel Picard approximation algorithm for general semilinear parabolic PDEs with gradient-dependent nonlinearities whose coefficient functions do not need to be constant. We also provide a full convergence and complexity analysis of our algorithm. To obtain our main results, we consider a particular stochastic fixed-point equation (SFPE) motivated by the Feynman-Kac representation and the Bismut-Elworthy-Li formula. We show that the PDE under consideration has a unique viscosity solution which coincides with the first component of the unique solution of the stochastic fixed-point equation. Moreover, the gradient of the unique viscosity solution of the PDE exists and coincides with the second component of the unique solution of the stochastic fixed-point equation. Furthermore, we also provide a numerical example in up to $300$ dimensions to demonstrate the practical applicability of our multilevel Picard algorithm.
In this study, we investigate an anisotropic weakly over-penalised symmetric interior penalty method for the Stokes equation {on convex domains}. Our approach is a simple discontinuous Galerkin method similar to the Crouzeix--Raviart finite element method. As our primary contribution, we show a new proof for the consistency term, which allows us to obtain an estimate of the anisotropic consistency error. The key idea of the proof is to apply the relation between the Raviart--Thomas finite element space and a discontinuous space. While inf-sup stable schemes of the discontinuous Galerkin method on shape-regular mesh partitions have been widely discussed, our results show that the Stokes element satisfies the inf-sup condition on anisotropic meshes. Furthermore, we also provide an error estimate in an energy norm on anisotropic meshes. In numerical experiments, we compare calculation results for standard and anisotropic mesh partitions.
In recent years, SPDEs have become a well-studied field in mathematics. With their increase in popularity, it becomes important to efficiently approximate their solutions. Thus, our goal is a contribution towards the development of efficient and practical time-stepping methods for SPDEs. Operator splitting schemes are a powerful tool for deterministic and stochastic differential equations. An example is given by domain decomposition schemes, where we split the domain into sub-domains. Instead of solving one expensive problem on the entire domain, we deal with cheaper problems on the sub-domains. This is particularly useful in modern computer architectures, as the sub-problems may often be solved in parallel. While splitting methods have already been used to study domain decomposition methods for deterministic PDEs, this is a new approach for SPDEs. We provide an abstract convergence analysis of a splitting scheme for stochastic evolution equations and state a domain decomposition scheme as an application of the setting. The theoretical results are verified through numerical experiments.
This paper addresses a factorization method for imaging the support of a wave-number-dependent source function from multi-frequency data measured at a finite pair of symmetric receivers in opposite directions. The source function is given by the inverse Fourier transform of a compactly supported time-dependent source whose initial moment or terminal moment for radiating is unknown. Using the multi-frequency far-field data at two opposite observation directions, we provide a computational criterion for characterizing the smallest strip containing the support and perpendicular to the directions. A new parameter is incorporated into the design of test functions for indicating the unknown moment. The data from a finite pair of opposite directions can be used to recover the $\Theta$-convex polygon of the support. Uniqueness in recovering the convex hull of the support is obtained as a by-product of our analysis using all observation directions. Similar results are also discussed with the multi-frequency near-field data from a finite pair of observation positions in three dimensions. We further comment on possible extensions to source functions with two disconnected supports. Extensive numerical tests in both two and three dimensions are implemented to show effectiveness and feasibility of the approach. The theoretical framework explored here should be seen as the frequency-domain analysis for inverse source problems in the time domain.
The hazard function represents one of the main quantities of interest in the analysis of survival data. We propose a general approach for parametrically 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 consider nonlinear solvers for the incompressible, steady (or at a fixed time step for unsteady) Navier-Stokes equations in the setting where partial measurement data of the solution is available. The measurement data is incorporated/assimilated into the solution through a nudging term addition to the the Picard iteration that penalized the difference between the coarse mesh interpolants of the true solution and solver solution, analogous to how continuous data assimilation (CDA) is implemented for time dependent PDEs. This was considered in the paper [Li et al. {\it CMAME} 2023], and we extend the methodology by improving the analysis to be in the $L^2$ norm instead of a weighted $H^1$ norm where the weight depended on the coarse mesh width, and to the case of noisy measurement data. For noisy measurement data, we prove that the CDA-Picard method is stable and convergent, up to the size of the noise. Numerical tests illustrate the results, and show that a very good strategy when using noisy data is to use CDA-Picard to generate an initial guess for the classical Newton iteration.
Covariance matrices of random vectors contain information that is crucial for modelling. Certain structures and patterns of the covariances (or correlations) may be used to justify parametric models, e.g., autoregressive models. Until now, there have been only few approaches for testing such covariance structures systematically and in a unified way. In the present paper, we propose such a unified testing procedure, and we will exemplify the approach with a large variety of covariance structure models. This includes common structures such as diagonal matrices, Toeplitz matrices, and compound symmetry but also the more involved autoregressive matrices. We propose hypothesis tests for these structures, and we use bootstrap techniques for better small-sample approximation. The structures of the proposed tests invite for adaptations to other covariance patterns by choosing the hypothesis matrix appropriately. We prove their correctness for large sample sizes. The proposed methods require only weak assumptions. With the help of a simulation study, we assess the small sample properties of the tests. We also analyze a real data set to illustrate the application of the procedure.
In this paper, we prove convergence for contractive time discretisation schemes for semi-linear stochastic evolution equations with irregular Lipschitz nonlinearities, initial values, and additive or multiplicative Gaussian noise on $2$-smooth Banach spaces $X$. The leading operator $A$ is assumed to generate a strongly continuous semigroup $S$ on $X$, and the focus is on non-parabolic problems. The main result concerns convergence of the uniform strong error $$E_{k}^{\infty} := \Big(\mathbb{E} \sup_{j\in \{0, \ldots, N_k\}} \|U(t_j) - U^j\|_X^p\Big)^{1/p} \to 0\quad (k \to 0),$$ where $p \in [2,\infty)$, $U$ is the mild solution, $U^j$ is obtained from a time discretisation scheme, $k$ is the step size, and $N_k = T/k$ for final time $T>0$. This generalises previous results to a larger class of admissible nonlinearities and noise, as well as rough initial data from the Hilbert space case to more general spaces. We present a proof based on a regularisation argument. Within this scope, we extend previous quantified convergence results for more regular nonlinearity and noise from Hilbert to $2$-smooth Banach spaces. The uniform strong error cannot be estimated in terms of the simpler pointwise strong error $$E_k := \bigg(\sup_{j\in \{0,\ldots,N_k\}}\mathbb{E} \|U(t_j) - U^{j}\|_X^p\bigg)^{1/p},$$ which most of the existing literature is concerned with. Our results are illustrated for a variant of the Schr\"odinger equation, for which previous convergence results were not applicable.