Partial differential equations are powerful tools for used to characterizing various physical systems. In practice, measurement errors are often present and probability models are employed to account for such uncertainties. In this paper we present a Monte Carlo scheme that yields unbiased estimators for expectations of random elliptic partial differential equations. This algorithm combines a multilevel Monte Carlo method (Giles (2008)) and a randomization scheme proposed by Rhee and Glynn (2012), (2013). Furthermore, to obtain an estimator with both finite variance and finite expected computational cost, we employ higher-order approximations.

We derive and analyze adaptive solvers for boundary value problems in which thedifferential operator depends affinely on a sequence of parameters. These methods convergeuniformly in the parameters and provide an upper bound for the maximal error. Numericalcomputations indicate that they are more efficient than similar methods that control theerror in a mean square sense.