Uncertainty quantification in computational fluid dynamics¶
In this tutorial we will learn
Polynomial chaos and Monte Carlo¶
This tutorial raises a very important question. Why bother using polynomials for estimating moments? What exactly is the advantage? Moreover, are we guaranteed that we will converge to the Monte Carlo solution? The answer is a resounding yes! Infact this is precisely what Dongbin Xiu and George Karniandakis showed in their seminal paper . As always we begin with some definitions: Parameter, Poly and Basis.
from equadratures import * import numpy as np
For our model problem, let’s consider Rosenbrock’s function
where we will assume that \(x_1\) and \(x_2\) are two uncertainties. We will assume that both parameters are Gaussians with \(\mu=1\) and \(\sigma=2\). Our objective is to compute the mean and variance in the output. We start by defining our computational model
def rosenbrock_fun(x): return (1 - x)**2 + 100*(x - x**2)**2
Next, we set the number of evaluation points in each direction. Lets opt for 7 points along each direction—more than sufficient to approximate the function exactly.
mu = 1 sigma = 2 variance = sigma**2 x1 = Parameter(distribution="Gaussian", shape_parameter_A=mu, shape_parameter_B=variance, order=4) x2 = Parameter(distribution="Gaussian", shape_parameter_A=mu, shape_parameter_B=variance, order=4) parameters = [x1, x2]
Now, we can set the problem up, compute the coefficients, and then ask Effective Quadratures to output the mean and the variance.
parameters = [x1, x2] basis = Basis('Tensor grid') uqProblem = Poly(parameters, basis, method='numerical-integration') uqProblem.computeCoefficients(rosenbrock_fun) uqProblem.set_model(rosenbrock_fun) mean, variance = uqProblem.get_mean_and_variance() print(mean, variance) >> 6804.000000000022 476659232.0000047
Now, we compare these results with Monte Carlo.
large_number = 1000000 s = sigma * np.random.randn(large_number,2) + mu f = np.zeros((large_number,1)) for i in range(0, large_number): f[i,0] = rosenbrock_fun([s[i,0], s[i,1]]) print np.mean(f), np.var(f) >> 6813.941920252046 483131338.5544447
The results are very close! In fact the polynomial approximation results are exact, because Rosenbrock’s function is a polynomial itself!
But what order should we use? This is a tough question to answer without any apriori knowledge of the function we wish to obtain statistical moments from. We defer this question to the later tutorials, but will explore the effect of the order on accuracy. The plots below show the convergence in mean and variance with different number of samples.