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
Solve a problem – Filter by language, license, keyword, owner, or search text to find code & info fast. 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