Senin, 02 November 2015

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

Mengubah File XPS menjadi PDF

Baru-baru ini saya dibingungkan karena tidak bisa membuka file .xps di rental komputer. Akhirnya setelah searching di google ternyata gampang kok mengubah file XPS jadi PDF, yang pastinya lebih familiar. Yaitu dengan convert online di https://xps2pdf.co.uk/. Selamat mencoba.

Source code MATLAB untuk Menyelesaikan Sistem dengan 3 Persamaan (Metode Adam-Bashforth-Moulton)

%Metode Adam-Bashforth-Moulton
clc;
clear all;
h=0.01;%lebar sub interval
t=0:h:2000;%batas interval
N=length(t)-1;
S=zeros(N+1,1);
I=zeros(N+1,1);
P=zeros(N+1,1);
%nilai awal
S(1)=50;
I(1)=20;
P(1)=30;
for i=1:3
    % Runge Kutta order 4 untuk mencari titik 1 sampai 3
    k1t=f1(t(i),S(i),I(i),P(i));
    l1t=f2(t(i),S(i),I(i),P(i));
    m1t=f3(t(i),S(i),I(i),P(i));
   
    k2t=f1(t(i)+h/2,S(i)+k1t*h/2,I(i)+l1t*h/2,P(i)+m1t*h/2);
    l2t=f2(t(i)+h/2,S(i)+k1t*h/2,I(i)+l1t*h/2,P(i)+m1t*h/2);
    m2t=f3(t(i)+h/2,S(i)+k1t*h/2,I(i)+l1t*h/2,P(i)+m1t*h/2);
   
    k3t=f1(t(i)+h/2,S(i)+k2t*h/2,I(i)+l2t*h/2,P(i)+m2t*h/2);
    l3t=f2(t(i)+h/2,S(i)+k2t*h/2,I(i)+l2t*h/2,P(i)+m2t*h/2);
    m3t=f3(t(i)+h/2,S(i)+k2t*h/2,I(i)+l2t*h/2,P(i)+m2t*h/2);
   
    k4t=f1(t(i)+h,S(i)+k3t*h,I(i)+l3t*h,P(i)+m3t*h);
    l4t=f2(t(i)+h,S(i)+k3t*h,I(i)+l3t*h,P(i)+m3t*h);
    m4t=f3(t(i)+h,S(i)+k3t*h,I(i)+l3t*h,P(i)+m3t*h);
   
    S(i+1)=S(i)+(k1t+2*k2t+2*k3t+k4t)*(h/6);
    I(i+1)=I(i)+(l1t+2*l2t+2*l3t+l4t)*(h/6);
    P(i+1)=P(i)+(m1t+2*m2t+2*m3t+m4t)*(h/6);
end
for i=4:N
%   Adam-Bashforth-Moulton
%   Prediktor
    S(i+1)=S(i)+(h/24)*(55*f1(t(i),S(i),I(i),P(i))-59*f1(t(i-1),S(i-1),I(i-1),P(i-1))...
        +37*f1(t(i-2),S(i-2),I(i-2),P(i-2))-9*f1(t(i-3),S(i-3),I(i-3),P(i-3)));
    I(i+1)=I(i)+(h/24)*(55*f2(t(i),S(i),I(i),P(i))-59*f2(t(i-1),S(i-1),I(i-1),P(i-1))...
        +37*f2(t(i-2),S(i-2),I(i-2),P(i-2))-9*f2(t(i-3),S(i-3),I(i-3),P(i-3)));
    P(i+1)=P(i)+(h/24)*(55*f3(t(i),S(i),I(i),P(i))-59*f3(t(i-1),S(i-1),I(i-1),P(i-1))...
        +37*f3(t(i-2),S(i-2),I(i-2),P(i-2))-9*f3(t(i-3),S(i-3),I(i-3),P(i-3)));
%   Korektor
    S(i+1)=S(i)+(h/24)*(9*f1(t(i+1),S(i+1),I(i+1),P(i+1))+19*f1(t(i),S(i),I(i),P(i))...
        -5*f1(t(i-1),S(i-1),I(i-1),P(i-1))+f1(t(i-2),S(i-2),I(i-2),P(i-2)));
    I(i+1)=I(i)+(h/24)*(9*f2(t(i+1),S(i+1),I(i+1),P(i+1))+19*f2(t(i),S(i),I(i),P(i))...
        -5*f2(t(i-1),S(i-1),I(i-1),P(i-1))+f2(t(i-2),S(i-2),I(i-2),P(i-2)));
    P(i+1)=P(i)+(h/24)*(9*f3(t(i+1),S(i+1),I(i+1),P(i+1))+19*f3(t(i),S(i),I(i),P(i))...
        -5*f3(t(i-1),S(i-1),I(i-1),P(i-1))+f3(t(i-2),S(i-2),I(i-2),P(i-2)));
end

tn=t(1:N+1);
Sn=S(1:N+1)';
In=I(1:N+1)';
Pn=P(1:N+1)';
%Menampilkan hasil iterasi
fprintf('%6.3f %12.6f %18.6f %23.6f\n',[tn; Sn; In; Pn]);
%Menampilkan grafik hasil iterasi
figure(3)
plot(t,S,'r',t,I,'b',t,P,'k');
legend('S','I','P')
title('Metode Adam Bashforth Moulton Orde 4');

grid on