Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
pecos_hermite_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// pecos_hermite_example
45//
46// usage:
47// pecos_hermite_example
48//
49// output:
50// prints the Hermite Polynomial Chaos Expansion of the simple function
51//
52// v = 1/(log(u)^2+1)
53//
54// where u = 1 + 0.4*H_1(x) + 0.06*H_2(x) + 0.002*H_3(x), x is a zero-mean
55// and unit-variance Gaussian random variable, and H_i(x) is the i-th
56// Hermite polynomial.
57//
58// Same as hermite_example, except uses Pecos to define the Hermite basis.
59
60#include "Stokhos.hpp"
61#include "HermiteOrthogPolynomial.hpp" // from Pecos
62
63int main(int argc, char **argv)
64{
65 try {
66
67 // Basis of dimension 3, order 5
68 const int d = 3;
69 const int p = 5;
70 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
71 for (int i=0; i<d; i++) {
72 bases[i] = Teuchos::rcp(new Stokhos::PecosOneDOrthogPolyBasis<int,double>(Teuchos::rcp(new Pecos::HermiteOrthogPolynomial), "Hermite", p));
73 }
74 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
75 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
76
77 // Quadrature method
78 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
79 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
80
81 // Triple product tensor
82 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
83 basis->computeTripleProductTensor();
84
85 // Expansion method
87
88 // Polynomial expansions
89 Stokhos::OrthogPolyApprox<int,double> u(basis), v(basis), w(basis);
90 u.term(0,0) = 1.0;
91 for (int i=0; i<d; i++) {
92 u.term(i,1) = 0.4 / d;
93 u.term(i,2) = 0.06 / d;
94 u.term(i,3) = 0.002 / d;
95 }
96
97 // Compute expansion
98 expn.log(v,u);
99 expn.times(w,v,v);
100 expn.plusEqual(w,1.0);
101 expn.divide(v,1.0,w);
102 //expn.times(v,u,u);
103
104 // Print u and v
105 std::cout << "v = 1.0 / (log(u)^2 + 1):" << std::endl;
106 std::cout << "\tu = ";
107 u.print(std::cout);
108 std::cout << "\tv = ";
109 v.print(std::cout);
110
111 // Compute moments
112 double mean = v.mean();
113 double std_dev = v.standard_deviation();
114
115 // Evaluate PCE and function at a point = 0.25 in each dimension
116 Teuchos::Array<double> pt(d);
117 for (int i=0; i<d; i++)
118 pt[i] = 0.25;
119 double up = u.evaluate(pt);
120 double vp = 1.0/(std::log(up)*std::log(up)+1.0);
121 double vp2 = v.evaluate(pt);
122
123 // Print results
124 std::cout << "\tv mean = " << mean << std::endl;
125 std::cout << "\tv std. dev. = " << std_dev << std::endl;
126 std::cout << "\tv(0.25) (true) = " << vp << std::endl;
127 std::cout << "\tv(0.25) (pce) = " << vp2 << std::endl;
128
129 // Check the answer
130 if (std::abs(vp - vp2) < 1e-2)
131 std::cout << "\nExample Passed!" << std::endl;
132 }
133 catch (std::exception& e) {
134 std::cout << e.what() << std::endl;
135 }
136}
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Class to store coefficients of a projection onto an orthogonal polynomial basis.
void plusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
Orthogonal polynomial expansions based on numerical quadrature.
void times(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
void divide(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
void log(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
int main(int argc, char **argv)