Error control in the numerical posterior distribution in the Bayesian UQ analysis of a semilinear evolution PDE

Error con trol in the n umerical p osterior distribution in the Ba y esian UQ analysis of a semilinear ev olution PDE Maria L. Daza-T orres ∗ ‡ J. Cricelio Mon tesinos-Lóp ez ∗ Marcos A. Capistrán ∗ J. Andrés Christen ∗ Heikki Haario † No vem ber 9, 2020 W e elab orate on results obtained in [ 1 ] for con trolling the n umerical posterior error for Ba yesian UQ problems, now considering forward maps arising from the solution of a semilinear ev olution partial differen tial equation. Results in [ 1 ] demand an estimate for the absolute global error (A GE) of the numeric forw ard map. Our con tribution is a numerical method for computing the A GE for semilinear ev olution PDEs and sho ws the potential applicability of [ 1 ] in this imp ortan t wide range family of PDEs. Numerical examples are giv en to illustrate the efficiency of the prop osed metho d, obtaining numerical posterior distributions for unknown parameters that are nearly iden tical to the corresp onding theoretical posterior, by k eeping their Bay es factor close to 1. 1 In tro duction A wide range of applications are concerned with the solution of an inv erse problem (IP) [ 2 , 3 , 4 , 5 , 6 , 7 , 8 ]: giv en some observ ations of the output, y = ( y 1 , . . . , y n ) , to determine the corresp onding inputs θ such that y i = F ( θ ) + error . W e refer to the ev aluation of F as solving the forward problem, and consequently , F is called the F orward Map (FM). In general, the FM is a complex non-linear map, with input parameters θ , defined b y an initial/b oundary v alue problem for a system of ordinary differential equations (ODEs) or partial differential equations (PDEs). Then, to ev aluate F ( θ ) , w e must solve an initial/b oundary v alue problem for a system of (O, P)DEs. IPs are typically ill-p osed: there may b e no solution, or the solution ma y not b e unique and ma y dep end sensitiv ely on y i [ 9 ]. A wa y to approac h these difficulties is to formulate the IP in ∗ Centro de Inv estigación en Matemáticas (CIMA T), Jalisco S/N, V alenciana, Guana juato, 36023, México. mdazatorr es, jose.montesinos, mar c os, jac at cimat.mx † Lappeenranta Universit y of T ec hnology , Department of computational and pro cess engineering, Lapp een- ranta, Finland and Finnish Meteorological Institute, Helsinki, Finland heikki.haario@lut.fi ‡ Corresponding author 1 the Bay esian framew ork. Stuart [ 10 ] studied conditions for the w ell-p osedness of the Bay esian form ulation of IPs. In this sc heme, a noise mo del is assumed for the observ ations, e.g., y i = F ( θ ) + ε i ; ε i ∼ N (0 , σ 2 ) . This observ ational mo del generates a probabilit y density giv en the parameter Φ = ( θ, σ ) , namely P Y | Φ ( y | θ, σ ) , for fixed data y , obtaining the lik eliho od function. Based on the av ailable infor- mation, a prior mo del P Φ ( · ) is stated for Φ , and a posterior distribution is obtained, P Φ | Y ( θ , σ | y ) = P Y | Φ ( y | θ , σ ) P Φ ( θ , σ ) P Y ( y ) . Explicit analytic forms are usually not av ailable for the posterior distributions, so sampling approac hes such as the Mon te Carlo Marko v Chain (MCMC) are required to c haracterize it. These metho ds in volv e repeated FM solutions used to define the likelihoo d function. Usually , we do not ha ve an analytical or computationally precise and straightforw ard imple- men tation of the FM. This necessarily inv olv es a n umerical appro ximation, F α ( n ) , where α ( n ) represen ts a discretization used to appro ximate the FM, leading to a numerical/appro ximate p osterior distribution. Th us, the numerical solution of the FM will in tro duce some numerical error in the p osterior distribution. At least theoretically , numerical errors in the FM can b e con trolled and reduced to an arbitrarily lo w lev el, through the use of finer discretizations, but what numerical error must b e tolerated in the FM to obtain a correct and acceptable n umerical p osterior distribution? Sev eral approaches start by building cheap computationally approximations of the FM and using these approximations as surrogates in the sampling pro cedure [ 11 , 12 , 13 , 14 ]. Although suc h approac hes can b e quite effective at reducing computation cost, there has been little analysis of posterior inference appro ximation. Recen tly , adaptiv e m ulti-fidelity techniques hav e been dev elop ed to con trol the n umerical posterior error for Ba yesian UQ [ 15 , 16 , 17 , 18 , 19 ]. In [ 19 ] prop osed an adaptive m ulti-fidelit y polynomial c haos (PC) MCMC algorithm to find a distribution that is “close” to the posterior in the sense of Kullback-Leibler divergence. Similar approac hes were proposed in [ 18 , 19 ] using an adaptive m ulti-fidelity PC based ensemble Kalman in version technique. Close in spirit to the w orks men tioned, in [ 20 ] prop oses the use of Ba yes F actors (BF; the o dds in fa vor) of the n umerical mo del vs the theoretical mo del. In an ODE framework, they sho w that the BF conv erge to 1, that is, b oth mo dels would b e equal, in the same order as the n umerical solver used. Later, this idea w as generalized in [ 1 ] to consider also PDEs and, more imp ortan tly , the use of the exp e cte d v alue of the BF s, b efore observing data. This results in more practical and w ork able guidelines in a more realistic m ultidimensional setting. The main result in [ 1 ] is a b ound for the exp ected BF. This b ound allows deciding what precision to run the solver, whic h could require less computational effort. Indeed, a reliable estimate of the error for the n umerical 2 metho d used is the central p oin t in the calculation of this b ound. Curren t efforts to estimate the discretization error fo cus on after-the-fact metho ds (i.e. a- p osteriori metho ds, we prefer to call then after-the-fact to av oid the obvious confusion with the Ba yesian jargon). These methods provide an error estimate only after the n umerical solution has been computed. They use the computed solution to the discrete equations, p ossibly with additional information supplied b y the equations, to estimate the error relative to the exact solution of the mathematical model [ 21 ]. Most of the previous w orks are based on higher error b ounds with asymptotic conv ergence when the mesh size tends to zero [ 22 , 23 , 24 , 25 ]. Unfortunately , these estimates imply “constants of stability” generally unkno wn and difficult to calculate. The resulting error estimation techniques, in practice, do not prov ide mathematically pro ven b ounds that, in general, can b e computed efficiently [ 26 ]. In this pap er, we derive an after-the-fact error estimate for a numerical approximation of the physical mo dels in volving a semi-linear evolution differential equation of the form: ∂ u ∂ t = D ∂ 2 u ∂ x 2 + F  u, ∂ u ∂ x , θ  , (1) defined on the region t ∈ [0 , τ ] , x ∈ [ a, b ] , with left and right b oundary conditions u ( a, t ) = g ( t ) and u ( b, t ) = h ( t ) , 0 ≤ t ≤ τ , (2) and initial condition u ( x, 0) = f ( x ) , a ≤ x ≤ b. (3) In Eq. ( 1 ), D is the diffusion co efficient, θ is a parameter (possibly a v ector) of in terest, and F is a non-linear op erator. This physical mo del arises in sev eral fields of science and engineering [ 27 , 28 ]. It is used to describ e many complex nonlinear settings in applications such as vibration and wa ve propaga- tion, fluid mechanics, plasma physics, quan tum mechanics, nonlinear optics, solid-state ph ysics, c hemical kinematics, ph ysical chemistry , p opulation dynamics, and many other areas of mathe- matical mo deling. The numerical solution for Eq. ( 1 ) is obtained b y discretizing first in the space with the finite difference (FD) method and solving the resulting system in time with the Runge-Kutta Cash-Karp (RKCK) metho d [ 29 ]. This scheme is widely used to solve numerically evolution partial differential equations [ 30 , 31 , 32 ]. How ever, numerical after-the-fact error estimates for these metho ds ha ve not y et b een derived. The idea behind our construction of the error estimates for the PDE in Eq. ( 1 ) is the a v ailable error estimates for the RKCK metho d. Our numerical metho d uses these error estimates, in time, for the resulting ODE system. The truncation error in tro duced for the appro ximation with finite differences is computed using the solutions in t wo differen t mesh sizes. In modern computers, the added computational effort can b e reduced to result equiv alent to solving the 3 PDE conv entionally (on a single mesh) since ev aluating the solution in tw o differen t meshes ma y b e easily parallelized. W e will incorporate this after-the-fact error estimate in the result of [ 1 ], for the solution of the Ba y esian Inv erse Problem (BIP) asso ciated with the PDE giv en in Eq. ( 1 ), to con trol the error in the p osterior distribution. Numerical examples are given to illustrate the efficiency of the prop osed metho d. W e obtain n umerical p osterior distributions, for unkno wn parameters, that are nearly iden tical to the corresp onding theoretical posterior, keeping their exp e cte d Bay es factor close to 1. The pap er is organized as follows. In Section 2 , we present a numerical metho d used for solving evolution partial differen tial equations numerically . In Section 3 , we derive our after-the- fact error estimate for semi-linear evolution differen tial equations. The accuracy of our error estimate is ev aluated for some classic examples. In Section 4 , w e prop ose an algorithm that incorp orates the after-the-fact error estimate to con trol the error in the posterior distribution. Numerical examples are given in Section 5 to illustrate the efficiency of the prop osed Algorithm. Finally , a conclusion is given in Section 6 . 2 Numerical Solution Here, w e introduce a common n umerical pro cedure for the solution of semilinear evolution partial differential equations. This pro cedure has b een widely used for solving evolution partial differen tial equations [ 30 , 31 , 32 ]. The basic idea of the metho d is to replace the spatial deriv ation in the PDE with an algebraic approximation in order to obtain an ODE system. The resulting system is then solv ed with a standard ODE solv er. W e discretize in the space with the FD metho d and solving the ODE system with the RKCK metho d. W e called this method FD- RK CK. F or simplicit y , w e denote ˙ u := ∂ u ∂ t , u 0 := ∂ u ∂ x , and F ( u, u 0 ) instead of F ( u, u 0 , θ ) . Moreo ver, without losing generality , w e can set D = 1 in Eq. ( 1 ). W e consider a one-dimensional uniform mesh, Ω h , on the region [ a, b ] , with no des x i , for i = 0 , 1 , . . . , N , where Ω h : a = x 0 < x 1 < · · · < x N = b, (4) and a constant step size h b et ween any tw o successive no des (i.e., h = x i − x i − 1 ). T o solv e the PDE in Eq. ( 1 ), we start b y linearizing F using the quasi-linearization metho d that w as in tro duced in [ 30 ] for solving nonlinear evolution partial differential equations. This metho d consists of separating the function F into a linear ( L ) and a nonlinear ( N ) comp onen t, and rewriting Eq. ( 1 ) in the form ˙ u = u 00 + L [ u, u 0 ] + N [ u, u 0 ] . (5) F or example, in Section 3 w e use the Fisher equation where F = ru (1 − u ) , thus L = r u and 4 N = − r u 2 . Afterwards, the nonlinear op erator N is approximated with a T aylor series, assuming that the difference u i +1 , · − u i, · and all its spatial deriv atives are small. Hence N [ u i +1 , · , u 0 i +1 , · ] ≈ N [ u i, · , u 0 i, · ] + φ 0 ,i [ u i, · , u 0 i, · ] · ( u i +1 , · − u i, · ) + φ 1 ,i [ u i, · , u 0 i, · ] · ( u 0 i +1 , · − u 0 i, · ) , (6) where u i, · := u ( x i , t ) is the solution of Eq. ( 1 ) ev aluated in ( x i , t ) , and φ k,i [ u i, · , u 0 i, · ] := ∂ N [ u i, · , u 0 i, · ] ∂ u ( k ) , k = 0 , 1 . F or simplicity , u ( i ) denotes the i-th deriv ative. Substituting Eq. ( 6 ) into Eq. ( 5 ), w e get ˙ u i +1 , · ≈ u 00 i +1 , · + L [ u i +1 , · , u 0 i +1 , · ]+ N [ u i, · , u 0 i, · ]+ φ 0 ,i [ u i, · , u 0 i, · ] · ( u i +1 , · − u i, · )+ φ 1 ,i [ u i, · , u 0 i, · ] · ( u 0 i +1 , · − u 0 i, · ) , (7) for i = 1 , ..., N − 2 . No w, the spatial partial deriv atives are approximated using the central difference form ula. F or simplicit y , w e use the simplest spatial deriv ativ e approximations here, while the analysis can be extended for other (e.g., five-point stencil) approximations as well, u 0 i, · ≈ u i +1 , · − u i − 1 , · 2 h , u 00 i ≈ u i +1 , · − 2 u i, · + u i − 1 , · h 2 , (8) for i = 1 , . . . , N − 1 , and u 0 0 , · ≈ u 1 , · − u 0 , · h . (9) Substituting Eqs. ( 8 )–( 9 ) in to Eq. ( 7 ), join t with the b oundary condions ( 2 ) and initial condition ( 3 ), we get the following semi-discrete differential equation: ˙ V h ( t ) = 1 h 2 A xx V h ( t ) + F ( t, V h ( t )) (10) V h (0) = U (0) (11) where V h ( t ) = ( v 1 , · , v 2 , · , . . . , v N − 1 , · ) T appro ximates U ( t ) = ( u 1 , · , u 2 , · , . . . , u N − 1 , · ) T , (12) U is the exact solution of the PDE ( 1 ) on the mesh Ω h , and F is the approximate op erator F in matrix form, see A for details. Remark 1. The semi-discr ete differ ential e quation ( 10 ) have a trunc ation err or O ( h p ) : (i) If F do es not have a nonline ar c omp onent, the quasi-line ar appr oximation ( 6 ) is not ne c es- sary. Thus, the trunc ation err or for the c entr al differ enc e formula is not affe cte d ( p = 2 ). 5 (ii) If N is non-line ar in u 0 , the quasi-line ar appr oximation ( 6 ) intr o duc es a trunc ation err or of first-or der, which is pr op agate d when u 0 is appr oximate d using the c entr al differ enc e formula. Thus, the or der of the trunc ation err or for ( 10 ) is less than 2 ( p < 2 ). (iii) If N is non-line ar in u and line ar in u 0 , the trunc ation err or intr o duc e d for the quasi-line ar appr oximation ( 6 ) is not pr op agate d as in c ase (ii). Then, the trunc ation err or in ( 10 ) is slightly affe cte d ( p ≈ 2 ). In order to solv e the resulting ODE system ( 10 )-( 11 ), with N − 2 equations, we use the RKCK metho d. This metho d uses six function ev aluations to calculate fourth and fifth-order accurate solutions. The difference b et ween these solutions is then tak en to b e the error (fourth-order) of the solution; see [ 33 ] for details. The av ailable error estimate is the reason to solv e the resulting ODE’s system with this Runge-Kutta (RK) method, and it will be used in turn, in Section 3 , for computing the after-the-fact error of the n umerical solution of Eq. ( 1 ). Setting G ( t, V h ( t )) = 1 h 2 A xx V h ( t ) + F ( t, V h ( t )) , a RK scheme applied to the ODE sys- tem ( 10 ), at a uniform time grid 0 = t 0 < t 1 < · · · < t n < · · · < t M − 1 < t M = τ ; t n +1 = t n + k , (13) is given by K 1 ,n = G ( t n , W · ,n ) , (14) K l,n = G   t n + c l k , W · ,n + k l − 1 X j =1 a lj K j,n   , l = 2 , 3 , . . . , 6 , (15) W · ,n +1 = W · ,n + k 6 X l =1 b l K l,n , n = 1 , . . . , M − 1 , where W · ,n +1 is the approximation for V h ( t n +1 ) , ( a lj ) are the Runge-Kutta coefficients, b = ( b 1 , b 2 , . . . , b 6 ) are the quadrature no des, and c = ( c 1 , c 2 , . . . , c 6 ) are the quadrature weigh ts of the RK scheme. k = ∆ t > 0 is the step size in time and define a uniform grid. In order to ha v e stable solutions in explicit schemes, the step size in time is related to the discretization through the Couran t-F riedrichs-Lewy (CFL) condition [ 34 ], which restricts the step size in time based on the eigenspectrum of the discretized spatial operator. The CFL condition for the FD-RKCK sc heme considering only the pure diffusion is ∆ t ∆ x 2 ≤ 1 4 b max , where b max = max i b i , and b 1 , b 2 , . . . , b 6 are the quadrature nodes for the RK metho d used, see B for details. 6 3 After-the-fact error estimates In this section, we prop ose a n umerical pro cedure to obtain an after-the-fact error estimate of the A GE, for the n umerical solution of Eq. ( 1 ). F or them, we use the error estimation in the time-stepping giv en for the RKCK metho d and estimate the leading term of the truncation error in space stepping. This sc heme can b e extended for differential equations of non-linear ev olution, but some additional considerations about the stabilit y of the solution must b e taken in to account. In Section 2 , w e obtained the semi-discrete differential equation ( 10 )–( 11 ), with a unique solution v ector, V h ( t ) , b eing a grid function on Ω h . This initial v alue problem solved with the RK CK metho d yields appro ximations W .,n to V h ( t n ) . The global error at the spatial mesh p oin ts at knot t n is defined by E h ( t n ) := W .,n − U ( t n ) , (16) where U is the exact solution of the PDE ( 1 ) on the mesh Ω h defined in ( 12 ). The v ector E h ma y also b e written as a combination of the ODE global error, this is defined as the error made b y the solv er, i.e., e h ( t n ) = W .,n − V h ( t n ) , (17) and the spatial discretization error defined by η h ( t n ) = V h ( t n ) − U ( t n ) . (18) The function η ( t ) represents the accum ulation of the spatial truncation error (TE) when we solv e ( 10 )-( 11 ), T E h ( t ) = G ( t, U ) − ˙ U ( t ) . (19) F rom Eqs. ( 17 )–( 18 ), the global error E h ( t n ) may b e written as the sum of the global time and spatial error, i.e., E h ( t n ) = e h ( t n ) + η h ( t n ) . (20) W e assume that u ( t, x ) is p -times differen tiable with respect to x and fourth-times contin- uously differentiable with resp ect to t . Then, it holds for the global space and time error that || η h || = O ( h p ) and || e h ( t n ) || = O ( k 4 ) , n = 1 , . . . , M , resp ectively . The ODE global error ( 17 ) is calculated using the error estimation of RKCK [ 30 ]. The spatial discretization error implementation based on ( 21 )-( 22 ) requires an estimation for the truncation error. The Richardson extrap olation [ 21 ] pro vides a suitable estimate of the truncation error. The idea is to calculate the solution using a one-step size h and then compute them again with half the space step ( h/ 2 ). The result obtained using tw o step s size is more accurate than using 7 the single-step size h. Their difference can b e used as an estimate of the truncation error, whic h is prop ortional to the p o wer of h . 3.1 Spatial discretization error W e can obtained an equation for the evolution of η ( t ) by adding terms to b oth sides of ( 10 ): ˙ V h ( t ) − ˙ U ( t n ) = G ( t, V h ) − G ( t, U ) + G ( t, U ) − ˙ U ( t n ) . F rom the initial condition ( 11 ) and using the definition η ( t ) in the ab o ve equation, the accumu- lation of the spatial discretization error is the solution to the initial v alue problem: ˙ η ( t ) = G ( t, V h ) − G ( t, U ) + T E h ( t ) , t ∈ (0 , τ ] (21) η (0) = 0 . (22) Assuming G to b e twice contin uously differentiable, we use the approximation: ∂ G ∂ V h ≈ G ( t, V h ) − G ( t, U ) V h − U . (23) Finally , we rewrite ( 21 )-( 22 ) to get ˙ η ( t ) = ∂ G ∂ V h η ( t ) + T E h ( t ) , t ∈ (0 , τ ] (24) η (0) = 0 . (25) The i n tegration of ( 24 )–( 25 ) is p erformed using M steps of size k of the RKCK metho d, as in the solution of the semi-discrete differential equation ( 10 )–( 11 ). In eac h RKCK step, ∂ G ∂ V h is appro ximated using the approximations W .,n and W .,n +1 to V h at time t n +1 , i.e., ∂ G ∂ V h ≈ G ( t n +1 , W .,n +1 ) − G ( t n +1 , W .,n ) W .,n +1 − W .,n . (26) 3.2 Spatial and time error The ODE global error ( 17 ) is computed by the error estimation giv en b y the RK CK metho d. This scheme uses an RK method with a fifth-order local truncation error to estimate the lo cal error in an RK metho d of fourth-order. Both with the same n umber of stages s = 6 , Runge Kutta matrix A , and w eights c , while their no des ˆ b and b , resp ectively , are differen t; see [ 30 ] for details. Let W · ,n +1 the n + 1 approximation of V h ( t n +1 ) of fourth-order, and let Y · ,n +1 b e obtained 8 b y the fifth-order metho d starting at W · ,n , namely W · ,n +1 = W · ,n + k s X i =1 b i K i,n and Y · ,n +1 = W · ,n + k s X i =1 ˆ b i K i,n , (27) The local truncation error ˆ τ · ,n +1 at node t n +1 of the RK metho d is defined as the error made in step n + 1 of the solver if starting at the exact v alue W · ,n . The estimation of ˆ τ · ,n +1 for the RK CK metho d is given b y ˆ τ · ,n +1 = Y · ,n +1 − W · ,n +1 = k s X i =1 ( ˆ b i − b i ) K i,n , and the global error at knot t n +1 is ˆ e .,n +1 = n +1 X j =1 ˜ τ · ,j . (28) In each RK iteration, w e solve the equation for the spatial discretization error ( 24 )–( 25 ), ˆ η n +1 = ˆ η n + k s X i =1 b i ˆ K i,n , (29) where ˆ K 1 ,n = H ( t n , ˆ η n ) ˆ K l,n = H   t n + c l k , η n + k l − 1 X j =1 a lj ˆ K j,n   , l = 2 , 3 , . . . , 6 . H is the right side of ( 24 ). F rom ( 28 )–( 29 ), an estimation for the global error E h ( 20 ) at knot t n is given by , ˆ E .,n +1 ≈ ˆ e .,n +1 + ˆ η n +1 . (30) Note that to solv e fully ( 24 )–( 25 ), we need an estimate for the truncation error. This estimation is done in parallel to be used in ( 29 ). Belo w w e give details for computing the truncation error. 3.3 Spatial truncation error The truncation error is the difference b et ween the discretized equations and the original partial differen tial equations. It con tains the errors due to the discretization of the PDE and the errors due to the grid. F or the finite difference scheme used to approximate the spatial op erator, we 9 ha ve that the truncation error at time t has rate order O ( h p ) , T E h ( t ) ≈ O ( h p ) . An efficient strategy to estimate the spatial truncation error b y Ric hardson extrap olation is prop osed in [ 21 ]. W e will adopt this approach to our setting. The actual mesh used to compute the numerical solution to the PDE is used as the fine mesh in the Richardson extrapolation pro cess. Supp ose w e are given a second semi-discretization of the PDE system ( 1 ), no w with doubled lo cal mesh sizes defined as follows, Ω 2 h := a = z 0 < z 1 < z 2 < . . . < z N/ 2 = b, z i = x 2 i , i = 0 , . . . , N . This mesh is called the coarse mesh. W e assume that the solution V 2 h ( t ) to the discretized PDE, on the coarse mesh 2 h , exists and is unique. The Ric hardson extrapolation giv es an estimation of the truncation error for the fine mesh at time t , d T E ( t ) ≈ V 2 h − R 2 h ( V h ) 2 p − 1 , (31) where R 2 h is the usual restriction operator defined b y R 2 h ( V h ) =  v 1 ,. , v 2 ,. , . . . , v ( N − 1) / 2 ,.  T , v i,. = v ( z i , t ) , z i ∈ Ω 2 h . Remark 2. In the c ompute of ( 31 ) , we have found in c omputational exp eriment that this term c ould b e appr oximate (on Ω h × (0 , τ ] ) at time t as d T E h ( t ) = h p τ . This estimation is valid for our settings, but is not applic able if another scheme for solving ( 1 ) is use d. Remark 3. The after-the-fact err or estimate of the AGE for the FD-RKCK solution of ( 1 ) , on Ω h × (0 , τ ] , is given by ˆ K = || ˆ E .,M || ∞ , wher e Ω h is define d in ( 4 ) , the discr etization grid for (0 , τ ] is define d in ( 13 ) , and ˆ E .,M is an estimate for the glob al err or ( 20 ) . In Algorithm 1 , w e describe the steps necessary to compute the numerical solution of Eq. ( 1 ), with the after-the-fact error estimation. W e call this algorithm DF-RK CK. T o test our algorithm, w e consider three classical semi-linear PDEs, of the form ( 1 ): Example 1 (Fisher equation), Example 2 (Fitzhugh-Nagumo equation), and Example 3 (Burgers-Fisher equation). The three examples used also hav e analytic solutions, allowing us to compute the actual n umerical error and compare it with our estimates. In Figure 1 , a graphical comparison 10 is sho wn b et ween our n umerical implementation approximations and the exact solution for the three examples. T able 1 shows the con vergences order of the solution obtained with the DF- RK CK Algorithm. It can be seen that the metho d achiev es full conv ergence for the error (order 2) for Examples 1 and 2 , but the order of conv ergence for Example 3 is 1 , and this is due to the non-linearit y of F in u 0 , as was mentioned b efore. Remark 4. T o c ompute the numeric al c onver genc e r ate, we use p = log 2  || u 4 h − u 2 h || ∞ || u 2 h − u h || ∞  . 0.0 0.2 0.4 0.6 0.8 1.0 t 0.2 0.4 0.6 0.8 1.0 u ( t , x ) Exac t N u me r i c 0.0 0.2 0.4 0.6 0.8 1.0 t 0.66 0.68 0.70 0.72 0.74 0.76 0.78 u ( t , x ) Exac t N u me r i c (a) (b) 0.0 0.2 0.4 0.6 0.8 1.0 t 0.2 0.4 0.6 0.8 1.0 u ( t , x ) Exac t N u me r i c (c) Figure (1) The analytical and appro ximate solution in x = 0 . 1 , 0 . 3 , 0 . 5 , 0 . 7 , with step sizes in space h = 0 . 0125 and in time k = 0 . 0001 , for (a) the Fisher’s equation, (b) the Fitzh ugh-Nagumo equation, and (c) the Burgers-Fisher equation. In Figure 2 , we show the maximum error b etw een the exact solution and the numerical solution for the three examples considered, comparing it to our error estimates. W e can see that the estimation prop osed for the absolute global error is an upper b ound for the exact error. The n umerical implemen tation has been performed in Python, using the scip y , n umpy , and 11 T able (1) Conv ergence Order h Example 1 Example 2 Example 3 0.0125 1.999092 1.992578 1.075049 0.0083 1.999617 1.995410 1.052625 0.00625 1.999789 1.996692 1.040450 0.005 1.999866 1.997418 1.032826 matplotlib pac k ages. F or the s ak e of repro ducibilit y , all code is av ailable in a Github rep ository [ 35 ]. 10 − 2 6 × 10 − 3 2 × 10 − 2 4 × 10 − 2 h 10 − 4 10 − 3 10 − 2 10 − 1 Er r or Examp l e 3 Examp l e 1 | u − u h | l ∞ Es t. Er r or 10 − 2 6 × 10 − 3 2 × 10 − 2 4 × 10 − 2 h 10 − 5 10 − 4 10 − 3 Er r or | u − u h | l ∞ Es t. Er r or (a) (b) Figure (2) The maxim um error b et ween the exact solution and our DF-RKCK method against the error estimation: (a) the Fisher and the Burgers-Fisher equation, (b) the Fitzhugh-Nagumo equation. Different step sizes ( h ) in space are taken and for time w e let k = αh 2 , with α = 3 / 4 . Example 1 (Fisher equation) . Fisher’s equation b elongs to the class of reaction-diffusion equation and is encoun tered in chemical kinetics and p opulation dynamics applications. The equation is giv en by ∂ u ∂ t = ∂ 2 u ∂ 2 x + r u (1 − u ) , (32) with b oundary and initial conditions u (0 , t ) = 1 (1 + e − 5 t ) 2 , 0 ≤ t ≤ τ , u (1 , t ) = 1 (1 + e 1 − 5 t ) 2 , 0 ≤ t ≤ τ , u ( x, 0) = 1 (1 + e x ) 2 , 0 ≤ x ≤ 1 . This PDE has the follo wing analytic close form solution u ( x, t ) = 1  1 + exp  p r 6 x − 5 r 6 t  2 , x ∈ [0 , 1] , and t ∈ [0 , τ ] , 12 where r is a parameter. The non-linear op erator is F ( u, u 0 , r ) = ru (1 − u ) ; hence the appropriate linear comp onen t is L = ru , and the non-linear comp onen t is N = − ru 2 ; we see that the op erator F do es not depend on u 0 so the method achiev es order 2, as can b e seen in T able 1 . F or the examples in Figs. ( 1 )–( 2 ) and T able 1 , we use r = 4 , and this parameter will b e tried to identify using synthetic data in section 5 . Example 2 (Fitzh ugh-Nagumo equation) . The Fitzhugh-Nagumo equation is giv en by ∂ u ∂ t = ∂ 2 u ∂ 2 x + u (1 − u ) ( u − a ) , 0 < a < 1 , (33) with b oundary and initial conditions u (0 , t ) = 1 2 (1 + a ) + 1 2 (1 − a ) tanh  1 − a 2  4 t ! , 0 ≤ t ≤ τ , u (1 , t ) = 1 2 (1 + a ) + 1 2 (1 − a ) tanh √ 2 (1 − a ) 1 4 +  1 − a 2  4 t ! , 0 ≤ t ≤ τ , u ( x, 0) = 1 2 (1 + a ) + 1 2 (1 − a ) tanh  √ 2 (1 − a ) x 4  , 0 ≤ x ≤ 1 . The analytic solution for this PDE is giv en by u ( x, t ) = 1 2 (1 + a )+ 1 2 (1 − a ) tanh √ 2 (1 − a ) x 4 +  1 − a 2  4 t ! , x ∈ [0 , 1] , and t ∈ [0 , τ ] , where a is a parameter. The non-linear op erator is F ( u, u 0 , a ) = u (1 − u ) ( u − a ) ; hence the appropriate linear comp onen t is L = − au , and the non-linear comp onen t is N = u 2 (1 − u + a ) ; w e see that the op erator F do es not depend on u 0 , so the metho d achiev es order 2, as can be seen in T able 1 . F or the examples in Figs. ( 1 )–( 2 ) and T able 1 , we use a = 0 . 3 and this parameter will b e tried to identify using syn thetic data in section 5 . Example 3 (Burgers-Fisher equation) . The Burgers-Fisher equation is given by ∂ u ∂ t = ∂ 2 u ∂ 2 x − r uu 0 + su (1 − u ) , (34) with the initial condition u ( x, 0) = 1 2 + 1 2 tanh  − r 4 x  , 0 ≤ x ≤ 1 , 13 and with b oundary and initial conditions u (0 , t ) = 1 2 + 1 2 tanh  r 2 8 + s 2  t  , 0 ≤ t ≤ τ , u (1 , t ) = 1 2 + 1 2 tanh  − r 4  1 −  r 2 + 2 s r  t  , 0 ≤ t ≤ τ , u ( x, 0) = 1 2 + 1 2 tanh  − r 4 x  , 0 ≤ x ≤ 1 . This problem also has an analytic solution giv en by u ( x, t ) = 1 2 + 1 2 tanh  − r 4  x −  r 2 + 2 s r  t  , x ∈ [0 , 1] , and t ∈ [0 , τ ] , where r and s are parameters. The non-linear op erator is F ( u, u 0 , a ) = − r uu 0 + su (1 − u ) ; hence the appropriate linear comp onent is L = su , and the non-linear comp onen t is N = − ruu 0 − su 2 . Of note is that the order of con vergence for the error is 1 b ecause F is nonlinear in u 0 . F or the examples in Figs. ( 1 )–( 2 ) and T able 1 , we use r = 4 . 5 and s = 5 . 5 , and these parameters will b e tried to identify using synthetic data in the next section. 4 Error Con trol in Bay esian UQ In this section, we discuss ho w to incorp orate the after-the-fact error estimate, prop osed in Section 3 , in the results of [ 1 ], to control the error in the p osterior distribution. W e follow the general setting of [ 1 ] for the statistical IP . Let Θ and V b e separable Banach spaces. Let F : Θ → V b e the FM (t ypically F ( θ ) , for all θ ∈ Θ , is the solution of a system of PDE’s) and H : V → A ⊆ R m the observ ation operator (e.g., H ( F ( θ )) is one particular state v ariable, for whic h w e hav e observ ations). The comp osition H ◦ F defines a mapping from the parameter space Θ to the data sample space in R m . Also, assume that f ( y | θ ) is a density for data y : f ( y | θ ) := f o ( y |H ( F ( θ ))); θ ∈ Θ , where f o ( y | η ( θ )) is a density function that in teracts with θ only through η ( θ ) ∈ R m . Let F α ( n ) b e a discretized version of the FM F , for some discretization α that dep ends on an in teger refinement n , e.g., a spatial step size in FD discretization. And, let f n ( y | θ ) := f o ( y |H ( F α ( n ) ( θ ))) b e the resulting discretized numerical likelihoo d. T o find reasonable guidelines, to c ho ose a discretization lev el, in [ 20 ] compare the n umeric p osterior with the theoretical posterior using Ba yesian model selection, namely Ba yes F actors (BF). Assuming an equal prior probabilit y π for b oth mo dels, the BF is the ratio of the nor- malization constants Z n ( y ) Z ( y ) , where Z ( y ) = Z f ( y | θ ) π ( θ ) dθ, 14 and Z n ( y ) is the corresp onding n umeric normalization constant. Later, in [ 1 ] try to con trol the BF b et ween the discretized mo del and the theoretical mo del, through the use of the Absolute BF (ABF), AB F := 1 2     Z n k ( y ) Z ( y ) − 1     . T o do that, they b ound the exp e cte d ABF (the EABF), E AB F = Z 1 2     Z n k ( y ) Z ( y ) − 1     Z ( y ) dy , in terms of estimates on the error in the numeric FM. In Theorem 6 , we state the main result of [ 1 ], and the follo wing are the assumptions required. Assume that we observe a pro cess y = ( y 1 , . . . , y m ) at lo cations z 1 , . . . , z m ∈ D . This is a general setting, to include PDEs and other IPs, in which the domain D may include, for example, space and time: z i = ( x i , t i ) . That is, z i is an observ ation at coordinates x i and at time t i . Assumption 4. Assume that, for al l y ∈ R m , the observation mo del f o ( y | η ) is uniformly Lipschitz c ontinuous on η , and for y ∈ R m , f o ( y | η ) is b ounde d. Mor e over, the FMs H ◦ F and H ◦ F α ( n ) ar e c ontinuous. Assumption 5. Assume a glob al err or c ontr ol of the numeric FM as ||H ( F ( θ )) − H ( F α ( n ) ( θ )) || ∞ < K α ( n ) . (35) Note that this is a glob al b ound, valid for al l θ ∈ Θ , and includes alr e ady the observational op er ator. That is, it is a glob al b ound, but is only a statement at the lo c ations H i ’s wher e e ach y i is observe d. Theorem 6. (Capistr án et al. [ 1 ]) With assumptions 4 – 5 , and assuming indep endent data, y , arising fr om a lo c ation-sc ale family, with sc ale p ar ameter σ 2 and lo c ation p ar ameter η = H ( F ( θ )) = ( H 1 ( F ( θ )) , . . . , H m ( F ( θ ))) T , namely f o ( y | η ) = m Y i =1 σ − 1 ρ  y i − η i σ  , (36) with ρ a b ounde d C 1 symmetric L eb esgue density in R , with R ∞ −∞ x 2 ρ ( x ) dx = 1 , then E AB F < ρ (0) K α ( n ) σ m. (37) Note that mo del ( 36 ) can b e written as y i = H i ( F ( θ )) + σ ε i , i = 1 , . . . , m, (38) 15 where each ε i has zero mean and unit v ariance, and its probability distribution function b elongs to the lo cation-scale family . 4.1 Cho osing a solv er discretization The bound obtained in Theorem 6 allo ws deciding what precision to run the solv er. The idea is to k eep the EABF b elo w a small threshold (e.g., 1 20 ) so that the BF is close to 1, and the difference b et w een the n umeric and the theoretical model is “not worth more than a bare mention” [ 36 , 37 ]. If we let the E AB F < b , we need the numerical error in the FM in ( 35 ) satisfies K α ( n ) < σ m b ρ (0) . (39) Note that, in practice, there is no need to establish the global b ound ( 35 ) theoretically , but rather b y a careful strategy for actual global error estimation. In most cases, the p osterior distribution is sampled using MCMC, whic h requires the appro ximated likelihoo d at each of man y iterations; an automatic pro cess of global error estimation and control will b e necessary to comply with ( 35 ). W e propose a MCMC algorithm with refinement to assure to comply the global b ound ( 35 ) for all θ in the parametric space of in terst. Assume we ha v e an algorithm to simulate from the posterior distribution. Algorithm 2 describ es a strategy for incorporating the b ound ( 39 ) and the after-the-fact error estimate, prop osed in Section 3 , to con trol the error in the p osterior distribution. The basic idea of Algorithm 2 is to start with a relativ ely large step size (e.g., h = 0 . 1 ), and the step size in time is established to keep the stabilit y condition k = 3 4 h 2 . At each iteration, θ i , of the MCMC, the FM, F h ( θ i ) , is computed, including the after-the-fact error estimate ˆ K h θ i , using Algorithm 1 . If the error in the FM do es not comply with the b ound in ( 39 ), then run the solv er again reducing the spatial step size by half. In the pro cess, we assure ( 35 ) for all θ ∈ Θ . 5 Numerical examples In this section, we use the three previous examples to show the p erformance of Algorithm 2 , in the solution of the corresponding BIP , using simulated data sets. W e sim ulate data as follows. The (synthetic) observ ations , y = ( y 1 , . . . , y m ) , are generated under an indep enden t Gaussian model f o ( y | η ) = m Y i =1 σ − 1 ρ  y i − η i σ  16 with ρ ( x ) = 1 √ 2 π e − x 2 2 , i.e., y i = H i ( F ( θ )) + σ ε i , (40) where the ε i ’s are indep endent and iden tically distributed as N (0 , 1) , θ ∈ Θ ⊂ R d is a vector of unkno wn parameters, and F ( θ ) represen ts the FM. In all our examples, we consider the v ariance, σ 2 , to b e known. The solution of ( 1 ) with its initial and b oundary conditions defines our FM, and we take H ( x ) = x as the observ ation op erator. W e consider the BIP to estimate the parameter θ giv en observ ations H i ( F ( θ )) = u ( x i , t 1 , θ ) at some p oints in space x i , for i = 1 , 2 , . . . , m and at a fixed time 0 < t 1 ≤ τ . W e let the system evolv e until time t 1 and then observ e it at the spacial lo cations x i ’s. The resulting observ ations are y i = u ( x i , t 1 , θ ) + σ ε i ; ε i ∼ N (0 , 1) . The IP will b e treated as a statistical inference problem under a Bay esian approach, setting a prior distribution on the unknown parameter. The IP will b e treated as a statistical inference problem under a Bay esian approach, setting a prior distribution, π Θ ( θ ) , on the unknown parameter, to obtain the p osterior distribution, π Θ | Y ( θ | y ) , from whic h all the required inferences are drawn [ 38 , 10 ]. The implemen tation was done using MCMC, through a generic MCMC algorithm, called the t-w alk [ 39 ]. Note that considering independent data with a Gaussian mo del, the first part of assumption 4 is right, and w e only require to verify that H ◦ F and H ◦ F α ( n ) are contin uous. Indeed, the latter is true if the observ ation operator is the identit y . With this scheme, all the necessary assumptions for Theorem 6 are satisfied. And, in this case, ρ (0) = 1 √ 2 π and the threshold B := σ m b ρ (0) in ( 39 ), for the n umerical error in the FM, is B = σ m √ 2 π 20 . (41) Example 7 (In verse Problem - Fisher’s equation) . W e consider the BIP to estimate θ = r in Fisher’s equation of Example 1 , given measuremen ts of H i ( F ( θ )) at time t 1 = 0 . 4 . The synthetic data are simulated with the error mo del ( 40 ), using the analytical solution for the FM, and the follo wing parameters: θ = 4 and σ = 0 . 007 , to maintain a 0.01 signal-to-noise ratio. The solution of ( 32 ) with its initial and b oundary conditions defines our FM. W e consider n = 8 observ ations at lo cations x i regularly spaced b et ween 0 and 1 . The data are plotted in Fig. 3 (a). Considering a tolerance b = 1 20 in ( 39 ) and with the standard error and sample size used, the error bound for the FM is B = 1 . 1 × 10 − 4 . W e require a prior distribution, π ( · ) , for the parameter θ ; it is assumed θ ∼ Gamma ( α 1 , β 1 ) with all known hyperparameters. Regarding the n umerical solver, we b egin with a (relatively) large step size, h = 0 . 05 , and the step size in time is established to k eep the stability condition k = 3 4 h 2 . Then, we start the Algorithm 2 . F or h = 0 . 01 , the b ound is achiev ed for all iterations. W e compare the p osterior distributions using the n umerical FM vs. the exact FM, with 200,000 iterations of the t-w alk; the histogram is rep orted with 150,000 samples since the first 17 (burn-in) 50,000 are discarded. The results are sho wn in Fig. 3 (b) and T able 2 . The differences observ ed in both results may b e attributed to the Monte Carlo sampling. 0.0 0.2 0.4 0.6 0.8 1.0 x 0.40 0.45 0.50 0.55 0.60 u ( x o b s , t o b s ) D ata T r u e mod e l 3.9 4.0 4.1 r 0 2 4 6 8 10 12 D e n s i ty P r i or T r u e N u me r i c Exac t (a) (b) Figure (3) (a) Fisher equation example data (blue p oin ts) and true mo del (black line), consid- ering r = 4 . (b) Comparison betw een n umerical (blue) and theoretical (magenta) posterior for parameter r . The green line represents the prior distribution. Example 8 (Fitzh ugh-Nagumo equation) . F or this example, the IP is to estimate θ = a in the Fitzh ugh-Nagumo equation of Example 2 , given measurements H i ( F ( θ )) at time t 1 = 0 . 3 . The syn thetic data are simulated with the error mo del ( 40 ), using the analytical solution for the FM, and the following parameters: θ = 0 . 3 and σ = 0 . 007 . The solution of ( 33 ) with its initial and b oundary conditions defines our FM. W e consider n = 8 observ ations at lo cations x i regularly spaced b etw een 0 and 1. The data are plotted in Fig. 4 (a). T o b e able to get the p osterior distributions, w e assume that θ ∼ Gamma ( α 2 , β 2 ) with all kno wn h yp erparameters. With the standard error and sample size used, and considering a tolerance b = 1 20 in ( 39 ), we hav e that the error b ound for the FM is B = 1 . 1 × 10 − 4 . Regarding the numerical solver, we b egin with a step size, h = 0 . 1 , and the step size in time k = 3 4 h 2 . Then, we start the Algorithm 2 . F or h = 0 . 0125 , the b ound is achiev ed for all iterations. W e compare the p osterior distributions using the n umerical FM vs. the exact FM, with 200,000 iterations of the t-w alk; the histogram is rep orted with 150,000 samples since the first (burn-in) 50,000 are discarded. The results are sho wn in Fig. 4 (b) and in T able 2 . The differences observed in b oth results may b e attributed to the Monte Carlo sampling. Example 9 (Burgers-Fisher equation) . F or this example, the IP is to estimate θ = ( r, s ) in the Burgers-Fisher equation of Example 7 , given measurements H i ( F ( θ )) at time t 1 = 0 . 2 . The syn thetic data are simulated with the error mo del ( 40 ), using the analytical solution for the FM, and the following parameters: θ = (4 . 5 , 5 . 5) and σ = 0 . 05 . The solution of ( 34 ) with its initial and boundary conditions defines our FM. W e consider n = 10 observ ations at lo cations x i regularly spaced b et ween 0 and 1. The data are plotted in Fig. 5 (a). T o get the p osterior distributions, we assume independent priors b et ween the parameters of 18 0.0 0.2 0.4 0.6 0.8 1.0 x 0.68 0.70 0.72 0.74 0.76 u ( x o b s , t o b s ) D ata T r u e mod e l 0.28 0.30 0.32 a 0 10 20 30 40 50 60 D e n s i ty T r u e P r i or N u me r i c Exac t (a) (b) Figure (4) (a) Fitzh ugh-Nagumo equation example data (blue p oints) and true model (blac k line), considering a = 0 . 3 . (b) Comparison b et ween numerical (blue) and theoretical (magenta) p osterior for parameter a . The green line represen ts the prior distribution. T able (2) Comparison of the p osterior mean (PM) of each parameter using the exact and the n umeric FM. Example 1 Example 2 Example 3 P arameter r a r s T rue 4 0.3 4.5 5.5 PM-Exact 3.9915 0.2988 4.1813 5.5476 PM-Numeric 3.9916 0.2989 4.1859 5.5274 the mo del. W e assume r ∼ Gamma ( α r , β r ) and s ∼ Gamma ( α s , β s ) with all known hyperpa- rameters. With the standard error and sample size used, and considering a tolerance b = 1 20 in ( 39 ), we hav e that the error bound for the FM is B = 6 × 10 − 4 . Regarding the numerical solv er, w e b egin with a step size h = 0 . 1 , and the step size in time k = 3 4 h 2 . Then, w e start the Algorithm 2 . F or h = 0 . 0017 , the b ound is achiev ed for all iterations. W e compare the p osterior distributions using the n umerical FM vs. the exact FM, with 200,000 iterations of the t-w alk; the histogram is rep orted with 150,000 samples since the first (burn-in) 50,000 are discarded. The results are shown in Fig. 5 (b)–(c) and in T able 2 . The differences observed in b oth results may b e attributed to the Monte Carlo sampling. As seen in Figs. 3 (b), 4 (b), and 5 (b)–(c), and T able 2 , the histograms and the posterior means obtained with the n umerical and the exact FM are practically iden tical. The small differences observ ed in b oth results may b e attributed to the effect of generating approximate samples from the p osterior distribution using MCMC metho ds. 19 0.0 0.2 0.4 0.6 0.8 1.0 x 0.4 0.5 0.6 0.7 0.8 0.9 u ( x o b s , t o b s ) D ata T r u e mod e l 0 2 4 6 8 r 0.0 0.1 0.2 0.3 0.4 D e n s i ty P r i or T r u e N u me r i c Exac t (a) (b) 0 2 4 6 8 s 0.0 0.1 0.2 0.3 0.4 0.5 D e n s i ty P r i or T r u e N u me r i c Exac t (c) Figure (5) (a) Burgers-Fisher equation example data (blue p oin ts) and true model (black line), considering θ = (4 . 5 , 5 . 5) . Histogram from the numerical (blue) and theoretical (magen ta) p osterior distribution for: (b) parameter r and (c) parameter s . The green line represents the prior distribution. 6 Conclusion This pap er propose d an error estimation for a class of partial differential equations motived b y its application in the uncertaint y quan tification area. Our error estimation allo ws us to apply the results obtained in [ 1 ] for con trolling the error in the respective numerical posterior for in verse problems that the forward mapping in volv es a semi-linear evolution PDE. W e presen ted three w orkout examples; in all cases, the n umerical error in the p osterior w as successfully con trolled, which led to a negligible increase in accuracy if the exact FM is considered. This, in turn, may result in CPU time sa ve, as cheaper/rougher solv ers are used. Although t wo numerical solutions are required for the error estimation, the added computa- tional effort can be reduced to result equiv alent to solving the PDE con ven tionally (on a single mesh) since ev aluating the solution in t wo different meshes ma y b e easily parallelized. F or future work, we plan to extend the method used for computing the error estimation 20 to nonlinear ev olution differen tial equations, but some consideration about the stabilit y of the solution and the conv ergence orders needs to b e added. 21 References [1] J. A. Christen, M. A. Capistrán, M. L. Daza-T orres, H. Flores-Argüedas, and J. Cricelio Mon tesinos-Lóp ez. Posterior distribution existence and error control in Banac h spaces in the Bay esian approach to UQ in inv erse problems. T echnical Rep ort 1712.03299, arXiv, Octob er 2018. [2] Jari P Kaipio and Colin F ox. The bay esian framew ork for inv erse problems in heat transfer. He at T r ansfer Engine ering , 32(9):718–753, 2011. [3] Caifang Cai, All Mohammad-Djafari, Sam uel Legoupil, and Thomas Ro det. Ba yesian data fusion and in v ersion in x-ray m ulti-energy computed tomography . In 2011 18th IEEE International Confer enc e on Image Pr o c essing , pages 1377–1380. IEEE, 2011. [4] Khosro w Chadan, David Colton, Lassi Päivärin ta, and William Rundell. An intr o duction to inverse sc attering and inverse sp e ctr al pr oblems . SIAM, 1997. [5] Da vid S Holder. Ele ctric al imp e danc e tomo gr aphy: metho ds, history and applic ations . CR C Press, 2004. [6] OR Burggraf. An exact solution of the inv erse problem in heat conduction theory and applications. Journal of He at tr ansfer , 86(3):373–380, 1964. [7] Ro el Snieder and Jeannot T ramp ert. Inv erse problems in geophysics. In W avefield inversion , pages 119–190. Springer, 1999. [8] Maria L Daza, Marcos A Capistrán, J Andrés Christen, and Lilí Guadarrama. Solution of the inv erse scattering problem from inhomogeneous media using affine inv ariant sampling. Mathematic al Metho ds in the Applie d Scienc es , 40(9):3311–3319, 2017. [9] Jari Kaipio and Erkki Somersalo. Statistic al and c omputational inverse pr oblems , v olume 160. Springer Science & Business Media, 2005. [10] Andrew M Stuart. In verse problems: a ba yesian p erspective. A cta numeric a , 19:451–559, 2010. [11] Y oussef M Marzouk, Habib N Na jm, and Larry A Rahn. Sto chastic spectral metho ds for ef- ficien t bay esian solution of inv erse problems. Journal of Computational Physics , 224(2):560– 586, 2007. [12] Da vid Galbally , Krzysztof Fidk owski, Karen Willcox, and Omar Ghattas. Non-linear mo del reduction for uncertaint y quantification in large-scale in verse problems. International jour- nal for numeric al metho ds in engine ering , 81(12):1581–1608, 2010. [13] Chad Lieberman, Karen Willcox, and Omar Ghattas. Parameter and state model reduc- tion for large-scale statistical inv erse problems. SIAM Journal on Scientific Computing , 32(5):2523–2542, 2010. 22 [14] Carl Edward Rasmussen, JM Bernardo, MJ Bay arri, JO Berger, AP Dawid, D Hec kerman, AFM Smith, and M W est. Gaussian pro cesses to speed up hybrid mon te carlo for expensive ba yesian integrals. In Bayesian Statistics 7 , pages 651–659, 2003. [15] Tiangang Cui, Y oussef M Marzouk, and Karen E Willcox. Data-driven mo del reduction for the bay esian solution of inv erse problems. International Journal for Numeric al Metho ds in Engine ering , 102(5):966–990, 2015. [16] Benjamin Peherstorfer, Karen Willcox, and Max Gunzburger. Surv ey of m ultifidelity meth- o ds in uncertain ty propagation, inference, and optimization. Siam R eview , 60(3):550–591, 2018. [17] Liang Y an and T ao Zhou. Adaptiv e m ulti-fidelity p olynomial chaos approach to ba y esian inference in inv erse problems. Journal of Computational Physics , 381:110–128, 2019. [18] Liang Y an and T ao Zhou. An adaptiv e m ultifidelity pc-based ensem ble k alman in version for inv erse problems. International Journal for Unc ertainty Quantific ation , 9(3), 2019. [19] Jinglai Li and Y oussef M Marzouk. Adaptiv e construction of surrogates for the ba yesian solution of in verse problems. SIAM Journal on Scientific Computing , 36(3):A1163–A1186, 2014. [20] M. Capistrán, J.A. Christen, and S. Donnet. Bay esian Analysis of ODE’s: solv er optimal accuracy and Bay es factors. Journal of Unc ertainty Quantific ation , 4(1):829–849, 2016. [21] Christopher Ro y . Review of discretization error estimators in scien tific computing. In 48th AIAA A er osp ac e Scienc es Me eting Including the New Horizons F orum and A er osp ac e Exp osition , page 126, 2010. [22] Iv o Babušk a and W erner C Rhein b oldt. A-p osteriori error estimates for the finite elemen t metho d. International Journal for Numeric al Metho ds in Engine ering , 12(10):1597–1615, 1978. [23] Marie E Rognes and Anders Logg. Automated goal-oriented error con trol i: Stationary v ariational problems. SIAM Journal on Scientific Computing , 35(3):C173–C193, 2013. [24] JP De SR Gago, DW Kelly , OC Zienkiewicz, and I Babusk a. A p osteriori error analysis and adaptiv e pro cesses in the finite elemen t metho d: Part ii—adaptive mesh refinemen t. International journal for numeric al metho ds in engine ering , 19(11):1621–1656, 1983. [25] Mark Ainsworth and J Tinsley Oden. A p osteriori err or estimation in finite element anal- ysis , volume 37. John Wiley & Sons, 2011. [26] Thomas Grätsc h and Klaus-Jürgen Bathe. A p osteriori error estimation tec hniques in practical finite element analysis. Computers & structur es , 83(4-5):235–265, 2005. 23 [27] Sp yros G T zafestas. Distribute d p ar ameter c ontr ol systems: The ory and applic ation , vol- ume 6. Elsevier, 2013. [28] T ok a Diagana. Semiline ar Evolution Equations and Their Applic ations . Springer, 2018. [29] Jeff R Cash and Alan H Karp. A v ariable order runge-kutta method for initial v alue prob- lems with rapidly v arying righ t-hand sides. ACM T r ansactions on Mathematic al Softwar e (TOMS) , 16(3):201–222, 1990. [30] PG Dlamini and M Khumalo. A new compact finite difference quasilinearization method for nonlinear evolution partial differential equations. Op en Mathematics , 15(1):1450–1462, 2017. [31] Mehdi Bastani and Dav o d Kho jasteh Salkuy eh. A highly accurate metho d to solv e fisher’s equation. Pr amana , 78(3):335–346, 2012. [32] Murat Sari and Gürhan Gürarslan. A sixth-order compact finite difference metho d for the one-dimensional sine-gordon equation. International Journal for Numeric al Metho ds in Biome dic al Engine ering , 27(7):1126–1138, 2011. [33] Ric hard L Burden and J Douglas F aires. Numerical analysis (ed.). Br o oks/Cole , 2000. [34] R Couran t, K F riedrichs, and H Lewy . On the partial difference equations op mathematical ph ysics. Mathematische Annalen , 100(1):32–74, 1928. [35] See supplemen tal material at h ttps://github.com/mdazatorres/Error_con trol , for the p ython co des. [36] RE Kass and AE Raftery . Ba yes factors. JOURNAL OF THE AMERICAN ST A TISTICAL ASSOCIA TION , 90:773–795, JUN 1995 1995. [37] H. Jeffreys. The ory of Pr ob ability . Oxford, Oxford, England, third edition, 1961. [38] C F ox, H Haario, and J A Christen. Bayesian The ory and Applic ations , chapter Inv erse problems Chapter 31. Oxford Universit y Press, 2013. [39] J Andrés Christen, Colin F ox, et al. A general purpose sampling algorithm for con tinuous distributions (the t-walk). Bayesian Analysis , 5(2):263–281, 2010. 24 A Details of the numerical solution F or the reader’s con venience, here w e describ e in detail the numerical metho d in tro duced in Section 2 . T o solv e the PDE in Eq. ( 1 ), we start by separating the function F into a linear ( L ) and a nonlinear ( N ) comp onen t and rewriting the Eq. ( 1 ) in the form ˙ u = u 00 + L [ u, u 0 ] + N [ u, u 0 ] . (42) The nonlinear op erator N is approximated with a T aylor series, assuming that the difference u i +1 , · − u i, · and all its spatial deriv atives are small, hence N [ u i +1 , · , u 0 i +1 , · ] ≈ N [ u i, · , u 0 i, · ] + φ 0 ,i [ u i, · , u 0 i, · ] · ( u i +1 , · − u i, · ) + φ 1 ,i [ u i, · , u 0 i, · ] · ( u 0 i +1 , · − u 0 i, · ) , (43) where u i, · := u ( x i , t ) is the solution of Eq. ( 1 ) ev aluated in ( x i , t ) and φ k,i [ u i, · , u 0 i, · ] := ∂ k N ∂ u ( k ) [ u i, · , u ( k ) i, · ] , k = 0 , 1 . Substituting Eq. ( 43 ) into Eq. ( 42 ), we get ˙ u i +1 , · = u 00 i +1 , · + L [ u i +1 , · , u 0 i +1 , · ]+ N [ u i, · , u 0 i, · ]+ φ 0 ,i [ u i, · , u 0 i, · ] · ( u i +1 , · − u i, · )+ φ 1 ,i [ u i, · , u 0 i, · ] · ( u 0 i +1 , · − u 0 i, · ) , (44) for i = 1 , ..., N − 2 . No w, the spatial partial deriv atives are approximated using the central difference formula ( 8 )-( 9 ). W e write in matrix form the central differences approximations V 0 h = 1 2 h A x V h + C x , (45) and V 00 h = 1 h 2 A xx V h + C xx , (46) where V h = ( v 1 , · , v 2 , · , . . . , v N − 1 , · ) T appro ximates U = ( u 1 , · , u 2 , · , . . . , u N − 1 , · ) T , A x =           0 1 0 . . . 0 − 1 0 1 . . . 0 0 − 1 0 . . . . . . . . . . . . . . . . . . 1 0 . . . 0 − 1 0           ( N − 2) × ( N − 2) , C x = 1 2 h          − u 0 , · 0 . . . 0 u N , ·          ( N − 2) × 1 , 25 A xx =           − 2 1 0 . . . 0 1 − 2 1 . . . . . . 0 1 − 2 · · · . . . . . . . . . . . . . . . 1 0 . . . 0 1 − 2           ( N − 2) × ( N − 2) , and C xx = 1 h 2          u 0 , · 0 . . . 0 u N , ·          ( N − 2) × 1 . No w, u 0 0 , · is approximated with the forw ard difference sc heme ( 9 ), whic h lea ves us b V 0 h = 1 2 h b A x V h + b C x , (47) b A x =           2 0 0 . . . 0 0 1 0 . . . . . . − 1 0 1 · · · . . . . . . . . . . . . . . . 1 0 . . . 0 1 0           ( N − 2) × ( N − 2) , and b C x = 1 2 h          − 2 u 0 , · − u 0 , · . . . 0 u N , ·          ( N − 2) × 1 . where b V h = ( v 0 , · , v 1 , · , . . . , v N − 2 , · ) T . Finally , substituting the appro ximate deriv atives in Eqs. ( 45 )–( 47 ) in to Eq. ( 7 ), join t with the b oundary conditions giv en in Eq. ( 2 ), we get the follo wing semi-discrete differential equation: ˙ V h = 1 h 2 A xx V h + C xx + L  V h , 1 2 h A x V h + C x  + N  b V h , 1 2 h b A x V h + b C  (48) + Φ 0  b V , 1 2 h b A x V h + b C  ·  V h − b V h  + Φ 1  b V h , 1 2 h b A x V + b C  ·  1 2 h ( A x − b A x ) V h + ( C x − b C x )  , where Φ k [ b V , b V 0 ] = ( φ k, 0 [ v 0 , · , v 0 0 , · ] , φ k, 1 [ v 1 , · , v 0 1 , · ] , . . . , φ k,N − 3 [ v N − 3 , · , v 0 N − 3 , · ]) T . Note that the right-hand side of equation ( 48 ) only depends on V and t , due v 0 , · and v N , · are known ( 2 ). Th us, we can write ( 48 ) in a compact form, ˙ V = 1 h 2 A xx V + F ( t, V ) (49) with 26 F ( t, V ) = C xx + L  V h , 1 2 h A x V h + C x  + N  b V h , 1 2 h b A x V h + b C  + Φ 0  b V , 1 2 h b A x V h + b C  ·  V h − b V h  + Φ 1  b V h , 1 2 h b A x V + b C  ·  1 2 h ( A x − b A x ) V h + ( C x − b C x )  . B Stabilit y Considerations W e briefly describe stability considerations for the DFRK metho d introduced in Section 2 . Let W · ,n +1 = W · ,n +  b 1 k h 2 A xx W · ,n + k b 1 F ( t n , W · ,n )  + 6 X i =2 b i k h 2 A xx   W · ,n + k i − 1 X j =1 a ij K j   + k 6 X i =2 b i F   t n + c i k , W · ,n + k i − 1 X j =1 a ij K j   the solution of ( 1 ) using the FD-RK CK metho d (see Section 2 ). T o determine the CFL condition, w e consider only the pure diffusion. Thus, the scheme is stable only if ρ ( A ) ≤ 1 [ 33 ], where A = b max k h 2 A xx b max = max i b i . The eigenv alues of A can b e shown to be µ i = − 4 λ  sin  iπ 2 N  2 , for i = 1 , 2 , . . . , N − 2 , where λ = b max k h 2 . So, the condition for stability consequently reduces to determining if ρ ( A ) = max 1 ≤ i ≤ N − 2      − 4 λ  sin  iπ 2 N  2      ≤ 1 , and this simplifies to 0 ≤ λ  sin  iπ 2 N  2 ≤ 1 4 , ∀ i = 1 , 2 , . . . , N − 2 . Stabilit y requires that this inequality condition hold as h → 0 , or, equiv alently , as N → ∞ , lim N →∞  sin  ( N − 1) π 2 N  2 = 1 . 27 Th us, stability o ccurs if only if 0 ≤ λ ≤ 1 4 . By definition λ = b max k h 2 , so this inequalit y requires that h and k b e chosen such that b max k h 2 ≤ 1 4 . The method conv erges to the solution with a rate of conv ergence O  h p + k 4  , provided b max k h 2 ≤ 1 4 . F or the numerical implementation, we take k = αh 2 , with α = 1 4 b max . 28 Algorithm 1 DF-RK CK Step 1: Initialization: • Spatial step size h . The step size in time is given for k eeping the stabilit y condition k = αh p • Initial conditions W h 0 and W 2 h 0 ; initial time t 0 ; parameter θ ; ˆ e ., 0 = 0 and d T E = 0 • The RK matrix A = ( a ij ) , the no des b and ˆ b , and the weigh ts c Step 2. Discretizing ( 1 ) with the FD metho d for h and 2 h , as is describ ed in Section 2 . Step 3. Solv e ( 10 ) with the Cash-Karp metho d for the step size h and 2 h : F or n = 1 , 2 , . . . , M : K h 1 ,n = G ( t n , W · ,n ) ; ˆ K 1 ,n = H ( t n , ˆ η n ) K 2 h 1 ,n = G ( t n , W · ,n ) Step 4. F or i = 2 , 3 , . . . , 6 : K h i,n = G   t n + c i k , W · ,n + k i − 1 X j =1 a ij K h j,n   ˆ K i,n = H   t n + c i k , ˆ η n + k i − 1 X j =1 a lj ˆ K j,n   K 2 h i,n = G   t n + c i k , W · ,n + k i − 1 X j =1 a ij K 2 h j,n   Step 5. Compute W h · ,n +1 = W h · ,n + k 6 X i =1 b i K h i ; ˆ η n +1 = ˆ η n + k 6 X i =1 b i ˆ K i,n W 2 h · ,n +1 = W 2 h · ,n + k 6 X i =1 b i K 2 h i ; d T E =   R 2 h ( W h · ,n +1 ) − W 2 h · ,n +1   ∞ h p (2 p − 1) ˆ e .,n +1 = ˆ e .,n + k 6 X i =1 ( b i − ˆ b i ) K h i b E .,n +1 = ˆ e .,n +1 + ˆ η n +1 Step 6. Compute the maximum absolute global error in the solution approximated W h , b K =    b E    ∞ Step 7: Output: W h , b K 29 Algorithm 2 Numerical refinement for the FM in the MCMC algorithm Step 1: Initialization: • Spatial step size h (large) • Standard error ( σ ) , sample size ( m ) , and ρ (0) • Calculate the error b ound B = σ m b ρ (0) , with a tolerance b , we suggest b = 1 20 • Initial v alue for the parameter, θ 0 • MCMC length M (num b er of simulations) Step 2 . F or i = 1 , 2 , . . . , M : Step 3. Compute the FM, F h ( θ i − 1 ) , and the error estimation ˆ K h θ i − 1 , using Algorithm 1 Step 4. If ˆ K h θ i − 1 > B - Set h = h/ 2 - Return to Step 3 Else - Simulate θ i with some MCMC algorithm Step 5: Output: ( θ 0 , θ 1 , . . . , θ M ) 30