Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_DiscretizedStieltjesBasisImp.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42template <typename ordinal_type, typename value_type>
44DiscretizedStieltjesBasis(const std::string& label,
45 const ordinal_type& p,
46 value_type (*weightFn)(const value_type&),
47 const value_type& leftEndPt,
48 const value_type& rightEndPt,
49 bool normalize,
50 Stokhos::GrowthPolicy growth) :
51 RecurrenceBasis<ordinal_type,value_type>(std::string("DiscretizedStieltjes -- ") + label, p, normalize, growth),
52 scaleFactor(1),
53 leftEndPt_(leftEndPt),
54 rightEndPt_(rightEndPt),
55 weightFn_(weightFn)
56{
57 // Set up quadrature points for discretized stieltjes procedure
58 Teuchos::RCP<const Stokhos::LegendreBasis<ordinal_type,value_type> > quadBasis =
59 Teuchos::rcp(new Stokhos::LegendreBasis<ordinal_type,value_type>(this->p));
60 quadBasis->getQuadPoints(200*this->p, quad_points, quad_weights, quad_values);
61
62 // Setup rest of recurrence basis
63 this->setup();
64}
65
66template <typename ordinal_type, typename value_type>
68DiscretizedStieltjesBasis(const ordinal_type& p,
69 const DiscretizedStieltjesBasis& basis) :
70 RecurrenceBasis<ordinal_type,value_type>(p, basis),
71 scaleFactor(basis.scaleFactor),
72 leftEndPt_(basis.leftEndPt_),
73 rightEndPt_(basis.rightEndPt_),
74 weightFn_(basis.weightFn_)
75{
76 // Set up quadrature points for discretized stieltjes procedure
77 Teuchos::RCP<const Stokhos::LegendreBasis<ordinal_type,value_type> > quadBasis =
78 Teuchos::rcp(new Stokhos::LegendreBasis<ordinal_type,value_type>(this->p));
79 quadBasis->getQuadPoints(200*this->p, quad_points, quad_weights, quad_values);
80
81 // Compute coefficients in 3-term recurrsion
82 computeRecurrenceCoefficients(p+1, this->alpha, this->beta, this->delta,
83 this->gamma);
84
85 // Setup rest of recurrence basis
86 this->setup();
87}
88
89template <typename ordinal_type, typename value_type>
92{
93}
94
95template <typename ordinal_type, typename value_type>
96bool
99 Teuchos::Array<value_type>& alpha,
100 Teuchos::Array<value_type>& beta,
101 Teuchos::Array<value_type>& delta,
102 Teuchos::Array<value_type>& gamma) const
103{
104 //The Discretized Stieltjes polynomials are defined by a recurrance relation,
105 //P_n+1 = \gamma_n+1[(x-\alpha_n) P_n - \beta_n P_n-1].
106 //The alpha and beta coefficients are generated first using the
107 //discritized stilges procidure described in "On the Calculation of DiscretizedStieltjes Polynomials and Quadratures",
108 //Robin P. Sagar, Vedene H. Smith. The gamma coefficients are then optionally set so that each
109 //polynomial has norm 1. If normalization is not enabled then the gammas are set to 1.
110
111 scaleFactor = 1;
112
113 //First renormalize the weight function so that it has measure 1.
114 value_type oneNorm = expectedValue_J_nsquared(0, alpha, beta);
115 //future evaluations of the weight function will scale it by this factor.
116 scaleFactor = 1/oneNorm;
117
118
119 value_type integral2;
120 //NOTE!! This evaluation of 'expectedValue_J_nsquared(0)' is different
121 //from the one above since we rescaled the weight. Don't combine
122 //the two!!!
123
124 value_type past_integral = expectedValue_J_nsquared(0, alpha, beta);
125 alpha[0] = expectedValue_tJ_nsquared(0, alpha, beta)/past_integral;
126 //beta[0] := \int_-c^c w(x) dx.
127 beta[0] = 1;
128 delta[0] = 1;
129 gamma[0] = 1;
130 //These formulas are from the above reference.
131 for (ordinal_type n = 1; n<k; n++){
132 integral2 = expectedValue_J_nsquared(n, alpha, beta);
133 alpha[n] = expectedValue_tJ_nsquared(n, alpha, beta)/integral2;
134 beta[n] = integral2/past_integral;
135 past_integral = integral2;
136 delta[n] = 1.0;
137 gamma[n] = 1.0;
138 }
139
140 return false;
141}
142
143template <typename ordinal_type, typename value_type>
144value_type
146evaluateWeight(const value_type& x) const
147{
148 return (x < leftEndPt_ || x > rightEndPt_) ? 0: scaleFactor*weightFn_(x);
149}
150
151template <typename ordinal_type, typename value_type>
152value_type
154expectedValue_tJ_nsquared(const ordinal_type& order,
155 const Teuchos::Array<value_type>& alpha,
156 const Teuchos::Array<value_type>& beta) const
157{
158 //Impliments a gaussian quadrature routine to evaluate the integral,
159 // \int_-c^c J_n(x)^2w(x)dx. This is needed to compute the recurrance coefficients.
160 value_type integral = 0;
161 for(ordinal_type quadIdx = 0;
162 quadIdx < static_cast<ordinal_type>(quad_points.size()); quadIdx++) {
163 value_type x = (rightEndPt_ - leftEndPt_)*.5*quad_points[quadIdx] +
164 (rightEndPt_ + leftEndPt_)*.5;
165 value_type val = evaluateRecurrence(x,order,alpha,beta);
166 integral += x*val*val*evaluateWeight(x)*quad_weights[quadIdx];
167 }
168
169 return integral*(rightEndPt_ - leftEndPt_);
170}
171
172template <typename ordinal_type, typename value_type>
173value_type
175expectedValue_J_nsquared(const ordinal_type& order,
176 const Teuchos::Array<value_type>& alpha,
177 const Teuchos::Array<value_type>& beta) const
178{
179 //Impliments a gaussian quadrature routineroutine to evaluate the integral,
180 // \int_-c^c J_n(x)^2w(x)dx. This is needed to compute the recurrance coefficients.
181 value_type integral = 0;
182 for(ordinal_type quadIdx = 0;
183 quadIdx < static_cast<ordinal_type>(quad_points.size()); quadIdx++){
184 value_type x = (rightEndPt_ - leftEndPt_)*.5*quad_points[quadIdx] +
185 (rightEndPt_ + leftEndPt_)*.5;
186 value_type val = evaluateRecurrence(x,order,alpha,beta);
187 integral += val*val*evaluateWeight(x)*quad_weights[quadIdx];
188 }
189
190 return integral*(rightEndPt_ - leftEndPt_);
191}
192
193template <typename ordinal_type, typename value_type>
194value_type
196eval_inner_product(const ordinal_type& order1, const ordinal_type& order2) const
197{
198 //Impliments a gaussian quadrature routine to evaluate the integral,
199 // \int_-c^c J_n(x)J_m w(x)dx. This method is intended to allow the user to
200 // test for orthogonality and proper normalization.
201 value_type integral = 0;
202 for(ordinal_type quadIdx = 0;
203 quadIdx < static_cast<ordinal_type>(quad_points.size()); quadIdx++){
204 value_type x = (rightEndPt_ - leftEndPt_)*.5*quad_points[quadIdx] +
205 (rightEndPt_ + leftEndPt_)*.5;
206 integral += this->evaluate(x,order1)*this->evaluate(x,order2)*evaluateWeight(x)*quad_weights[quadIdx];
207 }
208
209 return integral*(rightEndPt_ - leftEndPt_);
210}
211
212template <typename ordinal_type, typename value_type>
213value_type
215evaluateRecurrence(const value_type& x,
216 ordinal_type k,
217 const Teuchos::Array<value_type>& alpha,
218 const Teuchos::Array<value_type>& beta) const
219{
220 if (k == 0)
221 return value_type(1.0);
222 else if (k == 1)
223 return x-alpha[0];
224
225 value_type v0 = value_type(1.0);
226 value_type v1 = x-alpha[0]*v0;
227 value_type v2 = value_type(0.0);
228 for (ordinal_type i=2; i<=k; i++) {
229 v2 = (x-alpha[i-1])*v1 - beta[i-1]*v0;
230 v0 = v1;
231 v1 = v2;
232 }
233
234 return v2;
235}
236
237template <typename ordinal_type, typename value_type>
238Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
240cloneWithOrder(ordinal_type p) const
241{
242 return Teuchos::rcp(new Stokhos::DiscretizedStieltjesBasis<ordinal_type,value_type>(p,*this));
243}
expr val()
Generates three-term recurrence using the Discretized Stieltjes procedure.
value_type expectedValue_tJ_nsquared(const ordinal_type &order, const Teuchos::Array< value_type > &alpha, const Teuchos::Array< value_type > &beta) const
Approximates where = order.
Teuchos::Array< value_type > quad_weights
Quadrature points for discretized stieltjes procedure.
virtual bool computeRecurrenceCoefficients(ordinal_type n, Teuchos::Array< value_type > &alpha, Teuchos::Array< value_type > &beta, Teuchos::Array< value_type > &delta, Teuchos::Array< value_type > &gamma) const
Compute recurrence coefficients.
Teuchos::Array< Teuchos::Array< value_type > > quad_values
Quadrature values for discretized stieltjes procedure.
Teuchos::Array< value_type > quad_points
Quadrature points for discretized stieltjes procedure.
DiscretizedStieltjesBasis(const std::string &name, const ordinal_type &p, value_type(*weightFn)(const value_type &), const value_type &leftEndPt, const value_type &rightEndPt, bool normalize, GrowthPolicy growth=SLOW_GROWTH)
Constructor.
value_type eval_inner_product(const ordinal_type &order1, const ordinal_type &order2) const
Evaluate inner product of two basis functions to test orthogonality.
value_type evaluateWeight(const value_type &x) const
Evaluates the scaled weight function.
value_type evaluateRecurrence(const value_type &point, ordinal_type order, const Teuchos::Array< value_type > &alpha, const Teuchos::Array< value_type > &beta) const
Evaluate 3-term recurrence formula using supplied coefficients.
value_type expectedValue_J_nsquared(const ordinal_type &order, const Teuchos::Array< value_type > &alpha, const Teuchos::Array< value_type > &beta) const
Approximates where = order.
virtual Teuchos::RCP< OneDOrthogPolyBasis< ordinal_type, value_type > > cloneWithOrder(ordinal_type p) const
Clone this object with the option of building a higher order basis.
Legendre polynomial basis.
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship:
Teuchos::Array< value_type > alpha
Recurrence coefficients.
Teuchos::Array< value_type > beta
Recurrence coefficients.
ordinal_type p
Order of basis.
Teuchos::Array< value_type > gamma
Recurrence coefficients.
Teuchos::Array< value_type > delta
Recurrence coefficients.
virtual void setup()
Setup basis after computing recurrence coefficients.
GrowthPolicy
Enumerated type for determining Smolyak growth policies.