Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JQuadrature.cc
Go to the documentation of this file.
1
2#include <string>
3#include <iostream>
4#include <iomanip>
5#include <cmath>
6
8
9#include "JMath/JPolynome.hh"
10#include "JMath/JGauss.hh"
11
12#include "Jeep/JParser.hh"
13#include "Jeep/JMessage.hh"
14
15
16namespace GAUSS_LEGENDRE {
17
18 using namespace JPP;
19
20 const int numberOfPoints = 10;
21 const double epsilon = 1.0e-12;
22
23 const double xmin = -5.0;
24 const double xmax = +3.0;
25
26 /**
27 * Function.
28 */
29 const JPolynome f1(1.0, 2.0, 3.0);
30
31 /**
32 * Integral.
33 */
35}
36
37namespace GAUSS_LAGUERRE {
38
39 using namespace JPP;
40
41 const int numberOfPoints = 10;
42 const double a = 1.0;
43 const double epsilon = 1.0e-12;
44
45 const double xmin = 0.0;
46 const double xmax = +1.0e10; // approximation of infinity
47
48 /**
49 * Function.
50 */
51 inline double f1(const double x)
52 {
53 return sin(x);
54 }
55
56 /**
57 * Integral.
58 */
59 inline double F1(const double x)
60 {
61 // integral of sin(x) * x^1 * exp(-x)
62
63 return -0.5 * exp(-x) * (x*sin(x) + (x + 1.0)*cos(x));
64 }
65}
66
67namespace GAUSS_HERMITE {
68
69 using namespace JPP;
70
71 const int numberOfPoints = 10;
72 const double epsilon = 1.0e-12;
73
74 const double sigma[] = { 4.0, 0.5 };
75
76 const double xmin = -3.0 * sigma[0];
77 const double xmax = +3.0 * sigma[0];
78
79 /**
80 * Function.
81 */
82 const JGauss f1(0.0, sigma[0]);
83
84 /**
85 * Integral.
86 */
87 const JGauss F1(0.0, hypot(sigma[0], sigma[1]));
88}
89
90
91/**
92 * \file
93 *
94 * Example program to test JTOOLS::JQuadrature.
95 * \author mdejong
96 */
97int main(int argc, char **argv)
98{
99 using namespace std;
100 using namespace JPP;
101
102 double precision;
103 int debug;
104
105 try {
106
107 JParser<> zap("Example program to test JTOOLS::JQuadrature.");
108
109 zap['e'] = make_field(precision) = 1.0e-4;
110 zap['d'] = make_field(debug) = 3;
111
112 zap(argc, argv);
113 }
114 catch(const exception &error) {
115 FATAL(error.what() << endl);
116 }
117
118 {
119 using namespace GAUSS_LEGENDRE;
120
121 double V = F1(xmax) - F1(xmin);
122
123 JGaussLegendre engine(numberOfPoints, epsilon);
124
125 const double x0 = 0.5 * (xmax + xmin);
126 const double x1 = 0.5 * (xmax - xmin);
127
128 double W = 0.0;
129
130 for (JGaussLegendre::const_iterator i = engine.begin(); i != engine.end(); ++i) {
131
132 const double x = x0 + i->getX() * x1;
133
134 W += f1(x) * i->getY() * x1;
135 }
136
137 DEBUG("integral analytical " << FIXED(9,5) << V << endl);
138 DEBUG("integral numerical " << FIXED(9,5) << W << endl);
139
140 ASSERT(fabs(V - W) <= precision, "Test Gauss-Legendre integration.");
141 }
142 {
143 using namespace GAUSS_LAGUERRE;
144
145 const double V = F1(xmax) - F1(xmin);
146
147 JGaussLaguerre engine(numberOfPoints, a, epsilon);
148
149 double W = 0.0;
150
151 for (JGaussLaguerre::const_iterator i = engine.begin(); i != engine.end(); ++i) {
152
153 const double x = i->getX();
154
155 W += f1(x) * i->getY();
156 }
157
158 DEBUG("integral analytical " << FIXED(9,5) << V << endl);
159 DEBUG("integral numerical " << FIXED(9,5) << W << endl);
160
161 ASSERT(fabs(V - W) <= precision, "Test Gauss-Laguerre integration.");
162 }
163 {
164 using namespace GAUSS_HERMITE;
165
166 JGaussHermite engine(numberOfPoints, epsilon);
167
168 for (double x = xmin; x <= xmax; x += 0.05*(xmax - xmin)) {
169
170 const double V = F1(x);
171
172 double W = 0.0;
173
174 for (JGaussHermite::const_iterator i = engine.begin(); i != engine.end(); ++i) {
175
176 const double u = i->getX() * sigma[1] * sqrt(2.0);
177 const double v = i->getY() / sqrt(PI);
178
179 W += f1(x + u) * v;
180 }
181
182 DEBUG(FIXED(7,3) << x << ' ' << FIXED(9,5) << V << ' ' << FIXED(9,5) << W << endl);
183
184 ASSERT(fabs(V - W) <= precision * V, "Test Gauss-Hermite integration.");
185 }
186 }
187
188 return 0;
189}
190
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define ASSERT(A,...)
Assert macro.
Definition JMessage.hh:90
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
int main(int argc, char **argv)
Auxiliary classes for numerical integration.
int numberOfPoints
Definition JResultPDF.cc:22
Utility class to parse command line options.
Definition JParser.hh:1698
container_type::const_iterator const_iterator
Numerical integrator for .
Numerical integrator for .
Numerical integrator for .
const JGauss F1(0.0, hypot(sigma[0], sigma[1]))
Integral.
const double xmin
const double epsilon
const double sigma[]
const double xmax
const int numberOfPoints
const JGauss f1(0.0, sigma[0])
Function.
const int numberOfPoints
const double epsilon
const double xmin
const double a
double f1(const double x)
Function.
double F1(const double x)
Integral.
const double xmax
const int numberOfPoints
const JPolynome F1
Integral.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
const double xmax
const double xmin
const double epsilon
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Gauss function object.
Definition JMathlib.hh:1589
Recursive template class for polynomial function.
Definition JPolynome.hh:165
double getIntegral(const double x) const
Integral value.
Definition JPolynome.hh:276