Metode Numerik 1 Imam Fachruddin (Departemen Fisika, Universitas Indonesia)
Daftar Pustaka: •
P. L. DeVries, A First Course in Computational Physics (John Wiley & Sons, Inc., New York, 1994)
•
W. H. Press, et. al., Numerical Recipes in Fortran 77, 2nd Ed. (Cambridge University Press, New York, 1992) (online / free download: http://www.nrbook.com/a/bookfpdf.php)
Isi •
akar fungsi
•
solusi sistem persamaan linear
•
fitting dengan least square
•
interpolasi
•
integrasi
•
persamaan differensial
Akar Fungsi f(x) = 0
x=? akar fungsi f(x)
Contoh:
x-
1 =0 x
x2 = 1
3x2 = 6 - 7x
x = 1 dan -1
3x2 + 7x - 6 = (3x - 2)(x + 3) = 0
Pada dua contoh di atas akar fungsi dapat dicari secara analitik. Secara umum, tidak selalu begitu keadaannya.
x = 2/3 dan -3
Problem: Sebuah lampu dipasang di pinggir sebuah piringan berjari-jari 10 cm. Sebuah plat bercelah sempit diletakkan di dekat piringan itu. Tepat di belakang celah itu dipasang sebuah sensor cahaya yang menghadap tegak lurus ke celah. Piringan diputar konstan 1 rad/s dan plat beserta sensor digeser lurus konstan 10 cm/s. Saat ini posisi celah dan lampu seperti pada gambar 1. Kapan sensor cahaya menerima cahaya terbanyak? Sensor menerima cahaya terbanyak pada saat posisi lampu dan celah membentuk garis tegak lurus terhadap plat, seperti pada gambar 2. x = r cos (ωt) = vt r
lampu
gambar 1
gambar 2 celah
sensor
ω plat v
cos (t) = t
?
Plot cos(x) dan x: Grafik ini menunjukkan bahwa cos(x) = x pada x sedikit kurang dari 0.75.
Bisakah lebih akurat lagi?
Cari secara numerik akar fungsi dari f(x) = cos(x) - x
Bisection Prinsip: Kurung akar fungsi di antara dua batas, lalu paruh batas itu terus menerus sampai batas itu sedemikian sempit dan dengan demikian lokasi akar fungsi diketahui dengan keakuratan tertentu. Langkah: 1. Perkirakan akar fungsi (bisa dengan cara memplot fungsi). 2. Tentukan batas awal yang mengurung akar fungsi.
akar fungsi
a
df e c
b
3. Belah dua daerah berisi akar fungsi itu. 4. Tentukan daerah yang berisi akar fungsi. 5. Ulangi langkah 3 dan 4 sampai dianggap cukup. 6. Tentukan akar fungsi.
Batas e, f atau nilai di tengahnya bisa dipilih sebagai akar fungsi.
•
Menentukan daerah yang berisi akar fungsi: Jika z merupakan akar fungsi, maka f(x < z) dan f(x > z) saling berbeda tanda.
a
c
b
f(x)
f(a)*f(c) negatif, berarti di antara a & c ada akar fungsi.
z
f(b)*f(c) positif, berarti di antara b & c tidak ada akar fungsi •
Menentukan kapan proses pencarian akar fungsi berhenti: Proses pencarian akar fungsi dihentikan setelah keakuratan yang diinginkan dicapai, yang dapat diketahui dari kesalahan relatif semu.
kesalahan relatif semu =
perkiraan sebelum - perkiraan berikut perkiraan berikut
x
Kesalahan kesalahan mutlak = | perkiraan – nilai sebenarnya |
kesalahan relatif =
perkiraan – nilai sebenarnya nilai sebenarnya
Dalam perhitungan numerik, nilai sebenarnya justru sering tidak diketahui, yang didapat hanya perkiraan terbaik. Karena perkiraan langkah berikut dianggap lebih akurat, yaitu lebih mendekati nilai sebenarnya, maka kesalahan yang dihitung yaitu: kesalahan mutlak semu = | perkiraan sebelum – perkiraan berikut |
kesalahan relatif semu =
perkiraan sebelum - perkiraan berikut perkiraan berikut
Newton-Raphson Prinsip: Buat garis singgung kurva f(x) di titik di sekitar akar fungsi. Titik tempat garis singgung itu memotong garis nol ditentukan sebagai akar fungsi.
akar fungsi yang diperoleh
f(x)
x a akar fungsi sebenarnya
p(x) = garis singgung kurva f(x) di titik f(a)
f(x) c
x
a
Diperoleh:
p(x)
p(x) = f(a) + (x − a)f'(a) (f’(a) turunan pertama f(x) pada x = a)
p(c) = 0
c=a−
f(a) f'(a)
Langkah:
f(x)
1. Perkirakan akar fungsi.
c
2. Buat garis singgung pada titik sesuai akar fungsi yang diperkirakan itu, lalu cari titik potongnya dengan garis nol.
x
a akar fungsi sebenarnya
3. Titik potong itu merupakan perkiraan akar fungsi baru. 4. Ulangi langkah 2 dan 3 sampai dianggap cukup. 5. Titik potong garis nol dan garis singgung kurva yang terakhir dinyatakan sebagai akar fungsi.
f(x) a c
c=a−
f(a) f'(a)
x
1
f(x) Contoh perkiraan akar fungsi awal yang baik perkiraan akar fungsi makin mendekati akar fungsi sebenarnya.
x
2 1 f(x) Contoh perkiraan akar fungsi awal yang buruk perkiraan akar fungsi makin menjauhi akar fungsi sebenarnya.
x
2
Menghitung akar fungsi dengan metode Newton-Raphson:
xi+1 = xi −
f(xi ) f'(xi )
(i = 0, 1, 2, ...;
x0 = a )
kesalahan relatif semu:
∆ rel =
xi − xi+1 xi+1
Penghitungan dihentikan jika kesalahan relatif semu sudah mencapai / melampaui batas yang diinginkan.
Kecepatan Konvergensi Pencarian akar fungsi dimulai dengan perkiraan akar fungsi yang pertama, lalu diikuti oleh perkiraan berikutnya dan seterusnya sampai perkiraan yang terakhir, yang kemudian dinyatakan sebagai akar fungsi hasil perhitungan tersebut. Proses itu harus bersifat konvergen yaitu, selisih perkiraan sebelum dari yang setelahnya makin lama makin kecil. Setelah dianggap cukup, proses pencarian akar fungsi berhenti.
x2 − x1 > x3 − x2 > x4 − x3 ... xn − xn −1 ≤ ε (ε = bilangan kecil) Kecepatan konvergensi sebuah proses yaitu, kecepatan proses itu untuk sampai pada hasil akhir.
Contoh pencarian akar fungsi dengan metode Bisection: akar fungsi
x a
b
x1 x3
x2
x4
Jika εi ≡ xi+1 − xi , maka dari gambar diperoleh:
ε1 = x2 − x1 , ε2 = x3 − x2 , ε3 = x4 − x3 ε2 = 21 ε1 , ε3 = 21 ε2 Kecepatan konvergensi bersifat linear:
εi+1 = 21 εi
Pada metode Newton-Raphson:
xi+1 = xi −
f(xi ) f'(xi )
ekspansi deret Taylor:
εi ≡ xi − xi+1 =
f(xi ) f'(xi )
εi+1 =
f(xi+1 ) =? f'(xi+1 )
f(xi+1 ) = f(xi − εi ) = f(xi ) − εif'(xi ) + 21 εi2f''(xi ) − ... f'(xi+1 ) = f'(xi − εi ) = f'(xi ) − εif''(xi ) + ...
f(xi ) − εif'(xi ) + 21 εi2f''(xi ) − ... εi+1 = f'(xi ) − εif''(xi ) + ... f(xi ) − εif'(xi ) + 21 εi2f''(xi ) ≅ f'(xi ) ≅
f''(xi ) 2 εi 2f'(xi )
Kecepatan konvergensi pada metode Newton-Raphson bersifat kurang lebih kuadratik:
εi+1 ≅
f''(xi ) 2 εi 2f'(xi )
Dengan begitu, metode Newton-Raphson lebih cepat dari metode Bisection.
Contoh hasil pencarian akar fungsi untuk soal cos(x) = x:
metode
akar
f(akar)
jumlah langkah
Bisection
0.7390795
9.3692161E-06
12
Newton-Raphson
0.7390851
-7.7470244E-09
4
Keterangan:
•
Pencarian akar berhenti jika kesalahan relatif semu sama atau kurang dari 1.0E-05.
•
Batas awal kiri dan kanan untuk metode Bisection 0.72 dan 0.75.
•
Perkiraan akar fungsi pertama untuk metode Newton-Raphson 0.72.
Solusi Sistem Persamaan Linear Sistem persamaan linear:
a11x1
+ a12x2
+ a13x3 ⋯ + a1n xn
= b1
a21x1 a31x1
+ a22x2 + a32x2
+ a23x3 ⋯ + a2n xn + a33x3 ⋯ + a3n xn
= b2 = b3
⋮
⋮
an1x1
+ an2x2
⋮
⋱
⋮
+ an3x3 ⋯ + ann xn
⋮ = bn
n buah persamaan dengan n buah unknown
xj
aij dan bi
diketahui i, j = 1, 2, …, n
xj = ?
Soal:
Jawab:
2x − 3y + 2z = −6
(1)
− x + 2y − 3z = 2
(2)
x+ y− z =0
(3)
2x −
2x −
3y + 2z = −6
(1)
0.5y − 2z = −1
(2)
2.5y − 2z = 3
(3)
3y + 2z = −6
(1)
0.5y − 2z = −1
(2)
8z = 8
eliminasi x: pers. (2) + 0.5 pers. (1) pers. (3) – 0.5 pers. (1)
eliminasi y: pers. (3) – 5 pers. (2)
(3)
z =1 − 1 + 2z =2 0.5 − 6 + 3y − 2z = −1 x= 2 y=
3 persamaan dan 3 unknown
substitusi mundur: pers. (3) mencari z pers. (2) mencari y pers. (1) mencari x
Dalam bentuk matriks: Soal:
2 x − 6 2 −3 2 − 3 y = 2 −1 1 1 − 1 z 0
Jawab:
2 x − 6 2 −3 0 0.5 − 2 y = − 1 0 2.5 − 2 z 3 2 x − 6 2 −3 0 0.5 − 2 y = − 1 0 0 8 z 8
z =1 − 1 + 2z =2 y= 0.5 − 6 + 3y − 2z = −1 x= 2
Eliminasi Gauss Metode Eliminasi Gauss mencari solusi sebuah sistem persamaan linear dengan cara seperti ditunjukkan pada contoh sebelum ini:
a11 a21 a 31 ⋮ a n1
a12 a22 a32 ⋮ an2
a13 ⋯ a1n x1 b1 a23 ⋯ a2n x2 b2 a33 ⋯ a3n x3 = b3 ⋮ ⋱ ⋮ ⋮ ⋮ an3 ⋯ ann xn bn
aij(0) = aij , bi(0) = bi (k) ij
a
(k) i
b
(k -1) ij
aik(k-1) (k-1) − (k-1) akj akk
(k −1) i
aik(k-1) (k −1) − (k-1) bk akk
=a
=b
(k = 1, ..., n − 1;
a 0 0 ⋮ 0
(0) 11
(0) 12 (1) 22
a a 0 ⋮ 0
(0) 13 (1) 23 (2) 33
a a a ⋮ 0
⋯ a x1 b ⋯ a x2 b ⋯ a x3 = b ⋱ ⋮ ⋮ ⋮ (n-1) (n-1) x ⋯ ann n bn (0) 1n (1) 2n (2) 3n
(0) 1 (1) 2 (2) 3
i = k + 1, ..., n; j = k, ..., n)
aij(m) , bi(m) ≡ aij , bi pada langkah ke m halaman berikut
Substitusi mundur:
a11(0) 0 0 ⋮ 0
(0) a12 (1) a22 0
⋮ 0
(0) a13 ⋯ a1n(0) x1 b1(0) (1) (1) (1) a23 ⋯ a2n x2 b2 (2) (2) x3 = b3(2) a33 ⋯ a3n ⋮ ⋱ ⋮ ⋮ ⋮ (n-1) (n-1) x 0 ⋯ ann n bn
bn(n-1) xn = (n-1) ann (n- j-1) n-j
b xn − j =
n
∑a
−
k =n- j+1 (n- j-1) n - j,n - j
a
(n- j-1) n - j,k k
x
(j = 1, ..., n − 1)
A
X = B
atau
AX = B
Jadi, metode Eliminasi Gauss terdiri dari dua tahap: 1. triangulasi: mengubah matriks A menjadi matriks segitiga (matriks B dengan begitu juga berubah)
=
=
2. substitusi mundur: menghitung x mengikuti urutan terbalik, dari yang terakhir ( xn ) sampai yang pertama ( x1 )
Kasus Beberapa Sistem Persamaan Linear Pada kasus yang lebih umum bisa saja terdapat beberapa sistem persamaan linear dengan nilai B yang berlainan, namun memiliki nilai A yang sama. Dalam bentuk matriks sistem seperti ini dituliskan sebagai:
A
Keterangan: •
•
X
=
B
atau AX = B
A matriks n x n, X dan B matriks n x m, dengan m = jumlah sistem persamaan linear, n = jumlah persamaan / unknown dalam tiap sistem persamaan tersebut Tiap kolom matriks X merupakan solusi untuk kolom yang sama pada matriks B.
Langkah dan rumus pada metode Eliminasi Gauss berlaku sama untuk kasus ini. Hanya saja, di sini matriks X dan B terdiri dari beberapa kolom, bukan hanya satu.
Contoh dua sistem persamaan linear yang memiliki nilai A sama tapi B berbeda.
a11 a21 a 31 ⋮ a n1
a12 a22 a32 ⋮ an2
a13 a23 a33 ⋮ an3
⋯ a1n x11 b11 ⋯ a2n x21 b21 ⋯ a3n x31 = b31 ⋱ ⋮ ⋮ ⋮ ⋯ ann xn1 bn1
a11 a21 a 31 ⋮ a n1
a12 a22 a32 ⋮ an2
a11 a21 a 31 ⋮ a n1
a13 ⋯ a1n x11 a23 ⋯ a2n x21 a33 ⋯ a3n x31 ⋮ ⋱ ⋮ ⋮ an3 ⋯ ann xn1
a12 a22 a32 ⋮ an2
a13 a23 a33 ⋮ an3
⋯ a1n x12 b12 ⋯ a2n x22 b22 ⋯ a3n x32 = b32 ⋱ ⋮ ⋮ ⋮ ⋯ ann xn2 bn2
x12 b11 b12 x22 b21 b22 x32 = b31 b32 ⋮ ⋮ ⋮ xn2 bn1 bn2
Metode Eliminasi Gauss: •
rumus triangulasi:
aij(0) = aij , bir(0) = bir
aij(k) = aij(k-1) − bir(k) •
(k -1) ik (k -1) kk
a a
(i, j = 1, ..., n; r = 1, ..., m)
akj(k-1) (k = 1, ..., n − 1;i = k + 1, ..., n;
(k -1) a = bir(k −1) − ik(k-1) bkr(k −1) akk
aij(m) , bir(m) ≡ aij , bir pada langkah ke m
j = k, ..., n;r = 1, ..., m)
rumus substitusi mundur:
bnr(n-1) xnr = (n-1) ann (n- j-1) n- j,r
b xn − j,r =
(r = 1, ..., m) n
∑a
−
k =n - j+1 (n- j-1) n - j,n - j
a
(n- j-1) n- j,k kr
x
(j = 1, ..., n − 1; r = 1, ..., m)
Catatan: Dalam rumus-rumus metode Eliminasi Gauss terdapat pembagian oleh elemen diagonal matriks yaitu, oleh elemen diagonal matriks A. Jika secara kebetulan elemen diagonal itu nol, maka akan timbul error. Karena itu, pada setiap langkah dalam proses triangulasi matriks A perlu dilakukan pemeriksaan, apakah elemen matriks A yang bersangkutan sama dengan nol. Jika bernilai nol, maka baris berisi elemen diagonal nol itu harus ditukar dengan salah satu baris setelahnya, sehingga elemen diagonal menjadi bukan nol. Perubahan baris pada matriks A harus disertai perubahan baris yang sama pada matriks B.
Soal:
Jawab:
1 3 x1 2 2 −4 2 3 − 2 x2 2 −1 = 3 −4 1 2 x3 2 1 − 3 −1 5 x4 2
baris 2 ditukar dengan baris 3
1 3 x1 2 2 − 4 0 3.5 − 0.5 x2 3 0 = 0 2 − 0.5 − 2.5 x3 −1 0 − 1 − 1 .5 x 1 3 . 5 4
1 3 x1 2 2 − 4 2 − 0.5 − 2.5 x2 − 1 0 = 0 0 3.5 − 0.5 x3 3 0 0 0 2 x4 2
x4 = 1 x3 =
3 + 0.5x4 =1 3.5
1 3 x1 2 2 − 4 2 − 0.5 − 2.5 x2 − 1 0 = 0 0 3.5 − 0.5 x3 3 0 − 1 − 1 .5 x 1 3 . 5 4
1 3 x1 2 2 − 4 2 − 0.5 − 2.5 x2 − 1 0 = 0 0 3.5 − 0.5 x3 3 0 0 − 1.75 2.25 x4 0.5
− 1 + 0.5x3 + 2.5x4 =1 2 2 + 4x2 − x3 − 3x4 =1 x1 = 2 x2 =
Data Fitting dengan Metode Least Square Keterangan:
f(x) p(x)
x
•
f(xi ) mewakili data; i = 1, …, N; N = jumlah data
•
p(x) merupakan fungsi yang dicocokkan (fitted) terhadap data f( xi )
Sifat fitting: tidak selalu p(xi ) = f( xi ) untuk semua xi .
Prinsip penentuan fungsi p(x): •
p(x) merupakan polinomial orde m: m
p(x) = a0 + a1x + a2x + a3x + ... + amx = ∑ ajx j 2
3
m
j= 0
(Secara umum, p(x) juga bisa merupakan polinomial bentuk yang lain seperti, polinomial Legendre.) •
Selisih antara p(x) dan f(x) untuk titik data tertentu: m
∆i = f(xi ) − p(xi ) = f(xi ) − ∑ ajxij
(i = 1, ..., N )
j= 0
•
Jumlah kuadrat selisih antara p(x) dan f(x) untuk semua titik data: m 2 j 2 S = ∑ ∆i = ∑ (f(xi ) − p(xi ) ) = ∑ f(xi ) − ∑ ajxi i=1 i=1 i=1 j= 0 N
N
N
2
Fungsi p(x) ditentukan dengan mencari nilai aj (j = 0, …, m) yang membuat S bernilai minimum.
Titik Minimum g(x) g(a) merupakan titik minimum jika:
dg(x) dx x = a
a
2 dg(x) d g(x) = 0 dan >0 2 dx x = a dx x = a
x
Spesial: fungsi kuadratik g(x) = ax2 + bx + c dg(x) = 2ax + b dx 2 d g(x) = 2a dx 2
g(x) memiliki satu titik minimun jika a > 0 atau sebaliknya satu titik maksimum jika a < 0.
S merupakan fungsi kuadratik dalam aj (j = 0, …, m): 2
m N m j 2 2 2j S(a0 , ..., am ) = ∑ f(xi ) − ∑ ajxi = ∑ ∑ aj xi + ... + f (xi ) j= 0 i=1 j= 0 i=1 N
(
N m k ∂S(a0 , ..., am ) j = −2∑ f(xi ) − ∑ ajxi xi ∂ak i=1 j= 0
N ∂ 2S(a0 , ..., am ) = 2∑ xi2k > 0 2 ∂ak i=1
)
(k = 0, ..., m)
(k = 0, ..., m)
S memiliki satu titik minimum pada nilai aj (j = 0, …, m) tertentu.
Mencari aj (j = 0, …, m): N m k ∂S(a0 , ..., am ) j = −2∑ f(xi ) − ∑ ajxi xi = 0 ∂ak i=1 j= 0
N N j+k ∑ xi aj = ∑ f(xi )xik ∑ j= 0 i=1 i=1 m
Definisikan:
N
(k = 0, ..., m)
(k = 0, ..., m)
N
ckj ≡ ∑ x
bk ≡ ∑ f(xi )xik
j +k i
i=1
i=1 m
maka diperoleh sebuah sistem persamaan linear:
∑c
a = bk
kj j
(k = 0, ..., m)
j= 0
dalam bentuk matrik:
C
A = B
atau
CA = B
Jadi, aj (j = 0, …, m) diperoleh sebagai solusi persamaan linear CA = B.
Contoh: Terdapat tiga data f(x) yaitu, f(1) = 30, f(2) = 70 dan f(3) = 120. Cari fungsi p(x) yang dapat melukiskan data itu. Dari data itu jelas p(x) bukan fungsi linear. Jadi, dicoba fungsi kuadratik:
f(x)
p(x)
120
p(x) = a0 + a1x + a2x2 Sistem persamaan linier untuk mencari aj :
3 6 14 a0 220 6 14 36 a = 530 1 14 36 98 a 1390 2 3 6 14 a0 220 0 1 4 a1 = 45 0 0 1 a 5 2 Jadi,
p(x) = 5x(x + 5 )
70 30 1
2
3
a0 0 a1 = 25 a 5 2 Cek:
p(1) = 30, p(2) = 70, p(3) = 120
! K O
x
Contoh: Kuat medan listrik E di sekitar sebuah benda berbentuk lempeng diukur pada jarak 10 cm dari pusat massanya dan arah yang bervariasi. Arah dinyatakan dalam sudut θ terhadap sumbu y yang ditetapkan sebelum pengukuran. Diperoleh data sebagai berikut: θ [derajat]
E [V/cm]
10 15 20 25 30 35 40 45 50 55
0.01794775 0.03808997 0.05516225 0.05598281 0.04795629 0.04807485 0.06273566 0.07853982 0.07395442 0.04201338
Cari fungsi p(x) yang dapat melukiskan data itu.
y
θ
E
Dicoba beberapa polinomial dengan orde berbeda, diperoleh: a0 = 8.983713484853211E - 03 a1 = 1.324478388111303E - 03
m = 3: a2 = 3.487808787880805E - 05 a3 = - 8.085809790211842E - 07 S = 1.0339E - 03
a0 = - 3.557800654975570E - 02 a1 = 1.061996221844471E - 03 a2 = 8.802185976358352E - 04
m = 5: a3 = - 5.862332690401015E - 05 a4 = 1.362046192596346E - 06 a5 = - 1.063951754163944E - 08 S = 8.1573E - 05
a0 = - 1.757260839248139E - 02
m = 9:
a1 = 1.596300085173997E - 02 a2 = - 3.402768734407800E - 03
a0 = 1.864754537649403E - 01
a3 = 3.358961098305538E - 04
a2 = 4.007658091692495E - 03
a4 = - 1.368895999268855E - 05 a5 = 1.132254508386570E - 07 a6 = 8.262829873458547E - 09 a7 = - 2.741786330789355E - 10 a8 = 3.317446724324134E - 12 a9 = - 1.459511835946927E - 14 S = 1.7528E - 11
a1 = - 4.631839872868015E - 02 a3 = - 8.985715636865594E - 05
m = 7: a4 = - 3.230489224228010E - 06 a5 = 1.912806006890119E - 07 a6 = - 3.252863805243949E - 09 a7 = 1.876184315740421E - 11 S = 3.1629E - 07
Interpolasi Keterangan:
f(x) p(x)
x
•
f(xi ) mewakili data; i = 1, …, N; N = jumlah data
•
p(x) merupakan fungsi interpolasi berdasarkan data f(xi )
Sifat interpolasi: p(xi ) = f(xi ) untuk semua xi .
Interpolasi Lagrange Digunakan p(x), suatu polinomial berorde m = N – 1, dengan N = jumlah data:
p(x) = a0 + a1x + a2x2 + ... + aN-1x N-1 ≅ f(x) Nilai ai (i = 0, …, N-1) ditetukan dengan menetapkan bahwa untuk semua titik data: (i = 1, ..., N ) p(xi ) = f(xi ) Jadi, diperoleh persamaan linear:
p(x1 ) = a0 + a1 x1 + a2 x12 + ... + aN-1 x1N-1 = f(x1 ) p(x2 ) = a0 + a1x2 + a2 x22 + ... + aN-1x2N-1 = f(x2 ) p(x3 ) = a0 + a1x3 + a2x32 + ... + aN-1x3N-1 = f(x3 ) ... p(xN ) = a0 + a1xN + a2xN2 + ... + aN-1 xNN-1 = f(xN ) dan ai (i = 0, …, N-1) diperoleh sebagai solusi dari persamaan linear itu.
p(x1 ) = a0 + a1x1 = f(x1 )
N = 2:
p(x2 ) = a0 + a1 x2 = f(x2 ) a0 = −
x2f(x1 ) − x1f(x2 ) x1 − x2
a1 =
f(x1 ) − f(x2 ) x1 − x2
x − x2 x − x1 f(x1 ) + f(x2 ) p(x) = x1 − x2 x2 − x1
p(x1 ) = a0 + a1 x1 + a2 x12 = f(x1 )
N = 3:
p(x2 ) = a0 + a1x2 + a2 x22 = f(x2 ) p(x3 ) = a0 + a1x3 + a2x32 = f(x3 ) a0 =
(x2 − x3 )x2 x3f(x1 ) + (x3 − x1 )x3x1f(x2 ) + (x1 − x2 )x1x2f(x3 ) (x2 − x3 )x12 + (x3 − x1 )x22 + (x1 − x2 )x32
(x22 − x32 )f(x1 ) + (x32 − x12 )f(x2 ) + (x12 − x22 )f(x3 ) a1 = − (x2 − x3 )x12 + (x3 − x1 )x22 + (x1 − x2 )x32 a2 =
(x2 − x3 )f(x1 ) + (x3 − x1 )f(x2 ) + (x1 − x2 )f(x3 ) (x2 − x3 )x12 + (x3 − x1 )x22 + (x1 − x2 )x32
x − x1 x − x2 x − x2 x − x3 x − x1 x − x3 f(x1 ) + f(x2 ) + f(x3 ) p(x) = x − x x − x x − x x − x x − x x − x 2 1 3 3 2 1 2 1 2 3 1 3
N
Secara umum, untuk N data rumus interpolasi Lagrange:
p(x) = ∑ l(x, xi )f(xi ) i=1
x − xj l(x, xi ) = ∏ j≠i xi − xj
Untuk x = xk (k = 1, …, N):
xi − xj = 1, (i = k) ∏ xk − xj j≠i xi − xj = l(xk , xi ) = ∏ j≠i xi − xj ... xk − xk ... = 0, (i ≠ k) x −x j i l(xk , xi ) = δik
p(xk ) = f(xk )
Perlukah memakai semua N data yang ada? Pada bagian sebelum ini interpolasi menggunakan seluruh N data f( xi ) yang tersedia, yang berarti menggunakan polinomial p(x) berorde N-1. Kini, misal N = 4 dan x berada di sekitar x4 , maka diperoleh:
x − x2 x − x3 x − x4 l(x, x1 ) = x1 − x2 x1 − x3 x1 − x4
x − x1 x − x3 x − x4 l(x, x2 ) = x2 − x1 x2 − x3 x2 − x4
x − x1 x − x2 x − x4 l(x, x3 ) = x3 − x1 x3 − x2 x3 − x4
x − x1 x − x2 x − x3 l(x, x4 ) = x4 − x1 x4 − x2 x4 − x3
Dapat dilihat bahwa, l(x, x1 ) < l(x, x2 ) < l(x, x3 ) < l(x, x4 ) . Ini berarti, semakin jauh dari x pengaruh data f( xi ) semakin kecil dalam menentukan nilai p(x). Data yang penting yaitu yang berada di sekitar titik x. Karena itu, cukup data-data di sekitar titik x yang digunakan. Dengan kata lain, untuk interpolasi cukup digunakan polinomial p(x) berorde rendah, contoh berorde 3 (fungsi kubik).
Interpolasi Lagrange Kubik Interpolasi Lagrange Kubik menggunakan polinomial p(x) berorde 3 sebagai fungsi interpolasi:
p(x) = a0 + a1x + a2x2 + a3x3 ≅ f(x) Untuk mencari nilai aj (j = 0, 1, 2, 3) diperlukan 4 data f( xi ) di sekitar x:
f(x0 ), f(x1 ), f(x2 ), f(x3 )
(xi ≤ x ≤ xi+1 ;
x0 = xi-1 , x1 = xi , x2 = xi+1 , x3 = xi+2 )
untuk membentuk sistem persamaan linear:
a0 + a1xj + a2xj2 + a3xj3 = f(xj )
(j = 0, 1, 2, 3)
Langkah pertama dengan begitu, menentukan xj (j = 0, 1, 2, 3) dengan melihat posisi x di antara titik data xi (i = 1, …, N). Diperoleh 3
p(x) = ∑ l(x, xj )f(xj ) j= 0
x − xk l(x, xj ) = ∏ k ≠ j xj − xk
Interpolasi Multidimensi Jika data bergantung pada lebih dari satu variabel, maka dilakukan interpolasi multidimensi. Metode interpolasi yang telah disampaikan bisa dipakai untuk melakukan interpolasi multidimensi. Sebagai contoh di sini ditunjukkan interpolasi 2 dimensi. Untuk dimensi lebih tinggi berlaku cara yang sama.
n
m
i=1
j=1
p(x, y) = ∑ S(x, xi )∑ S(y, yj )f(xi , yj )
Pada contoh di atas, interpolasi menggunakan (n x m) data f(x,y). Interpolasi dilakukan per dimensi: Untuk satu titik data x tertentu dilakukan interpolasi di sepanjang sumbu y, hal yang sama dilakukan untuk semua titik data x yang lain. Prinsip yang sama berlaku untuk interpolasi berdimensi lebih tinggi.
Contoh, interpolasi Lagrange kubik:
3
3
i= 0
j= 0
p(x, y) = ∑ l(x, xi )∑ l(y, yj )f(xi , yj ) x − xk l(x, xi ) = ∏ k ≠i xi − xk y − ys l(y, yj ) = ∏ s ≠ j yj − ys
Kembali ke contoh problem least square: Kuat medan listrik E di sekitar sebuah benda berbentuk lempeng diukur pada jarak 10 cm dari pusat massanya dan arah yang bervariasi. Arah dinyatakan dalam sudut θ terhadap sumbu y yang ditetapkan sebelum pengukuran. Diperoleh data sebagai berikut:
θ [derajat]
E [V/cm]
10 15 20 25 30 35 40 45 50 55
0.01794775 0.03808997 0.05516225 0.05598281 0.04795629 0.04807485 0.06273566 0.07853982 0.07395442 0.04201338
y
θ
Dengan interpolasi, cari nilai p(x) di sepanjang titik data.
E
Integrasi Menghitung luas daerah di bawah kurva: f(x)
f(x) analitik
numerik
b
∑ w f(x )
∫ f(x) dx
i
i
a
a
b b
N
a
i=1
I = ∫ f(x) dx ≅ ∑ wif(xi )
i
x
a
b
x
Integral numerik sering disebut juga sebagai quadrature; integrasi numerik disebut sebagai integrasi dgn menjumlah quadrature.
Meski tidak terlihat pada rumus akhir, pada integrasi numerik integrand f(x) diinterpolasi dengan suatu polinomial:
b
N
a
i=1
I = ∫ f(x) dx ≅ ∑ wif(xi )
f(x) ≅ p(x)
Akan dibahas: •
quadrature trapezoid
•
quadrature Simpson
polinomial
Quadrature Trapezoid Kurva integrand f(x) diinterpolasi dengan sebuah garis lurus (f(x) diinterpolasi dengan fungsi linier / polinomial orde 1): b
b
N
I = ∫ f(x) dx ≅ ∫ p(x) dx = ∑ wip(xi ), a
a
p(x) = r + sx
i=1
f(x) Untuk menarik garis lurus diperlukan minimal 2 titik, dipilih titik f(a) dan f(b):
p(x)
p(a) = f(a), p(b) = f(b) b
∫ p(x) dx a
a
b
x
Dengan diketahui hanya p(a) dan p(b) (r dan s tidak dicari), maka integrasi numerik dikerjakan untuk N = 2: b
2
∫ p(x) dx = ∑ w p(x ) = w p(x ) + w p(x ) = w p(a) + w p(b) i
a
1
1
2
2
1
w1 , w2 = ?
2
i=1
Mencari w1 dan w2 :
p(x) = r + sx
i
b
∫ (r + sx) dx = w (r + sa) + w (r + sb) 1
2
a
1 r(b - a) + s(b2 − a2 ) = r(w1 + w2 ) + s(aw1 + bw2 ) 2 w1 + w2 = b - a 1 aw1 + bw2 = (b2 − a2 ) 2 b
Rumus quadrature trapezoid:
I = ∫ f(x) dx ≅ a
h (f(a) + f(b) ) 2
w1 = w2 =
1 (b − a) 2
(h = b − a)
luas trapezoid (lihat gambar)
Quadrature Simpson & Boole Cara yang sama seperti pada quadrature trapezoid bisa dipakai untuk polinomial p(x) orde lebih tinggi. Contoh, quadrature Simpson memakai p(x) fungsi kuadratik / polinomial orde 2 untuk menginterpolasi integrand f(x): c
c
N
I = ∫ f(x) dx ≅ ∫ p(x) dx = ∑ wip(xi ), a
a
p(x) = r + sx + tx2
i=1
f(x) Untuk membuat kurva kuadratik diperlukan minimal 3 titik, dipilih titik f(a), f(b) dan f(c):
p(x)
p(a) = f(a), p(b) = f(b),
b
∫ p(x) dx
p(c) = f(c)
a
dengan
b=
a+c 2
a
b
c
x
Integrasi numerik dikerjakan untuk N = 3: c
3
∫ p(x) dx = ∑ w p(x ) = w p(a) + w p(b) + w p(c) i
a
i
1
2
w1 , w2 , w3 = ?
3
i=1
Mencari w1 , w2 , w3:
p(x) = r + sx + tx2 c
2 2 2 (r + sx + tx ) dx = w (r + sa + ta ) + w (r + sb + tb ) 1 2 ∫ a
r(c - a) +
1 1 s(c2 − a2 ) + t(c3 − a3 ) = 2 3
w1 + w2 + w3 = c - a 1 aw1 + bw2 + cw3 = (c2 − a2 ) 2 1 a2w1 + b2w2 + c2w3 = (c3 − a3 ) 3
+ w3 (r + sc + tc2 ) r(w1 + w2 + w3 ) + s(aw1 + bw2 + cw3 ) + t(a2w1 + b2w2 + c2w3 )
w1 = w3 = w2 =
1 (c − a) 6
2 (c − a) 3
c
Diperoleh Rumus quadrature Simpson:
I = ∫ f(x) dx ≅ a
dengan h =
h (f(a) + 4f(b) + f(c) ) 3
c−a yaitu jarak antar titik xi tempat f(x) dihitung: h = b − a = c − b 2
Dengan cara yang sama, menggunakan p(x) polinomial orde 3 diperoleh rumus quadrature Simpson 3 8 : d
I = ∫ f(x) dx ≅ a
3h (f(a) + 3f(b) + 3f(c) + f(d) ) 8
d-a = b − a = c − b = d − c h = 3
dan dengan p(x) polinomial orde 4 rumus quadrature Boole: e
2h (7f(a) + 32f(b) + 12f(c) + 32f(d) + 7f(e) ) I = ∫ f(x) dx ≅ 45 a
e-a h = 4
= b− a = c −b = d−c = e − d
Integrasi Komposit Polinomial orde rendah memadai untuk menginterpolasi sebuah fungsi dalam daerah yang sempit. Untuk daerah yang lebar diperlukan orde yang lebih tinggi. Alternatif lain yaitu, membagi daerah fungsi yang lebar itu dalam beberapa daerah yang sempit, lalu di tiap daerah yang sempit itu digunakan polinomial orde rendah untuk interpolasi. Quadrature trapezoid dan Simpson pada dasarnya memadai untuk daerah integrasi yang sempit, namun dengan membagi daerah integrasi dalam beberapa daerah yang sempit, maka quadrature trapezoid dan Simpson bisa dipakai juga untuk daerah integrasi yang lebar. Integral total merupakan jumlah semua integral untuk daerah yang sempit. Integrasi seperti ini disebut integrasi komposit. Bergantung pada integrand f(x), daerah integrasi yang lebar bisa dibagi dalam beberapa daerah sempit yang sama atau berbeda panjang. Juga, semua integral untuk daerah yang sempit bisa dihitung menurut rumus quadrature yang sama, misal semuanya trapezoid, atau berbeda-beda, sesuai kurva di tiap daerah sempit itu. Kasus sederhana yaitu, bila daerah integrasi dibagi sama panjang dan untuk tiap daerah digunakan rumus quadrature yang sama.
Contoh, daerah integrasi [a,b] dibagi dalam N bagian sama panjang. b
a +d
a +2d
b -d
b
a
a
a +d
b-2d
b -d
I = ∫ f(x) dx =
•
∫ f(x) dx +
∫ f(x) dx + ... +
∫ f(x) dx + ∫ f(x) dx
b−a d = N
integrasi komposit menggunakan quadrature trapezoid b
I = ∫ f(x) dx ≅ h[21 (f0 + fN ) + f1 + f2 + ... + fN −1 ] a
h=
•
b−a , fi = f(a + ih), i = 0, ..., N N
integrasi komposit menggunakan quadrature Simpson b
I = ∫ f(x) dx ≅ a
2h 1 [2 (f0 + f2N ) + 2(f1 + f3 + ... + f2N−1 ) + f2 + f4 + ... + f2N−2 ] 3 h=
b−a , fi = f(a + ih), i = 0, ..., 2N 2N
Integrasi komposit trapezoid untuk daerah integrasi [a,b] yang dibagi 8 sama panjang: b
I = ∫ f(x) dx ≅ h[21 (f0 + f8 ) + f1 + f2 + f3 + f4 + f5 + f6 + f7 ] a
f(x)
a
h
b
x
Integrasi komposit yang menggunakan quadrature trapezoid dan Simpson; daerah integrasi [a,b] yang dibagi 3: b
I = ∫ f(x) dx ≅ a
h1 (fa + 2fa +h1 + fc ) + h2 (fc + 4fc+h2 + fb ) 2 3
f(x) Simpson
trapezoid
a
b
c h1
h1
2h2
x
Integrasi Monte Carlo Mungkin saja cara-cara integrasi numerik yang sudah disampaikan sulit atau tidak bisa diterapkan untuk mengevaluasi suatu integral. Pada keadaan ini, integrasi Monte Carlo dapat dipilih. Integrasi Monte Carlo tidak menggunakan interpolasi seperti pada cara-cara integrasi numerik sebelum ini. Integral dianggap sebagai satu persegi panjang, dengan lebar daerah integrasi dan tinggi nilai rata-rata integrand f(x), yang diperoleh melalui statistik dengan memanfaatkan bilangan acak: f(x)
1 n < f(x) >= ∑ f(xi ) n i=1 xi = bilangan acak : a ≤ xi ≤ b
b
1 n I = ∫ f(x)dx ≅ (b - a) ∑ f(xi ) n i=1 a
(b-a) x a
b
Persamaan Differensial Persamaan differensial (PD) yang dimaksud yaitu persamaan differensial biasa, bukan persamaan differensial parsial, untuk orde 1 dan 2. Dua masalah yang akan dibahas yaitu: •
PD dengan syarat awal
•
PD dengan syarat batas
PD dengan Syarat Awal dy = f(x, y) dx
Bentuk umum PD orde 1:
y'=
Diketahui:
y(x0 ) = y0 y
Integrasi:
y(x) = ?
x
Masalah persamaan differensial berubah menjadi masalah persamaan integral.
∫ dy = ∫ f(x, y)dx y0
x0
x
y(x) = y0 + ∫ f(x, y)dx x0
x0 +h
Dicari y(x) pada titik x = x0 + h :
y(x0 + h) = y0 +
∫ f(x, y)dx
x0
Setelah y(x0 + h) didapat, selanjutnya dicari y(x0 + 2h) . Demikian seterusnya.
Metode Euler Menurut metode Euler: f(x,y)
f(x0 , y0 )
f(x,y) dianggap konstan dan dihitung pada x = x0.
x
x0 x0 + h Diperoleh:
y(x0 + h) yg diperoleh
y(x) x0 +h
y(x0 + h) ≅ y0 + f(x0 , y0 )
∫ dx
x0
≅ y0 + hf(x0 , y0 )
y(x0 + h) sebenarnya
y0 x
x0 x0 + h
Metode Euler yang Dimodifikasi f(x,y)
f(x0 + 21 h, y(x0 + 21 h))
Modifikasi dilakukan dalam memilih nilai f(x,y) yang dianggap konstan. Dipilih f(x,y) pada titik x = x0 + 21 h : 1 2
x
1 2
f(x0 + h, y(x0 + h))
x0 x0 + h x0 + 21 h
dengan y(x0 + 21 h) dihitung memakai metode Euler:
y(x0 + 21 h) ≅ y0 + 21 hf(x0 , y0 )
y(x0 + h) yg diperoleh
y(x)
Diperoleh:
y(x0 + h) ≅ y0 + hf(x0 + 21 h, y(x0 + 21 h)) ≅ y0 + hf(x0 + 21 h, y0 + 21 hf(x0 , y0 ))
y0 x
x0 x0 + h x0 + 21 h
y(x0 + h) sebenarnya
Metode Euler yang Lebih Baik (Improved) Kali ini dipakai nilai f(x,y) yang merupakan rata-rata dari dua nilai f(x,y), masing-masing pada titik x0 dan x0 + h : 1 2
[f(x0 , y0 ) + f(x0 + h, y(x0 + h))]
f(x,y)
Ini sama dengan menggunakan quadrature trapezoid untuk mengevaluasi integral: x0 +h
∫ f(x, y)dx ≅ h[f(x , y 1 2
0
0
) + f(x0 + h, y(x0 + h))]
x
x0 x0 + h
x0
dengan y(x0 + h) dihitung memakai metode Euler:
y(x0 + h) ≅ y0 + hf(x0 , y0 ) Diperoleh:
y(x0 + h) ≅ y0 + 21 h[f(x0 , y0 ) + f(x0 + h, y(x0 + h))] ≅ y0 + 21 h[f(x0 , y0 ) + f(x0 + h, y0 + hf(x0 , y0 ))]
PD Orde 2
Bentuk umum PD orde 2:
d2 y y''= 2 = f(x, y, y') dx
Diketahui:
y(x0 ) = y0 , y'(x0 ) = y'0
Definisikan fungsi baru u:
u = y' u0 = y'0
Masalah PD orde 2 berubah menjadi masalah PD orde 1.
y(x) = ?
y'= u(x, y) u'= f(x, y, u)
Contoh penyelesaian dengan metode Euler yang lebih baik (improved):
u'= f(x, y, u)
y'= u(x, y)
u(x0 + h) = u0 + 21 h(f0 + f1 ) f0 = f(x0 , y0 , u0 )
y(x0 + h) = y0 + 21 h(u0 + u1 ) u0 = y'0
f1 = f(x0 + h, y0 + hu0 , u1 )
u1 = u0 + hf0
Alur perhitungan:
y0 , u0
f0
u1
f1 , y(x0 + h), u(x0 + h)
x0 + h → x0 , u(x0 + h) → u0 , y(x0 + h) → y0
PD dengan Syarat Batas Contoh, gelombang yang merambat di sepanjang tali bisa digambarkan dengan PD orde 2. Jika ujung-ujung tali itu diikat sehingga tidak bisa bergerak, maka kita temui kasus PD dengan syarat batas.
terikat
terikat
Bentuk umum PD orde 1 & 2 linear:
Diketahui:
(1) y'= f(x, y) = d(x) − e(x)y (2) y''= g(x, y, y') = a(x) − b(x)y − c(x)y'
x0 ≤ x ≤ xn y(x0 ) = y0
y(x) = ?
y(xn ) = yn Dicari yi = y(xi ) pada titik xi = x0 + ih (i = 1, ..., n − 1) dengan h =
xn − x0 . n
Metode Finite Differences (1) y'+e(x)y = d(x)
yi+1 − yi−1 2h y − 2yi + yi−1 y''i ≅ i+1 h2 y'i ≅
(2) y''+c(x)y'+b(x)y = a(x) (1) y'i +ei yi = di (2) y''i +ci y'i +bi yi = ai
yi+1 − yi−1 + ei yi ≅ di 2h y − 2yi + yi−1 yi+1 − yi−1 (2) i+1 + c + bi yi ≅ ai i 2 h 2h
(1)
Jadi, pada akhirnya ditemui masalah sistem persamaan linear:
(1) − yi−1 + 2eihyi + yi+1 ≅ 2dih ch ch (2) 1 − i yi−1 − 2 − bih2 yi + 1 + i yi+1 ≅ aih2 2 2
(
)
yang dapat diselesaikan menggunakan metode, contoh, Eliminasi Gauss. Namun, sistem persamaan linear juga dapat diselesaikan dengan cara lain yaitu, iterasi. Contoh: iterasi Jacobi dan iterasi Gauss-Siedel.
Iterasi Jacobi sistem persamaan linear:
n
∑a x ij
j
= bi (i = 1, ..., n)
j=1
n 1 solusi: xi = bi − ∑ aijxj aii j≠i
Pencarian solusi dimulai dengan nilai awal xi(0) (i = 1, …, n) hasil perkiraan / tebakan. Dengan nilai tebak awal ini diperoleh nilai perkiraan berikut xi(1) melalui: (1) i
x
n 1 (0) = bi − ∑ aijxj (i = 1, ..., n) aii j≠i
Demikian seterusnya berulang-ulang, nilai perkiraan pada langkah ke k diperoleh dari nilai perkiraan pada langkah ke k-1: (k) i
x
n 1 = bi − ∑ aijxj(k-1) (i = 1, ..., n) aii j≠i
Pencarian dihentikan setelah didapat nilai xi yang konvergen yaitu, yang tidak atau sedikit berubah dari nilai yang diperoleh pada langkah sebelumnya:
xi(k-1) 1 − (k) < ε, xi
ε = bilangan kecil
Iterasi Gauss-Siedel
Rumus iterasi Jacobi dapat ditulis:
(k) i
x
1 (k -1) (k -1) = bi − ∑ aijxj − ∑ aijxj aii ji
Jika pada tiap langkah pencarian dilakukan dengan urutan i yang makin besar, (k) (k) maka semua xji (k) sudah diperoleh ketika mencari xi . (k) (k) Karena itu, nilai xji itu bisa langsung dipakai untuk mencari xi(k), sehingga iterasi mencapai nilai konvergen menjadi lebih cepat:
1 (k) (k -1) b − a x − a x ∑ ∑ i ij j ij j aii ji
(i = 1, 2, ..., n)
1 (k -1) (k) = bi − ∑ aijxj − ∑ aijxj aii ji
(i = n, ..., 2, 1)
xi(k) = (k) i
x
Iterasi seperti ini disebut iterasi Gauss-Siedel.
Sebelum ini dikenal metode Eliminasi Gauss dan LU Decomposition untuk mencari solusi sebuah sistem persamaan linear. Pada metode ini terdapat substitusi mundur dan maju. Pada substitusi mundur (maju), nilai xi dihitung dari nilai xj>i ( xji ( xj
Aplikasi Iterasi Jacobi dan Gauss-Siedel pada PD dengan Syarat Batas PD orde 1: (1) y'= d(x) − e(x)y
− yi−1 + 2eihyi + yi+1 ≅ 2dih
PD orde 2: (2) y''= a(x) − b(x)y − c(x)y'
cih ch 1 − yi−1 − (2 − bih2 )yi + 1 + i yi+1 ≅ aih2 2 2
Iterasi Jacobi:
(1) yi(k) ≅
1 2dih + yi(k-1-1) − yi(k+1-1) 2eih
(2) yi(k) ≅
1 cih (k-1) cih (k-1) 2 − a h + (2 − bih2 ) i 1 − 2 yi-1 + 1 + 2 yi+1
(
)
Iterasi Gauss-Siedel (contoh untuk i membesar, i = 1, …, n-1):
(1) yi(k) ≅
1 (k -1) 2dih + yi(k) -1 − yi+1 2eih
(2) yi(k) ≅
1 cih (k) cih (k-1) 2 − a h + yi+1 1 − yi-1 + 1 + i 2 (2 − bih ) 2 2
(
Catatan, sesuai syarat batas:
y0(k) = y0
)
yn(k) = yn
