source: trunk/source/global/HEPNumerics/test/G4IntegratorTest.cc@ 1340

Last change on this file since 1340 was 1315, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 7.9 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4IntegratorTest.cc,v 1.10 2006/08/21 12:24:32 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
29//
30// Test program for G4Integrator class. The function std::exp(-x)*std::cos(x) is
31// integrated between zero and two pi. The exact result is 0.499066278634
32//
33// History:
34//
35// 05.04.97 V.Grichine, first implementation
36// 04.09.99 V.Grichine, integration of member function from a class and main,
37// as well as integration of global scope functions
38
39#include "G4ios.hh"
40#include "globals.hh"
41#include "G4SimpleIntegration.hh"
42#include "G4Integrator.hh"
43
44G4double GlobalFunction( G4double x ){ return std::exp(-x)*std::cos(x) ; }
45
46G4double GlobalCos( G4double x ){ return std::cos(x) ; }
47
48G4double GlobalHermite(G4double x){ return x*x*std::cos(x) ; }
49
50
51G4double fY;
52G4int fN = 100;
53
54G4double X1 = -3.25/5.;
55G4double X2 = 3.25/5.;
56
57G4double Y1 = -7.5/5.;
58G4double Y2 = 7.5/5.;
59
60G4double Harp(G4double x)
61{
62
63 G4double tmp = std::sqrt(1 + x*x + fY*fY);
64 tmp *= 1 + x*x + fY*fY ;
65 return 1/tmp;
66}
67
68
69G4double HarpX(G4double y)
70{
71 fY = y;
72 G4SimpleIntegration myIntegrand(Harp);
73 return myIntegrand.Simpson(X1,X2,fN);
74 // return myIntegrand.Gauss(X1,X2,fN);
75}
76
77
78
79
80
81
82G4double HarpY()
83{
84 // G4int i;
85 G4SimpleIntegration myIntegrand(HarpX);
86 // G4double sum =0.;
87 // for(i = 0; i < fN; i++)
88 G4double result = myIntegrand.Simpson(Y1,Y2,fN);
89 return result/twopi;
90
91}
92
93class B
94{
95 typedef G4double (B::* PBmem)(G4double);
96
97public:
98
99 B(){;}
100
101 ~B(){;}
102
103 G4double TestFunction(G4double x){ return std::exp(-x)*std::cos(x) ; }
104
105 G4double CosFunction(G4double x) { return std::cos(x) ; }
106
107 G4double TestHermite(G4double x){ return x*x*std::cos(x) ; }
108
109 G4double HarpX(G4double);
110 G4double HarpY(G4double);
111
112
113 void Integrand() ;
114
115};
116
117
118void B::Integrand()
119{
120 G4int i, n ;
121 G4double pTolerance;
122 G4double simpson1=0., simpson2=0. ;
123 G4double legendre1=0., legendre2=0. ;
124 G4double chebyshev1=0., chebyshev2=0. ;
125 G4double adaptg1=0., adaptg2=0.;
126 G4double a = 0.0 ;
127 G4double b = twopi ;
128
129 G4SimpleIntegration myIntegrand(GlobalFunction) ;
130
131 G4Integrator<B,PBmem> integral;
132
133 B bbb ;
134
135 G4cout<<"Iter"
136 <<"\t"<<"Simpson "<<"\t"<<"Simpson "
137 <<"\t"<<"Legendre"<<"\t"<<"Legendre"
138 <<"\t"<<"Chebyshev"<<"\t"<<"Chebyshev"<<G4endl ;
139
140 for(i=1;i<=20;i++)
141 {
142 n = 2*i ;
143
144 simpson1 = integral.Simpson(this,&B::TestFunction,a,b,n) ;
145 legendre1 = integral.Legendre(this,&B::TestFunction,a,b,n) ;
146 simpson2 = integral.Simpson(bbb,&B::TestFunction,a,b,n) ;
147 legendre2 = integral.Legendre(bbb,&B::TestFunction,a,b,n) ;
148 chebyshev1 = integral.Chebyshev(this,&B::TestFunction,a,b,n) ;
149 chebyshev2 = integral.Chebyshev(bbb,&B::TestFunction,a,b,n) ;
150
151 G4cout<<n
152 <<"\t"<<simpson1<<"\t"<<simpson2
153 <<"\t"<<legendre1<<"\t"<<legendre2
154 <<"\t"<<chebyshev1<<"\t"<<chebyshev2<<G4endl ;
155 }
156 G4cout<<G4endl ;
157
158 for(i=0;i<8;i++)
159 {
160 pTolerance = std::pow(10.0,-i) ;
161 adaptg1 = integral.AdaptiveGauss(bbb,&B::TestFunction,a,b,pTolerance) ;
162 adaptg2 = integral.AdaptiveGauss(this,&B::TestFunction,a,b,pTolerance) ;
163 G4cout<<pTolerance<<"\t"<<adaptg1<<"\t"<<adaptg2<<G4endl;
164 }
165 for(i=1;i<20;i++)
166 {
167 n = 1*i ;
168 G4double laguerre1=0., laguerre2=0.;
169 laguerre1 = integral.Laguerre(bbb,&B::CosFunction,0.0,n) ;
170 laguerre2 = integral.Laguerre(this,&B::CosFunction,0.0,n) ;
171 G4cout<<"n = "<<n<<"\t"<<"exact = 0.5 "
172 <<" and n-point GaussLaguerre = "
173 <<laguerre1<<"\t"<<laguerre2<<G4endl ;
174 }
175 for(i=1;i<20;i++)
176 {
177 n = 1*i ;
178 G4double hermite1=0., hermite2=0;
179 G4double exactH = 2*0.125*std::sqrt(pi)*std::exp(-0.25) ;
180 hermite1 = integral.Hermite(bbb,&B::TestHermite,n) ;
181 hermite2 = integral.Hermite(this,&B::TestHermite,n) ;
182 G4cout<<"n = "<<n<<"\t"<<"exact = "<<exactH
183 <<" and n-point GaussHermite = "
184 <<hermite1<<"\t"<<hermite2<<G4endl ;
185 }
186 G4double exactJ = pi*0.4400505857 ;
187
188 for(i=1;i<20;i++)
189 {
190 n = 1*i ;
191 G4double jacobi1=0., jacobi2=0.;
192 jacobi1 = integral.Jacobi(bbb,&B::CosFunction,0.5,0.5,n) ;
193 jacobi2 = integral.Jacobi(this,&B::CosFunction,0.5,0.5,n) ;
194 G4cout<<"n = "<<n<<"\t"<<"exact = "<<exactJ
195 <<" and n-point Gauss-Jacobi = "
196 <<jacobi1<<"\t"<<jacobi2<<G4endl ;
197 }
198}
199
200
201int main()
202{
203 B myIntegration ;
204
205 myIntegration.Integrand() ;
206
207 G4Integrator<B,function> iii ;
208 G4int i, n ;
209 G4double a = 0.0 ;
210 G4double b = twopi ;
211 G4double simpson3,legendre,legendre10,legendre96,chebyshev ;
212
213 G4cout<<"Global function integration"<<G4endl ;
214 G4cout<<"n = "<<"\t"<<"Simpson"<<"\t"
215 <<"\t"<<"Legendre""\t"<<"Chebyshev"<<G4endl ;
216 for(i=1;i<=30;i++)
217 {
218 n = 2*i ;
219 simpson3 = iii.Simpson(GlobalFunction,a,b,n) ;
220 legendre = iii.Legendre(GlobalFunction,a,b,n) ;
221 chebyshev = iii.Chebyshev(GlobalFunction,a,b,n) ;
222 G4cout<<n<<"\t"<<simpson3<<"\t"<<legendre<<"\t"<<chebyshev<<G4endl ;
223 }
224 legendre10 = iii.Legendre10(GlobalFunction,a,b) ;
225 legendre96 = iii.Legendre96(GlobalFunction,a,b) ;
226 G4cout<<"Legendre 10 points = "<<"\t"<<legendre10<<G4endl ;
227 G4cout<<"Legendre 96 points = "<<"\t"<<legendre96<<G4endl ;
228
229 for(i=0;i<8;i++)
230 {
231 G4double pTolerance = std::pow(10.0,-i) ;
232 G4double adaptg = iii.AdaptiveGauss(GlobalFunction,a,b,pTolerance) ;
233 G4cout<<pTolerance<<"\t"<<adaptg<<G4endl;
234 }
235 for(i=1;i<20;i++)
236 {
237 n = 1*i ;
238 G4double laguerre = iii.Laguerre(GlobalCos,0.0,n) ;
239 G4cout<<"n = "<<n<<"\t"<<"exact = 0.5 "
240 <<" and n-point Laguerre = "
241 <<laguerre<<G4endl ;
242 }
243 for(i=1;i<20;i++)
244 {
245 n = 1*i ;
246 G4double exactH = 2*0.125*std::sqrt(pi)*std::exp(-0.25) ;
247 G4double hermite = iii.Hermite(GlobalHermite,n) ;
248 G4cout<<"n = "<<n<<"\t"<<"exact = "<<exactH
249 <<" and n-point Hermite = "
250 <<hermite<<G4endl ;
251 }
252 G4double exactJ = pi*0.4400505857 ;
253
254 for(i=1;i<20;i++)
255 {
256 n = 1*i ;
257 G4double jacobi = iii.Jacobi(GlobalCos,0.5,0.5,n) ;
258 G4cout<<"n = "<<n<<"\t"<<"exact = "<<exactJ
259 <<" and n-point Jacobi = "
260 <<jacobi<<G4endl ;
261 }
262
263
264
265 G4cout<<"Ivanchenko integral = "<<HarpY()<<G4endl;
266
267
268 return 0;
269}
Note: See TracBrowser for help on using the repository browser.