PRAKTIKUM METODA NUMERIK MATERI: Aproksimasi Integral
Tujuan:
Mengimplementasikan berbagai metoda aproksimasi integral dengan menggunakan MATLAB. 1. Menghitung aproksimasi integral dengan formula kuadratur dasar 2. Memahami order konvergensi aproksimasi integral dan interprestasinya
Kompetensi:
3. Menghitung aproksimasi integral dengan formula kuadratur bersusun 4. Menghitung integral data tabular 5. Membandingkan akurasi dan konvergensi beberapa metoda
RANGKUMAN TEORI Diberikan integral tertentu
ˆ
b
I :=
f (x)dx. a
Integral ini diaproksimasi oleh formula kuadratur yang berbentuk
Q(f ) :=
n X
wi f (xi ).
k=1 Kuadratur Dasar (Midpoint, Trapesium dan Simpson):
b−a M (f ) = (b − a)f , T (f ) = (f (a) + f (b)) , 2 b−a a+b S(f ) = f (a) + 4f + f (b) . 6 2 a+b 2
Formula Error:
EM (f ) =
f 00 (ξ) f 00 (ξ) f (4) (ξ) (b − a)3 , ET (f ) = − (b − a)3 , ES (f ) = − (b − a)5 . 24 12 2880
Kuadratur Bersusun:
Pn/2 M (f ) = 2h k=1 f (a + (2k − 1)h). Pn−1 • Trapesium: T (f ) = h2 f (a) + f (b) + 2 k=1 f (a + kh) .
•
Midpoint:
•
Simpson:
S(f ) =
h 3
f (a) + f (b) + 4
P n2
k=1
f (a + (2k − 1)h) + 2
P n2 −1 k=1
f (a + 2kh) .
Formula error untuk kuadratur bersusun diberikan sebagai berikut
ˆ I(f ) − M (f ) :=
b
f (x)dx − M (f ) = a
1
(b − a)h2 00 · f (ξ). 24
ˆ
b
ˆ
a b
f (x)dx − S(f ) =
I(f ) − S(f ) := a Di sini mesh
Proyek 1:
h
(b − a)h2 00 · f (ξ). 12
f (x)dx − T (f ) =
I(f ) − T (f ) :=
h :=
dihitung berdasarkan rumus
(b − a)h4 (4) · f (ξ). 180
b−a n .
Implementasi metoda kuadratur dasar.
´1
1. Aproksimasi integral
0
2
xe−x dx
dihitung manual melalui command window
dengan menggunakan formula kuadratur dasar. 2. Menghitung kesalahan nyata masing-masing aproksimasi (tentunya harus tahu dulu eksaknya). 3. Menyusun m- le untuk implementasi metoda kuadratur dasar dengan input: batas integral
a
dan
b,
output nilai aproksimasi
M, T
dan
S.
Sebagai pengantar, lakukan manual melalui command window untuk menghitung aproksimasi
f 0 (0.5)
dengan menggunakan formula (1), (2) dan (3). De nisikan dulu fungsi
f
sebagai m- le
function y = fun55(x) y = x.*exp(-x.^2); Simpanlah m- le ini dengan nama barkan pada domain
[0, 1]
fun55.m.
Untuk melihat pro l fungsi ini, coba gam-
sbb
>>x=0:0.01:1; >>y = fun55(x); >>plot(x,y) Untuk melihat daerah yang relevan dengan nilai integral yang dimaksud, berikan perintah berikut
>>area(x,y) 1
0.8
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
1
Figure 1: Luas sebagai integral Ketiga formula tersebut dihitung dengan perintah berikut:
>>M=(1-0)*feval('fun55',(0+1)/2); >>T=(1-0)/2*(feval('fun55',0)+feval('fun55',1)); >>S=(1-0)/6*(feval('fun55',0)+4*feval('fun55',(0+1)/2)+feval('fun55',1)); 2
Untuk menghitung error nyata dihitung dulu nilai eksak integral di atas, yaitu
ˆ
1
I=
xe
−x2
ˆ 1 1 2 d(e−x ) − 2 0 i 1h 2 1 − e−x 2 0 1 −1 − e −1 2 0.31606.
dx =
0
= = = Selanjutnya gunakan rumus error nyata
Nyata
EM
= |M − I|,
dan seterusnya juga untuk
trapesium dan Simpson. Secara umum kita dapat menyusun m- le sebagai berikut.
function [M,T,S]=kuad_dasar(a,b) % input: batas integrasi a dan b % output: ketiga metoda kuadratur dasar, M untuk midpoint, % T untuk trapesium dan S untuk Simpson. M = (b-a)*feval('fun',(a+b)/2); T = (b-a)/2*(feval('fun',a)+feval('fun',b)); S = (b-a)/6*(feval('fun',a)+4*feval('fun',(a+b)/2)+feval('fun',b)); %fungsi yang diintegralkan function y = fun(x) y = x*exp(-x^2); Nilai eksak integral fungsi ini dengan batas
I=−
a
dan
b
diberikan oleh
1 exp(−b2 ) − exp(−a2 ) . 2
Gunakan m- le ini untuk menghitung nilai eksak, aproksimasi dan error untuk beberapa domain berikut.
Midpoint Trapesium Simpson Interval Eksak Aproksimasi Errornya Aproksimasi Errornya Aproksimasi Errornya [0, 0.25] [0, 0.5] [0, 1] [0, 2] [0, 2.2] Pada interval
[0, 1]
diperoleh
>> [M,T,S]=kuad_dasar(0,1) M = 0.38940 T = 0.18394 S = 0.32091 >> I = -0.5*(exp(-1^2)-exp(-0^2)) I = 0.31606 >> E_M=abs(M-I);E_T=abs(T-I);E_S=abs(S-I); Silakan angkanya disikan di kolom yang sesuai pada tabel di atas. Lakukan cara yang sama untuk melengkapi tabel ini.
Proyek 2:
Pada proyek ini disusun m- le dengan input batas integral
a, b
dan banyak
subinterval dan output Aproksimasi dengan menggunakan kuadratur bersusun. Masing-masing metoda disusun satu m- le.
3
Formula midpoint dan Simpson di atas mensyaratkan banyak node ganjil atau sedangkan trapesium boleh interval
[a, b]
menjadi
n
n
n
genap,
ganjil. Untuk menyusun m- le yang dimaksud kita pecah
subinterval seragam dengan menetapkan partisi berikut
a =: x0 < x1 < x2 < · · · < xn := b. Karena
a, b
dan
n
diambil sebagai variabel masukan maka ditetapkan
tisi ini dituliskan pada MATLAB dalam bentuk array deks pada MATLAB dimulai dari dengan
x(1)
1,
indeks
0
x = a:h:b.
tidak dikenal. Jadi
pada MATLAB. Himpunan node dengan indeks ganjil
{x(2), x(4), x(6), · · · } pada MATLAB. {x2 , x4 , x6 , · · · } menjadi{x(3), x(5), x(7), · · · }. esuaikan menjadi
h :=
b−a n . Par-
Perlu diketahui in-
a := x0 bersesuaian {x1 , x3 , x5 , · · · } dis-
Begitu juga node indeks genap
Pada midpoint kita cukup menghitung nilai fungsi pada semua node dengan indeks
2h. Perhatikan node dengan x(2:2:end-1). Perintah end merujuk
ganjil, kemudian di jumlahkan dan hasilnya dikalikan dengan indeks ganjil diperoleh dengan menuliskan array
pada indeks terakhir pada array x. Berikut m- le untuk midpoint diterapkan pada fungsi 2
f (x) = xex . function M=midpoint(a,b,n) %masukkan n genap h = (b-a)/n; x = a:h:b;%mendefinisikan partisi seragam if(rem(n,2)~=0)%cek n genap. error('Anda memasukkan nilai n ganjil, seharusnya genap') end node_ganjil = x(2:2:end-1);%mengumpulkan node dengan indeks ganjil. M = 2*h*sum(feval('fun',node_ganjil)); %fungsi yang diintegralkan function y = fun(x) y = x.*exp(x.^2); Pada metoda Simpson kita fokuskan pada titik node interior
x1 , x2 , · · · xn−1 .
Hitung
nilai fungsi pada semua node dengan indeks ganjil, dijumlahkan dan dikalikan dengan
4.
Juga, hitung nilai fungsi pada semua node dengan indeks genap, dijumlahkan dan
dikalikan dengan
2.
Selanjutnya, terapkan rumus Simpson di atas.
function S=simpson(a,b,n) %masukkan n genap h = (b-a)/n; x = a:h:b; if(rem(n,2)~=0) error('Anda memasukkan nilai n ganjil, seharusnya genap') end node_ganjil = x(2:2:end-1); node_genap = x(3:2:end-2); S = (h/3)*(feval('fun',a)+feval('fun',b)+... 4*sum(feval('fun',node_ganjil))+2*sum(feval('fun',node_genap))); %fungsi yang diintegralkan function y = fun(x) y = x.*exp(x.^2); Pada metoda trapesium kita cukup menghitung semua nilai fungsi pada titik interior
x1 , x2 , · · · , xn−1 ,
dijumlahkan dan hasilnya dikalikan dengan
rumus trapesium di atas.
function T=trapes(a,b,n) %masukkan n bulat positif (boleh genap atau ganjil) h = (b-a)/n; x = a:h:b; node_interior = x(2:end-1); 4
2.
Selanjutnya terapkan
T = (h/2)*(feval('fun',a)+feval('fun',b)+... 2*sum(feval('fun',node_interior))); %fungsi integral function y = fun(x) y = x.*exp(x.^2); Melalui command window diberikan perintah berikut.
>> M=midpoint(0,2,20) M = 26.0054 >> S=simpson(0,2,20) S = 26.8058 >> T=trapes(0,2,20) T = 27.2060 Hasil di atas menggunakan
n = 20.
Untuk
n = 15,
diperoleh hasil berikut
>> T=trapes(0,2,15) T = 27.5201 >> M=midpoint(0,2,15) ??? Error using ==> midpoint Anda memasukkan nilai n ganjil, seharusnya genap Ternyata midpoint tidak dapat dieksekusi karena seharusnya n genap, tetapi dimasukkan ganjil.
Proyek 3:
Integral fungsi dalam bentuk tabular.
Pada bidang terapan umumnya rumus fungsi tidak diketahui secara eksplisit, namun hanya tersedia data tabular yang biasanya diperoleh dari hasil pengukuran dalam suatu eksperimen. Dengan asumsi jarak antara dua pengukuran tidak sama dan terdapat fungsi
f
sehingga
yk = f (xk ), k = 0, 1, · · · , n
walaupun
f
tidak diketahui secara eksplisit.
Table 1: Data tabular hasil pengukuran
x0 y0
x1 y1
x2 y2
··· ···
xk yk
··· ···
xn−1 yn−1
xn yn
Untuk mengaproksimasi fungsi yang tersajikan dalam bentuk tabular seperti ini umumnya dilakukan dengan menggunakan teknik interpolasi. Paling sering digunakan interpolasi linear pada setiap pasangan titik berdekatan.
Ilustrasinya diberikan pada
gambar 2. (xk,yk) (xn-1,yn-1)
(x1,y1) (x2,y2)
(xn,yn)
(x0,y0)
a
b
Figure 2: Interpolasi linear data tabular
5
x1 x0 x1 − x0
x2 x1 x − xn
··· ··· ···
xn xn−1 xn − xn−1
y0 y1 y0 + y1
··· ··· ···
y1 y2 y1 + y2
yn−1 yn yn−1 + yn
Konsekuensi dengan teknik ini maka metoda yang cocok adalah trapesium takseragam
n X hk
T1 (f ) =
k=1 di mana
hk := xk − xk−1
2
(f (xk−1 ) + f (xk ))
adalah lebar subinterval
Ik := [xk−1 , xk ]
dan
f (xk ) = yk .
Lagi, cara cerdas menyusun m- le MATLAB adalah dengan menggunakan operasi array. Himpunan para lebar subinterval
hk
dan nilai
yk−1 + yk
dapat diperoleh dengan mudah
sebagai berikut: Bila atas
h := {h1 , h2 , · · · , hn } dan yy := {y0 +y1 , y1 +y2 , · · · , yn−1 +yn } maka formula di dapat
ditulis
T1=sum((h/2).*yy);
dengan
operasi
array
MATLAB
berikut
sangat sederhana. Notasi .* digunakan untuk perkalian antar
komponen. Berikut m- le yang dimaksud.
function T1 = int_tab(datax,datay) %input: datax adalah array para x_k, datay adalah array para y_k a1=datax(2:end);a2=datax(1:end-1);h=a1-a2; %mendefinisikan array h b1=datay(1:end-1);b2=datay(2:end);yy=b1+b2; % mendefinisikan array yy T1 = sum((h/2).*yy); Untuk memastikan m- le ini benar kita terapkan pada kasus seragam.
Sebelumnya
15 subinterval pada interval [0, 2] sebanyak 16 titik partisi seragam pada
sudah dicobakan trapesium seragam dengan mengambil dengan nilai 27.5201. Secara ekuivalen diambil
[0, 2].
Berikut perintah pada command window.
>> >> >> T1
datax=linspace(0,2,16); datay=datax.*exp(datax.^2); T1 = int_tab(datax,datay) = 27.5201
Hasilnya persis sama dengan sebelumnya. Selanjutnya andaikan kita mempunyai pasangan data berikut
>>datax=[0 0.2000 1.4000 1.5000 >>datay=[0 0.2082 9.9391 14.2316
0.5000 1.6000 0.6420 20.6973
0.6000 0.7000 0.9000 1.2000 ... 1.7500 1.8500 1.9000 2.0000]; 0.8600 1.1426 2.0231 5.0648 ... 37.4166 56.6950 70.2355 109.1963];
>> T1 = int_tab(datax,datay) T1 = 27.3677
Sesungguhnya data ini diperoleh dari fungsi secara tidak seragam pada interval dengan ketika diambil seragam.
[0, 2].
y = xex
2
yang dicuplik sebanyak
14
titik
Ternyata nilai integralnya cukup mendekati
Berdasarkan hasil eksperimen ini maka dapat dipa-
hami bahwa pengambilan partisi seragam hanya dimaksudkan untuk kepraktisan, bukan pertimbangan akurasi. Bahkan partisi takseragam dapat memberikan hasil yang lebih akurat.
Pertanyaan untuk tugas responsi Anda diberikan masalah untuk menghitung integral berikut
ˆ
1
exp(x) cos(2πx)dx. 0
6
1. Menyelidiki pola error berdasarkan pengurangan nilai mesh banyak partisi
n
n.
h
atau penambahan
Lengkapilah tabel berikut
Midpoint Trapesium Simpson Eksak Aproksimasi Errornya Aproksimasi Errornya Aproksimasi Errornya
2 4 8 16 32
2. Gambarkan gra k
x=
1 n versus
y = En
untuk masing-masing metoda. Informasi
apa yang Anda dapatkan. Bagaimana kaitannya dengan order konvergensi. 3. Untuk masing-masing
n,
tentukan batas atas kesalahan yang mungkin terjadi dan
bandingkan dengan kesalahan nyatanya. Anda dapat menggunakan bantuan gra k untuk menentukan nilai terkecil
M
untuk setiap
ξ ∈ (0, 1).
K
dan
M
sehingga
|f 00 (ξ)| ≤ K
dan
|f (4) (ξ)| ≤
Susunlah pada sebuah tabel, kemudian gambarkan
gra knya. Apakah hasilnya sesuai dengan latar belakang teoritis yang ada.
7
