License New BSD license
Lines 55
Keywords
booles rule (1) newton cotes (1) numerical computing (8) numerical integration (1) numerics (5) simpsons rule (1)
Permissions
Owner: Stou S.
Viewable by Everyone
Editable by All Siafoo Users
Hide
Easily highlight source code for your blog with our Syntax Highlighter. Join Siafoo Now or Learn More

Numerical Integration with Newton-Cotes Atom Feed 0

In Brief To use it to integrate a function like , create a functor like:... more
# 's
 1/*
2 * NewtonCotes.hpp
3 *
4 * Created on: Oct 15, 2009
5 * Author: stou@icapsid.net
6 * License: New BSD
7 */
8
9#ifndef HPP_NEWTON_COTES
10#define HPP_NEWTON_COTES
11
12/*
13 * Note: My tests indicate that gcc -O3 should inline this function into the specific algorithms below
14 * resulting in identical code with the algorithms having the formula built in.
15 */
16template<typename T, typename F> T NewtonCotes(F f, T a, T b, unsigned int segments, unsigned int degree, unsigned int divisor, const T coeff[]){
17
18 T seg_step = (b - a)/segments;
19 T f_val = T(0);
20
21 for(unsigned int s = 0; s < segments; ++s){
22
23 T seg_a = a + seg_step*s;
24 T seg_b = seg_a + seg_step;
25 T f_subval = T(0);
26
27 for(unsigned int i = 0; i < degree + 1; ++i){
28 f_subval += coeff[i] * f(seg_a + i*(seg_b - seg_a)/degree);
29 }
30
31 f_val += ((seg_b - seg_a)/T(divisor)) * f_subval;
32 }
33
34 return f_val;
35}
36
37/*
38 * Simpsons Rule
39 */
40template<typename T, typename F> T Simpsons(F f, T a, T b, unsigned int segments){
41
42 const unsigned int degree = 2;
43 const unsigned int divisor = 6;
44 const T coeff[] = {1., 4., 1.};
45
46 return NewtonCotes<T, F>(f, a, b, segments, degree, divisor, coeff);
47}
48
49/*
50 * Simpsons 3/8 Rule
51 */
52template<typename T, typename F> T Simpsons3_8(F f, T a, T b, unsigned int segments){
53
54 const unsigned int degree = 3;
55 const unsigned int divisor = 8;
56 const T coeff[] = {1., 3., 3., 1.};
57
58 return NewtonCotes<T, F>(f, a, b, segments, degree, divisor, coeff);
59}
60
61/*
62 * Booles Rule
63 */
64template<typename T, typename F> T Booles(F f, T a, T b, unsigned int segments){
65
66 const unsigned int degree = 4;
67 const unsigned int divisor = 90;
68 const T coeff[] = {7., 32., 12., 32., 7.};
69
70 return NewtonCotes<T, F>(f, a, b, segments, degree, divisor, coeff);
71}
72
73#endif // HPP_NEWTON_COTES

To use it to integrate a function like , create a functor like:

# 's
1struct PythagorianTest{
2 float operator()(float t) {
3 return sin(t)*sin(t) + cos(t)*cos(t);
4 }
5};

and then call one of the integration methods like:

# 's
1Simpsons<float, PythagorianTest>(PythagorianTest(), 0, Pi, 10)

A full working example:

# 's
 1#include "NewtonCotes.hpp"
2#include <iostream>
3#include <cmath>
4
5#define Pi 3.14159265
6
7using namespace std;
8
9struct ConstTest{
10 float operator()(float t) { return 1;}
11};
12
13struct PythagorianTest{
14 float operator()(float t) {
15 return sin(t)*sin(t) + cos(t)*cos(t);
16 }
17};
18
19int main(){
20 cout<<"Simpsons: "<<Simpsons<float, PythagorianTest>(PythagorianTest(), 0, Pi, 10)<<endl;
21 cout<<"Simpsons3_8: "<<Simpsons3_8<float, PythagorianTest>(PythagorianTest(), 0, Pi, 10)<<endl;
22 cout<<"Booles: "<<Booles<float, PythagorianTest>(PythagorianTest(), 0, Pi, 10)<<endl;
23 return 0;
24}

This code isn't extensively tested and I haven't used it in my project yet... so drop me a line if you find a problem or an improvement (or just edit this)

Newton-Cotes on MathWorld or Wikipedia