MATLAB: Simpson’s 1/3 Rule Implementation Issue (Code Provided)

code providedMATLABnumerical methodssimpsons rulesimsons1/3

I have written the following code to Implement Simpson's 1/3 Rule, but it does not seem to give the correct answer. I am unsure what I have done wrong in my code.
MATLAB CODE:
function Integration(a,b,x,n,true)
a=0; %Lower Integral Bound
b=1; %Upper Integral Bound
x=[a,.25,.5,.75,b];
n=4; %Even Number of Segments
f=zeros(n+1,1);
true=.602298; %True Value of Equation
h=(b-a)/n; %width of each segment
f(1)=(x(1)^0.1)*(1.2-x(1))*(1-exp(20*(x(1)-1))); %Value For Initial Bound
for i=2:n
f(i)=2*(x(i)^0.1)*(1.2-x(i))*(1-exp(20*(x(i)-1))); %Values for In-Between Bounds
end
f(n+1)=(x(n+1)^0.1)*(1.2-x(n+1))*(1-exp(20*(x(n+1)-1))); %Value for Final Bound
q=sum(f,'all');
I=(h/2)*(q); %Width Times Average Height
disp('USING TRAPEZOIDAL')
fprintf('I=')
disp(I)
et=abs((true-I)/true)*100;
fprintf('True Error=')
disp(et)
disp("USING SIMPSON'S 1/3")
f2=zeros(n,1);
f2(1)=(x(1)^0.1)*(1.2-x(1))*(1-exp(20*(x(1)-1)));
for j=2:2:n-2
f2(j)=2*(x(j)^0.1)*(1.2-x(j))*(1-exp(20*(x(j)-1)));
end
for j=3:2:n-1
f2(j)=4*(x(j)^0.1)*(1.2-x(j))*(1-exp(20*(x(j)-1)));
end
f2(n)=(x(n)^0.1)*(1.2-x(n))*(1-exp(20*(x(n)-1)));
v=sum(f2,'all');
I2=(h/3)*v;
fprintf('I=')
disp(I2)
%ea=abs((true-I2)/true)*100;
%printf('True Error=')
%disp(ea)
%Nobody Cares, Work Harder
%Keep Hammering
end

Best Answer

  • The following function might help.
    function o = simp(a,b,f,n)
    x=linspace(b,a,n+1);
    o=0;
    for k=1:n
    xx=[x(k);(x(k)+x(k+1))/2;x(k+1)];
    y=f(xx);
    c=([xx.^2,xx,ones(3,1)]\y)';
    p=polyint(c);
    o=o+polyval(p,x(k))-polyval(p,x(k+1));
    end
    Evaluate the function using your equation desiring integration.
    f=@(x)x.^.1.*(1.2-x).*(1-exp(20*(x-1)));
    output=simp(0,1,f,4);