Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
recurrence_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// recurrence_example
45//
46// usage:
47// recurrence_example
48//
49// output:
50// Prints the recurrence coefficients for the first 5 normalized polynomials
51// orthogonal wrt the given weight. Follows up by printing the computed
52// norms and outputting a 11 point gaussian quadrature rule.
53// Demonstrate orthogonality by outputting the maximum computed
54// |<psi_i, psi_j>| for j != i.
55
56#include <iostream>
57#include <iomanip>
58
59#include "Stokhos.hpp"
60
61double weightFunction(const double& x){
62 if(2*fabs(x+2) < 1){
63 return exp(-1/(1-pow(2*(x+2),2)));
64 }else if(2*fabs(x-2)<1){
65 return exp(-1/(1-pow(2*(x-2),2)));
66 }else{
67 return 0;
68 }
69}
70
71int main(int argc, char **argv)
72{
73 try {
74 const int p = 5;
75 const double leftEndPt = -3;
76 const double rightEndPt = 3;
77
78 //Generate a basis up to order p with the support of the weight = [-5,5] using normalization.
79 Teuchos::RCP<const Stokhos::DiscretizedStieltjesBasis<int,double> > basis =
80 Teuchos::rcp(new Stokhos::DiscretizedStieltjesBasis<int,double>("Beta",p,&weightFunction,leftEndPt,rightEndPt,true));
81 Teuchos::Array<double> alpha;
82 Teuchos::Array<double> beta;
83 Teuchos::Array<double> delta;
84 Teuchos::Array<double> gamma;
85 Teuchos::Array<double> norm_sq;
86 norm_sq = basis->norm_squared();
87 basis->getRecurrenceCoefficients(alpha, beta, delta, gamma);
88
89 //Fetch alpha, beta, gamma and the computed norms and print.
90 for(int i = 0; i<p; i++){
91 std::cout << "alpha[" << i <<"]= " << alpha[i] << " beta[" << i <<"]= " << beta[i] << " gamma[" << i <<"]= " << gamma[i] << "\n";
92 }
93
94 for(int i = 0; i<=p; i++){
95 std::cout << "E(P_"<<i<<"^2) = "<< norm_sq[i] <<"\n";
96 }
97
98 Teuchos::Array<double> quad_points;
99 Teuchos::Array<double> quad_weights;
100 Teuchos::Array<Teuchos::Array< double > > quad_values;
101 basis->getQuadPoints(20, quad_points, quad_weights, quad_values);
102 for(int i = 0; i<quad_points.size(); i++){
103 std::cout << "x_i = "<<quad_points[i]<<" w_i = "<< quad_weights[i] <<" " << i << " / " << quad_points.size()<< "\n";
104 }
105
106 double maxOffDiag = 0;
107 double currentInnerProduct;
108 for(int i = 0; i<=p; i++){
109 for(int j = 0; j<i; j++){
110 currentInnerProduct = fabs(basis->eval_inner_product(i, j));
111 if(currentInnerProduct > maxOffDiag){
112 maxOffDiag = currentInnerProduct;
113 }
114 }
115 }
116 std::cout<<"Maximum Off Diagonal Inner Product is: " << maxOffDiag << "\n";
117 }
118
119 catch (std::exception& e) {
120 std::cout << e.what() << std::endl;
121 }
122
123
124}
Generates three-term recurrence using the Discretized Stieltjes procedure.
int main(int argc, char **argv)
double weightFunction(const double &x)