Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_KL_OneDExponentialCovarianceFunctionImp.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
42#include "Teuchos_Assert.hpp"
43
44template <typename value_type>
47 const value_type& a,
48 const value_type& b,
49 const value_type& L_,
50 const int dim_name,
51 Teuchos::ParameterList& solverParams) :
52 L(L_),
53 eig_pair(M)
54{
55 // Get parameters with default values
56 magnitude_type eps = solverParams.get("Bound Perturbation Size", 1e-6);
57 magnitude_type tol = solverParams.get("Nonlinear Solver Tolerance", 1e-10);
58 int max_it = solverParams.get("Maximum Nonlinear Solver Iterations", 100);
59
60 value_type aa, alpha, omega, lambda;
61 int i=0;
62 double pi = 4.0*std::atan(1.0);
63 int idx = 0;
64
65 aa = (b-a)/2.0;
66 while (i < M-1) {
67 alpha = aa/L;
68 omega = bisection(EigFuncCos(alpha), idx*pi, idx*pi+pi/2.0-eps,
69 tol, max_it) / aa;
70 lambda = 2.0*L/(L*L*omega*omega + 1.0);
71 eig_pair[i].eig_val = lambda;
74 i++;
75
76 omega = bisection(EigFuncSin(alpha), idx*pi+pi/2.0+eps, (idx+1)*pi,
77 tol, max_it) / aa;
78 lambda = 2.0*L/(L*L*omega*omega + 1.0);
79 eig_pair[i].eig_val = lambda;
82 i++;
83
84 idx++;
85 }
86 if (i < M) {
87 omega = bisection(EigFuncCos(alpha), idx*pi, idx*pi+pi/2.0-eps,
88 tol, max_it) / aa;
89 lambda = 2.0*L/(L*L*omega*omega + 1.0);
90 eig_pair[i].eig_val = lambda;
93 }
94}
95
96template <typename value_type>
97template <class Func>
98value_type
100newton(const Func& func, const value_type& a, const value_type& b,
101 magnitude_type tol, int max_num_its)
102{
103 value_type u = (a+b)/2.0;
104 value_type f = func.eval(u);
105 int nit = 0;
106 while (Teuchos::ScalarTraits<value_type>::magnitude(f) > tol &&
107 nit < max_num_its) {
108 std::cout << "u = " << u << " f = " << f << std::endl;
109 value_type dfdu = func.deriv(u);
110 u -= f / dfdu;
111 f = func.eval(u);
112 ++nit;
113 }
114 TEUCHOS_TEST_FOR_EXCEPTION(nit >= max_num_its, std::logic_error,
115 "Nonlinear solver did not converge!" << std::endl);
116
117 return u;
118}
119
120template <typename value_type>
121template <class Func>
122value_type
124bisection(const Func& func, const value_type& a, const value_type& b,
125 magnitude_type tol, int max_num_its)
126{
127 value_type low, hi;
128 value_type fa = func.eval(a);
129 value_type fb = func.eval(b);
130 TEUCHOS_TEST_FOR_EXCEPTION(fa*fb > value_type(0.0), std::logic_error,
131 "Bounds [" << a << "," << b << "] must bracket the root!" << std::endl <<
132 "f(a) = " << fa << ", f(b) = " << fb << std::endl)
133
134 if (fa <= 0.0) {
135 low = a;
136 hi = b;
137 }
138 else {
139 low = b;
140 hi = a;
141 }
142
143 int nit = 0;
144 value_type u = low + (hi - low)/2.0;
145 value_type f = func.eval(u);
146 while ((Teuchos::ScalarTraits<value_type>::magnitude(hi - low) > 2.0*tol ||
147 Teuchos::ScalarTraits<value_type>::magnitude(f) > tol) &&
148 nit < max_num_its) {
149 //std::cout << "u = " << u << " f = " << f << std::endl;
150 if (f <= 0.0)
151 low = u;
152 else
153 hi = u;
154 u = low + (hi - low)/2.0;
155 f = func.eval(u);
156 ++nit;
157 }
158 TEUCHOS_TEST_FOR_EXCEPTION(nit >= max_num_its, std::logic_error,
159 "Nonlinear solver did not converge!" << std::endl);
160
161 return u;
162}
One-dimensional eigenfunction for exponential covariance function.
value_type bisection(const Func &func, const value_type &a, const value_type &b, magnitude_type tol, int max_num_its)
A basic root finder based on bisection.
Teuchos::ScalarTraits< value_type >::magnitudeType magnitude_type
value_type newton(const Func &func, const value_type &a, const value_type &b, magnitude_type tol, int max_num_its)
A basic root finder based on Newton's method.
OneDExponentialCovarianceFunction(int M, const value_type &a, const value_type &b, const value_type &L, const int dim_name, Teuchos::ParameterList &solverParams)
Constructor.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
Nonlinear function whose roots define eigenvalues for cos() eigenfunction.
Nonlinear function whose roots define eigenvalues for sin() eigenfunction.