Source Code Program (MATLAB) Adaptif Kuadratur Simpson

function [I,err,iflg]=adpsim(a,b,tol,fun)

% implementasi adaptif kuadratur Simpson

% Masukkan integrand,fun.

% a,b :Batas integrasi,tol:toleransi eror absolute

% errest: estimasi error

% iflg: Modus pengembalian , memberikan jumlah subinterval di mana

% jumlah maksimum ( levmax = 10 ) dari terbagi dua diperlukan dan

% nilai diterima secara default. Semakin besar iflg, kepercayaan semakin

% berkurang

% harus memiliki nilai yang dihitung, y.

% nofun: jumlah fungsi yang dievaluasi

% inisialisasi

I=0;iflg=0;jflg=0;err=0;levmax=20;

fsave=zeros(levmax,3);xsave=zeros(levmax,3);simp=zeros(levmax);

a=input('a=');

b=input('b=');

tol=input('toleransi awal=');

%fun=@(x)sqrt(x); %integrand

fun=@(x)(sin(x))^2; %integrand

tol2=tol+10*eps;

tol1=tol2*15/(b-a);

x=a:(b-a)/4:b;

for j=1:5

    f(j)=feval(fun,x(j));

end

level=1;

%level=0 berarti seluruh interval tertutup , maka selesai

while level>0

for k=1:3

    fsave(level,k)=f(k+2);

    xsave(level,k)=x(k+2);

end

h=(x(5)-x(1))/4;

simp(level)=(h/3)*(f(3)+4*f(4)+f(5));

if jflg<=0

    s1=2*(h/3)*(f(1)+4*f(3)+f(5));

end

sl=(h/3)*(f(1)+4*f(2)+f(3));

s2=sl+simp(level);

d=abs(s1-s2);

if d<=tol1*4*h

    level=level-1;

    jflg=0;

    I=I+s2;

    err=err+d/15;

    if level<=0

        fprintf('nilai integral=  %.10f \n',I)

        return

    end

    for j=1:3

        jj=2*j-1;

        f(jj)=fsave(level,j);

        x(jj)=xsave(level,j);

    end

   

else

    level=level+1;

    s1=sl;

    if level <= levmax

        jflg=1;

        f(5)=f(3);f(3)=f(2);

        x(5)=x(3);x(3)=x(2);

    else

        iflg=iflg+1;

        level=level-1;

        jflg=0;

        I=I+s2;

        err=err+d/15;

    end

end

    for k=1:2

        kk=2*k;

        x(kk)=.5*(x(kk+1)+x(kk-1));

        f(kk)=feval(fun,x(kk));

    end

end