//=========================================================================== // GoTools Core - SINTEF Geometry Tools Core library, version 2.0.1 // // Copyright (C) 2000-2007, 2010 SINTEF ICT, Applied Mathematics, Norway. // // This program is free software; you can redistribute it and/or // modify it under the terms of the GNU General Public License // as published by the Free Software Foundation version 2 of the License. // // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with this program; if not, write to the Free Software // Foundation, Inc., // 59 Temple Place - Suite 330, // Boston, MA 02111-1307, USA. // // Contact information: E-mail: tor.dokken@sintef.no // SINTEF ICT, Department of Applied Mathematics, // P.O. Box 124 Blindern, // 0314 Oslo, Norway. // // Other licenses are also available for this software, notably licenses // for: // - Building commercial software. // - Building software whose source code you wish to keep private. //=========================================================================== 00015 #ifndef _INTEGRATION_H_ 00016 #define _INTEGRATION_H_ 00017 00018 00019 #include <math.h> 00020 #include <functional> 00021 #include <numeric> 00022 #include <limits> 00023 #include <stdexcept> 00024 #include "GoTools/utils/errormacros.h" 00025 00026 namespace Go { 00027 00028 00034 template <typename Functor> 00035 void trapezoidal(Functor& f, double a, double b, double& s, int n) 00036 { 00037 DEBUG_ERROR_IF(n < 1, "Bad function argument to trapezoidal()."); 00038 if (n == 1) { 00039 s = 0.5 * (b - a) * (f(a) + f(b)); 00040 } else { 00041 int num_int_pts = 1 << (n-2); 00042 //int num_int_pts = std::power(2, n-2); 00043 double spacing = (b - a) / num_int_pts; 00044 double x = a + 0.5 * spacing; 00045 double sum = 0.0; 00046 for (int i = 0; i < num_int_pts; ++i) { 00047 sum += f(x); 00048 x += spacing; 00049 } 00050 s += (b - a) * sum / num_int_pts; 00051 s *= 0.5; 00052 } 00053 } 00054 00055 00058 template <typename Functor> 00059 double simpsons_rule(Functor& f, double a, double b, 00060 const double eps = 1.0e-6, const int max_iter = 20) 00061 { 00062 const double one_third = double(1) / double(3); 00063 double result = 0; 00064 double tz = 0; 00065 double last_tz = std::numeric_limits<double>::min(); 00066 double last_result = std::numeric_limits<double>::min(); 00067 00068 for (int j = 1; j <= max_iter; ++j) { 00069 trapezoidal(f, a, b, tz, j); 00070 result = (4.0 * tz - last_tz) * one_third; 00071 if ((fabs(result - last_result) < eps * fabs(last_result)) || 00072 (fabs(result) < eps && fabs(last_result) < eps && j > 6)) { 00073 return result; 00074 } 00075 last_result = result; 00076 last_tz = tz; 00077 } 00078 MESSAGE("Too many steps in simpsons_rule."); 00079 throw std::runtime_error("Too many steps in simpsons_rule."); 00080 } // 00081 00084 template <typename Functor> 00085 double gaussian_quadrature(Functor& f, double a, double b) 00086 { 00087 static double x[] = { 0.1488743389, 00088 0.4333953941, 00089 0.6794095682, 00090 0.8650633666, 00091 0.9739065285 }; 00092 static double weight[] = { 0.2955242247, 00093 0.2692667193, 00094 0.2190863625, 00095 0.1494513491, 00096 0.0666713443 }; 00097 00098 double midpt = 0.5 * (b + a); 00099 double half_length = 0.5 * (b - a); 00100 double scaled_result = double(0); 00101 double step = double(0); 00102 for (int i = 0; i < 5; ++i) { 00103 step = half_length * x[i]; 00104 scaled_result += weight[i] * (f(midpt + step) + f(midpt - step)); 00105 } 00106 // rescale result to correspond to actual interval size 00107 return scaled_result * half_length; 00108 } 00109 00110 00111 00112 // Help functor to simpsons_rule2D and gaussian_quadrature2D 00114 template <typename Functor> 00115 class Integrate2ndFunctor { 00116 public: 00117 Integrate2ndFunctor(Functor& f, double a, double b) 00118 : f_(f), a_(a), b_(b) { } 00119 00120 double operator() (double t) const 00121 { 00122 std::binder1st<Functor> fu(f_, t); 00123 return gaussian_quadrature(fu, a_, b_); 00124 } 00125 00126 private: 00127 Functor f_; 00128 double a_, b_; 00129 00130 }; 00132 00135 template <typename Functor2D> 00136 double simpsons_rule2D(Functor2D& f, 00137 double ax, double bx, double ay, double by, 00138 const double eps = 1.0e-6, const int max_iter = 20) 00139 { 00140 Integrate2ndFunctor<Functor2D> fu(f, ay, by); 00141 return simpsons_rule(fu, ax, bx, eps, max_iter); 00142 } 00143 00144 00147 template <typename Functor2D> 00148 double gaussian_quadrature2D(Functor2D& f, 00149 double ax, double bx, double ay, double by) 00150 { 00151 Integrate2ndFunctor<Functor2D> fu(f, ay, by); 00152 return gaussian_quadrature(fu, ax, bx); 00153 } 00154 00155 00156 } // namespace GO 00157 00158 00159 #endif // _INTEGRATION_H 00160
Generated on Tue Sep 21 15:44:17 2010 for GoTools Core by  doxygen 1.6.3