{'fraction_non_alphanumeric': 0.03861743352189212, 'fraction_numerical': 0.009442786512850207, 'mean_word_length': 4.545634185500803, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Probability Measures for Numerical Solutions of Differential Equations

In this paper, we present a formal quantification of epistemic uncertainty induced by numerical solutions of ordinary and partial differential equation models. Numerical solutions of differential equations contain inherent uncertainties due to the finite dimensional approximation of an unknown and implicitly defined function. When statistically analysing models based on differential equations describing physical, or other naturally occurring, phenomena, it is therefore important to explicitly account for the uncertainty introduced by the numerical method. This enables objective determination of its importance relative to other uncertainties, such as those caused by data contaminated with noise or model error induced by missing physical or inadequate descriptors. To this end we show that a wide variety of existing solvers can be randomised, inducing a probability measure over the solutions of such differential equations. These measures exhibit contraction to a Dirac measure around the true unknown solution, where the rates of convergence are consistent with the underlying deterministic numerical method. Ordinary differential equations and elliptic partial differential equations are used to illustrate the approach to quantifying uncertainty in both the statistical analysis of the forward and inverse problems.

Abstract: In this paper, we present a formal quantification of epistemic uncertainty induced by numerical solutions of ordinary and partial differential equation models. Numerical solutions of differential equations contain inherent uncertainties due to the finite dimensional approximation of an unknown and implicitly defined function. When statistically analysing models based on differential equations describing physical, or other naturally occurring, phenomena, it is therefore important to explicitly account for the uncertainty introduced by the numerical method. This enables objective determination of its importance relative to other uncertainties, such as those caused by data contaminated with noise or model error induced by missing physical or inadequate descriptors. To this end we show that a wide variety of existing solvers can be randomised, inducing a probability measure over the solutions of such differential equations. These measures exhibit contraction to a Dirac measure around the true unknown solution, where the rates of convergence are consistent with the underlying deterministic numerical method. Ordinary differential equations and elliptic partial differential equations are used to illustrate the approach to quantifying uncertainty in

Introduction and Motivation
The numerical analysis literature has developed a large range of efficient algorithms for solving ordinary and partial differential equations, which are typically designed to solve a single problem as efficiently as possible [12,8]. When classical numerical methods are placed within statistical analysis, however, we argue that significant difficulties can arise. Statistical models are commonly constructed by naively replacing the analytical solution with a numerical method; these models are then typically analysed by computational methods that sample the solver with many slightly varying configurations. This is problematic, not only because the numerical solutions are merely approximations, but also because the assumptions made by the numerical method induce correlations across multiple numerical solutions. While the distributions of interest commonly do converge asymptotically as the solver mesh becomes dense (e.g., in statistical inverse problems [6]), we argue that at a finite resolution the analysis may be vastly overconfident as a result of these unmodelled errors.
The purpose of this paper is to address these issues by the construction and rigorous analysis of novel probabilistic integration methods for both ordinary and partial differential equations. The approach in both cases is similar: we identify the key discretisation assumptions and introduce a local random field, in particular a Gaussian field, to reflect our uncertainty in those assumptions. The probabilistic solver may then be sampled repeatedly to interrogate the uncertainty in the solution. For a wide variety of commonly used numerical methods, our construction is straightforward to apply and preserves the formal order of convergence of the original method.
Furthermore, we demonstrate the value of these probabilistic solvers in statistical inference settings. Analytic and numerical examples show that using a classical non-probabilistic solver with inadequate discretisation when performing inference can lead to inappropriate and misleading posterior concentration in a Bayesian setting. In contrast, the probabilistic solver reveals the structure of uncertainty in the solution, naturally limiting posterior concentration as appropriate.
Our strategy fits within the emerging field known as Probabilistic Numerics [14], a perspective on computational methods pioneered by Diaconis [7], and subsequently Skilling [23]. Under this framework, solving differential equations is recast as the statistical inference problem of finding the latent, unobserved function that solves the differential equation, based on observable constraints. The output from such an inference problem is a probability measure over functions that satisfy the constraints imposed by the specific differential equation. This measure formally quantifies the uncertainty in candidate solution(s) of the differential equation and allows us to coherently propagate this uncertainty when carrying out statistical analyses, for example, in uncertainty quantification [25] or Bayesian inverse problems [6].

Review of existing work
Earlier work in the numerical analysis literature including randomisation in the approximate integration of ordinary differential equations (ODEs) includes [5,24]. However, the connections between probability and numerical analysis were elucidated in the foundational paper [7], which explains this connection in the context of Bayesian statistics and sets out an agenda for Probabilistic Numerics. It was in [23] that a concrete demonstration of how numerical integration of ordinary differential equations could be interpreted as a Bayesian inference method, and this idea was converted into an operational methodology, with consistency proofs given in [4].
The methodology detailed in [23,4] (explored in parallel in [13]) builds a probabilistic structure which contains two sources of randomness: an overarching Gaussian process structure and an internal randomisation of the mean process, resulting in different draws of the mean for each realisation. However, careful analysis reveals that, because of the conjugacy inherent in linear Gaussian models, the posterior variance of the outer Gaussian process is independent of the observed constraints, and thus does not reveal any information about the uncertainty in the particular problem being solved. The randomisation of the mean process is more useful, and is related to the strategy we present here. Those initial efforts were limited to first order Euler schemes [12], and were later extended to Runge-Kutta methods [22], although the incomplete characterisation of uncertainty remains. The existing literature has focussed only on ordinary differential equations and a formal probabilistic construction for partial differential equations (PDEs) that is generally applicable remains an open issue. All of these open issues are resolved by this work. Unlike [23,4,22], which are exclusively Bayesian by construction, our probabilistic solvers return a probability measure which can then be employed in either frequentist or Bayesian inference frameworks.
Our motivation for enhancing inference problems with models of discretisation error is similar to the more general concept of model error, as developed by [18]. However, our focus on errors arising from the discretisation of differential equations leads to more specialised methods. Existing strategies for discretisation error include empirically fitted Gaussian models for PDE errors [17] and randomly perturbed ODEs [2]; the latter partially coincides with our construction, but our motivation and analysis is distinct. Recent work [3] uses Bayes factors to analyse the impact of discretisation error on posterior approximation quality. Probabilistic models have also been used to study error propagation due to rounding error; see [11].
The remainder of the paper has the following structure: Section 2 introduces and formally analyses the proposed probabilistic solvers for ODEs. Section 3 explores the characteristics of random solvers employed in the statistical analysis of both forward and inverse problems. Then, we turn to elliptic PDEs in Section 4, where several key steps of the construction of probabilistic solvers and their analysis have intuitive analogues in the ODE context. Finally, an illustrative example of PDE inference problem is presented in Section 5. 1

Probability Measures via Probabilistic Time Integrators
Consider the following ordinary differential equation (ODE): where u(·) is a continuous function taking values in R n . 2 We let Φ t denote the flow-map for Equation (2.1), so that u(t) = Φ t u(0) . The conditions ensuring that this solution exists will be formalised in Assumption 2, below. Deterministic numerical methods for the integration of this equation on timeinterval [0, T ] will produce an approximation to the equation on a mesh of points {t k = kh} K k=0 , with Kh = T , (for simplicity we assume a fixed mesh). Let u k = u(t k ) denote the exact solution of (2.1) on the mesh and U k ≈ u k denote the approximation computed by the numerical method, based on evaluating f (and possibly higher derivatives) at a finite set of points that are generated during numerical integration. Typically these methods output a single discrete solution {U k } K k=0 , possibly augmented with some type of error indicator, but do not statistically quantify the uncertainty remaining in the path.
Let X a,b denote the Banach space C([a, b]; R n ). The exact solution of (2.1) on the time-interval [0, T ] may be viewed as a Dirac measure δ u on X 0,T at the element u that solves the ODE. We will construct a probability measure µ h on X 0,T , that is straightforward to sample from both on and off the mesh, for which h quantifies the size of the discretisation step employed, and whose distribution reflects the uncertainty resulting from the solution of the ODE. Convergence of the numerical method is then related to the contraction of µ h to δ u .
We briefly summarise the construction of the numerical method. Let Ψ h : R n → R n denote a classical deterministic one-step numerical integrator over time-step h, a class including all Runge-Kutta methods and Taylor methods for ODE numerical integration [12]. Our numerical methods will have the property that, on the mesh, they take the form where ξ k (h) are suitably scaled, i.i.d. Gaussian random variables. That is, the random solution iteratively takes the standard step, Ψ h , followed by perturbation with a random draw, ξ k (h), modelling uncertainty that accumulates between mesh points. The discrete path {U k } K k=0 is straightforward to sample and in general is not a Gaussian process. Furthermore, the discrete trajectory can be extended into a continuous time approximation of the ODE, which we define as a draw from the measure µ h . The remainder of this section develops these solvers in detail and proves strong convergence of the random solutions to the exact solution, implying that µ h → δ u in an appropriate sense. Finally, we establish a close relationship between our random solver and a stochastic differential equation (SDE) with small mesh-dependent noise.

Probabilistic time-integrators: general formulation
The integral form of Equation (2.1) is (2. 3) The solutions on the mesh satisfy and may be interpolated between mesh points by means of the expression We may then write where g(s) = f u(s) is an unknown function of time. In the algorithmic setting we have approximate knowledge about g(s) through an underlying numerical method. A variety of traditional numerical algorithms may be derived based on approximation of g(s) by various simple deterministic functions g h (s). Perhaps the simplest such numerical method arises from invoking the Euler approximation that In particular, if we take t = t k+1 and apply this method inductively the corresponding numerical scheme arising from making such an approximation to g(s) in (2.6) is U k+1 = U k + hf (U k ). Now consider the more general one-step numerical method U k+1 = Ψ h (U k ). This may be derived by approximating g(s) in (2.6) by We note that all consistent (in the sense of numerical analysis) one-step methods will satisfy d dτ The approach based on the approximation (2.8) leads to a deterministic numerical method which is defined as a continuous function of time. Specifically we have U (s) = Ψ s−t k (U k ), s ∈ [t k , t k+1 ). Consider again the Euler approximation, for which Ψ τ (U ) = U + τ f (U ), and whose continuous time interpolant is then given by U (s) = U k + (s − t k )f (U k ), s ∈ [t k , t k+1 ). Note that this produces a continuous function, namely an element of X 0,T , when extended to s ∈ [0, T ]. The preceding development of a numerical integrator does not acknowledge the uncertainty that arises from lack of knowledge about g(s) in the interval s ∈ [t k , t k+1 ). We propose to approximate g stochastically in order to represent this uncertainty, taking where the {χ k } form an i.i.d. sequence of Gaussian random functions defined We will choose C h to shrink to zero with h at a prescribed rate, and also to ensure that χ k ∈ X 0,h almost surely. The functions {χ k } represent our uncertainty about the function g. The corresponding numerical scheme arising from such an approximation is given by where the i.i.d. sequence of functions {ξ k } lies in X 0,h and is given by Note that the numerical solution is now naturally defined between grid points, via the expression In the case of the Euler method, for example, we have (2.12) and, between grid points, This method is illustrated in Figure 1.  The random integrator in (b) outputs the path in red; we overlay the standard Euler step constructed at each step, before it is perturbed (blue).
While we argue that the choice of modelling local uncertainty in the flow-map as a Gaussian process is natural and analytically favourable, it is not unique. It is possible to construct examples where the Gaussian assumption is invalid; for example, when a highly inadequate time step is used, a systemic bias may be introduced. However, in regimes where the underlying deterministic method performs well, the centred Gaussian assumption is a reasonable prior.
Finally, we comment on sampling off-grid. Imagine that we wish to sample the continuous time solution at point s ∈ (t k , t k+1 ). This might be performed in two different ways, given the solution generated on-grid, up to time t k , producing sample U k . In the first method, we draw the sample U (s) = Ψ s−t k (U k ) + ξ k (s − t k ). Having done so we must be careful to draw . This is not complicated, because all random variables are Gaussian given available information, but needs to be respected via conditioning on the end points of the grid. An alternative is to draw the sample , which the conditional Gaussian structure facilitates.

Strong convergence result
To prove the strong convergence of our probabilistic numerical solver, we first need two assumptions quantifying properties of the random noise and of the underlying deterministic integrator, respectively. In what follows we use ·, · and | · | to denote the Euclidean inner-product and norm on R n . We denote the Frobenius norm on R n×n by | · | F , and E h denotes expectation with respect to the i.i.d. sequence {χ k }.
Here, and in the sequel, K is a constant independent of h, but possibly changing from line to line. The covariance kernel C h is not uniquely defined, but the uncertainty at the beginning of the step must be zero and increase with step length, scaling exactly as h 2p+1 at the end of the interval. Implementation requires selection of the constant matrix Q; in our numerical examples we assume a scaled identity matrix, Q = σI, and we discuss one possible strategy for choosing σ in Section 3.1.
Assumption 2. The function f and a sufficient number of its derivatives are bounded uniformly in R n in order to ensure that f is globally Lipschitz and that the numerical flow-map Ψ h has uniform local truncation error of order q + 1: Remark 2.1. We assume globally Lipschitz f , and bounded derivatives, in order to highlight the key probabilistic ideas, whilst simplifying the numerical analysis. Future work will address the non-trivial issue of extending of analyses to weaken these assumptions. In this paper, we provide numerical results indicating that a weakening of the assumptions is indeed possible.
This theorem implies that every probabilistic solution is a good approximation of the exact solution in both a discrete and continuous sense. Choosing p ≥ q is natural if we want to preserve the strong order of accuracy of the underlying deterministic integrator; we proceed with the choice p = q, introducing the maximum amount of noise consistent with this constraint.
Proof. We first derive the convergence result on the grid, and then in continuous time. From (2.9) we have Subtracting equation (2.16) from (2.15) and defining e k = u k − U k , we get Taking the Euclidean norm and taking expectations gives, using Assumption 1 and the independence of the ξ k , Noting that Φ h is globally Lipschitz with constant bounded by 1 + Lh under Assumption 2, we then obtain Using Cauchy-Schwarz on the inner-product, and the fact that Φ h is Lipschitz with constant bounded independently of h, we get Application of the Gronwall inequality gives the desired result. Now we turn to continuous time. We note that, for s ∈ [t k , t k+1 ), Let F t denote the σ-algebra of events generated by the {ξ k } up to time t.
Subtracting we obtain, using Assumptions 1 and 2 and the fact that Φ s−t k has Lipschitz constant of the form 1 + O(h),

Now taking expectations we obtain
Using the on-grid error bound gives the desired result, after noting that the constants appearing are uniform in 0 ≤ kh ≤ T.

Examples of probabilistic time-integrators
The canonical illustration of a probabilistic time-integrator is the probabilistic Euler method already described. In this section we describe two methods which showcase the generality of the approach. The first is the classical Runge-Kutta method which defines a one-step numerical integrator as follows: The method has local truncation error in the form of Assumption 2 with q = 4. It may be used as the basis of a probabilistic numerical method (2.11), and hence (2.9) at the grid-points. Thus, provided that we choose to perturb this integrator with a random process χ k satisfying Assumption 1 with p ≥ 4, Theorem 2.2 shows that the error between the probabilistic integrator based on the classical Runge-Kutta method is, in the mean square sense, of the same order of accuracy as the deterministic classical Runge-Kutta integrator. The second is an integrated Ornstein-Uhlenbeck process, derived as follows. Define, on the interval s ∈ [t k , t k+1 ), the pair of equations Here W is a standard Brownian motion and Λ and Σ are invertible matrices, possibly depending on h. The approximating function g h (s) is thus defined by V (s), an Ornstein-Uhlenbeck process. Integrating (2.18b) we obtain Note that the h-dependence of C h comes through the time-interval on which χ k is defined, and through Λ and Σ. Integrating (2.18a), using (2.19), we obtain The numerical method (2.20) may be written in the form (2.11), and hence (2.9) at the grid-points, with the definition This integrator is first order accurate and satisfies Assumptions 2 with p = 1.
Choosing to scale Σ with h so that q ≥ 1 in Assumptions 1 leads to convergence of the numerical method with order 1. Had we carried out the above analysis in the case Λ = 0 we would have obtained the probabilistic Euler method (2.13), and hence (2.12) at grid points, used as our canonical example in the earlier developments.

Backward error analysis
The idea of backward error analysis is to identify a modified equation which is solved by the numerical method either exactly, or at least to a higher degree of accuracy than the numerical method solves the original equation. In the context of differential equations this modified equation will involve the step-size h. In the setting of ordinary differential equations, and the random integrators introduced in this section, we will show that the modified equation is a stochastic differential equation (SDE) in which only the matrix Q from Assumption 1 enters; the details of the random processes used in our construction do not enter the modified equation. This universality property underpins the methodology we introduce as it shows that many different choices of random processes all lead to the same effective behaviour of the numerical method.
We introduce the operators L and L h defined so that, for all φ ∈ C ∞ (R n , R), Thus L := f · ∇ and e hL h is the kernel for the Markov chain generated by the probabilistic integrator (2.2). In fact we never need to work with L h itself in what follows, only with e hL h , so that questions involving the operator logarithm do not need to be discussed. We now introduce a modified ODE and a modified SDE which will be needed in the analysis that follows. The modified ODE is whilst the modified SDE has the form The precise choice of f h is detailed below. Letting E denote expectation with respect to W , we introduce the operators L h and L h so that, for all φ ∈ C ∞ (R n , R), where : denotes the inner product on R n×n which induces the Frobenius norm, that is, The fact that the deterministic numerical integrator has uniform local truncation error of order q + 1 (Assumption 2) implies that, since φ ∈ C ∞ , (2.27) The theory of modified equations for classical one-step numerical integration schemes for ODEs [12] establishes that it is possible to find f h in the form We work with this choice of f h in what follows. Now for our stochastic numerical method we have Furthermore the last term has mean of size O(|ξ k (h)| 4 ). From Assumption 1 we From this it follows that Finally we note that (2.26) implies that Thus we have Balancing these terms, in what follows we make the choice l = 2p − q. If l < 0 we adopt the convention that the drift f h is simply f. With this choice of q we obtain In this case, the second part of Theorem 2.2 gives us the first weak order result in Theorem 2.4. Additionally, Assumption 3 is stronger than we need, but allows us to highlight probabilistic ideas whilst keeping overly technical aspects of the numerical analysis to a minimum. More sophisticated, but structurally similar, analysis would be required for weaker assumptions on f . Similar considerations apply to the assumptions on φ.
Theorem 2.4. Consider the numerical method (2.9) and assume that Assumptions 1 and 3 are satisfied. Then, for φ ∈ C ∞ function with all derivatives bounded uniformly on R n , we have that where u andũ solve (2.1) and (2.24) respectively.
Proof. We prove the second bound first. Let w k = E φ(ũ(t k ))|ũ(0) = u and W k = E h φ(U k )|U 0 = u). Then let δ k = sup u∈R n |W k − w k |. It follows from the Markov property that Using (2.34) and Assumption 3 we obtain Iterating and employing the Gronwall inequality gives the second error bound. Now we turn to the first error bound, comparing with the solution u of the original equation (2.1). From (2.31) and then (2.27) we see that This gives the first weak error estimate, after using the stability estimate on e hL from Assumption 3.
Example 2.5. Consider the probabilistic integrator derived from the Euler method in dimension n = 1. We thus have q = 1, and we hence set p = 1.
The results in [10] allow us to calculate f h with l = 1. The preceding theory then leads to strong order of convergence 1, measured relative to the true ODE (2.1), and weak order 3 relative to the SDE With these results now available, the following section provides an empirical study of our probabilistic integrators.

Statistical Inference and Numerics
This section explores applications of the randomised ODE solvers developed in Section 2. First, we study forward uncertainty propagation and propose a method for calibrating the problem-dependent scaling constant, σ, against classic error indicators. 4 Secondly, we employ the calibrated measure within Bayesian inference problems, demonstrating that the resulting posteriors exhibit more consistent behaviour over varying step-sizes than with naive use of a deterministic integrators. Throughout this section we use the FitzHugh-Nagumo model to illustrate ideas [21]. This is a two-state non-linear oscillator, with states (V, R) and parameters (a, b, c), governed by the equations This particular example does not satisfy the stringent Assumptions 2 and 3 and the numerical results shown demonstrate that, as indicated in Remarks 2.1 and 2.3, our theory will extend to weaker assumptions on f , something we will address in future work.

Calibrating forward uncertainty propagation
Consider Equation (3.1) with fixed initial conditions V (0) = −1, R(0) = 1, and parameter values (.2, .2, 3). Figure 2 shows draws of the V species trajectories from the measure associated with the probabilistic Euler solver with p = q = 1, for various values of the step-size and fixed σ = 0.1. The contraction of the measure towards the reference solution, as h shrinks, is clearly evident. Furthermore, the uncertainty exhibits interesting, non-Gaussian structure where trajectories disagree whether to begin the steep phase of the oscillation.
Although the rate of contraction is governed by the underlying deterministic method, the scale parameter, σ, completely controls the apparent uncertainty in the solver. To illustrate this, Figure 3 fixes the step-size, h = 0.1, and shows that rescaling the noise can create any apparent level of certainty desired, including high confidence in an incorrect solution. This tuning problem exists in general, since σ is problem dependent and cannot obviously be computed analytically. Therefore, we propose to calibrate σ to replicate the amount of error suggested by classical error indicators. In the following discussion, we often explicitly denote the dependence on h and σ with superscripts, hence the probabilistic solver is U h,σ and the corresponding deterministic solver is U h,0 . Define the true error as e(t) = u(t) − U h,0 (t). Then we assume there is some computable error indicator E(t) ≈ e(t), defining E k = E(t k ). The simplest error indicators might compare differing step-sizes, E(t) = U h,0 (t) − U 2h,0 (t), or differing order methods, as in a Runge-Kutta 4-5 scheme.
We proceed by constructing a probability distribution π(σ) that is maximised when the overall amount of error produced by our probabilistic integrator matches that suggested by the error indicator. We perform this scale matching by: (i) using a Gaussian approximation of our random solver at each step k,μ h,σ k = N (E(U h,σ k ), V(U h,σ k )); and (ii) by constructing a second Gaussian measure from the deterministic solver, U h,0 k , and the available error indicator, . Thus (ii) is a natural Gaussian representation of the uncertainty based upon information from deterministic methods. We construct π(σ) by penalising the distance between these two normal distributions at every step: π(σ) ∝ k exp −d(μ h,σ k , ν σ k ) . We find that the Bhattacharyya distance (closely related to the Hellinger metric) works well [16], since it diverges quickly if either the mean or variance differs between the two inputs. Computing π(σ) requires the mean and variance of the process U h,σ k , which can be estimated using Monte Carlo. If the ODE state is a vector, we take the product of the univariate Bhattacharyya distances. Note that this calibration depends not only on the problem of interest, but also on the initial conditions and any parameters of the ODE.
Returning to the FitzHugh-Nagumo model, we apply this strategy for choosing the noise scale. We sampled from π(σ) using pseudomarginal MCMC [1]. The sampled densities are shown in Figure 4; the similarity in the inferred σ suggests that the asymptotic convergence rate is nearly exact in this case. Since this distribution is well peaked, MCMC mixes very quickly and is well represented by the MAP value, hence we proceed using σ * = arg max π(σ). Next, we examine the quality of the scale matching by plotting the magnitudes of the random variation against the error indicator in Figure 5, for several different step-sizes, observing good agreement of the marginal variances. Although the marginal variances are empirically matched to the error indicator, our measure is still reveals non-Gaussian structure and correlations in time not revealed by the deterministic analysis. As described, this procedure requires fixed inputs to the ODE, but it is straightforward to marginalise out a prior distribution over input parameters.

Constructing Bayesian posterior inference problems
Given the calibrated probabilistic ODE solvers described above, let us consider how to incorporate them into inference problems. Assume we are interested in inferring parameters of the ODE given noisy observations of the state. Specifically, we wish to infer parameters θ ∈ R d for the differential equationu = f (u, θ), with fixed initial conditions u(t = 0) = u 0 (a straightforward modification may include inference on initial conditions). Assume we are provided with data d ∈ R m , d j = u(τ j ) + η j at some collection of times τ j , corrupted by i.i.d. noise, η j ∼ N (0, Γ). If we have prior Q(θ), the posterior we wish to explore is, P(θ | d) ∝ Q(θ)L(d, u(θ)), where density L compactly summarises this likelihood model.  The red curves depict fifty samples of the magnitude of the difference between a standard Euler solver for several step-sizes and the equivalent randomised variant, using σ * , maximising π(σ).
The standard computational strategy is to simply replace the unavailable trajectory u with a numerical approximation, inducing approximate posterior P h,0 (θ | d) ∝ Q(θ)L(d, U h,0 (θ)). Informally, this approximation will be accurate when the error in the numerical solver is small compared to Γ and often converges formally to P(θ | d) as h → 0 [6]. However, the error in U h,0 might be non-Gaussian and highly correlated, which could make it surprisingly easy for errors at finite h to have substantial impact.
In this work, we are concerned about the undue optimism in the predicted variance; without a rigorous treatment of the solver uncertainty, the posterior can concentrate around an arbitrary parameter value even when the deterministic solver is inaccurate, and is merely able to reproduce the data by coincidence. The more conventional concern is that any error in the solver will be transferred into posterior bias. Practitioners commonly alleviate both concerns by tuning solvers until they are nearly perfect, however, we note that this may be computationally prohibitive in many contemporary statistical applications.
We can construct a different posterior that includes the uncertainty in the solver by taking an expectation over random solutions to the ODE Intuitively, this construction favours parameters that exhibit agreement with the entire family of uncertain trajectories, not just the one produced by the deterministic integrator. The typical effect of this expectation is to increase the posterior uncertainty on θ, preventing the inappropriate posterior collapse we are concerned about. Indeed, if the integrator cannot resolve the underlying dynamics, h p+1/2 σ will be large. Then U h,σ (θ, ξ) is independent of θ, hence the prior is recovered, Notice that as h → 0, both the measures P h,0 and P h,σ collapse to the posterior using the analytic solution, P, hence both methods are correct. We do not expect the bias of P h,σ to be improved, since all of the averaged trajectories are of the same quality as the deterministic solver in P h,0 . We now construct an analytic inference problem demonstrating these behaviours.
Example 3.1. Consider inferring the initial condition, u 0 ∈ R, of the scalar linear differential equation,u = λu, with λ > 0. We apply a numerical method to produce the approximation U k ≈ u(kh). We observe the state at some times t = kh, with additive noise η k ∼ N (0, γ 2 ): d k = U k + η k . If we use a deterministic Euler solver, the model predicts U k = (1 + hλ) k u 0 . These model predictions coincide with the slightly perturbed problem hence error increases with time. However, the assumed observational model does not allow for this, as the observation variance is γ 2 at all times. In contrast, our proposed probabilistic Euler solver predicts where we have made the natural choice p = q, where σ is the problem dependent scaling of the noise and the ξ k are i.i. d. N (0, 1). For a single observation, η k and every ξ k are independent, so we may rearrange the equation to consider the perturbation as part of the observation operator. Hence, a single observation at k has effective variance Thus, late time observations are modelled as being increasingly inaccurate. Consider inferring u 0 , given a single observation d k at time k. If a Gaussian prior N (m 0 , ζ 2 0 ) is specified for u 0 , then the posterior is N (m, ζ 2 ), where The observation precision is scaled by (1 + hλ) 2k because late time data contain increasing information. Assume that the data are d k = e λkh u † 0 + γη † , for some given true initial condition u † 0 and noise realisation η † . Consider now the asymptotic regime where h is fixed and k → ∞. For the standard Euler method, where γ h = γ, we see that ζ 2 → 0, whilst m (1 + hλ) −1 e hλ k u † 0 . Thus the inference scheme becomes increasingly certain of the wrong answer: the variance tends to zero and the mean tends to infinity.
In contrast, with a randomised integrator, the fixed h, large k asymptotics are Thus, the mean blows up at a modified rate, but the variance remains positive.
We take an empirical Bayes approach to choosing σ, that is, using a constant, fixed value σ * = arg max π(σ), chosen before the data is observed. Joint inference of the parameters and the noise scale suffer from well-known MCMC mixing issues in Bayesian hierarchic models.
We now return to the FitzHugh-Nagumo model; given fixed initial conditions, we attempt to recover parameters θ = (a, b, c) from observations of both species at times τ = 1, 2, . . . , 40. The priors are log-normal, centred on the true value with unit variance, and with observational noise Γ = 0.001. The data is generated from a high quality solution, and we perform inference using Euler integrators with various step-sizes, h ∈ {0.005, 0.01, 0.02, 0.05, 0.1}, spanning a range of accurate and inaccurate integrators.
We first perform the inferences with naive use of deterministic Euler integrators. We simulate from each posterior using delayed rejection MCMC [9] for 100,000 steps, discarding the first 10% as burn in. The posteriors are shown in Figure 6. Observe the undesirable concentration of every posterior, even those with poor solvers; the posteriors are almost mutually singular, hence clearly the posterior widths are meaningless.
Secondly, we repeat the experiment using our probabilistic Euler integrators, with results shown in Figure 7. We use a noisy pseudomarginal MCMC method, whose fast mixing is helpful for these initial experiments [20]. These posteriors are significantly improved, exhibiting greater mutual agreement and obvious increasing concentration with improving solver quality. The posteriors are not perfectly nested, possible evidence that our choice of scale parameter is imperfect, or that the assumption of locally Gaussian error deteriorates. As expected, the bias of these posteriors is essentially unchanged.

Probabilistic Solvers for Partial Differential Equations
We now turn to present a framework for probabilistic solutions to partial differential equations, working within the finite element setting. Our approach to PDE problems is analogous to that for ODEs, except that now we randomly perturb the finite element basis functions. The assumptions required, and the strong convergence results, closely resemble their ODE counterparts.

Probabilistic finite element method for variational problems
Let V be a Hilbert space of real valued functions defined on a bounded polygonal domain D ⊂ R d . Consider a weak formulation of a linear PDE specified via a symmetric bilinear form a : V × V −→ R, and a linear form r : V −→ R to give the problem of finding u ∈ V : a(u, v) = r(v), ∀v ∈ V. This problem can be approximated by specifying a finite dimensional subspace V h ⊂ V and seeking a solution in V h instead. This leads to a finite dimensional problem to be solved for the approximation U : This is known as the Galerkin method. We will work in the setting of finite element methods, assuming that V h = span{φ j } J j=1 where φ j is locally supported on a grid of points {x j } J j=1 . The parameter h is introduced to measure the diameter of the finite elements. We will also assume that Any element U ∈ V h can then be written as from which it follows that U (x k ) = U k . The Galerkin method then gives the following equation for U = (U 1 , . . . , U J ) T : where A jk = a(φ j , φ k ) and r k = r(φ k ). In order to account for uncertainty introduced by the numerical method, we will assume that each basis function φ j can be split into the sum of a systematic part φ s j and random part φ r j , where both φ j and φ s j satisfy the nodal property (4.2), hence φ r j (x k ) = 0. Furthermore we assume that each φ r j shares the same compact support as the corresponding φ s j , preserving the sparsity structure of the underlying deterministic method.

Strong convergence result
As in the ODE case, we begin our convergence analysis with assumptions constraining the random perturbations and the underlying deterministic approximation. The bilinear form a(·, ·) is assumed to induce an inner-product, and then norm via · 2 a = a(·, ·); furthermore we assume that this norm is equivalent to the norm on V. Throughout E h denotes expectation with respect to the random basis functions.
a ≤ Ch 2 min{p,q} . As for ODEs, the solver accuracy is limited by either the amount of noise injected or the convergence rate of the underlying deterministic method; the choice p = q is natural since it introduces the maximum amount of noise that does not disturb the deterministic rate of convergence.
Proof. Recall the Galerkin orthogonality property which follows from subtracting the approximate variational principle from the true variational principle: it states that, for e = u − U , a(e, v) = 0, ∀v ∈ V h . (4.5) From this it follows that To see this note that, for any v ∈ V h , the orthogonality property (4.5) gives a(e, e) = a(e, e + U − v) = a(e, u − v). Thus, by Cauchy-Schwarz, e 2 a ≤ e a u − v a , ∀v ∈ V h implying (4.6). We now set, for v ∈ V, By the mean-zero and independence properties of the random basis functions we deduce that The result follows from Assumptions 4 and 5.

Poisson solver in two dimensions
Consider a Poisson equation with Dirichlet boundary conditions in dimension d = 2, namely We set V = H 1 0 (D) and H to be the space L 2 (D) with inner-product ·, · and resulting norm | · | 2 = ·, · . The weak formulation of the problem has the form (4.1) with Now consider piecewise linear finite elements satisfying the assumptions of Section 4.2 in [15] and take these to comprise the set {φ s j } J j=1 . Then h measures the width of the triangulation of the finite element mesh. Assuming that f ∈ H it follows that u ∈ H 2 (D) and that Thus q = 1. We choose random basis members {φ r j } J j=1 so that Assumptions 4 hold with p = 1. Theorem 4.1 then shows that, for e = u − U , E h e 2 a ≤ Ch 2 . Note, however, that in the deterministic case we expect an improved rate of convergence in the function space H. We show that such a result also holds in our random setting.
Note that, under Assumptions 4, if we introduce γ jk that is 1 when two basis functions have overlapping support, and 0 otherwise, then γ jk is symmetric and there is constant C, independent of j and J, such that J k=1 γ jk ≤ C. Now let ϕ solve the equation a(ϕ, v) = e, v , ∀v ∈ V. Then ϕ H 2 ≤ C|e|. We define ϕ s and ϕ r in analogy with the definitions of v s and v r . Following the usual arguments for application of the Aubin-Nitsche trick [15], we have |e| 2 = a(e, ϕ) = a(e, ϕ − ϕ s − ϕ r ). Thus where, by Sobolev embedding (d = 2 here), a j := ϕ(x j )/ ϕ H 2 satisfies max 1≤j≤J |a j | ≤ C. Note, however, that the a j are random and correlated with all of the random basis functions. Using this, together with (4.8), in (4.9), we obtain We see that |e| ≤ C e a h 2 + J j=1 J k=1 a j a k a(φ r j , φ r k ) Taking expectations, using that p = q = 1, we find, using Assumptions 4, that E h |e| ≤ Ch E h e 2 a 1 2 ≤ Ch 2 as desired. Thus we recover the extra order of convergence over the rate 1 in the · a norm (although the improved rate is in L 1 (Ω; H) whilst the lower rate of convergence is in L 2 (Ω; V).)

PDE Inference and Numerics
We now perform numerical experiments using probabilistic solvers for elliptic PDEs. Specifically, we perform inference in a 1D elliptic PDE, ∇·(κ(x)∇u(x)) = 4x for x ∈ [0, 1], given boundary conditions u(0) = 0, u(1) = 2. We represent log κ as piecewise constant over ten equal-sized intervals; the first, on x ∈ [0, .1) is fixed to be one to avoid non-identifiability issues, and the other nine are given a prior θ i = log κ i ∼ N (0, 1). Observations of the field u are provided at x = (0.1, 0.2, . . . 0.9), with i.i.d. Gaussian error, N (0, 10 −5 ); the simulated observations were generated using a fine grid and quadratic finite elements, then perturbed with error from this distribution.
Again we investigate the posterior produced at various grid sizes, using both deterministic and randomised solvers. The randomised basis functions are draws from a Brownian bridge conditioned to be zero at the nodal points, implemented in practice with a truncated Karhunen-Loève expansion. The covariance operator may be viewed as a fractional Laplacian, as discussed in [19]. The scaling σ is again determined by maximising the distribution described in Section 3.1, where the error indicator compares linear to quadratic basis functions, and we marginalise out the prior over the κ i values.
The posteriors are depicted in Figures 8 and 9. As in the ODE examples, the deterministic solvers lead to incompatible posteriors for varying grid sizes. In contrast, the randomised solvers suggest increasing confidence as the grid is refined, as desired. The coarsest grid size uses an obviously inadequate ten elements, but this is only apparent in the randomised posterior.

Conclusions
We have presented theory and methods for probabilistic approaches to the numerical solution of both ordinary and partial differential equations. These methods give rise to a probability measure whose qualitative and quantitative properties account for and quantify the uncertainty induced by finite dimensional approximation. We provide a theoretical analysis of the properties of these probabilistic integrators and demonstrate that they induce more coherent inferences in illustrative empirical examples. This work opens the door for statistical analysis that explicitly incorporates numerical uncertainty in many important classes of contemporary statistical problems across the sciences, including engineering, climatology, and biology. Drawing parallels to model error [18], we can consider our intrusive modifications to be a highly-specialised discrepancy model, designed using our intimate knowledge of the structure and properties of numerical methods. Further study is required to compare our approach to existing methods, and we hope to develop other settings where modifying the internal structure of numerical methods can produce computationally and analytically tractable measures of uncertainty. Developing robust, practical solvers, and efficiently performing computations with the additional uncertainty we propose remains a challenge. Finally, future work may be able to exploit our coherent uncertainty propagation to optimise the tradeoff between solver cost and statistical variance.