Lesson of 25 June 2008¶
A simple C-like library for numerical integration
The declaration file
File intrule.h
1/*
2// A simple procedural implementation of some numerical
3// integration rules.
4//
5*/
6
7#include <stdio.h>
8
9typedef double valueType ;
10typedef int indexType ;
11
12typedef valueType (*functionType) ( valueType x ) ;
13
14// declaration of the prototype of the trapezoidal integration rule
15valueType
16trap( functionType fun, // pointer to the function to integrate
17 valueType a, // left boundary
18 valueType b, // right boundary
19 indexType N ) ; // number of interval
20
21// declaration of the prototype of the midpoint integration rule
22valueType
23midpoint( functionType fun, // pointer to the function to integrate
24 valueType a, // left boundary
25 valueType b, // right boundary
26 indexType N ) ; // number of interval
27
28// declaration of the prototype of the Simpson integration rule
29valueType
30simpson( functionType fun, // pointer to the function to integrate
31 valueType a, // left boundary
32 valueType b, // right boundary
33 indexType N ) ; // number of interval
34
The definition file file
File intrule.cc
:open:
1/*
2// A simple procedural implementation of some numerical
3// integration rules.
4//
5*/
6
7#include "intrule.h"
8
9// definition of the prototype of the trapezoidal integration rule
10valueType
11trap( functionType fun,
12 valueType a,
13 valueType b,
14 indexType N ) {
15 valueType h = (b-a)/N ;
16 valueType res = h*((*fun)(a) + fun(b))/2 ;
17 // ^ ^
18 // | + a shortcut or syntactic glue from C
19 // + the "correct way"
20 for ( indexType i = 1 ; i < N ; ++i ) {
21 valueType xi = a + h*i ;
22 res += h*fun(xi) ;
23 }
24 return res ;
25}
26
27// definition of the prototype of the midpoint integration rule
28valueType
29midpoint( functionType fun,
30 valueType a,
31 valueType b,
32 indexType N ) {
33 valueType h = (b-a)/N ;
34 valueType res = 0 ;
35 for ( indexType i = 0 ; i < N ; ++i ) {
36 valueType xi_plus_half = a + h*(i+0.5) ;
37 res += h*fun(xi_plus_half) ;
38 }
39 return res ;
40}
41
42// definition of the prototype of the Simpson integration rule
43valueType
44simpson( functionType fun,
45 valueType a,
46 valueType b,
47 indexType N ) { // N must be even!
48 valueType h = (b-a)/N ;
49 valueType res = h*((*fun)(a) + fun(b))/3 ;
50 for ( indexType i = 1 ; i < N ; i += 2 ) {
51 valueType xi = a + h*i ;
52 res += (4*h/3)*fun(xi) ;
53 }
54 for ( indexType i = 2 ; i < N-1 ; i += 2 ) {
55 valueType xi = a + h*i ;
56 res += (2*h/3)*fun(xi) ;
57 }
58 return res ;
59}
A driver program
File main.cc
:open:
1/*
2// A simple procedural implementation of some numerical
3// integration rules.
4//
5*/
6
7#include "intrule.h"
8#include <math.h>
9
10static
11valueType
12testFun( valueType x ) {
13 return x*exp(-x) ;
14}
15
16int
17main() {
18
19 valueType a = 0 ;
20 valueType b = 10 ;
21
22 for ( indexType N = 10 ; N <= 1000 ; N *= 2 ) {
23
24 valueType t = trap ( testFun, a, b, N ) ; // call trapeziodal integrator
25 valueType m = midpoint ( testFun, a, b, N ) ; // call midpoint integrator
26 valueType s = simpson ( testFun, a, b, N ) ; // call Simpson integrator
27
28 // print the result
29 printf( "Integration with %d intervals\n", N);
30 printf( "Trapezoidal rule %10lf\n", t);
31 printf( "Midpoint rule %10lf\n", m);
32 printf( "Simspon rule %10lf\n", s);
33 }
34 printf( "All done!\n");
35 return 0;
36}
An Object Oriented version of the library
The declaration file
File intrule.h
:open:
1/*
2// A simple object oriented implementation of some numerical
3// integration rules.
4//
5*/
6
7typedef double valueType ;
8typedef int indexType ;
9
10typedef valueType (*functionType) ( valueType x ) ;
11
12class numericalIntegrator {
13protected:
14 functionType fun ;
15public:
16 // default initialilizer
17 numericalIntegrator()
18 : fun(0)
19 {}
20
21 numericalIntegrator(functionType ffun)
22 : fun(ffun)
23 {}
24
25 // set up a funtion to integrate
26 void
27 setup( functionType ffun ) {
28 fun = ffun ;
29 }
30
31 //
32 // To integrate a function if INT is of type numericalIntegrator
33 // INT . setup(fun);
34 // res = INT(a,b,N) ;
35 //
36 valueType
37 operator () ( valueType a, valueType b, indexType N ) ;
38
39} ;
40
41class trap : public numericalIntegrator {
42public:
43 // default initialilizer
44 trap()
45 : numericalIntegrator()
46 {}
47
48 trap(functionType ffun )
49 : numericalIntegrator(ffun)
50 {}
51
52 //
53 // To integrate a function if INT is of type numericalIntegrator
54 // INT . setup(fun);
55 // res = INT(a,b,N) ;
56 //
57 valueType
58 operator () ( valueType a, valueType b, indexType N ) ;
59
60} ;
61
62class midpoint : public numericalIntegrator {
63public:
64 // default initialilizer
65 midpoint()
66 : numericalIntegrator()
67 {}
68
69 midpoint(functionType ffun )
70 : numericalIntegrator(ffun)
71 {}
72
73 //
74 // To integrate a function if INT is of type numericalIntegrator
75 // INT . setup(fun);
76 // res = INT(a,b,N) ;
77 //
78 valueType
79 operator () ( valueType a, valueType b, indexType N ) ;
80
81} ;
82
83class simpson : public numericalIntegrator {
84public:
85 // default initialilizer
86 simpson()
87 : numericalIntegrator()
88 {}
89
90 simpson(functionType ffun )
91 : numericalIntegrator(ffun)
92 {}
93
94 //
95 // To integrate a function if INT is of type numericalIntegrator
96 // INT . setup(fun);
97 // res = INT(a,b,N) ;
98 //
99 valueType
100 operator () ( valueType a, valueType b, indexType N ) ;
101
102} ;
The definition file file
File intrule.cc
:open:
1/*
2// A simple object oriented implementation of some numerical
3// integration rules.
4//
5*/
6
7#include "intrule.h"
8
9valueType
10trap::operator () ( valueType a, valueType b, indexType N ) {
11 valueType h = (b-a)/N ;
12 valueType res = h*(fun(a) + fun(b))/2 ;
13 for ( indexType i = 1 ; i < N ; ++i ) {
14 valueType xi = a + h*i ;
15 res += h*fun(xi) ;
16 }
17 return res ;
18}
19
20valueType
21midpoint::operator () ( valueType a, valueType b, indexType N ) {
22 valueType h = (b-a)/N ;
23 valueType res = 0 ;
24 for ( indexType i = 0 ; i < N ; ++i ) {
25 valueType xi_plus_half = a + h*(i+0.5) ;
26 res += h*fun(xi_plus_half) ;
27 }
28 return res ;
29}
30
31valueType
32simpson::operator () ( valueType a, valueType b, indexType N ) {
33 valueType h = (b-a)/N ;
34 valueType res = h*((*fun)(a) + fun(b))/3 ;
35 for ( indexType i = 1 ; i < N ; i += 2 ) {
36 valueType xi = a + h*i ;
37 res += (4*h/3)*fun(xi) ;
38 }
39 for ( indexType i = 2 ; i < N-1 ; i += 2 ) {
40 valueType xi = a + h*i ;
41 res += (2*h/3)*fun(xi) ;
42 }
43 return res ;
44}
A driver program
File main.cc
:open:
1/*
2// A simple object oriented implementation of some numerical
3// integration rules.
4//
5*/
6
7#include "intrule.h"
8#include <cmath>
9#include <iostream>
10
11using namespace std ;
12
13static
14valueType
15testFun( valueType x ) {
16 return x*exp(-x) ;
17}
18
19int
20main() {
21
22 trap tInstance( testFun ) ;
23 midpoint mInstance ;
24 simpson sInstance ;
25
26 mInstance . setup( testFun ) ;
27 sInstance . setup( testFun ) ;
28
29 valueType a = 0 ;
30 valueType b = 10 ;
31
32 for ( indexType N = 10 ; N <= 1000 ; N *= 2 ) {
33
34 valueType t = tInstance( a, b, N ) ; // call trapeziodal integrator
35 valueType m = mInstance( a, b, N ) ; // call midpoint integrator
36 valueType s = sInstance( a, b, N ) ; // call Simpson integrator
37
38 // print the result
39 cout << "Integration with " << N << " intervals\n"
40 << "Trapezoidal rule " << t << "\n"
41 << "Midpoint rule " << m << "\n"
42 << "Simspon rule " << s << "\n" ;
43 }
44 cout << "All done!\n" ;
45 return 0;
46}
An Object Oriented version of the library with virtual interface
The declaration file
File intrule.h
:open:
1/*
2// A simple object oriented implementation of some numerical
3// integration rules.
4//
5*/
6
7#include <iostream>
8
9typedef double valueType ;
10typedef int indexType ;
11
12typedef valueType (*functionType) ( valueType x ) ;
13
14class numericalIntegrator {
15protected:
16 functionType fun ;
17public:
18 // default initialilizer
19 numericalIntegrator()
20 : fun(0)
21 {}
22
23 numericalIntegrator(functionType ffun)
24 : fun(ffun)
25 {}
26
27 virtual
28 ~numericalIntegrator()
29 {}
30
31 // set up a funtion to integrate
32 void
33 setup( functionType ffun ) {
34 fun = ffun ;
35 }
36
37 //
38 // To integrate a function if INT is of type numericalIntegrator
39 // INT . setup(fun);
40 // res = INT(a,b,N) ;
41 // declaring the operator PURE VIRTUAL
42 virtual
43 valueType
44 operator () ( valueType a, valueType b, indexType N ) const = 0 ;
45
46} ;
47
48class trap : public numericalIntegrator {
49public:
50 // default initialilizer
51 trap()
52 : numericalIntegrator()
53 {}
54
55 trap(functionType ffun )
56 : numericalIntegrator(ffun)
57 {}
58
59 virtual
60 ~trap()
61 {}
62
63 //
64 // To integrate a function if INT is of type numericalIntegrator
65 // INT . setup(fun);
66 // res = INT(a,b,N) ;
67 //
68 virtual
69 valueType
70 operator () ( valueType a, valueType b, indexType N ) const ;
71
72} ;
73
74class midpoint : public numericalIntegrator {
75public:
76 // default initialilizer
77 midpoint()
78 : numericalIntegrator()
79 {}
80
81 midpoint(functionType ffun )
82 : numericalIntegrator(ffun)
83 {}
84
85 virtual
86 ~midpoint()
87 {}
88
89 //
90 // To integrate a function if INT is of type numericalIntegrator
91 // INT . setup(fun);
92 // res = INT(a,b,N) ;
93 //
94 virtual
95 valueType
96 operator () ( valueType a, valueType b, indexType N ) const ;
97
98} ;
99
100class simpson : public numericalIntegrator {
101public:
102 // default initialilizer
103 simpson()
104 : numericalIntegrator()
105 {}
106
107 simpson(functionType ffun )
108 : numericalIntegrator(ffun)
109 {}
110
111 virtual
112 ~simpson()
113 {}
114
115 //
116 // To integrate a function if INT is of type numericalIntegrator
117 // INT . setup(fun);
118 // res = INT(a,b,N) ;
119 //
120 virtual
121 valueType
122 operator () ( valueType a, valueType b, indexType N ) const ;
123
124} ;
The definition file
File intrule.cc
:open:
1/*
2// A simple object oriented implementation of some numerical
3// integration rules.
4//
5*/
6
7#include "intrule.h"
8
9valueType
10trap::operator () ( valueType a, valueType b, indexType N ) const {
11 valueType h = (b-a)/N ;
12 valueType res = h*(fun(a) + fun(b))/2 ;
13 for ( indexType i = 1 ; i < N ; ++i ) {
14 valueType xi = a + h*i ;
15 res += h*fun(xi) ;
16 }
17 return res ;
18}
19
20valueType
21midpoint::operator () ( valueType a, valueType b, indexType N ) const {
22 valueType h = (b-a)/N ;
23 valueType res = 0 ;
24 for ( indexType i = 0 ; i < N ; ++i ) {
25 valueType xi_plus_half = a + h*(i+0.5) ;
26 res += h*fun(xi_plus_half) ;
27 }
28 return res ;
29}
30
31valueType
32simpson::operator () ( valueType a, valueType b, indexType N ) const {
33 valueType h = (b-a)/N ;
34 valueType res = h*((*fun)(a) + fun(b))/3 ;
35 for ( indexType i = 1 ; i < N ; i += 2 ) {
36 valueType xi = a + h*i ;
37 res += (4*h/3)*fun(xi) ;
38 }
39 for ( indexType i = 2 ; i < N-1 ; i += 2 ) {
40 valueType xi = a + h*i ;
41 res += (2*h/3)*fun(xi) ;
42 }
43 return res ;
44}
A driver program
File main.cc
:open:
1/*
2// A simple object oriented implementation of some numerical
3// integration rules.
4//
5*/
6
7#include "intrule.h"
8#include <cmath>
9#include <iostream>
10
11using namespace std ;
12
13static
14valueType
15testFun( valueType x ) {
16 return x*exp(-x) ;
17}
18
19int
20main() {
21
22 numericalIntegrator * tInstance = new trap( testFun ) ;
23 numericalIntegrator * mInstance = new midpoint( testFun ) ;
24 numericalIntegrator * sInstance = new simpson( testFun ) ;
25
26 valueType a = 0 ;
27 valueType b = 10 ;
28
29 for ( indexType N = 10 ; N <= 1000 ; N *= 2 ) {
30
31 valueType t = (*tInstance)( a, b, N ) ; // call trapeziodal integrator
32 valueType m = (*mInstance)( a, b, N ) ; // call midpoint integrator
33 valueType s = (*sInstance)( a, b, N ) ; // call Simpson integrator
34
35 // print the result
36 cout << "Integration with " << N << " intervals\n"
37 << "Trapezoidal rule " << t << "\n"
38 << "Midpoint rule " << m << "\n"
39 << "Simspon rule " << s << "\n" ;
40 }
41
42 delete tInstance ;
43 delete mInstance ;
44 delete sInstance ;
45
46 cout << "All done!\n" ;
47
48 return 0;
49}
The lesson in a zip file
The zip file lesson4.zip