.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/surrogates/pce/pce_euler_UQ.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_surrogates_pce_pce_euler_UQ.py: Physics-Informed PCE: Uncertainty Quantification of Euler Beam ====================================================================== In this example, we use Physics-informed Polynomial Chaos Expansion for approximation and UQ of Euler beam equation. .. GENERATED FROM PYTHON SOURCE LINES 10-11 Import necessary libraries. .. GENERATED FROM PYTHON SOURCE LINES 14-27 .. code-block:: default from UQpy.distributions import Uniform, JointIndependent from UQpy.surrogates import * # load PC^2 from UQpy.surrogates.polynomial_chaos.polynomials.TotalDegreeBasis import TotalDegreeBasis from UQpy.surrogates.polynomial_chaos.physics_informed.ConstrainedPCE import ConstrainedPCE from UQpy.surrogates.polynomial_chaos.physics_informed.PdeData import PdeData from UQpy.surrogates.polynomial_chaos.physics_informed.PdePCE import PdePCE from UQpy.surrogates.polynomial_chaos.physics_informed.Utilities import * from UQpy.surrogates.polynomial_chaos.physics_informed.ReducedPCE import ReducedPCE .. GENERATED FROM PYTHON SOURCE LINES 28-29 We then define functions used for construction of Physically Constrained PCE (PC :math:`^2`) .. GENERATED FROM PYTHON SOURCE LINES 32-95 .. code-block:: default # Definition of PDE/ODE def pde_func(standardized_sample, pce): der_order = 4 deriv_pce = derivative_basis(standardized_sample, pce, derivative_order=der_order, leading_variable=0) * ((2 / 1) ** der_order) pde_basis = deriv_pce return pde_basis # Definition of the source term def pde_res(standardized_sample): load_s = load(standardized_sample) return -load_s[:, 0] def load(standardized_sample): return const_load(standardized_sample) def const_load(standardized_sample): l = (1 + (1 + standardized_sample[:, 1]) / 2).reshape(-1, 1) return l # Definition of the function for sampling of boundary conditions def bc_sampling(nsim=1000): # BC sampling nsim_half = round(nsim / 2) sample = np.zeros((nsim, 2)) real_ogrid_1d = ortho_grid(nsim_half, 1, 0.0, 1.0)[:, 0] sample[:nsim_half, 0] = np.zeros(nsim_half) sample[:nsim_half, 1] = real_ogrid_1d sample[nsim_half:, 0] = np.ones(nsim_half) sample[nsim_half:, 1] = real_ogrid_1d return sample # define sampling and evaluation of BC for estimation of error def bc_res(nsim, pce): physical_sample= np.zeros((2, 2)) physical_sample[1, 0] = 1 standardized_sample = polynomial_chaos.Polynomials.standardize_sample(physical_sample, pce.polynomial_basis.distributions) der_order = 2 deriv_pce = np.sum( derivative_basis(standardized_sample, pce, derivative_order=der_order, leading_variable=0) * ((2 / 1) ** der_order) * np.array(pce.coefficients).T, axis=1) return deriv_pce # Definition of the reference solution for an error estimation def ref_sol(physical_coordinate, q): return (q + 1) * (-(physical_coordinate ** 4) / 24 + physical_coordinate ** 3 / 12 - physical_coordinate / 24) .. GENERATED FROM PYTHON SOURCE LINES 96-98 Beam equation is parametrized by one deterministic input variable (coordinate :math:`x`) and random load intensity :math:`q`, both modeled as Uniform random variables. .. GENERATED FROM PYTHON SOURCE LINES 99-111 .. code-block:: default # number of input variables nrand = 1 nvar = 1 + nrand # definition of a joint probability distribution dist1 = Uniform(loc=0, scale=1) dist2 = Uniform(loc=0, scale=1) marg = [dist1, dist2] joint = JointIndependent(marginals=marg) .. GENERATED FROM PYTHON SOURCE LINES 112-113 The physical domain is defined by :math:`x \in [0,1]` .. GENERATED FROM PYTHON SOURCE LINES 114-120 .. code-block:: default # define geometry of physical space geometry_xmin = [0] geometry_xmax = [1] .. GENERATED FROM PYTHON SOURCE LINES 121-124 Boundary conditions are prescribed as .. math:: \frac{d^2 u(0,q)}{d x^2}=0, \qquad \frac{d^2 u(1,q)}{d x^2}=0, \qquad u(0,q)=0, \qquad u(1,q)=0 .. GENERATED FROM PYTHON SOURCE LINES 125-141 .. code-block:: default # number of BC samples nbc = 2 * 10 # derivation orders of prescribed BCs der_orders = [0, 2] # normals associated to precribed BCs bc_normals = [0, 0] # sampling of BC points bc_xtotal = bc_sampling(20) bc_ytotal = np.zeros(len(bc_xtotal)) bc_x = [bc_xtotal, bc_xtotal] bc_y = [bc_ytotal, bc_ytotal] .. GENERATED FROM PYTHON SOURCE LINES 142-143 Here we construct an object containing general physical information (geometry, boundary conditions) .. GENERATED FROM PYTHON SOURCE LINES 144-146 .. code-block:: default pde_data = PdeData(geometry_xmax, geometry_xmin, der_orders, bc_normals, bc_x, bc_y) .. GENERATED FROM PYTHON SOURCE LINES 147-148 Further we construct an object containing PDE physical data and PC :math:`^2` definitions of PDE .. GENERATED FROM PYTHON SOURCE LINES 149-152 .. code-block:: default pde_pce = PdePCE(pde_data, pde_func, pde_source=pde_res, boundary_conditions_evaluate=bc_res) .. GENERATED FROM PYTHON SOURCE LINES 153-154 Get dirichlet boundary conditions from PDE data object that are further used for construction of initial PCE object. .. GENERATED FROM PYTHON SOURCE LINES 155-160 .. code-block:: default dirichlet_bc = pde_data.dirichlet x_train = dirichlet_bc[:, :-1] y_train = dirichlet_bc[:, -1] .. GENERATED FROM PYTHON SOURCE LINES 161-163 Finally, PC :math:`^2` is constructed using Karush-Kuhn-Tucker normal equations solved by ordinary least squares and Least Angle Regression. .. GENERATED FROM PYTHON SOURCE LINES 164-197 .. code-block:: default least_squares = LeastSquareRegression() p = 9 PCEorder = p polynomial_basis = TotalDegreeBasis(joint, PCEorder, hyperbolic=1) # create initial PCE object containing basis, regression method and dirichlet BC initpce = PolynomialChaosExpansion(polynomial_basis=polynomial_basis, regression_method=least_squares) initpce.set_data(x_train, y_train) # construct a PC^2 object combining pde_data, pde_pce and initial PCE objects pcpc = ConstrainedPCE(pde_data, pde_pce, initpce) # get coefficients of PC^2 by least angle regression pcpc.lar() # get coefficients of PC^2 by ordinary least squares pcpc.ols() # evaluate errors of approximations real_ogrid = ortho_grid(100, nvar, 0.0, 1.0) yy_val_pce = pcpc.lar_pce.predict(real_ogrid).flatten() yy_val_pce_ols = pcpc.initial_pce.predict(real_ogrid).flatten() yy_val_true = ref_sol(real_ogrid[:, 0], real_ogrid[:, 1]).flatten() err = np.abs(yy_val_pce - yy_val_true) tot_err = np.sum(err) print('\nTotal approximation error by PC^2-LAR: ', tot_err) err_ols = np.abs(yy_val_pce_ols - yy_val_true) tot_err_ols = np.sum(err_ols) print('Total approximation error by PC^2-OLS: ', tot_err_ols) .. GENERATED FROM PYTHON SOURCE LINES 198-199 Once the PC :math:`^2` is constructed, we use ReducedPce class to filter out influence of deterministic variable in UQ .. GENERATED FROM PYTHON SOURCE LINES 200-203 .. code-block:: default reduced_pce = ReducedPCE(pcpc.lar_pce, n_deterministic=1) .. GENERATED FROM PYTHON SOURCE LINES 204-205 ReducedPce is used for estimation of local statistical moments and quantiles :math:`\pm2\sigma` .. GENERATED FROM PYTHON SOURCE LINES 206-250 .. code-block:: default coeff_res = [] var_res = [] mean_res = [] vartot_res = [] lower_quantiles_modes = [] upper_quantiles_modes = [] n_derivations = 4 sigma_mult = 2 beam_x = np.arange(0, 101) / 100 for x in beam_x: mean = np.zeros(n_derivations + 1) var = np.zeros(n_derivations + 1) variances = np.zeros((n_derivations + 1, nrand)) lq = np.zeros((1 + n_derivations, nrand)) uq = np.zeros((1 + n_derivations, nrand)) for d in range(1 + n_derivations): if d == 0: coeff = (reduced_pce.evaluate_coordinate(np.array(x), return_coefficients=True)) else: coeff = reduced_pce.derive_coordinate(np.array(x), derivative_order=d, leading_variable=0, return_coefficients=True, derivative_multiplier=transformation_multiplier(pde_data, 0, d)) mean[d] = coeff[0] var[d] = np.sum(coeff[1:] ** 2) variances[d, :] = reduced_pce.variance_contributions(coeff) for e in range(nrand): lq[d, e] = mean[d] + sigma_mult * np.sqrt(np.sum(variances[d, :e + 1])) uq[d, e] = mean[d] - sigma_mult * np.sqrt(np.sum(variances[d, :e + 1])) lower_quantiles_modes.append(lq) upper_quantiles_modes.append(uq) var_res.append(variances) mean_res.append(mean) vartot_res.append(var) mean_res = np.array(mean_res) vartot_res = np.array(vartot_res) var_res = np.array(var_res) lower_quantiles_modes = np.array(lower_quantiles_modes) upper_quantiles_modes = np.array(upper_quantiles_modes) .. GENERATED FROM PYTHON SOURCE LINES 251-252 PC :math:`^2` and corresponding Reduced PCE can be further used for local UQ. Obtained results are depicted in figure. .. GENERATED FROM PYTHON SOURCE LINES 253-294 .. code-block:: default # plot results import matplotlib.pyplot as plt fig, ax = plt.subplots(1, 4, figsize=(14.5, 3)) colors = ['black', 'blue', 'green', 'red'] for d in range(2, 1 + n_derivations): ax[d - 1].plot(beam_x, mean_res[:, d], color=colors[d - 1]) ax[d - 1].plot(beam_x, mean_res[:, d] + sigma_mult * np.sqrt(vartot_res[:, d]), color=colors[d - 1]) ax[d - 1].plot(beam_x, mean_res[:, d] - sigma_mult * np.sqrt(vartot_res[:, d]), color=colors[d - 1]) ax[d - 1].plot(beam_x, lower_quantiles_modes[:, d, 0], '--', alpha=0.7, color=colors[d - 1]) ax[d - 1].plot(beam_x, upper_quantiles_modes[:, d, 0], '--', alpha=0.7, color=colors[d - 1]) ax[d - 1].fill_between(beam_x, lower_quantiles_modes[:, d, 0], upper_quantiles_modes[:, d, 0], facecolor=colors[d - 1], alpha=0.05) ax[d - 1].plot(beam_x, np.zeros(len(beam_x)), color='black') ax[0].plot(beam_x, mean_res[:, 0], color=colors[0]) ax[0].plot(beam_x, mean_res[:, 0] + sigma_mult * np.sqrt(vartot_res[:, 0]), color=colors[0]) ax[0].plot(beam_x, mean_res[:, 0] - sigma_mult * np.sqrt(vartot_res[:, 0]), color=colors[0]) ax[0].plot(beam_x, lower_quantiles_modes[:, 0, 0], '--', alpha=0.7, color=colors[0]) ax[0].plot(beam_x, upper_quantiles_modes[:, 0, 0], '--', alpha=0.7, color=colors[0]) ax[0].fill_between(beam_x, lower_quantiles_modes[:, 0, 0], upper_quantiles_modes[:, 0, 0], facecolor=colors[0], alpha=0.05) ax[0].plot(beam_x, np.zeros(len(beam_x)), color='black') ax[3].invert_yaxis() ax[1].invert_yaxis() ax[3].set_title(r'Load $\frac{\partial^4 w}{\partial x^4}$', y=1.04) ax[2].set_title(r'Shear Force $\frac{\partial^3 w}{\partial x^3}$', y=1.04) ax[1].set_title(r'Bending Moment $\frac{\partial^2 w}{\partial x^2}$', y=1.04) ax[0].set_title(r'Deflection $w$', y=1.04) for axi in ax.flatten(): axi.set_xlabel(r'$x$') plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.000 seconds) .. _sphx_glr_download_auto_examples_surrogates_pce_pce_euler_UQ.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/SURGroup/UQpy/master?urlpath=lab/tree/notebooks/auto_examples/surrogates/pce/pce_euler_UQ.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: pce_euler_UQ.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: pce_euler_UQ.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_