Simpson's 1/3 Code in Matlab

Simpson's 1/3 rules Code implementation in MatLab

Simpson's 1/3 rule

From Book of Steven Chopra -

Aside from applying the trapezoidal rule with finer segmentation, another way to obtain a more accurate estimate of an integral is to use higher-order polynomials to connect the points. 

For example, 

If there is an extra point midway between f(a) and f(b), the three points can be connected with a parabola.

 If there are two points equally spaced between f(a) and f(b), the four points can be connected with a third-order polynomial. The formulas that result from taking the integrals under these polynomials are called Simpson’s rules.

Simpson’s 1/3 Rule to find integral is :

Simpson's 1/3 rules Code implementation in MatLab

Question is:

f (x) = 0.2 + 25x - 200x2 + 675x3 - 900x4 + 400x5

from a 5 0 to b 5 0.8. Recall that the exact integral is 1.640533. Find the Integral and exact error

Solution in normal numerical mathematics of Simpson's 1/3 rules:

f(0) =  0.2 

f(0.4) = 2.456 

f(0.8) = 0.232

Therefore,above figure equation can be used to compute

I = (0.8 - 0) * (0.2 + 4*(2.456) + 0.232)/6

   = 1.367467

which represents an exact error of

E = 1.640533 - 1.367467

   =  0.2730667 

Et =  16.6%

Solution in matlab coding of Simpson's 1/3 rules:

% Author : Maniruzzaman Akash
% Code   : Simpsons 1/3 Rules Code in Matlab Implementation

% Function 
f = @ (x)  .2+ 25*x - 200*x^2+675*x^3-900*x^4+400*x^5;

disp('Press 1 for set values ');
disp('Press 2 for see results ');
disp('Press 0 for exit the program ');
disp('---------------------------------');

isInputGiven = false;
while(1)
    choose = input('Give your choose option: ');
    if(choose == 1)
        a = input('Enter lower limit : ');

        b = input('Enter upper limit : ');
        
        exact_integral = input('Enter exact integral : ');

        x0 = a; 
        x1 = (b-a)/2; 
        x2 = b; 


        fx0 = feval(f,x0);
        fx1 = feval(f,x1);
        fx2 = feval(f,x2);

        integral = ((fx0 + (4*fx1) + fx2)/6 * (b-a));
        isInputGiven = true;
    end;
    
    if(choose == 2)
        if(isInputGiven == true)
            integral_output = sprintf('Integral is = %0.5f',integral);
            disp(integral_output);
            
            exact_error = ((exact_integral - integral)/exact_integral) * 100 ;
            exact_error_output = sprintf('Exact Error = %0.5f',exact_error);
            disp(exact_error_output);
        end;
    end;
    
    if(choose == 0)
        disp('Exiting the code..');
        exit(0);
    end;
end;

Output Of the above code will look like:

Simpson's 1/3 rules Code implementation in MatLab

Simpson 3/8 is just as the problem. Just Equation is the different for that math and change in equation in code also.
Thanks.