Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
exp_moment_example.cpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
44#include <iostream>
45#include <iomanip>
46
47#include "Stokhos.hpp"
48
50
51// This example uses PC expansions for computing moments of
52//
53// u = exp(x1 + ... + xd)
54//
55// where x1, ..., xd are uniform random variables on [-1,1]. The methods
56// are compared to computing the "true" value via very high-order quadrature.
57// Because of the structure of the exponential, the moments can easily
58// be computed in one dimension.
59
60int main(int argc, char **argv)
61{
62 try {
63
64 // Compute "true" 1-D mean, std. dev using quadrature
65 const unsigned int true_quad_order = 200;
66 basis_type tmp_basis(1);
67 Teuchos::Array<double> true_quad_points, true_quad_weights;
68 Teuchos::Array< Teuchos::Array<double> > true_quad_values;
69 tmp_basis.getQuadPoints(true_quad_order, true_quad_points,
70 true_quad_weights, true_quad_values);
71 double mean_1d = 0.0;
72 double sd_1d = 0.0;
73 for (unsigned int qp=0; qp<true_quad_points.size(); qp++) {
74 double t = std::exp(true_quad_points[qp]);
75 mean_1d += t*true_quad_weights[qp];
76 sd_1d += t*t*true_quad_weights[qp];
77 }
78
79 const unsigned int dmin = 1;
80 const unsigned int dmax = 4;
81 const unsigned int pmin = 1;
82 const unsigned int pmax = 5;
83
84 // Loop over dimensions
85 for (unsigned int d=dmin; d<=dmax; d++) {
86
87 // compute "true" values
88 double true_mean = std::pow(mean_1d,static_cast<double>(d));
89 double true_sd = std::pow(sd_1d,static_cast<double>(d)) -
90 true_mean*true_mean;
91 true_sd = std::sqrt(true_sd);
92 std::cout.precision(12);
93 std::cout << "true mean = " << true_mean << "\t true std. dev. = "
94 << true_sd << std::endl;
95
96 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
97
98 // Loop over orders
99 for (unsigned int p=pmin; p<=pmax; p++) {
100
101 // Create product basis
102 for (unsigned int i=0; i<d; i++)
103 bases[i] = Teuchos::rcp(new basis_type(p));
104 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
105 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
106
107 // Create approximation
108 int sz = basis->size();
109 Stokhos::OrthogPolyApprox<int,double> x(basis), u(basis);
110 for (unsigned int i=0; i<d; i++) {
111 x.term(i,1) = 1.0;
112 }
113
114 // Tensor product quadrature
115 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
116 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
117
118 // Triple product tensor
119 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
120 basis->computeTripleProductTensor();
121
122 // Quadrature expansion
124 quad);
125
126 // Compute PCE via quadrature expansion
127 quad_exp.exp(u,x);
128 double mean = u.mean();
129 double sd = u.standard_deviation();
130
131 std::cout.precision(4);
132 std::cout.setf(std::ios::scientific);
133 std::cout << "d = " << d << " p = " << p
134 << " sz = " << sz
135 << "\tmean err = "
136 << std::fabs(true_mean-mean) << "\tstd. dev. err = "
137 << std::fabs(true_sd-sd) << std::endl;
138 }
139
140 }
141 }
142 catch (std::exception& e) {
143 std::cout << e.what() << std::endl;
144 }
145}
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
value_type mean() const
Compute mean of expansion.
value_type standard_deviation() const
Compute standard deviation of expansion.
Orthogonal polynomial expansions based on numerical quadrature.
void exp(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
virtual void getQuadPoints(ordinal_type quad_order, Teuchos::Array< value_type > &points, Teuchos::Array< value_type > &weights, Teuchos::Array< Teuchos::Array< value_type > > &values) const
Compute quadrature points, weights, and values of basis polynomials at given set of points points.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
int main(int argc, char **argv)
Stokhos::LegendreBasis< int, double > basis_type