All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Pages
gaussian_quadratur.h
1 /* Copyright (C) 2004-2015 MBSim Development Team
2  *
3  * This library is free software; you can redistribute it and/or
4  * modify it under the terms of the GNU Lesser General Public
5  * License as published by the Free Software Foundation; either
6  * version 2.1 of the License, or (at your option) any later version.
7  *
8  * This library is distributed in the hope that it will be useful,
9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11  * Lesser General Public License for more details.
12  *
13  * You should have received a copy of the GNU Lesser General Public
14  * License along with this library; if not, write to the Free Software
15  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
16  *
17  * Contact: martin.o.foerg@googlemail.com
18  */
19 
20 #ifndef GAUSSIAN_QUADRATUR_H_
21 #define GAUSSIAN_QUADRATUR_H_
22 
23 #include <fmatvec/fmatvec.h>
24 #include <mbsim/functions/function.h>
25 
26 #include <mbsim/mbsim_event.h>
27 
28 namespace MBSim {
29 
30  static double transformVariablePosition(double x, double l, double u) {
31  return (u - l) / 2. * x + (l + u) / 2.;
32  }
33 
39  static double pW[15][2] = {{0., 2.}, {-0.57735026919, 1.}, {0.57735026919, 1}, {-0.774596669241, 5. / 9.}, {0., 8. / 9}, {0.774596669241, 5. / 9.}, {-0.861136311594053, 0.347854845137454}, {-0.339981043584856, 0.652145154862546}, {0.339981043584856, 0.652145154862546}, {0.861136311594053, 0.347854845137454}, {-0.906179845938664, 0.236926885056189}, {-0.538469310105683, 0.478628670499366}, {0, 0.568888888888889}, {0.538469310105683, 0.478628670499366}, {0.906179845938664, 0.236926885056189}};
40 
41  template <class Ret>
43  public:
44  GaussLegendreQuadrature(fmatvec::Function<Ret(double)> * function, int dimension, int subIntegrals = 1, int Nb = 3) :
45  function(function), dim(dimension), subIntegrals(subIntegrals), Nb(Nb) {
46  }
48  }
49 
56  const Ret integrate(double l, double u) {
57  double currLow = l;
58  double step = (u - l) / subIntegrals;
59 
60  Ret result(dim, fmatvec::INIT, 0.);
61 
62  for (size_t i = 0; i < subIntegrals; i++) {
63  result += integrateSingle(currLow, currLow + step);
64  currLow += step;
65  }
66 
67  return result;
68  }
69 
70  /* GETTER/SETTER*/
71  void setSubIntegrals(int subInts_) {
72  subIntegrals = subInts_;
73  }
74 
79 
83  int dim;
84 
88  size_t subIntegrals;
89 
93  size_t Nb;
94 
95  protected:
96  const Ret integrateSingle(double l, double u) {
97  if (Nb > 5)
98  throw MBSim::MBSimError("It is not possible to use more than 5 Gauss-Points atm!");
99 
100  int posInpW = 0;
101 
102  for (unsigned int i = 1; i < Nb; i++) {
103  posInpW += i;
104  }
105 
106  Ret result(dim, fmatvec::INIT, 0.);
107  for (unsigned int i = 0; i < Nb; i++) {
108  int index = i + posInpW;
109  double p = transformVariablePosition(pW[index][0], l, u);
110  double w = pW[index][1];
111  result += w * (*function)(p);
112  }
113 
114  return (u - l) / 2. * result;
115  }
116  };
117 
118  namespace GaussTschebyschowQuadrature {
128  template <class Ret>
129  const Ret integrate(fmatvec::Function<Ret(double)> * function, double l, double u, size_t Nb = 3) {
130  Ret result((*function)(l).size(), fmatvec::INIT, 0.);
131  for (size_t i = 0; i < Nb; i++) {
132  double x = cos((2. * i - 1) / (2. * Nb) * M_PI);
133  double p = transformVariablePosition(x, l, u);
134  result += sqrt(1 - x * x) * (*function)(p);
135  }
136 
137  return (u - l) / 2. * M_PI / Nb * result;
138  }
139  }
140 }
141 
142 #endif /* GAUSSIAN_QUADRATUR_H_ */
int dim
dimension of the result
Definition: gaussian_quadratur.h:83
Definition: gaussian_quadratur.h:42
size_t Nb
number of gauss points that should be used
Definition: gaussian_quadratur.h:93
const Ret integrate(double l, double u)
function to integrate a MBSim-Function numerically using the gaussian quadrature
Definition: gaussian_quadratur.h:56
size_t subIntegrals
number of integrals the original one should be split to
Definition: gaussian_quadratur.h:88
basic error class for mbsim
Definition: mbsim_event.h:38

Impressum / Disclaimer / Datenschutz Generated by doxygen 1.8.5 Valid HTML