Numerical Integration

The Learning Point‎ > ‎Mathematics‎ > ‎Numerical Analysis‎ > ‎

Introduction

Need prerequisite knowledge of Newtons Forward Interpretation formula.

Forward interpretation formula is a statistical method of finding the value function f(x) which lies near the beginning of the table( we have data of y=f(x) for x₀, x₁, x₂, .. x_iwhich is equally spaced).

Suppose I is an integration of f(x) in the limit (a,b). Let us divide the interval (a,b) into n equal parts of width h so that

[a= x₀, x₁] ,[ x₁, x₂ ] , [x₂, x₃], ... [x_n-1, x_n = b]

by using Newtons Forward Interpretation formula in integration of f(x) we get the equation called Newtons General Quadrature.

Weddle's Rule

Description

  Integration of f(x) in range (a,b), first we divide the interval into n equally spaced intervals just like Airthematic Progression.
                        x₀ = a
                        x₁ = x₀ + h
                        x₂ = x₀ + 2h or x₁ + h
                        ...
                        ...
                        x_n = b or x₀ + nh or x_n-1 + h
Values of f(x) can be shown as.

x	f(x)
x₀	y₀
x₁	y₁
x₂	y₂
x₃	y₃
...	...
x_n	y_n

Now using the formula of Newtons General Quadrature Method, if we integrate the function as below interval size.

[x₀ , x₀+ 6h]

    [x₀ + 6h, x₀ + 12h]
    ....
    ....

[x₀ + (n-6)h, x₀ + nh] and by adding all terms we get the integration of f(x) in range (a,b).

Finally, after adding formula of Weddle's Rule is for n=6 looks like:

I = 3h/10[ y₀ + 5*y₁ + y₂ + 6*y₃ + y₄ + 5*y₅ + y₆]

Programming implementation of Weddle's Rule in C++

#include<iostream>
#include<math.h>
using namespace std;

double f(double x){
    return 1/(1 + x*x);
}

int main(){
    // Define your function
    double a,b,h;
    cout<<"Give Limit of integration (a,b):";
    cin>>a>>b;
    int n;
    //do{
        cout<<"Define value of h:";
        cin>>h;
        n = ceil((b-a)/h);
    //}while(n%6 != 0);

    double integration=0;
    for(int j=0; j<=n;j++){
        if(j%2 == 0)
            integration += f(a + j*h);
        else if(j%3 == 0)
            integration += 6*f(a + j*h);
        else
            integration += 5*f(a + j*h);
    }

    integration = (3*h*integration)/10;
    cout<<"The value of Integration of f(x) limit (a,b) is:"<<integration<<endl;
    return 0;
}

Code execution looks like:

Give Limit of integration (a,b):0 1
Define value of h:0.1666667
The value of Integration of f(x) limit (a,b) is:0.7854

Note: Method is applied when the number of subintervals is multiple of 6 or 6, the above code is only for n=6, try out code as per your requirements, and do experiments.

Simpson's 3/8th Rule

Description

x	f(x)
x₀	y₀
x₁	y₁
x₂	y₂
x₃	y₃
...	...
x_n	y_n

Now using the formula of Newtons General Quadrature Method, if we integrate the function as below interval sizes.

[x₀ , x₀ + 3h]

    [x₀ + 3h, x₀ + 6h]
    ....
    ....

[x₀ + (n-3)h, x₀ + nh] and by adding them all we get the intergration of f(x) in range (a,b).

Finally, after adding formula of Simpson's 3/8^th Rule looks like:

I = 3h/8[ y₀ + 2 (y₃ + y₆ + y₉ + ... ) + 3(y₂ + y₄ + y₅ + y₇ + ...) + y_n]

Programming implementation of Simpson's3/8th Rule in C++

#include<iostream>
#include<math.h>
using namespace std;

double f(double x){
    return 1/(1 + x*x);
}

int main(){
    // Define your function
    double a,b,h;
    cout<<"Give Limit of integration (a,b):";
    cin>>a>>b;
    int n;
    //do{
        cout<<"Define value of h:";
        cin>>h;
        n = ceil((b-a)/h);
    //}while(n%3 != 0);

    double integration=0;
    for(int j=0; j<=n;j++){
        if(j==0 || j==n)
            integration += f(a + j*h);
        else if(j%3 == 0)
            integration += 2*f(a + j*h);
        else
            integration += 3*f(a + j*h);
    }

    integration = (3*h*integration)/8;
    cout<<"The value of Integration of f(x) limit (a,b) is:"<<integration<<endl;

    return 0;
}

Code execution looks like:

Give Limit of integration (a,b):0 1
Define value of h:0.1666667
The value of Integration of f(x) limit (a,b) is:0.785396

Note: Simpson's 3/8^th rule is only applied when the number of subintervals is multiple of 3.
Recommend to try out the code on other functions with different interval size h, test the accuracy of the method.

Simpson's 1/3rd Rule

Description

x	f(x)
x₀	y₀
x₁	y₁
x₂	y₂
x₃	y₃
...	...
x_n	y_n

Now using the formula of Newtons General Quadrature Method, if we integrate the function as below interval sizes.

[x₀ , x₀ + 2h]

    [x₀ + 2h, x₀ + 4h]
    ....
    ....

[x₀ + (n-2)h, x₀ + nh] and by adding them all we get the intergration of f(x) in range (a,b).

Finally, after adding formula of Simpson's 1/3^rd Rule is looks like:

I = h/3[ y₀+ 4(y₁ + y₃+ y₅ + .. + y_n-1) + 2( y₂ + y₄ + ... + y_n-2) + y_n]

Programming implementation of Simpson's 1/3 rd Rule:

#include<iostream>
#include<math.h>
using namespace std;

double f(double x){
    return 1/(1 + x*x);
}

int main(){

    // Define your function
    double a,b,h;
    cout<<"Give Limit of integration (a,b):";
    cin>>a>>b;
    cout<<"Define value of h:";
    cin>>h;

    int n = ceil((b-a)/h);

    double integration=0;
    for(int j=0; j<=n;j++){
        if(j==0 || j==n)
            integration += f(a + j*h);
        else if(j%2 == 1)
            integration += 4*f(a + j*h);
        else
            integration += 2*f(a + j*h);
    }

    integration = (h*integration)/3;
    cout<<"The value of Integration of f(x) limit (a,b) is:"<<integration<<endl;

    return 0;
}

Execution of above code:

Give Limit of integration (a,b):0 1
Define value of h:0.25
The value of Integration of f(x) limit (a,b) is:0.785392

Note: Simpson's 1/3^rd rule applied only when the number of subintervals is even. Recommend to try out the code on other function with different interval size h, test the accuracy of the method.

Trapezoidal Rule

Description

Integration of f(x) in range (a,b), first we divide the interval into n equally spaced intervals just like Airthematic Progression with difference h.

                        x₀ = a
                        x₁ = x₀ + h
                        x₂ = x₀ + 2h or x₁ + h
                        ...
                        ...
                        x_n = b or x₀ + nh or x_n-1 + h
Values of f(x) can be shown as.

x	f(x)
x₀	y₀
x₁	y₁
x₂	y₂
x₃	y₃
...	...
x_n	y_n

Now using the formula of Newtons General Quadrature Method, if we integrate the function as below interval sizes.

[x₀ , x₀ + h]

    [x₀ + h, x₀+ 2h]
    ....
    ....

[x₀ + (n-1)h, x₀ + nh] and by adding them all, we get the integration of f(x) in range (a,b).

Finally after adding, formula of Trapezoidal Rule is looks like:

I = h/2[ y₀ + 2 (y₁ + y₂ + y₃ + ... + y_n-1) + y_n ]

Programming implementation of Trapezoidal Rule in C++

#include<iostream>
using namespace std;

double f(double x){
    return 1/(1 + x*x);
}

int main(){

    // Define your function
    double a,b,h;
    cout<<"Give Limit of integration (a,b):";
    cin>>a>>b;
    cout<<"Define value of h:";
    cin>>h;

    double integration=0;
    for(double i=a; i<=b; i+= h){
        if(i==a || i==b)
            integration += f(i);
        else
            integration += 2*f(i);
    }

    integration = (h*integration)/2;
    cout<<"The value of Integration of f(x) limit (a,b) is:"<<integration<<endl;

    return 0;
}

Execution of above code:

Give Limit of integration (a,b):0 1
Define value of h:0.2
The value of Integration of f(x) limit (a,b) is:0.783732

Note: Trapezoidal Rule is very simple and good approximate method to find integration with good accuracy( increasing n or decreasing h some time have a bad effect on result accuracy).

Recommend to try out the code on other function with different interval size h, test the accuracy of the method.

The Learning Point

Navigation

Navigation