A modified perturbation method for mathematical models with randomness: An analysis through the steady‐state solution to Burgers' partial differential equation

The variability of the data and the incomplete knowledge of the true physics require the incorporation of randomness into the formulation of mathematical models. In this setting, the deterministic numerical methods cannot capture the propagation of the uncertainty from the inputs to the model output. For some problems, such as the Burgers' equation (simplification to understand properties of the Navier–Stokes equations), a small variation in the parameters causes nonnegligible changes in the output. Thus, suitable techniques for uncertainty quantification must be used. The generalized polynomial chaos (gPC) method has been successfully applied to compute the location of the transition layer of the steady‐state solution, when a small uncertainty is incorporated into the boundary. On the contrary, the classical perturbation method does not give reliable results, due to the uncertainty magnitude of the output. We propose a modification of the perturbation method that converges and is comparable with the gPC approach in terms of efficiency and rate of convergence. The method is even applicable when the input random parameters are dependent random variables.


INTRODUCTION
The numerical simulation of a complex system leads to errors that are usually grouped into three distinct families: model errors, data errors and numerical errors. 1 The first two types of error reflect the inherent uncertainty in the model predictions. Hence, randomness must be incorporated from the beginning into the model formulation to predict the true physics correctly. This is specially true in mathematical models that show high sensitivity to the input parameters. The Burgers' partial differential equation (PDE) is a simplification of the incompressible Navier-Stokes equations for the flow of a Newtonian fluid, { ∇ · u = 0, in D × (0, T), u t + (u · ∇)u + ∇p − ∇ 2 u = 0, in D × (0, T).
Here, D ⊆ R d , T > 0, is the constant fluid density, p is the fluid pressure, is the fluid viscosity, and u is the fluid velocity represented as a vector field. 2 If we consider the first component of the second equation in Equation (1), call u the first component of u in the spatial domain, let d = 1 and D = (a, b), and the pressure is neglected, we arrive at the viscous Burgers' PDE: u t + uu x = u xx , x ∈ (a, b), t ∈ (0, T), where = ∕ is a new constant, the kinematic viscosity. It is the simplest equation that combines nonlinear propagation effects uu x and diffusive effects u xx . The presence of the viscosity > 0 smooths out the shock discontinuity that would develop in the inviscid case. There is thus a transition layer z, defined as the zero of the solution u(t, z) = 0 and which is a region of rapid variation, whose location at the steady state is supersensitive to small variations in and the boundary conditions. 3 Consider the specific Burgers' PDE where D = (a, b) = (−1, 1). There is a small perturbation > 0 in the left boundary. In Physics, this may be due to error measurements. The steady-state solution (u t = 0) satisfies the ordinary differential equation (ODE) problem The location of the transition layer is the zero u(z) = 0. When > 0 is random and the solution u(x) becomes a stochastic process, appropriate stochastic numerical methods must be carried out to capture this location. In Xiu and Karniadakis, 4 several techniques were compared: brute-force Monte Carlo simulation, generalized polynomial chaos (gPC) expansions together with the stochastic Galerkin projection technique, and the perturbation method. 5 The Monte Carlo simulation is a popular statistical method based on obtaining realizations of the solution u(x), by generating random realizations of the input random parameter. 6 The statistics of the solution are estimated from the sample statistics. Although it is a robust and easy to implement method, it may become inefficient due to its nonmonotonic and slow convergence. On the other hand, the perturbation method expands u(x) as a power series in the centered inputs; satisfactory approximations using finite-term series are subject to small variations in the inputs and the model output. [7][8][9] In fact, for the Burgers' Equation (2), no convergence was perceived due to its supersensitivity nature. The gPC approach showed the main advantage for Equation (2). The solution is expanded in terms of orthogonal polynomials in the input parameters, and combined with a stochastic Galerkin projection technique, 10-12 rapid approximations of the solution statistics are obtained.
In this paper, we propose a modification of the classical perturbation method to obtain faithful representations of the model solution. The transition layer of the stochastic solution to Equation (2) will be accurately determined in the random scenario. The efficiency and accuracy will be comparable with the gPC approach. We will work, first, with a random left boundary value 1 + , and second, with a random viscosity being nonindependent to . This is of high mathematical interest, as most of the methods in uncertainty quantification rely on parameterizations based on independent inputs. 5

A MODIFIED PERTURBATION METHOD
The classical perturbation method is mainly applied to random differential equations. If there is one random input parameter , the solution u(x) = u(x, ) is formally expanded as a Taylor-based series around the mean value E[ ] 5,7 : where u i (x) are unknown deterministic functions, given by i u(x, )| =E[ ] ∕i!. The leading coefficient u 0 (x) represents the unperturbed part. The expansion is formally substituted into the governing equation, and by matching the terms according to the powers of − E[ ], deterministic differential equations for u 0 (x), u 1 (x), and so forth are obtained recursively. Due to the complexity of the derivation, in practice the expansion (3) is typically stopped at i ≈ 2. Once the expansion is truncated, the statistics of u(x) are approximated using the corresponding statistics of the finite-term series. Good approximations are provided only when the variability of and u(x) is of small magnitude. Indeed, consider the random Lebesgue space The random space L 2 is the correct framework for uncertainty quantification, as it corresponds to the set of random variables with finite variance and its topology preserves the convergence of the expectation and the variance. Suppose that the analytic expression (3) holds in L 2 as ] is the 2i-centered absolute moment. Two conditions are necessary to derive satisfactory approximations from a truncated perturbed expansion (3): 1. m 2i is not large. Vaguely, this holds if the dispersion and the tails of the distribution of are moderate. 2. u(x, ·) is not too sensitive to small perturbations in , so that its derivatives with respect to are moderate. In the case of the Burgers' equation, the solution presents high variations and is extremely sensitive to small changes in the boundary conditions, thus the reason of the perturbation method performing poor.
Clearly, the expansion (3) is equivalent to the representation (4) we loose the concept of expansion in terms of the random perturbation − E[ ] and we do not have a Taylor series proper, it might be notationally simpler to use in formal manipulations in x.
In order to avoid poor approximations of u(x) when the amount of variability of and u(x) is high, we propose a different approach. Truncate (4). We seek an approximation of the form Notice that we are making the deterministic coefficients u i,N (x) that depend on the truncation level N. Now, we use the ideas of gPC expansions, despite no orthogonality relations holding here. We multiply by k , k = 1, … , N, and apply expectation E: The problem in this derivation is that f depends precisely on the unknown solution u(x). We approximate each f k using a Gaussian quadrature rule for integration: Here,  is the support of , f is the probability density function of , Q is the quadrature degree, y (j) are the nodes defined as the zeros of the Q-th degree orthogonal polynomial with respect to the weight function f , and (j) > 0 are the weights. For each y (j) , one solves the deterministic governing equation to determine u(x, y (j) ). Gauss quadratures are exact for polynomials of degree 2Q − 1; they yield the highest degree of accuracy. Letf Q = (̃Q k ) 0≤k≤N . We solve the linear system . The final approximation becomes The statistics of u(x) can be approximated from the statistical moments of . For instance, the expectation and the variance of u(x) are estimated as where In the terminology of uncertainty quantification, this method may be encompassed under the class of pseudospectral methods. 5 It is a spectral method, in the sense that it reconstructs the functional dependence of the solution on the input random parameter. Therefore, it is expected to inherit the spectral convergence with N of this type of methods. But it is also a collocation method, as it is based on nodes locations y (j) and an ensemble of deterministic solutions {u(x, ( ) )} Q =1 to achieve the final approximation. The good thing about collocation methods is the easy implementation. Unlike the stochastic Galerkin projection method, which requires new deterministic solvers and adaptation of existing codes for the Galerkin system of coupled differential equations (intrusive approach), the modified perturbation method only requires Q applications of the deterministic solver for the original governing equation (non-intrusive approach), which may even be carried out in parallel. The construction of the Galerkin system may not even be clear in some cases, for example, when dealing with a simple linear wave equation subject to random transport velocity. 13 Obvious changes allow adapting the modified perturbation method to stochastic systems with several input random parameters = ( 1 , … , s ). This is useful for complex models, where there may be different sources of uncertainty. The solution u(x) = u(x, ) is approximated as where i = (i 1 , … , i s ) represents a multi-index of modulus |i| = i 1 + … +i s , and i = i 1 1 … i s s . All of the previous formulas for the 1D case are easily adapted to this multivariate setting: In this case, when s is low or moderate, a cubature rule for the integration f k based on a tensor product construction from Q-th degree univariate Gaussian quadratures can be applied: k ≈̃Q k . When has bounded support, a simple tensor Gauss-Legendre quadrature may solve the integration problem (in this case f is not considered as the weight function, but as part of the integrand). When s is large, the integration becomes more challenging and schemes based on sparse grids shall be chosen. Notice that no independence assumption for the input parameters has been required. The final approximation becomes The modified perturbation method suffers from the so-called curse of dimensionality. When the number of input random parameters s is large or the quadrature degree Q must be increased, the efficiency is penalized, as lots of deterministic realizations of the governing equation must be solved (Q s ; exponential growth with s). The convergence of the integration schemes depends on the smoothness of the corresponding integrand: from exponential convergence if it is analytic, to maybe sub-algebraic convergence if it is merely continuous. If the integrand is smooth with large derivatives (sharp peaks), then the convergence rapidity deteriorates. When the integrand is discontinuous, convergence can still be achieved as long as the quadrature nodes do not coincide with the discontinuity points.
Another issue that we have not commented yet is the possible ill-conditioning of the matrix G. For example, if s = 1 and ∼ Uniform(0, 1), then G is the well-known Hilbert matrix, which is the typical example of ill-conditioned matrix. In practice, due to the fast convergence of Equation (6) with N, we expect a small or moderate truncation order N, so G may not pose serious numerical challenges. Finally, Equation (5) may entail difficulties due to catastrophic cancellation (loss of precision when subtracting two nearly equal numbers, which might be the case for G ij and E i E j ).
To conclude this section, we state and prove a theoretical result to validate the modified perturbation method. The core idea of the method is that it constructs the best polynomial approximation in L 2 . = ( 1 , … , s ). Let u(x) = u(x, ) be the solution. If has bounded support, and u(x, ·) is continuous on the support of , for each x, then

Proposition 2.1. Consider a stochastic model with input random parameters
where the limit on N is regarded in the random Lebesgue space L 2 , and the limit on Q is pointwise. Moreover, the convergence is faster depending on the smoothness of u(x, ·) (this kind of convergence is referred to as spectral).
Proof. The continuity of u(x, ·) gives the convergence of the quadrature rule for integration: k = lim Q→∞̃Q k . Then u N (x) = lim Q→∞ũ N,Q (x), that is,

By definition of the coefficients
Then, considering the Hilbert space L 2 with the inner product E[·], and the closed vector subspace P N ( ) of s-variate polynomials of degree less than or equal to N, we deduce that S N u(x) is the best approximation to u(x) in P N ( ) ⊆ L 2 . On the other hand, as u(x, ·) is continuous, for each x the Weierstrass approximation theorem gives a sequence of polynomials {q N ( )} ∞ N=1 , of degrees N, that converges uniformly in to u(x, ). As a consequence, Remark 1. It is a standard result of deterministic ODEs that, if the map of the governing equation is of class C q (0 ≤ q ≤ ∞) with respect to x, u and the parameters , then the solution u(x, ) is C q in x and . 14 Thus, the condition on the continuity of u(x, ·) in the proposition is actually mild.

Random boundary condition
Consider the steady-state problem (2). In Xiu and Karniadakis, 4 a comparative study between several stochastic methods was performed when the left boundary perturbation follows a Uniform(0, 0.1) distribution and the viscosity is 0.05. Although the variability of this input is small, the propagation to the model output may yield a variability as high as 40%.
We are interested in the location of the transition layer z of the mean steady-state solution. Given the stochastic solution u(x), its mean is x  →ū(x) = E[u(x)]. We seek z such thatū(z) = 0. We also want to determine the standard deviation of the solution at z, (z) = √ V[u(z)]. In Tables 1 and 2, we show our computations of z and (z) for the Monte Carlo simulation, the gPC-based Galerkin method with Legendre-chaos (we have chosen Legendre polynomials to perform the gPC expansion because has a Uniform distribution 5 ), and the classical and modified perturbation methods. The deterministic problem (2) has been solved using the trapezoidal rule for the associated first-order ODE (1-step Adams-Moulton method). Due to its A-stability, it works well for this stiff problem, compared with explicit methods (such as the classical Runge-Kutta scheme). 15 , Ch. 4 The transition layer is computed using the Newton's method for roots. The Monte Carlo simulation is time-consuming, specially for this problem characterized by supersensitivity and the challenge of computing accurately   its numerical solution a lot of times (once per realization). The Legendre-chaos gives converged results up to several significant digits for small orders of basis. The classical perturbation method does not provide converged results due to the large variability of the model output with (see item 2 at the beginning of the previous section); we tabulate the results up to a fourth-order expansion. Finally, the estimates of the modified perturbation method are also reported, for truncation level N and quadrature degree Q. A univariate Gauss-Legendre quadrature rule for integration for the expansion coefficients has been utilized. We observe fast convergence with N, as the gPC approach: for N = 8 and Q = 8 convergence up to five significant digits is obtained. In terms of efficiency, the gPC-based Galerkin method of degree N requires the numerical solution of N + 1 coupled differential equations, which needs additional efforts in adapting code, whereas the modified perturbation method requires solving Q uncoupled governing differential equations. Figure 1 depicts the average steady-state solutionū(x) and the standard deviation (x) of the steady-state solution, computed using the modified perturbation method for (N, Q) = (8,8), together with the extreme deterministic solutions for = 0 and = 0.1. This figure is analogous to that presented in 4 employing the gPC-Galerkin method. We observe that the highest dispersion of the solution is produced near the transition layer, as it is the region of fastest variation.

Dependent random boundary condition and viscosity
In Equation (2), let ( , ) = (0.1Z 1 , 0.01Z 2 + 0.04), where (Z 1 , Z 2 ) ∼ Dirichlet(9, 6, 6): This is a complex case where the input random parameters are nonindependent. The boundary perturbation varies between 0 and 0.1, whereas the viscosity lies in (0.04, 0.05). These small perturbations produce significant changes in the model output u. As before, we are interested in the location of the transition layer z of the average steady-state solution, u(z) = 0. We use the modified perturbation method with respect to = (Z 1 , Z 2 ). The deterministic coefficients of the expansion are computed using tensor Gauss-Legendre quadratures, each one of degree Q, on [0, 1] 2 . In this example,    the curse of dimensionality becomes evident, as the number of quadrature nodes is Q 2 (Q nodes per dimension) and the number of summands in the modified perturbation expansion is (N + 2)(N + 1)∕2 (number of multi-indices i = (i 1 , i 2 ) such that i 1 , i 2 ≥ 0 and |i| = i 1 + i 2 ≤ N). Although the efficiency is penalized, accurate results are still obtained. The numerical results are reported in Table 3. We show the approximations of z and (z) of the modified perturbation method, as well as the Monte Carlo estimates for validation. The numerical details on the resolution are the same as the previous section: trapezoidal rule for each realization in the governing ODE and Newton's method for the transition layer. The estimates of the modified perturbation method agree up to five significant digits for N = 4 and Q = 25. It was checked that the convergence with N is very fast, as it usually occurs with the spectral methods, but Q has to be taken sufficiently large, otherwise we make significant errors in the final approximations, specially for the standard deviation.

CONCLUSION
In this paper, we have proposed a modification of the classical perturbation method for uncertainty quantification. Given a stochastic system, the solution is approximated in the mean square sense via a finite-term power series on the input random parameters. The coefficients of the power series, instead of being obtained by matching terms according to the powers, are determined by applying inner products in random space and then solving a linear system. The inner products involving the unknown model solution require the use of numerical integration schemes.
These new approximations present spectral convergence to the model solution in the mean square sense. However, the method suffers from the curse of dimensionality and the ill-conditioning of the linear system matrix.
Burgers' equation is a benchmark for the analysis of new stochastic techniques, due to its simplicity in formulation and its stiffness and supersensitivity with respect to the parameters. Our modified perturbation method has been applied to determine the steady-state solution to Burgers' equation and the position of the transition layer, when the left boundary value and the viscosity are random variables, not necessarily independent.
Albeit our method has been efficient for this problem, its success and advantage is problem-dependent. Further analysis shall be carried out in the future. We remark that, if the dimension of the random space were larger, the computational cost would be prohibitively large, and other methods shall be sought.
Alternative perturbation methods, of different nature, were applied to some deterministic equations to derive solutions in the form of series expansions. [16][17][18][19] These methods have not yet been randomized and could become new efficient approaches for uncertainty quantification.