DIKTAT KULIAH (3 sks) MX 211: Metode Numerik (Revisi Terakhir: Juni 2009 ) 2009 )
Oleh: Didit Budi Nugroho, M.Si.
Program Studi Matematika Fakultas Sains dan Matematika Universitas Kristen Satya Wacana
KATA PENGANTAR Naskah ini ditulis ketika penulis mengajar Metode Numerik (MX 211) di Universitas Kristen Satya Wacana pada Semester 2 tahun 2008-2009 dan juga Trimester 2 tahun 2008-2009 di Soe, Nusa Tenggara Timur. Catatan ini membentuk naskah dasar untuk kuliah "Metode Numerik". Tujuan yang ingin dicapai dari mata kuliah Metode Numerik adalah memahami konsep dasar metode numeris, mengetahui kekurangan dan kelebihan dari setiap metode, dan juga mampu untuk untuk mengaplik mengaplikasik asikann annya. ya. Karena Karena itu, naskah naskah ini difokusk difokuskan an pada pemahaman konsep matematis dasar dan penyelesaian masalah menggunakan metode numerik numerik dengan bantuan program program MatLab. MatLab. Beberapa Beberapa metode disajikan disajikan dengan dengan algo algo-ritma dalam pseudocode sehingga pseudocode sehingga pembaca atau pengguna bisa mengimplementasikan dalam bahasanya sendiri. Di sini, analisis (seperti untuk galat) tidak diberikan secara mendalam. Naskah ini memerlukan masukan dan saran dari pembaca demi perbaikan dan pengembangan secara terus menerus. Harapannya adalah bahwa naskah ini memberikan manfaat yang lebih dalam pengajaran Metode Numerik. Salatiga, Juli 2009 Didit B. Nugroho
i
KATA PENGANTAR Naskah ini ditulis ketika penulis mengajar Metode Numerik (MX 211) di Universitas Kristen Satya Wacana pada Semester 2 tahun 2008-2009 dan juga Trimester 2 tahun 2008-2009 di Soe, Nusa Tenggara Timur. Catatan ini membentuk naskah dasar untuk kuliah "Metode Numerik". Tujuan yang ingin dicapai dari mata kuliah Metode Numerik adalah memahami konsep dasar metode numeris, mengetahui kekurangan dan kelebihan dari setiap metode, dan juga mampu untuk untuk mengaplik mengaplikasik asikann annya. ya. Karena Karena itu, naskah naskah ini difokusk difokuskan an pada pemahaman konsep matematis dasar dan penyelesaian masalah menggunakan metode numerik numerik dengan bantuan program program MatLab. MatLab. Beberapa Beberapa metode disajikan disajikan dengan dengan algo algo-ritma dalam pseudocode sehingga pseudocode sehingga pembaca atau pengguna bisa mengimplementasikan dalam bahasanya sendiri. Di sini, analisis (seperti untuk galat) tidak diberikan secara mendalam. Naskah ini memerlukan masukan dan saran dari pembaca demi perbaikan dan pengembangan secara terus menerus. Harapannya adalah bahwa naskah ini memberikan manfaat yang lebih dalam pengajaran Metode Numerik. Salatiga, Juli 2009 Didit B. Nugroho
i
DAFTAR ISI KATA PENGANTAR
i
DAFTAR ISI
ii
DAFTAR GAMBAR
iv
DAFTAR TABEL
v
1 Pendahuluan
1
2 Galat 2.1 Pengertian Galat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Penghitungan Galat . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3 3 4
3 Pencarian Akar 3.1 Akar-akar Persamaan . . . . 3.2 Metode Bagi Dua . . . . . . 3.3 Metode Posisi Palsu . . . . 3.4 Metode Iterasi Titik Tetap 3.5 Metode Newton-Raphson . 3.6 Metode Garis Potong . . . .
. . . . . .
6 6 7 10 11 14 17
4 Metode Iterasi 4.1 Iterasi Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Iterasi Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Iterasi SOR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20 20 23 25
5 Interpolasi Polinomial 5.1 Interpolasi Linear . . . 5.2 Interpolasi Kuadratik 5.3 Interpolasi Newton . . 5.4 Interpolasi Lagrange .
. . . .
28 28 28 29 32
6 Interpolasi Spline 6.1 Spline Linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Spline Kuadratik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Spline Kubik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
36 36 39 42
. . . .
. . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
ii
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
. . . . . .
. . . .
DAFTAR ISI
iii
7 Regresi Kuadrat Terkecil 7.1 Regresi Linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Regresi Polinomial Deraja ajat Dua . . . . . . . . . . . . . . . . . . . . . . . 7.3 Linearisasi Fungsi Tak Linear . . . . . . . . . . . . . . . . . . . . . . . .
50 50 52 55
8 Diferensiasi Numerik 8.1 Turunan Pertama . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.2 Ekstrapol polasi Richardson . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.3 Turunan Kedua . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57 57 60 61
9 Persamaan 9.1 Metode 9.2 Metode 9.3 Metode 9.4 Metode
. . . .
63 63 66 68 71
10 Integrasi Numerik 10.1 At Aturan Trap esium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.2 Atu Aturan-aturan Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . .
74 74 76
DAFTAR PUSTAKA
80
Diferensial Biasa Euler . . . . . . . . Heun . . . . . . . . Titik Tengah . . . . Runge-Kutta . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
DAFTAR GAMBAR 3.1 Ilustrasi gra…s untuk akar hampiran dalam metode bagi dua. . . . . . . 3.2 Ilustrasi gra…s untuk akar hampiran dalam metode posisi palsu. . . . . . 3.3 Ilustrasi gra…s untuk akar hampiran dalam metode iterasi titik tetap y = x , y = g (x). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Ilustrasi gra…s untuk akar hampiran dalam metode Newton-Raphson. . 3.5 Ilustrasi gra…s untuk akar hampiran dalam metode Secant. . . . . . . .
7 10
5.1 Kurva polinomial Newton untuk Contoh 5.3. . . . . . . . . . . . . . . .
33
6.1 6.2 6.3 6.4
38 41 42 48
Spline kuadratik untuk Contoh 6.2. . . . . . Spline kuadratik untuk Contoh 6.3. . . . . . Pendekatan dengan polinomial spline kubik. Spline kuadratik untuk Contoh 6.4. . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
7.1 Kurva linear dengan metode kuadrat terkecil untuk titik-titik data dalam Contoh 7.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Kurva polinomial derajat dua dengan metode kuadrat terkecil untuk titik-titik data dalam Contoh 7.2. . . . . . . . . . . . . . . . . . . . . . .
11 14 17
53 55
9.1 Ilustrasi dari penurunan metode Euler. . . . . . . . . . . . . . . . . . . . 9.2 Ilustrasi dari penurunan metode Heun. . . . . . . . . . . . . . . . . . . . 9.3 Ilustrasi penurunan metode titik tengah. . . . . . . . . . . . . . . . . . .
64 66 68
10.1 Ilustrasi dari aturan trapesium. . . . . . . . . . . . . . . . . . . . . . . .
75
iv
DAFTAR TABEL 5.1
Tabel beda terbagi hingga untuk orde 3. . . . . . . . . . . . . . . . . . .
30
7.1
Linearisasi dari fungsi tak linear dengan transformasi data. . . . . . . .
56
8.1 Tabel ekstrapolasi Richardson untuk beda maju sampai hampiran orde 4. 60 8.2 Tabel ekstrapolasi Richardson untuk beda pusat sampai hampiran orde 8. 61 9.1 Perbandingan hasil dari metode-metode penyelesaian untuk persamaan diferensial y 0 (t) = y (t) t2 + 1, 0 t 2, y (0) = 0:5. . . . . . . . . . .
v
73
Bab 1
Pendahuluan Tujuan Pembelajaran:
Mengetahui apa yang dimaksud dengan metode numerik. Mengetahui kenapa metode numerik perlu dipelajari. Mengetahui langkah-langkah penyelesaian persoalan numerik. Metode numerik merupakan teknik untuk menyelesaikan masalah matematika dengan pengoperasian aritmatika (hitungan). Beberapa alasan mengapa kita harus mempelajari metode numerik: 1. Metode numerik merupakan alat bantu pemecahan masalah matematika yang sangat ampuh. Metode numerik mampu menangani sistem persamaan besar, ketidaklinearan, dan geometri yang rumit yang dalam praktek rekayasa seringkali tidak mungkin dipecahkan secara analitik. 2. Di pasaran banyak tersedia program aplikasi numerik komersil. Penggunaan aplikasi tersebut menjadi lebih berarti bila kita memiliki pengetahuan metode numerik agar kita dapat memahami cara paket tersebut menyelesaikan persoalan. 3. Kita dapat membuat sendiri program komputer tanpa harus membeli paket programnya. Seringkali beberapa persoalan matematika tidak selalu dapat diselesaikan oleh program aplikasi. Sebagai contoh, terdapat program aplikasi tertentu yang tidak dapat dipakai untuk menghitung integrasi lipat dua, atau lipat tiga. Mau tidak mau, kita harus menulis sendiri programnya. Untuk itu, kita harus mempelajari cara pemecahan integral lipat dua atau lebih dengan metode numerik. 4. Metode numerik menyediakan sarana untuk memperkuat kembali pemahaman matematika, karena metode numerik ditemukan dengan cara menyederhanakan matematika yang lebih tinggi menjadi operasi matematika yang mendasar. Langkah-langkah penyelesaian persoalan numerik: 1. Identi…kasi masalah. 2. Memodelkan masalah secara matematis. 3. Identi…kasi metode numerik yang diperlukan untuk menyelesaikan masalah. 1
Bab 1. Pendahuluan
2
4. Implementasi metode dalam komputer. 5. Analisis hasil akhir: implementasi, metode, model dan masalah. Jenis-jenis persoalan matematika yang akan diselesaikan secara numerik dalam naskah ini: 1. Pencarian akar-akar persamaan tak linear. 2. Metode iteratif untuk penyelesaian sistem persamaan linear 3. Interpolasi linear, kuadrat, Newton, dan spline. 4. Regresi kuadrat terkecil. 5. Diferensiasi numerik. 6. Persamaan diferensial biasa. 7. Integrasi numerik.
Bab 2
Galat Tujuan Pembelajaran:
Mengetahui de…nisi dan jenis-jenis galat. Mengetahui bagaimana menghitung galat. 2.1
Pengertian Galat
Penyelesaian secara numerik dari suatu persamaan matematik hanya memberikan nilai perkiraan (hampiran) yang mendekati nilai eksak (yang sebenarnya) dari penyelesaian analitis. Berarti dalam penyelesaian numerik tersebut terdapat galat terhadap nilai eksak. Ada tiga macam galat: 1. Galat bawaan, terjadi karena kekeliruan dalam menyalin data, salah membaca skala, atau karena kurangnya pengertian mengenai hukum-hukum …sik dari data yang diukur. 2. Galat pembulatan (round-o¤ error ), terjadi karena tidak diperhitungkannya beberapa angka terakhir dari suatu bilangan. Sebagai contoh, 3.1415926 dapat dibulatkan menjadi 3.14. 3. Galat pemotongan (truncation error ), terjadi karena tidak dilakukannya hitungan sesuai dengan prosedur matematik yang benar. Sebagai contoh, turunan pertama dari V (t) terhadap t dihitung dengan prosedur dV V V (ti+1 ) = = t dt ti+1
V (t ) : t i
i
Contoh lain yaitu pengambilan beberapa suku awal dari deret Taylor: f (xi+1) = f ( xi ) + f 0 (xi ) (xi+1
+ dengan Rn =
f 000 (xi )
3!
(xi+1
f (n+1) (c) (xi+1 (n + 1)!
00
x ) + f 2!(x ) (x +1 x )2 i
i
i
f (n) (xi ) (xi+1 xi ) + ::: + n! 3
x ) +1 dimana c 2 [x ; x +1]. i
n
i
3
i
i
x ) i
n
+ Rn
Bab 2. Galat
4
Dari deret Taylor di atas, dipunyai hampiran orde-0: f ( xi+1 )
f (x ) ; i
hampiran orde-1: f (xi+1 )
f (x ) + f 0 (x ) (x +1 x ) ; i
i
i
i
hampiran orde-2: 00
f ( xi+1 )
f (x ) + f 0 (x ) (x +1 x ) + f 2!(x ) (x +1 x )2 : i
i
i
i
i
i
i
Contoh 2.1 Diberikan fungsi f ( x) =
0:1x4 0:15x3 0:5x2 0:25x + 1:2:
Turunan pertama dan kedua dari f (x) berturut-turut yaitu f 0 (x) = f 00 (x) =
0:4x3 0:45x2 x 0:25; 1:2x2 0:9x 1;
Dimulai dari x = 0, diperoleh hampiran orde-1 untuk f (1): f (1)
f (0) + f 0 (0)(1 0) = 1:2 0:25 = 0:95;
dan hampiran orde-2 untuk f (1): f (1)
2.2
0
f (0) + f (0)(1 0) +
f 00 (0)
2!
(1
0)2 = 1 :2 0:25 0:5 = 0:9:
Penghitungan Galat
Untuk galat pembulatan dan pemotongan, hubungan antara hasil yang eksak dengan hampirannya dapat dirumuskan oleh nilai eksak = hampiran + galat. Dengan menyusun kembali persamaan di atas, diperoleh E s = galat = nilai eksak - hampiran
dimana subskrip s menunjukkan bahwa galat adalah galat sejati. Kelemahan dari de…nisi di atas adalah bahwa tingkat besaran dari nilai yang diperiksa sama sekali tidak diperhatikan. Sebagai contoh, galat satu sentimeter jauh lebih berarti jika yang diukur adalah paku ketimbang jembatan. Salah satu cara untuk memperhitungkan besarnya besaran yang sedang dievaluasi adalah dengan menormalkan galat terhadap nilai eksak, yaitu galat relatif =
nilai eksak - hampiran . nilai eksak
Bab 2. Galat
5
Galat relatif dapat juga dikalikan dengan 100% agar dapat dinyatakan sebagai s = persen galat relatif =
nilai eksak - hampiran nilai eksak
100%.
Dicatat bahwa untuk metode numerik, nilai eksak hanya akan diketahui jika fungsi yang ditangani dapat diselesaikan secara eksak. Jika tidak demikian, maka alternatifnya adalah menormalkan galat dengan menggunakan hampiran terbaik yang tersedia dari nilai eksak, yaitu terhadap hampiran itu sendiri, seperti yang dirumuskan oleh h
= =
galat hampiran 100% hampiran hampiran sekarang hampiran sebelumnya hampiran sekarang
100%
dengan subskrip h menunjukkan bahwa galat dinormalkan terhadap nilai hampiran.
Bab 3
Pencarian Akar Tujuan Pembelajaran:
Mengetahui metode-metode pencarian akar dari persamaan tak linear. Mengetahui keunggulan dan kelemahan dari setiap metode. Mengaplikasikan metode-metode pencarian akar untuk persamaan yang sama. 3.1
Akar-akar Persamaan
Akar atau pembuat nol dari suatu fungsi adalah nilai-nilai dari variabel bebas yang membuat fungsi bernilai nol. Sebagai contoh, penyelesaian analitik untuk fungsi kuadratik f (x) = ax 2 + bx + c = 0 diberikan oleh
p b b2 4ac x = : 2a
Contoh lain, kita tidak mungkin mendapatkan suatu penyelesaian analitik untuk fungsi f (x) = e x x = 0.
Secara umum, terdapat dua kelompok metode untuk pencarian akar dari persamaan tak linear: 1. Metode Pengurung. Sesuai dengan namanya, tebakan akar dalam metode ini selalu berada "dalam kurung" atau berada pada kedua sisi dari nilai akar. Karena itu, di sini diperlukan dua tebakan awal untuk akar. Metode ini mempunyai suatu keunggulan yaitu konvergen (makin lama makin mendekati nilai sebenarnya), dan mempunyai kelemahan yaitu konvergensinya relatif lambat. Contoh dari metode pengurung yaitu metode bagi dua (bisection ) dan metode posisi palsu ( false position ). 2. Metode Terbuka. Dalam metode ini, pencarian dimulai dari suatu nilai tunggal variabel bebas, atau dua nilai yang tidak perlu mengurung akar. Metode ini mempunyai suatu kelemahan yaitu tidak selalu konvergen, tetapi mempunyai keunggulam yaitu jika konvergen maka konvergensinya lebih cepat daripada metode pengurung. Contoh dari metode terbuka yaitu metode iterasi titik tetap ( …xed-point iteration ), metode Newton-Raphson, dan metode garis potong ( secant ).
6
Bab 3. Pencarian Akar
7
Penggunaan metode pengurung didasarkan pada teorema berikut ini. Teorema 3.1 Diberikan f : [a; b] Jika f ( a) f ( b) < 0, maka terdapat c
! R adalah kontinu, dimana a; b 2 R dan a < b. 2 (a; b) sedemikian sehingga f (c) = 0.
3.2
Metode Bagi Dua
Metode ini mengasumsikan bahwa fungsi f (x) adalah kontinu pada interval [a1 ; b1 ], serta f (a1 ) dan f ( b1 ) mempunyai tanda berlawanan, artinya f ( a1 ) f ( b1 ) < 0. Karena itu terdapat minimal satu akar pada interval [ a1 ; b1 ].
Idenya adalah interval selalu dibagi dua sama lebar. Jika fungsi berubah tanda sepan jang suatu subinterval, maka letak akarnya kemudian ditentukan ada di tengah-tengah subinterval. Proses ini diulangi untuk memperoleh hampiran yang diperhalus. (Lihat Gambar 3.1.) y
a1
x2
x3
x1
b1
x
y = f ( x)
Gambar 3.1: Ilustrasi gra…s untuk akar hampiran dalam metode bagi dua.
Dicatat bahwa terdapat beberapa kriteria penghentian pencarian akar jika diberikan suatu toleransi keakuratan , yaitu s < , h < , atau nmaks = N
2 N:
Proses untuk metode bagi dua diberikan seperti dalam Algoritma 1.
Contoh 3.2 Selesaikan persamaan x2 metode bagi dua sampai lima iterasi.
3
= 0 dalam interval [1; 2] menggunakan
Penyelesaian. Proses metode bagi dua adalah seperti berikut ini. Iterasi 1: a1
= 1 =
b1
= 2;
x1
=
)
a1 + b1
2
f ( a1 ) =
=
2;
1+2 = 1:5 = 2
)
f ( x1 ) =
0:75:
Bab 3. Pencarian Akar
8
Algorithm 1 Algoritma Metode Bagi Dua Masukan: fungsi kontinu: f ( x) interval yang mengurung akar: [ a1 ; b1 ] maksimum iterasi: N N toleransi keakuratan: "; misalnya " = 105 Penghitungan Inti: Ketika n N dan h , an + bn Hitung: xn = . 2 Tentukan subinterval mana yang akan mengurung akar: a) Jika f ( an ) f (xn ) < 0, maka an+1 = a n , bn+1 = x n . b) Jika f ( an ) f ( xn ) > 0, maka an+1 = x n , bn+1 = b n . c) Jika f (an ) f ( xn ) = 0, maka diperoleh akar sama dengan xn . Berhenti.
2
Hitung: h =
xn
x 1 100%, n 1 x n
n
Hasil akhir: akar xn sedemikian sehingga f (xn )
0:
Iterasi 2: Diamati bahwa f ( a1 ) f ( x1 ) > 0, maka
a2
= x1 = 1:5 =
b2
= b1 = 2; 1:5 + 2 a2 + b2 = = = 1 :75 = f ( x2 ) = 0 :0625 ; 2 2 1:75 1:5 x2 x1 = 100% = 100% = 14 :29%: x2 1:75
x2 2
)
f ( a2 ) =
0:75;
)
Iterasi 3: Diamati bahwa f ( a2 ) f ( x2 ) < 0, maka a3 b3 x3 3
=)
= a2 = 1:5 f ( a3 ) = 0:75; = x2 = 1:75; 1:5 + 1 :75 a3 + b3 = = = 1 :625 = 2 2 1:625 1:75 x3 x2 = 100% = 1:625 x3
)
Iterasi 4: Diamati bahwa f ( a3 ) f ( x3 ) > 0, maka a4 b4 x4 4
=)
f ( x3 ) =
0:3594;
100% = 7 :69%:
= x3 = 1:625 f ( a4 ) = 0:3594 ; = b3 = 1:75; 1:625 + 1 :75 a4 + b4 = = = 1 :6875 = f ( x4 ) = 0:1523; 2 2 1:6875 1:625 x4 x2 = 100% = 100% = 3:7%: 1:6875 x4
)
Bab 3. Pencarian Akar
9
Iterasi 5: Diamati bahwa f ( a4 ) f ( x4 ) > 0, maka
x4 = 1:6875 =)
a5
=
b5
= b4 = 1:75; 1:6875 + 1 :75 a5 + b5 = = = 1 :7187 = f (x5 ) = 0:0459 ; 2 2 x5 x4 1:7187 1:6875 = 100% = 100% = 1:82%: 1:7187 x5
x5 5
f ( a5 ) =
0:1523;
)
Jadi pada iterasi ke-5 diperoleh akar hampiran x = 1:7187.
H
Dalam MatLab, Algoritma 1 untuk metode bagi dua diimplementasikan dalam fungsi bagidua() berikut ini. function [x,galat] = BagiDua(f,X,N,tol) % BagiDua Menyelesaikan persamaan f(x) = 0 menggunakan metode bagi dua. % % Input: f = fungsi dari x, gunakan fungsi inline(’ekspresi’,’x’) % X = [a b] = vektor titik-titik ujung interval dengan a < b % N = maksimum iterasi % tol = toleransi keakuratan % % Output: x = akar hampiran yang memenuhi kriteria % galat = persen galat relatif % ---PENGHITUNGAN INTI: if nargin < 4, tol = 1e-3; end %dinamakan BagiDua(f,X,N) if nargin < 3, N = 100; end %dinamakan BagiDua(f,X) a = X(1); % batas kiri interval b = X(2); % batas kanan interval x = (a+b)/2; % hampiran awal n = 1; galat = 1; % supaya iterasi pertama dapat dikerjakan while ( n <= N & galat > tol ) % kriteria penghentian if f(a)*f(x)<0 b = x; elseif f(a)*f(x)>0 a = x; else break end xnew = (a+b)/2; % titik tengah interval galat = abs((xnew - x)/xnew)*100; % galat relatif x = xnew; n = n+1; end
Dimisalkan x s adalah akar sejati dari f ( x) = 0. Dicatat bahwa kita dapat menurunkan batas galat untuk metode bagi dua seperti berikut ini.
j xs
xn+1
j j
bn+1
=
j
n
j 1 2
an+1
1 = bn 2
b1
a1j
an =
j j
2
j 1 2
bn1
a 1j = ::: n
Bab 3. Pencarian Akar
10
Ini berarti bahwa telah ditunjukkan xs xn+1 k xs xn , dengan k = 21 . Lebih lanjut, ini dinamakan konvergensi linear dan k dinamakan konstanta galat asimtotik.
j
j j j
Dalam contoh di atas, banyaknya langkah yang diperlukan untuk menjamin bahwa galat kurang dari 10 3 dihitung seperti berikut: n
jj 1 2
3.3
2
1
10
=)
3
=
n ln(2)
1 n 103 = 2n 2 3ln(10) = n 10:
)
) )
103
Metode Posisi Palsu
Metode posisi palsu adalah metode pencarian akar persamaan dengan memanfaatkan kemiringan dan selisih tinggi dari dua titik batas interval yang mengurung akar. Metode ini merupakan salah satu alternatif untuk mempercepat konvergensi. y
x3 x2
x1
a1
b1
x
y = f ( x)
Gambar 3.2: Ilustrasi gra…s untuk akar hampiran dalam metode posisi palsu.
Idenya adalah menghitung akar (yang merupakan titik ujung interval baru) yang merupakan absis untuk titik potong antara sumbu x dengan garis lurus yang melalui kedua titik yang absisnya adalah titik-titik ujung interval lama. (Lihat Gambar 3.2) Rumus untuk mencari akar adalah sebagai berikut. Diasumsikan bahwa fungsi f (x) adalah kontinu pada interval [ an ; bn ], dan f ( an ) f ( bn ) < 0. Garis yang melalui titik ( an ; f ( an )) dan ( bn ; f ( bn )) mempunyai persamaan
y
f (a ) = f (bb ) f a (a n
n)
n
n
n
(x
n) :
a
Garis memotong sumbu x jika y = 0, sehingga diperoleh titik absis sebagai hampiran akar yaitu xn = a n
f (bb ) af (a n
n
n
n)
f ( an ) :
(3.1)
Proses untuk metode posisi palsu adalah seperti metode bagi dua tetapi penghitungan xn menggunakan rumus (3.1).
Bab 3. Pencarian Akar
11
Contoh 3.3 Selesaikan persamaan x2 metode posisi palsu sampai lima iterasi.
3
= 0 dalam interval [1; 2] menggunakan
Penyelesaian. Berikut ini adalah tabel penyelesaian sampai lima iterasi. n
an
bn
xn
f ( xn )
1 2 3 4 5 6
1 1:66667 1:72727 1:73170 1:73202 1:73204
2 2 2 2 2 2
1:6667 1 :7273 1 :7317 1 :7320 1 :7321 1 :7321
0:22222222 0:01652892 0:00118976 0:00008543 0:00000613 0:00000044
n =
xn xn1 xn
%
3 :50877192 0 :25608194 0 :01840773 0 :00132172 0 :00009489
Jadi pada iterasi ke-6 diperoleh akar hampiran x = 1:7320 dengan f ( x) = H
3.4
0:000006.
Metode Iterasi Titik Tetap
Metode iterasi titik tetap adalah metode yang memisahkan x sedemikian sehingga f (x) = 0 ekivalen dengan x = g (x). Selanjutnya p adalah suatu akar dari f ( x) jika hanya jika p adalah suatu titik tetap dari g (x). Kita mencoba untuk menyelesaikan x = g (x) dengan menghitung xn = g (xn1 ) ; n = 1; 2;:::
dengan menggunakan tebakan awal x 0 . Ilustrasi gra…s untuk penyelesaian x n diberikan oleh Gambar 3.3 y y = x
y = g ( x)
x0 x2 x4
x3 x1
x
Gambar 3.3: Ilustrasi gra…s untuk akar hampiran dalam metode iterasi titik tetap y = x , y = g (x).
Proses untuk metode bagi dua diberikan seperti dalam Algoritma 2.
Contoh 3.4 Gunakan metode iterasi titik tetap sampai lima iterasi untuk menyelesaikan x2 3 = 0 dengan tebakan awal x0 = 1:5.
Bab 3. Pencarian Akar
12
Algorithm 2 Algoritma Iterasi Titik Tetap Masukan: fungsi kontinu: f ( x), g (x) tebakan awal: x0 maksimum iterasi: N N toleransi keakuratan: "; misalnya " = 105 Penghitungan Inti: Ketika n N dan h , Hitung: xn = g (xn1 ).
2
Hitung: h =
x 1 100% x
xn
n
n
Hasil akhir: akar xn sedemikian sehingga f (xn )
0:
Penyelesaian. Persamaan dapat diubah ke beberapa bentuk: Kasus 1:
x =
Kasus 2:
3
)
x
x = x
Kasus 3:
= g 1 (x) =
x = x
x2
xn =
3 xn1
) x = x 1 x2 1 3 x2 1 3 =) x = x 1 2
3 = g 2 (x) =
x2
2 3 = g3 (x)
n
n
n
n
n
n
Penyelesaian untuk Kasus 3: Iterasi 1: x1 1
= =
Iterasi 2: x2 2
1:875
3 1 :8752 3 = 1 :875 = 1 :6172 ; x1 2 2 1:6172 1:875 100% = 15:94%:
=
=
1:6172
2 2 x2 2 3 = 1 :6172 1 :61722 3 = 1 :8095; 1:8095 1:6172 100% = 10:63%:
= x2
3
=
x4
= x3
Iterasi 4:
4
3 1 :52 3 x0 = 1 :5 = 1 :875; 2 2 1:875 1:5 100% = 20%:
x21
Iterasi 3: x3
x20
=
1:8095
2 2 x3 2 3 = 1 :8095 1 :80952 3 = 1 :6723; 1:6723 1:8095 100% = 8:20%:
1:6723
Bab 3. Pencarian Akar
13
Iterasi 5: 2 2 x4 2 3 = 1 :6723 1 :67232 3 = 1 :7740; 1:7740 1:6723 100% = 5:73%:
= x4
x5
=
4
1:7740
Berikut ini adalah penyelesaian sampai lima iterasi untuk ketiga kasus. n
xn
Kasus 1 1:5 2 1:5 2 1:5 2
0 1 2 3 4 5
Kasus 2 1:5000 2:2500 0:1875 3:1523 3:7849 15:1106
Kasus 3 1 :5 1 :875 1 :6172 1 :8095 1 :6723 1 :7740
Terlihat bahwa Kasus 1 dan Kasus 2 adalah divergen, tetapi Kasus 3 adalah konvergen. H
Berikut ini diberikan kriteria kekonvergenan dari metode iterasi titik tetap. Teorema 3.5 Diambil k = max g0 (x) . Metode iterasi titik tetap adalah konvergen axb
jika hanya jika k < 1.
j
j
Bukti. Dimisalkan xs adalah akar sejati dari f ( x) = 0, maka
jx x j = jg (x ) g (x 1)j = s
n
s
n
g0 ( ) (xs
x 1) k jx x 1j n
berdasarkan Teorema Nilai Rata-rata. Karena itu
s
n
jx x j k jx x 1j k2 jx x 2j k jx x0j : s
n
s
n
s
n
n
s
Kita mencatat: 1. Telah ditunjukkan bahwa xs xn k xs titik tetap adalah konvergen secara linear.
j j j x 1j, yang berarti metode iterasi n
2. Karena k g 0 (xs ) , kita akan memilih fungsi iterasi g (xs ) sedemikian sehingga g0 (xs ) adalah kecil.
j
j
j
j
Dari contoh sebelumnya dipunyai xs = Kasus 1: Kasus 2: Kasus 3:
p 3 = 1:73205, maka
x32 =) g0 (x ) = 1 : divergen; g 0 (x) = 1 2x =) g0 (x ) = 2:4641 : divergen; g 0 (x) = 1 x =) g 0 (x ) = 0:73025 : konvergen : g 0 (x) =
s
s
s
Bab 3. Pencarian Akar
14
Dalam MatLab, Algoritma 2 untuk metode iterasi titik tetap diimplementasikan dalam fungsi TitikTetap() berikut ini. function [x,galat] = TitikTetap(g,x0,N,tol) % TitikTetap Menyelesaikan persamaan x = g(x), ekivalen dengan f(x)=0, % menggunakan metode iterasi titik tetap. % % Input: g = fungsi dari x, gunakan fungsi inline(’ekspresi’,’x’) % x0 = tebakan awal % N = maksimum iterasi % tol = toleransi keakuratan % % Output: x = akar hampiran yang memenuhi kriteria % galat = persen galat relatif % ---PENGHITUNGAN INTI: if nargin < 4, tol = 1e-3; end if nargin < 3, N = 100; end n = 1; galat = 1; while ( n <= N & galat > tol ) x = g(x0); galat = abs((x - x0)/x)*100; x0 = x; n = n+1; end
3.5
Metode Newton-Raphson
Metode Newton-Raphson adalah metode pendekatan yang menggunakan satu titik awal dan mendekatinya dengan memperhatikan kemiringan kurva pada titik tersebut. Pen jelasan gra…s mengenai metode ini adalah seperti dalam Gambar 3.4. y
x0
x1
x
x2
y = f ( x)
Gambar 3.4: Ilustrasi gra…s untuk akar hampiran dalam metode Newton-Raphson.
Diasumsikan bahwa fungsi f (x) adalah kontinu. Idenya adalah menghitung akar yang merupakan titik potong antara sumbu x dengan garis singgung pada kurva di titik (xn1 ; f ( xn1 )). Kemiringan kurva di titik tersebut adalah f 0 (xn1 ), sehingga garis singgung mempunyai persamaan y
f (x 1) = f 0 (x 1) (x x 1) : n
n
n
Bab 3. Pencarian Akar
15
Karena itu diperoleh akar hampiran dengan mengambil y = 0, yaitu xn = x n1
f f 0 ((xx 11)) : n
(3.2)
n
Algorithm 3 Algoritma Metode Newton-Raphson Masukan: fungsi kontinu: f ( x), f 0 (x) tebakan awal: x0 maksimum iterasi: N N toleransi keakuratan: "; misalnya " = 105 Penghitungan Inti: Ketika n N dan h , f (xn1 ) Hitung: xn = x n1 . f 0 (xn1 )
2
Hitung: h =
xn
x 1 100% x
n
n
Hasil akhir: akar xn sedemikian sehingga f (xn )
0:
Contoh 3.6 Gunakan metode Newton-Raphson untuk menyelesaikan x2 ngan tebakan awal x0 = 1:5.
3 = 0 de-
Penyelesaian. Karena f ( x) = x 2 3, f 0 (x) = 2x dan dibentuk rumus pencarian akar:
xn = x n1
x2n1
3 ; n = 1; 2;:::
2xn1
Iterasi 1: x1 1
=
Iterasi 2: x2 2
=
3
= =
1:75
3 1 :752 3 : 2x1 = 1 75 2 1:75 = 1 :73214; 1:73214 1:75 100% = 1:0309278% : 1:73214
= x1
Iterasi 3: x3
2 2 x02x0 3 = 1 :5 12:5 1:53 = 1 :75; 1:75 1:5 100% = 14:2857142% :
= x0
x21
x22
3 1 :73214 2 3 = 1 :73214 = 1 :73205; x2 2x2 2 1:73214 1:73205 1:73214 100% = 0:0053143% : 1:73205
Bab 3. Pencarian Akar
16
Iterasi 4: x4 4
2 23 x32x3 3 = 1 :73205 12:73205 1:73205 = 1 :73205; 1:73205 1:73205 100% = 0%:
= x3 =
1:73205
Jadi pada iterasi ke-4 diperoleh akar hampiran x = 1:73205.
H
Metode Newton-Raphson merupakan suatu contoh dari metode iterasi titik tetap, x n = g (xn1 ), dimana fungsi iterasinya adalah g (x) = x
f f 0 ((xx)) :
Dari sini diperoleh bahwa 0
g (xs ) = 1
[ f 0 (xs )]2 f ( xs ) f 00 (xs ) =0 [f 0 (xs )]2
dengan mengingat bahwa f (xs ) = 0 dan f 0 (xs ) = 0. Ini mengakibatkan bahwa metode Newton adalah konvergen lebih cepat daripada secara linear; pada kenyataannya dipunyai
6
jx x j C jx x 1j2 s
n
s
n
yaitu konvergen kuadratik. Dalam MatLab, Algoritma 3 untuk metode Newton-Raphson diimplementasikan dalam fungsi NewRap() berikut ini. function [x,galat] = NewRap(f,f1,x0,N,tol) % NewtonRaphson Menyelesaikan persamaan f(x) = 0 menggunakan metode Newton-Raphson. % % Input: f = fungsi dari x, gunakan fungsi inline(’ekspresi’,’x’) % f1 = turunan pertama dari f(x), cari dengan fungsi \QTR{bf}{diff} % x0 = tebakan awal % N = maksimum iterasi % tol = toleransi keakuratan % % Output: x = akar hampiran yang memenuhi kriteria % galat = persen galat relatif % ---PENGHITUNGAN INTI: if nargin < 4, tol = 1e-3; end if nargin < 3, N = 100; end n = 1; galat = 1; while ( n <= N & galat > tol ) x = x0-f(x0)/f1(x0); % persamaan (3.2) galat = abs((x - x0)/x)*100; x0 = x; n = n+1; end
Bab 3. Pencarian Akar
3.6
17
Metode Garis Potong
Masalah yang ada dalam metode Newton adalah terkadang sulit untuk mendapatkan turunan pertama f 0 (x). Alternatifnya adalah turunan f 0 (xn1 ), kemiringan garis di (xn1 ; f ( xn1 )), dihampiri oleh kemiringan garis potong yang melalui ( xn2 ; f ( xn2 )) dan ( xn1 ; f ( xn1 )): f (xn2 ) f ( xn1 ) f 0 (xn1 ) :
x 2 x 1
n
Jadi iterasi xn yaitu xn = x n1
n
f (f x (x1)2()xf 2 (x x1)1) : n
n
n
n
(3.3)
n
Penjelasan gra…s mengenai metode ini adalah seperti dalam Gambar 3.5. y
x4 x0
x2
x1
x-1
x3
x
y = f ( x)
Gambar 3.5: Ilustrasi gra…s untuk akar hampiran dalam metode Secant.
Dicatat: 1. Persamaan di atas memerlukan dua tebakan awal x tetapi tidak memperhatikan perubahan tanda dari f (x). 2. Kemudian dapat ditunjukkan bahwa
jx x j C jx x 1j1 6 ; s
n
s
n
:
sehingga metode garis potong lebih cepat daripada metode iterasi titik tetap, tetapi lebih lambat daripada metode Newton-Raphson. Proses untuk metode garis potong diberikan seperti dalam Algoritma 4.
Contoh 3.7 Gunakan metode garis potong untuk menyelesaikan x2 tebakan awal x1 = 1 dan x0 = 2.
3 = 0 dengan
Bab 3. Pencarian Akar
18
Algorithm 4 Algoritma Garis Potong Masukan: fungsi kontinu: f ( x) dua tebakan awal: x1 , x0 maksimum iterasi: N N toleransi keakuratan: "; misalnya " = 105 Penghitungan Inti: Ketika n N dan h , f ( xn1 ) (xn2 xn1 ) Hitung: xn = x n1 . f ( xn2 ) f (xn1 )
2
Hitung: h =
xn
x 1 100% x
n
n
Hasil akhir: akar xn sedemikian sehingga f (xn )
0:
Penyelesaian. Dituliskan rumus untuk metode garis potong : xn
x2n1
3 (x 2 x 1) x 1 x2 2 3 x2 1 3 x2 1 3 (x 2 x 1 ) x 1 ; n = 1; 2;::: x2 2 x2 1
=
n
=
n
n
n
n
n
n
n
n
n
n
Iterasi 1:
Iterasi 2: x2
2
Iterasi 3: x3
3
4
x1
= x0
1
=
3 (x1
x21
1:66667 2 1:66667
x21
= x1
x20
3 (x0
x20
x22
x1 )
x21
3 (x1
x21
= x3
3 (x2
x22
=2
22
3 (1 2) = 1:66667 ; 12 22
x23
= 1:732051 ; 1:732051 1:73214 = 1:732051
= 1:66667
1:666672 22
3 (2 1:66667) 1:666672
100% = 3:50877% :
x2 )
x22
= 1:73214; 1:73214 1:72727 = 1:73214
x23
x0 )
100% = 20% :
= 1:72727 ; 1:72727 1:66667 = 1:72727
= x2
Iterasi 4: x4
x20
= 1 :72727
1:72727 2 3 (1:66667 1:72727) 1:666672 1:727272
100% = 0:28116% :
x3 )
= 1 :73214
1:73214 2 3 (1:72727 1:73214) 1:727272 1:732142
100% = 0:00532% :
Bab 3. Pencarian Akar Jadi pada iterasi ke-4 ke- 4 diperoleh akar hampiran x = 1:732051. 732051 .
19 H
Dalam MatLab, Algoritma Algoritma 4 untuk untuk metode garis potong diimplementa diimplementasik sikan an dalam dalam fungsi Secant() berikut ini. function function [x,galat] [x,galat] = Secant(f,x Secant(f,x_1,x0 _1,x0,N,t ,N,tol) ol) % Secant Secant Menyel Menyelesa esaika ikan n persam persamaan aan f(x) = 0 menggu menggunak nakan an metode metode % garis potong (Secant). % % Input: Input: f = fungsi fungsi dari x, gunaka gunakan n fungsi fungsi inline(’ inline(’eks ekspre presi’ si’,’x ,’x’) ’) % x_1 = tebakan awal pertama % x0 = tebakan awal kedua % N = maksimum iterasi % tol = toleransi keakuratan % % Output Output: : x = akar akar hampir hampiran an yang yang memenu memenuhi hi kriter kriteria ia % galat = persen galat relatif % ---PENGHI ---PENGHITUNGA TUNGAN N INTI: if nargi nargin n < 4, tol = 1e-3 1e-3; ; end end if nargi nargin n < 3, N = 100; 100; end n = 1; galat = 1; while ( n <= N & galat > tol ) x = x0-f(x0)* x0-f(x0)*(x_1(x_1-x0)/ x0)/(f(x_ (f(x_1)-f( 1)-f(x0)) x0)); ; % persamaan persamaan (3.3) (3.3) galat galat = abs((x abs((x - x0)/x) x0)/x)*10 *100; 0; x_1 x_1 = x0; x0; x0 = x; n = n+1; end
Bab 4
Metode Iterasi Tujuan Pembelajaran:
Mengetahui metode-metode penyelesaian sistem persamaan linear secara iterasi. Mengaplikasikan iterasi Jacobi, Gauss-Seidel, dan SOR. Untuk menyelesaikan suatu sistem persamaan linier Ax = b, dipunyai pilihan antara metode metode langsu langsung ng atau atau metode metode iteras iterasi. i. Contoh Contoh dari metode metode langsu langsung ng yaitu yaitu metode invers, invers, eliminasi eliminasi Gauss, Gauss, dan dekomposis dekomposisii LU . Metode iterasi iterasi mempuny mempunyai ai kelebihan kelebihan dalam hal kemurahan memori dan waktu CPU . Metode ini dimulai dari penentuan nilai awal awal vektor vektor x0 sebagai sebagai suatu penyelesa penyelesaian ian awal awal untuk untuk x. Terdapa erdapatt beberapa beberapa metode iterasi seperti iterasi iterasi Jacobi, Jacobi, iterasi iterasi Gauss-Sei Gauss-Seidel, del, dan iterasi iterasi SOR. Suatu klas besar dari metode iterasi dapat dide…nisikan sebagai berikut. Sistem Ax = b dituliskan menjadi Qx = (Q A) x + b untuk suatu matriks Q yang berorde sama dengan A. Selanjutnya dide…nisikan suatu skema iterasi ke- k:
Qx(k) = (Q
A) x( 1) + b; k
k = 1; 2; :::; :::;
(4.1)
untuk x(0) adalah adalah suatu tebakan tebakan awal. Diasumsik Diasumsikan an bahwa Q adalah non-singular , sehingga skema iterasi menghasilkan suatu barisan tunggal vektor-vektor x(k) .
4.1 4.1
Iter Iteras asii Jaco Jacobi bi
Pertama kali dicatat bahwa matriks A dapat dituliskan sebagai A = L + D + U , dengan L adalah matriks segitiga bawah, D adalah matriks diagonal, dan U adalah matriks segitiga segitiga atas. Iterasi Iterasi Jacobi Jacobi memilih memilih Q = D . Karena itu (k) x
= D1 (D = D
1
n b
A) x(k1) + D1 b
(A
D) x(k1)
o
;
k = 1; 2; :::; :::;
dengan asumsi bahwa masukan-masukan diagonal dari A tidak sama dengan nol (jika tidak maka dilakukan penukaran baris-baris dan kolom-kolom untuk mendapatkan su(k) atu sistem yang yang ekivalen) ekivalen).. Untuk Untuk langk langkah ke- k, komponen-komponen xi dinyatakan oleh n 1 (k) (k1) (4.2) xi = bi aij x j ; i = 1; 2;:::;n aii
0 X @
1 A
j =1; j 6 =i
20
Bab 4. Metode Iterasi
21
Contoh 4.1 Diketahui sistem persamaan linear: 10x1
x2 + 2x3 x1 + 11 x2 x3 + 3x4 2x1 x2 + 10 x3 x4 3x2 x3 + 8 x4
= 6; = 25; =
11;
= 15:
Sistem persamaan tersebut diubah susunannya menjadi x1
=
x2
=
x3
=
x4
=
1 (6 + x2 2x3 ) ; 10 1 (25 + x1 + x3 3x4 ) ; 11 1 ( 11 2x1 + x2 + x4 ) ; 10 1 (15 3x2 + x3 ) : 8
Kita bisa menyatakan bahwa nilai x1 , x2 , x3, dan x4 yang berada di ruas kiri sebagai Sementar ara a itu, nilai nilai x1 , x2 , x3, dan x4 yang berada di ruas kanan sebagai x(baru) . Sement lama) ( . Secara umum, sistem persamaan tersebut dapat dituliskan menjadi seperti x (k) x1 (k) x2 (k) x3 (k) x4
= = = =
1 (k1) (k1) 2x3 +6 ; x2 10 1 (k1) (k1) (k1) + x3 3x4 + 25 ; x1 11 1 (k1) (k1) (k1) 2x1 + x2 + x4 11 ; 10 1 (k1) (k1) 3x2 + x3 + 15 : 8
Dimisalkan bahwa nilai-nilai awal x(0) = 0. Pada Pada langkah langkah k = 1, diperoleh nilai-nilai (1) untuk x : 6 (1) 25 (1) 11 (1) 15 (1) x1 = ; x2 = ; x3 = ;x = : 10 11 10 4 8
Setelah Setelah nilai-nilai nilai-nilai x(1) diper diperoleh, oleh, penghit penghitungan ungan diulangi diulangi kembali kembali dengan nilai k = 2 (2) untuk memperoleh nilai-nilai x : (2)
(2)
(2)
x1 = 1 :0473 ; x2 = 1:7159; x3 =
0:8052; x(2) 4 = 0 :8852 :
Proses diulang lagi untuk nilai-nilai k beriku berikutny tnya. a. Berikut Berikut ini diberi diberikan kan suatu suatu hasil hasil sampai langkah ke-4 ke-4. k (k)
x1
(k)
x2
(k)
x3
(k) x4
0 0 0 0 0
1 0:6000 2:2727 1:1000 1:8852
2 1 :0473 1 :7159 0:8052 0 :8852
3 0:9326 2:0530 1:0493 1:1309
4 1 :0152 1 :9537 0:9681 0 :9738
Bab 4. Metode Iterasi
22
Untuk kriteria penghentian iterasi, kita bisa menggunakan suatu norma dari vektor: l1 =
Sebagai contoh, vektor
kxk1 = maks jx j : 1 i
i n
x = (3;
5) memiliki norma 2; 8; 5) memiliki l1 = kxk1 = maks fj3j ; j2j ; j8j ; j5jg = 8 : 1 i n
Sekarang kita ambil kriteria penghentian dalam metode iterasi: l1 =
(k) x
x
(k1)
1
(k)
= maks xi 1in
(k1)
x
i
:
Berdasarkan kriteria tersebut, metode iterasi Jacobi dapat dinyatakan seperti dalam Algoritma 5, sedangkan implementasi dalam MatLab diberikan oleh fungsi jacobi() . Algorithm 5 Algoritma 5 Algoritma Iterasi Jacobi Masukan: matriks koe…sien sistem: A Rnn vektor konstanta sistem: b Rn vektor penyelesaian awal: x(0) Rn maksimum iterasi: N N toleransi keakuratan: "; misalnya " = 106 Penghitungan Inti: Dibentuk matriks: P = D 1 (D A) dan Q = D 1 b. Ketika k N dan norm ", Hitung: xn = Q P x(n1) norm = maks 1in xi Hasil akhir: vektor x sedemikian sehingga Ax b:
2 2
2 2
2
j j
function function X = jacobi(A, jacobi(A,b,X0 b,X0,N,to ,N,tol) l) % jacobi jacobi Menyel Menyelesa esaika ikan n SPL AX = b menggu menggunak nakan an iterasi iterasi Jacobi. Jacobi. % % Inpu Input: t: A = matr matrik iks s koef koefis isie ien n dari dari sistem sistem % b = vektor konstan dari sistem % X0 = penyelesaian awal % N = maksimum iterasi % tol = toleransi keakuratan % % Output Output: : X = penyel penyelesa esaian ian sistem sistem if nargi nargin n < 5, tol = 1e-6 1e-6; ; end end if nargi nargin n < 4, N = 1000 1000; ; end end if nargin nargin < 3, X0 = zeros( zeros(siz size(b e(b)); )); end n = size(A size(A,1) ,1); ; X = X0; P = zeros( zeros(n,n n,n); ); for i = 1:n for j = 1:n if j ~= i
Bab 4. Metode Iterasi
23
P(i,j) = A(i,j)/A(i,i); % elemen dari matriks inv(D)*(D-A) end end Q(i) = b(i)/A(i,i); % elemen dari matriks inv(D)*b end while k <= N & norma > tol X = Q’ - P*X; % persamaan (4.2) norma = max(abs(X-X0)); X0 = X; end
4.2
Iterasi Gauss-Seidel
Dari persamaan (4.2) dicatat bahwa komponen-komponen dari x(k) diketahui, tetapi tidak digunakan, ketika penghitungan komponen-komponen sisanya. Metode GaussSeidel merupakan suatu modi…kasi dari metode Jacobi, yaitu semua komponen-komponen terakhir yang dihitung dipergunakan. Prosedur dalam metode Gauss-Seidel diperoleh dengan memilih Q = D + L. Karena itu, dari (4.1) didapatkan bentuk matriks (D + L) x(k) = (D + L Dx(k) (k) x
A) x( 1) + b Lx( ) U x( 1) + b D1 b Lx( ) U x( 1)
= =
Jadi, skemanya adalah (k)
xi
=
0 X @ i1
1 aii
bi
k
k
k
k
k
n
(k)
aij x j
j =1
(k1)
X
aij x j
j =i+1
1 A
;
;
k = 1; 2; ::::
i = 1; 2;:::;n
(k)
dengan asumsi bahwa untuk langkah k , komponen-komponen x j , 1 diketahui.
(4.3)
j i 1, sudah
Contoh 4.2 Diperhatikan kembali sistem persamaan linear pada contoh sebelumnya. Sistem persamaan diubah susunannya menjadi (k ) x1
=
(k ) x2
=
(k )
=
x3
(k ) x4
=
1 (k1) x2 10
(k1) 2x3
+6 ;
1 (k) (k1) (k1) 3x4 + 25 ; x1 + x3 11 1 (k) (k) (k1) 2x1 + x2 + x4 11 ; 10 1 (k) (k) 3x2 + x3 + 15 : 8
Pada langkah k = 1, diperoleh nilai-nilai untuk x(1) : (1)
(1)
(1)
x1 = 0 :6000 ; x2 = 2:3272; x3 =
0:9873; x(1) 4 = 0 :8789 :
Setelah nilai-nilai x(1) diperoleh, penghitungan diulangi kembali dengan nilai k = 2
Bab 4. Metode Iterasi
24
untuk memperoleh nilai-nilai x(2): (2)
(2)
(2)
x1 = 1 :0300 ; x2 = 2:0369; x3 =
1:0145; x(2) 4 = 0 :9843 :
Proses diulang lagi untuk nilai-nilai k berikutnya. Berikut ini diberikan suatu hasil sampai langkah ke-4. k (k) x1 (k) x2 (k) x3 (k) x4
0 0 0 0 0
1 0:6000 2:3272 0; 9873 0:8789
2 1 :0300 2 :0369 1:0145 0 :9843
3 1:0066 2:0036 1:0025 0:9984
4 1 :0009 2 :0003 1:0003 0 :9998
Untuk keperluan komputasi, iterasi (4.3) dinyatakan secara khusus untuk i = 1 dan i = n berturut-turut:
0 X @ n
(k)
x1
=
b1
(k1)
a1 j x j
j =2
1 A
; x(nk) =
0 X @
n1
1 ann
bn
(k)
anj x j
j =1
1 A
:
(4.4)
Sekarang kita bisa menyatakan skema untuk iterasi Gauss-Seidel seperti dalam Algoritma 6 dan implementasi dalam MatLab diberikan oleh fungsi gauseid() . Algorithm 6 Algoritma Iterasi Gauss-Seidel Masukan: matriks koe…sien sistem: A Rnn vektor konstanta sistem: b Rn vektor penyelesaian awal: x(0) Rn maksimum iterasi: N N toleransi keakuratan: "; misalnya " = 106 Penghitungan Inti: Ketika 1 k N dan norm ", (k) Hitung: x1 menggunakan (4.4) untuk i = 2 : 1 : n 1 (k) Hitung: xi menggunakan (4.3) (k) xn menggunakan (4.4) Hitung: norm = maks 1in xi Hasil akhir: vektor x sedemikian sehingga Ax b:
2 2
2
2
j j
function X = gauseid(A,b,X0,N,tol) % gauseid Menyelesaikan SPL AX = b menggunakan iterasi Gauss-Seidel. % % Input: A = matriks koefisien dari sistem % b = vektor kolom untuk nilai konstanta dari sistem % X0 = penyelesaian awal % N = maksimum iterasi % tol = toleransi keakuratan % % Output: X = penyelesaian sistem
Bab 4. Metode Iterasi
25
if nargin < 5, tol = 1e-6; end if nargin < 4, N = 1000; end if nargin < 3, X0 = zeros(size(b)); end n = size(A,1); X = X0; while k <= N & norma > tol X(1,:) = (b(1)-A(1,2:n)*X(2:n,:))/A(1,1); % persamaan (4.4) for i = 2:n-1 tmp = b(i,:)-A(i,1:i-1)*X(1:i-1,:)-A(i,i+1:n)*X(i+1:n,:); X(i,:) = tmp/A(i,i); % persamaan (4.3) end X(n,:) = (b(n,:)-A(n,1:n-1)*X(1:n-1,:))/A(n,n); % persamaan (4.4) norma = max(abs(X-X0)); X0 = X; end
4.3
Iterasi SOR
Berikutnya diperhatikan suatu metode untuk mempercepat konvergensi dari metode iterasi. Dipilih 1 Q = D + L !
dengan ! adalah suatu faktor skala, maka iterasi (4.1) menjadi
1 !
1 !
(k)
D + L
x
=
D + L
(k) x
=
Dx(k)
=
(k) x
=
1 !
(k) x
1
!
D + L
1
!
(k1) x
A
1 D
U
1
Lx(k) +
!
=
!D
1
(k1) x
1 D
!D 1 Lx(k) + 1
(k1) x
+b
+b
U
(k1) x
!
!D 1 U
Lx(k)
(k1)
+ Dx
+b
(k1) x
+ !D 1 b
+ U x(k1)
untuk k = 1; 2;:::. Secara jelas, ini diproses dengan cara: (k)
xi
(k1)
= xi
= (1
a!
ii
0X @
(k1)
!) x
i
(k)
a ji xi
(k1)
+ aii xi
i
! + aii
X X
+
(k1)
a ji xi
i>j
0 X @ i1
bi
j =1
n
(k)
aij x j
j =i+1
b
1 A 1 A
(k1)
aij x j
bi
(4.5)
untuk i = 1; 2;:::;n, dan diasumsikan bahwa untuk langkah ke- k, komponen-komponen (k) x j , 1 j i 1, sudah diketahui.
Untuk ! = 1, iterasi (4.5) memberikan metode Gauss-Seidel. Untuk 0 < ! < 1, prosedurnya dinamakan metode under-relaxation dan dapat digunakan untuk memperoleh konvergensi dari beberapa sistem yang tidak konvergen oleh metode Gauss-Seidel. Untuk ! > 1, prosedurnya dinamakan metode overrelaxation , yang digunakan untuk mempercepat konvergensi bagi sistem yang konvergen oleh teknik Gauss-Seidel.
Bab 4. Metode Iterasi
26
Metode-metode tersebut disingkat SOR untuk Successive Overrelaxation dan digunakan untuk penyelesaian sistem linier yang muncul dalam penyelesaian numeris dari persamaan diferensial parsial tertentu. Contoh 4.3 Diketahui sistem persamaan linear: 4x1 + 3 x2 = 24; 3x1 + 4 x2
x3
x2 + 4x3
= 30; =
24:
Persamaan untuk metode relaksasi dengan ! = 1:25: (k) x1
= =
(k)
x2
= =
(k)
x3
= =
1 :25 (k1) (1 + 24 3x2 4 (k1) (k1) 0:25x1 0:9375x2 + 7:5; (k1) 1 :25 (k) (k1) (1 1:25) x2 + 30 3x1 + x3 4 (k) (k1) (k1) 0:9375 x1 0:25x2 + 0:3125x3 + 9:375; (k1) 1 :25 (k) (1 1:25) x3 + 24 + x2 4 (k) (k1) 0:3125x2 0:25x3 7:5:
(k1) 1:25) x1
Tabel berikut ini menampilkan hasil penghitungan sampai langkah ke-4 menggunakan penyelesaian awal x(0)= 0. k (k) x1 (k) x2 (k) x3
0 1 1 1
1 7:5000 2:3438 6:7676
2 3 :4277 3 :4607 4:7266
3 3:3987 3:8465 5:1163
4 3 :0442 3 :9606 4:9832
Skema iterasi (4.5) diimplementasikan dengan fungsi MatLab SOR() berikut ini. function X = SOR(A,b,X0,N,tol,w) % SOR Menyelesaikan SPL Ax = b menggunakan iterasi SOR % % Input: A = matriks koefisien dari sistem % b = vektor kolom untuk nilai konstanta dari sistem % X0 = penyelesaian awal % w = faktor skala % tol = toleransi keakuratan % % Output: X = penyelesaian sistem if nargin < 6, w = 1.25; end if nargin < 5, tol = 1e-6; end if nargin < 4, N = 1000; end if nargin < 3, X0 = zeros(size(b)); end n = size(A,1); X = X0; for k = 1:N X(1) = (1-w)*X(1)+w*(b(1)-A(1,2:n)*X(2:n,:))/A(1,1);
Bab 4. Metode Iterasi for i = 2:n-1 tmp = b(i,:)-A(i,1:i-1)*X(1:i-1,:)-A(i,i+1:n)*X(i+1:n,:); X(i) = (1-w)*X(i)+w*tmp/A(i,i); end X(n) = (1-w)*X(n)+w*(b(n,:)-A(n,1:n-1)*X(1:n-1,:))/A(n,n); galat = max(abs(X-X0)); if galat < tol, break; end X0 = X; end
27
Bab 5
Interpolasi Polinomial Tujuan Pembelajaran:
Mengetahui metode penentuan suatu polinomial yang melewati semua titik data yang diberikan. Mengaplikasikan interpolasi polinomial untuk menaksir nilai antara titik-titik data. Interpolasi menghubungkan titik-titik data diskret dalam suatu cara yang masuk akal sehingga dapat diperoleh taksiran layak dari titik-titik data di antara titik-titik yang diketahui. Dicatat bahwa kurva interpolasi melalui semua titik data.
5.1
Interpolasi Linear
Interpolasi linear merupakan interpolasi paling sederhana dengan mengasumsikan bahwa hubungan titik-titik antara dua titik data adalah linear. Karena itu digunakan pendekatan fungsi linear antara dua titik data, misalnya (x0 ; y0 ) dan ( x1 ; y1 ). Persamaan yang menghubungkan kedua titik tersebut yaitu y = y 0 +
y1 x1
y0 (x x0) : x0
Contoh 5.1 Hitung taksiran y untuk x = 2 dengan menggunakan interpolasi linear untuk data: (1; 0) dan (4; 1:386294). Penyelesaian. Taksiran y untuk x = 2 yaitu y
= y0 +
y1 x1
y0 (2 x0) = 0 + 1 :386294 0 (2 1) 41 x0
= 0:462098 :
5.2
Interpolasi Kuadratik
Interpolasi kuadratik menentukan titik-titik antara tiga titik data dengan menggunakan pendekatan fungsi kuadrat. Dibentuk persamaan kuadrat yang melalui tiga titik data, misalnya ( x0 ; y0 ), ( x1 ; y1 ), dan ( x2 ; y2 ), yaitu y = a 0 + a1 (x
x0) + a2 (x x0) (x x1) : 28
Bab 5. Interpolasi Polinomial
29
Kita akan menentukan b 0 , b1 , b 2 sedemikian sehingga persamaan di atas melalui ketiga titik data yang diberikan. x
=
x
=
x
= =
) =) =) =
)
x0 ; y = y 0 =
) a0 = y0; y1 y0 x1 ; y = y 1 =) y0 + b1 (x1 x0 ) = y 1 =) a1 = ; x1 x0 y1 y0 (x2 x0 ) + a2 (x2 x0 ) (x2 x1 ) = y 2 x2 ; y = y 2 =) y0 + x1 x0 y2 y0 y1 y0 a2 = (x2 x0 ) (x2 x1 ) (x1 x0 ) (x2 x1 ) y2 y1 y1 y0 y1 y0 a2 = + (x2 x0 ) (x2 x1 ) (x2 x0 ) (x2 x1 ) (x1 x0 ) (x2 x1 ) y2 y1 y1 y0 a2 = (x2 x0 ) (x2 x1 ) (x2 x0 ) (x1 x0 ) a2 = : x2 x0 y2 y1 x2 x1
y1 y0 x1 x0
Contoh 5.2 Hitung taksiran y untuk x = 2 dengan menggunakan interpolasi kuadratik untuk data: (1; 0), (4; 1:386294), (6; 1:791759). Penyelesaian. Koe…sien-koe…sien dari persamaan interpolasi kuadratik: a0 a1 a2
= 1; 1:386294 = 4 1
=
0 = 0 :4620981 ;
1:7917591:386294 64
6
0 1 386294 41 = 0:0518731 ; :
1
sehingga interpolasi kuadratik untuk data yang diberikan yaitu y = 1 + 0 :4620981 ( x
1) 0:0518731 2 (x 1) (x 4) :
Karena itu, nilai y untuk x = 2 yaitu y = 1 + 0 :4620981 (2
5.3
1) 0:0518731 2 (2 1)(2 4) = 0:5658444 :
Interpolasi Newton
Secara umum, n + 1 titik data, misalnya ( x0 ; y0 ), ( x1 ; y1 ), ::: , ( xn ; yn ) dapat dicocokkan dengan suatu polinomial berderajat n yang mempunyai bentuk y
= f n (x) = a 0 + a1 (x
x0) + a2 (x x0) (x x1) + +a (x x0 ) (x x1 ) (x x 1 ) : n
n
(5.1)
Bab 5. Interpolasi Polinomial
30
Persamaan-persamaan yang digunakan untuk menghitung koe…sien-koe…sien yaitu a0
= y0 ;
a1
=
a2
=
y2 y1 x2 x1
.. . a3
=
y0 = y [x1; x0] ; x0 = y [x2; x1] y [x1; x0] = y [x ; x ; x ] ; 2 1 0 x2 x0 x2 x0
y1 x1
y1 y0 x1 x0
y [xn ; xn1 ;:::;x2 ; x1 ] y [xn1 ; xn2 ;:::;x1 ; x0 ] = y [xn ; xn1 ;:::;x1 ; x0 ] : xn x0
Nilai fungsi berkurung siku dinamakan beda terbagi hingga dan dide…nisikan sebagai y [xi ; xi1 ] = y [xi2 ; xi1 ; xi ] = y [xn ; xn1 ;:::;x1 ; x0 ] =
yi yi1 (beda terbagi hingga pertama), xi xi1 y [xi2 ; xi1 ] y [xi1 ; xi ] (beda terbagi hingga kedua), xi2 xi y [xn ;:::;x1 ] y [xn1 ;:::;x0 ] (beda terbagi hingga ke-(5.2) n). xn x0
Persamaan-persamaan di atas adalah rekursif, yaitu beda orde lebih tinggi dihitung dengan mengambil beda dari orde lebih rendah. Sebagai contoh, penghitungan koe…sien-koe…sien dalam polinomial berderajat tiga dapat diperoleh secara berturut-turut mulai dari baris kedua dalam Tabel 5.1. Tabel 5.1: Tabel beda terbagi hingga untuk orde 3.
n
0 1 2 3
xn x0
yn b0 = y 0
x1 x2 x3
y1 y2 y3
pertama y y0 b1 = x11 x0 y2 y1 = b11 x2 x1 y y2 b12 = x33 x2
kedua b1 b2 = bx11 2 x0 b11 b21 = bx123 x1
ketiga b2 b3 = bx21 3 x0
Contoh 5.3 Hitung taksiran y untuk x = 2 dengan menggunakan interpolasi polinomial Newton untuk empat titik data: (1; 0) ; (4; 1:386294) ; (5; 1:609438) ; dan (6; 1:791759) : Penyelesaian. Disusun tabel beda terbagi hingga: n
xn
0 1 2 3
1 4 5 6
yn a0 = 0
1:386294 1:609438 1:791759
pertama a1 = 0:4621 a11 = 0:2231 a12 = 0:1824
kedua a2 = 0:0598 a21 = 0:0203
ketiga a3 = 0:0079
Bab 5. Interpolasi Polinomial
31
Diperoleh polinomial Newton orde ketiga: y = 0:4621 (x
1) 0:0598 (x 1) (x 4) + 0:0079 (x 1) (x 4) (x 5) :
(5.3)
Karena itu, nilai y untuk x = 2 yaitu y
= 0 + 0 :4621 (2 = 0:6289:
1) 0:0598 (2 1)(2 4) + 0:0079 (2 1)(2 4)(2 5)
Dicatat bahwa e…siensi komputasi yang lebih baik untuk menuliskan polinomial Newton (5.1) adalah dalam bentuk perkalian bersarang seperti y = ((
(a (x x 1) + a 1) (x x 2) + ) + a1) (x x0) + a0; n
n
n
n
(5.4)
untuk mendapatkan koe…sien-koe…sien dari suku-suku xn , xn1 , ..., x2 , x, dan x0 . Dalam MatLab ini dapat dilakukan dengan X; % X = [x_0 x_1 x_2 ... x_n] a; % a = [a_0 a_1 a_2 ... a_n] K = koef(n+1); % koefisien dari x^n for i = n:-1:1 % koefisien dari x^n, ..., x^0 K = [K a(i)] - [0 K*X(i)]; % a(i)*(x - x_(i - 1))+a_(i - 1) end
Fungsi MatLab poliNewton() menyusun Tabel Beda Terbagi Hingga seperti Tabel 5.1, mengkonstruksi polinomial Newton dan selanjutnya menghitung nilai y untuk suatu nilai x yang diberikan. Fungsi ini juga menampilkan plot titik-titik data dan kurva polinomial Newton. Di sini dicatat bahwa koe…sien-koe…sien dari polinomial dinyatakan sebagai vektor yang disusun dalam derajat yang menurun. function [T,K,yi]=poliNewton(X,Y,xi) % poliNewton Menyusun Tabel Beda Terbagi Hingga untuk interpolasi polinomial Newton % dan menghitung nilai fungsi interpolasi untuk suatu nilai x = xi. % % Input: X = data x % Y = data y yang berkorespondensi dengan x % xi = suatu nilai di antara titik-titik data x % % Output: T = Tabel Beda Terbagi Hingga mulai pertama % K = vektor yang memuat koefisien dari x^n, x^(n-1), ..., x, x^0 % yi = f(xi) % % ---PENGHITUNGAN INTI: N = length(X)-1; % banyaknya pasangan data % konstruksi tabel beda terbagi hasil = zeros(N+1,N+3); hasil(:,1) = [0:N]’; % iterasi hasil(:,2) = X’; % nilai x hasil(:,3) = Y’; % nilai y
Bab 5. Interpolasi Polinomial
32
for j=4:N+3 for i=1:N+3-(j-1) % beda terbagi orde 1,2, ..., n-1 hasil(i,j) = (hasil(i+1,j-1)-hasil(i,j-1))/(hasil(j+i-3,2)-hasil(i,2)); end end Tabel = hasil; % % T a
---OUTPUT: tabel beda terbagi hingga: = Tabel(:,4:end); = hasil(1,3:end); % a = [a_0 a_1 a_2 ... a_N]
% koef x^n, ..., x^0: K = a(N+1); % koefisien dari x^N for i = N:-1:1 K = [K a(i)] - [0 K*X(i)]; %a(i)*(x - x_(i - 1))+a_(i - 1) end % nilai y untuk suatu nilai xi: yi = sum(K’.*xi.^([N:-1:0]’)); % plot titik-titik data dan kurva polinomial xsim = X(1):0.1:X(end); for k=1:length(xsim) ysim(k) = sum(K’.*xsim(k).^([N:-1:0]’)); end plot(X,Y,’bo’,xsim,ysim,’r-’)
Sebagai contoh, dengan menjalankan fungsi MatLab poliNewton untuk data dalam 5.3 akan diperoleh fungsi polinomial Newton orde tiga yaitu y = 0:0079x3
0:1386x2 + 0:9894x 0:8587
yang adalah sama dengan hasil (5.3). Kurva dari fungsi tersebut diberikan dalam Gambar 5.1.
5.4
Interpolasi Lagrange
Interpolasi polinomial Lagrange hanyalah perumusan ulang dari polinomial Newton yang menghindari penghitungan beda terbagi hingga. Berikut ini adalah penurunan bentuk Lagrange secara langsung dari interpolasi polinomial Newton. Untuk kasus orde pertama dipunyai f 1 (x) = y 0 + ( x
x0) y [x1; x0] :
(5.5)
(5.6)
Beda terbagi hingga pertama y [x1 ; x0 ] =
y1 x1
y0 x0
dapat dirumuskan ulang sebagai y [x1 ; x0 ] =
y1 x1
x0
+
y0 x0
x1 :
Bab 5. Interpolasi Polinomial
33
Gambar 5.1: Kurva polinomial Newton untuk Contoh 5.3.
Karena itu, dengan mensubstitusikan (5.6) ke (5.5) akan dihasilkan f 1 (x) = y 0 +
x x1
x0 y1 + x x0 y0: x0 x0 x1
Selanjutnya dengan pengelompokkan suku-suku yang serupa dan penyederhanaan akan diperoleh bentuk Lagrange: f 1 (x) =
x x0
x1 y0 + x x0 y1: x1 x1 x0
Dengan cara yang serupa, yaitu beda-beda terbagi hingga dirumuskan ulang, akan dihasilkan polinomial Lagrange derajat n untuk n + 1 titik data, misalnya (x0 ; y0 ), (x1 ; y1 ), :::, ( xn ; yn ): n
y = f n (x) =
X
Li (x) yi
(5.7)
x : x
(5.8)
i=0
dengan
n
Li (x) =
Y
j =0;j 6 =i
x xi
j
j
Dalam hal ini, L i (x) dinamakan sebagai polinomial koe…sien Lagrange dan mempunyai sifat: 0 ; i = k Li (xk ) = : 1 ; i = k
6
Sebagai contoh, dari rumus di atas, versi orde kedua untuk polinomial Lagrange mem-
Bab 5. Interpolasi Polinomial
34
punyai bentuk yaitu (x
x1) (x x2) y0 + (x x0) (x x2) y1 (x0 x1 ) (x0 x2 ) (x1 x0 ) (x1 x2 ) (x x0 ) (x x1 ) + y2 : (x2 x0 ) (x2 x1 )
f 2 (x) =
Contoh 5.4 Diperhatikan kembali data dalam Contoh 5.3. Polinomial-polinomial koe …sien Lagrange: L0 (x) =
= L1 (x) =
= L2 (x) =
= L3 (x) =
=
(x
x1) (x x2) (x x3) = (x 4) (x 5) (x 6) (x0 x1 ) (x0 x2 ) (x0 x3 ) (1 4)(1 5)(1 6) 601 (x 4) (x 5) (x 6) = 601 () ; (x x0 ) (x x2 ) (x x3 ) (x 1) (x 5) (x 6) = (x1 x0 ) (x1 x2 ) (x1 x3 ) (4 1)(4 5)(4 6) 1 ( x 1) (x 5) (x 6) ; 6 (x x0 ) (x x1 ) (x x3 ) (x 1) (x 4) (x 6) = (x2 x0 ) (x2 x1 ) (x2 x3 ) (5 1)(5 4)(5 6) 14 (x 1) (x 4) (x 6) ; (x x0 ) (x x1 ) (x x2 ) (x 1) (x 4) (x 5) = (x3 x0 ) (x3 x1 ) (x3 x2 ) (6 1)(6 4)(6 5) 1 (x 1) (x 4) (x 5) : 10
Jadi, polinomial Lagrange yang sesuai dengan 4 titik data yang diberikan: y
=
=
601 (x 4) (x 5) (x 6) 0 + 16 (x 1) (x 5) (x 6) 1:386294 14 (x 1) (x 4) (x 6) 1:609438 + 101 (x 1) (x 4) (x 5) 1:791759 0:231049( x 1) (x 5) (x 6) 0:4023595 ( x 1) (x 4) (x 6) +0:1791759 ( x 1) (x 4) (x 5) :
Koe…sien-koe…sien dari suku-suku xn , xn1 , ..., x, x0 dalam polinomial Lagrange orde n dapat dibentuk melalui perkalian dua polinomial: (bn xn +
n nx
+ b1x + b0) (c
+
+ c1x + c0) = d2
2n
nx
+
d1x + d0
dengan min(k;n )
dk =
X
m=maks(0;k n)
bkm cm ,
untuk k = 2n; 2n
1; :::; 1; 0:
Operasi ini dapat dilakukan dengan menggunakan perintah conv() yang sudah diberikan oleh MatLab seperti diilustrasikan berikut ini.
Bab 5. Interpolasi Polinomial
35
>> b = [1 -1]; c = [1 1 1]; >> d = conv(a,b) d = 1 0 0 -1 % artinya bahwa (x-1)*(x^2+x+1)=x^3+0*x^2+0*x-1
Fungsi MatLab poliLagrange() mencari koe…sien-koe…sien dari polinomial Lagrange (5.7) bersama-sama dengan setiap polinomial koe…sien Lagrange (5.8). Fungsi ini juga menampilkan plot titik-titik data dan kurva polinomial Lagrange. Di sini dicatat bahwa koe…sien-koe…sien dari polinomial dinyatakan sebagai vektor yang disusun dalam dera jat yang menurun. function [K,L,yi]=poliLagrange(X,Y,xi) % poliLagrange Mencari koefisien dari polinomial Lagrange dan menghitung % nilai fungsi interpolasi untuk suatu nilai x = xi. % % Input: X = data x % Y = data y yang berkorespondensi dengan x % xi = suatu nilai di antara titik-titik data x % % Output: K = vektor yang memuat koefisien dari x^n, x^(n-1), ..., x, x^0 % L = vektor yang memuat polinomial koefisien Lagrange % yi = f(xi) % % ---PENGHITUNGAN INTI: N = length(X)-1; % banyaknya pasangan data L = zeros(N+1,N+1); % membentuk polinomial koefisien Lagrange for i=1:N+1 P = 1; for j=1:N+1 if i~=j P = conv(P,[1 -X(j)])/(X(i)-X(j)); end end L(i,:) = P; end % menentukan koefisien dari polinomial Lagrange K = Y*L; % menghitung nilai y = f(xi): yi = sum(K’.*xi.^([N:-1:0]’)); % plot titik-titik data dan kurva polinomial xsim = X(1):0.1:X(end); for k=1:length(xsim) ysim(k) = sum(K’.*xsim(k).^([N:-1:0]’)); end plot(X,Y,’bo’,xsim,ysim,’r-’)
Bab 6
Interpolasi Spline Tujuan Pembelajaran:
Mengetahui interpolasi spline untuk mendapatkan kurva mulus di antara dua titik data. Mengaplikasikan interpolasi spline untuk sekumpulan data. Pada bab sebelumnya telah dibahas mengenai interpolasi titik-titik data ( x0 ; y0 ) sampai (xn ; yn ) menggunakan suatu polinomial berderajat n. Namun terdapat kasus dimana fungsi-fungsi ini memberikan hasil yang salah. Pendekatan alternatifnya adalah menerapkan polinomial-polinomial berderajat lebih rendah pada sebagian titik data. Polinomial penghubung tersebut dinamakan fungsi-fungsi spline. De…nisi 6.1 Suatu fungsi f (x) dinamakan suatu spline berderajat k jika 1. Domain dari S adalah suatu interval [a; b]. 2. S; S 0 ;:::;S (k1) kontinu pada [a; b]. 3. Terdapat titik-titik xi sedemikian sehingga a = x0 < x1 < ::: < xn = b dan juga S adalah suatu polinomial berderajat k pada setiap [xi ; xi+1 ]. Dengan kata lain, spline adalah potongan-potongan fungsi polinomial dengan turunanturunan memenuhi kendala-kendala kekontinuan tertentu. Ketika k = 1, spline dinamakan spline linear. Ketika k = 2, spline dinamakan spline kuadratik. Ketika k = 3, spline dinamakan spline kubik.
6.1
Spline Linear
Kita mencoba mencari suatu fungsi spline linear S (x) sedemikian sehingga S (xi ) = y i untuk 0 i n. Diambil
S (x) =
8>> < >>:
S 0 (x) S 1 (x)
x x1 x x2
; x0 ; x1
.. .
.. .
S n1 (x) ; xn1
x x
dimana setiap S i (x) adalah linear. 36
n
Bab 6. Interpolasi Spline
37
Diperhatikan fungsi linear S i (x). Garis ini melalui titik (xi ; yi ) dan (xi+1 ; yi+1 ), sehingga kemiringan dari S i (x) yaitu mi =
yi+1 xi+1
y : x i
i
Kita dapat juga mengatakan bahwa garis tersebut melalui titik (xi ; yi ) dan (x; S (x)) untuk sembarang x [xi ; xi+1 ], sehingga
2
mi =
S i (x) yi ; x xi
yang memberikan S i (x) = yi + mi (x xi ) yi+1 yi = yi + (x xi+1 xi
x ) i
(6.1)
Contoh 6.2 Konstruksi spline kuadratik untuk data berikut: x y
0:0 1:3
0:1 4:5
0:4 2:0
0:5 2:1
0:75 5:0
1:0 3:0
Penyelesaian. 4 :5 0:1 2 :0 S 1 (x) = 4:5 + 0:4 2 :1 S 2 (x) = 2:0 + 0:5 5:0 S 3 (x) = 2:1 + 0:75 3:0 S 4 (x) = 5:0 + 1:0
1:3 (x 0) = 1:3 + 32x; 0 4:5 (x 0:1) = 16 25 x; 3 3 0:1 2:0 (x 0:4) = 1:6 + x; 0:4 2:1 (x 0:5) = 3:7 + 11:6x; 0:5 5:0 (x 0:75) = 11 8x: 0:75
[0:0; 0:1] : S 0 (x) = 1:3 + [0:1; 0:4] : [0:4; 0:5] : [0:5; 0:75] : [0:75; 1:0] :
Jadi spline adalah potongan linear, yaitu linear di antara setiap titik data. Persamaan (6.1) dapat dituliskan kembali sebagai S i (x) = a i x + bi ;
dengan ai =
yi+1 xi+1
y x
i i
i = 0; 1;:::;n
dan bi = y i
1
(6.2)
a x : i i
Fungsi MatLab spline1() mengkonstruksi spline linear dengan menghitung koe…sienkoe…sien spline linear (6.2) untuk koordinat-koordinat (x; y) dari titik-titik data. Fungsi ini juga menghitung taksiran nilai y untuk suatu nilai x di antara titik-titik data. Selain itu, kurva spline linear akan ditampilkan oleh fungsi tersebut.
Bab 6. Interpolasi Spline function [S,yi] = spline1(X,Y,xi) % spline2 Mengkonstruksi fungsi spline linear untuk titk-titik data (x,y). % % Input: X = [x0,x1,...,xn] % Y = [y0,y1,...,yn] % xi = suatu nilai di antara x0 dan xn % % Output: S = [ai,bi], koefisien-koefisien spline linear y = ai*x + bi % yi = taksiran nilai pada x = xi % % ---PENGHITUNGAN INTI: if nargin < 3, xi = X(1); end n = length(X); S = []; % koefisien-koefisien spline linear for i=1:n-1 a = (Y(i+1)-Y(i))/(X(i+1)-X(i)); b = Y(i)-a*X(i); if xi>X(i) & xi
Gambar 6.1: Spline kuadratik untuk Contoh 6.2.
38
Bab 6. Interpolasi Spline
39
Dari contoh di atas telah ditunjukkan bahwa kekurangan utama spline-spline linear adalah ketidakmulusannya. Artinya, pada titik-titik data di mana dua spline bertemu, kemiringannya berubah secara mendadak. Secara formal ini berarti bahwa turunan pertama dari fungsi tidak kontinu pada titik-titik tersebut. Kelemahan ini diatasi oleh penggunaan polinomial spline orde yang lebih tinggi.
6.2
Spline Kuadratik
Tidak seperti spline linear, spline kuadratik tidak dide…nisikan sepenuhnya oleh nilainilai di xi . Berikut ini kita perhatikan alasannya. Spline kuadratik dide…nisikan oleh S i (x) = a i x2 + bi x + ci :
Jadi terdapat 3 n parameter untuk mende…nisikan S (x). Diperhatikan titik-titik data: x0 y0
x1 y1
x2 y2
xn yn
Syarat-syarat untuk menentukan 3 n parameter dijelaskan seperti berikut ini. 1. Setiap subinterval [xi ; xi+1 ], untuk i = 0; 1; 2;:::;n samaan berkaitan dengan S i (x), yaitu S i (xi ) = y i
dan
1, memberikan dua per-
S i (xi+1) = y i+1 :
Jadi dari sini dipunyai 2 n persamaan. 2. Syarat pada kontinuitas dari S 0 (x) memberikan suatu persamaan tunggal untuk setiap titik dalam xi , i = 1; 2;:::;n 1, yaitu
S i01 (xi ) = S i0 (xi ) :
Jadi dari sini dipunyai n 1 persamaan. Sekarang totalnya terdapat 3 n 1 persamaan, tetapi karena terdapat 3 n parameter yang tidak diketahui maka sistem mempunyai kekurangan ketentuan.
3. Pilihan-pilihan yang mungkin untuk melengkapi kekurangan ketentuan yaitu S 0 (x0 ) = 0 atau
S 00 (x0 ) = 0:
Sekarang dimisalkan zi = S i0 (xi ). Karena S i (xi ) = yi , S i0 (xi ) = zi , dan S i0 (xi+1 ) = zi+1 , maka kita dapat mende…nisikan S i (x) =
zi+1 zi (x 2 (xi+1 xi )
x )2 + z (x x ) + y : i
i
i
i
(6.3)
Bab 6. Interpolasi Spline
40
Selanjutnya, dengan pengambilan x = x i+1 diperoleh yi+1
y y +1 y yi+1 i
i
=
i
=
zi+1 zi (xi+1 2 (xi+1 xi )
i
i
i
i
i
i
x )2 + z (x +1 x ) + y z +1 z (x +1 x ) + z (x +1 x ) 2 z +1 + z (x +1 x ) : 2
= S i (xi+1 ) = i
i
i
i
i
i
i
i
i
i
Jadi, kita dapat menentukan zi+1 dari zi : zi+1 = 2
yi+1 xi+1
y z : x i
i
(6.4)
i
Contoh 6.3 Konstruksi spline kuadratik untuk data berikut: x y
0:0 1:3
0:1 4:5
0:4 2:0
0:5 2:1
dengan penetapan z0 = 0. Penyelesaian. Pertama kali dihitung nilai-nilai zi : z1 z2 z3
y1 x1 y2 = 2 x2 y3 = 2 x3
= 2
y0 z0 = 2 4:5 1:3 0 = 64; 0:1 0:0 x0 y1 z1 = 2 2:0 4:5 64 = 242 ; 0:4 0:1 3 x1 y2 z2 = 2 2:1 2:0 + 242 = 248 ; 0:5 0:4 3 3 x2
Jadi, fungsi spline kuadratik S (x): S 0 (x) =
x0)2 + z0 (x x0) + y0
x1)2 + z1 (x x1) + y1
z1 z0 ( x 2 (x1 x0 )
= 320x2 + 1:3 = 320 x2 + 1:3; untuk 0 :0 S 1 (x) =
= = S 2 (x) =
z2 z1 ( x 2 (x2 x1 )
(x 0:1)2 + 64 ( x 0:1) + 4:5 2170 9 1010 194 2170 x2 + x ; untuk 0 :1 x 0:4 9 9 45 z3 z2 ( x x2 )2 + z2 (x x2 ) + y2 2 (x3 x2 )
2450 242 (x 0:4)2 (x 0:4) + 2 3 3 2450 2 2202 4948 = x x+ ; untuk 0 :4 3 3 30 =
x 0:1
x 0:5
Persamaan (6.3) dapat dituliskan kembali sebagai S i (x) = a i x2 + bi x + ci ;
i = 0; 1;:::;n
1
(6.5)
Bab 6. Interpolasi Spline dengan ai =
zi+1 zi ; bi = z i 2 (xi+1 xi )
41
2a x ; c = a x2 z x + y : i i
i
i i
i i
i
Fungsi MatLab spline2() mengkonstruksi (6.5) untuk koordinat-koordinat titik data (x; y ) dengan suatu syarat batas S 0 (x0 ) = 0. Fungsi ini juga menghitung taksiran nilai y untuk suatu nilai x di antara titik-titik data. Selain itu, kurva dari spline kuadratik untuk titik-titik data yang diberikan juga akan ditampilkan.
Gambar 6.2: Spline kuadratik untuk Contoh 6.3.
function [S,yi] = spline2(X,Y,xi) % spline2 Mengkonstruksi fungsi spline kuadrat untuk titik-titik data (x,y). % % Input: X = [x0,x1,...,xn] % Y = [y0,y1,...,yn] % xi = suatu nilai di antara x0 dan xn % % Output: S = [ai,bi,ci], koefisien-koefisien spline kuadrat y = ai*x^2+bi*x+ci % yi = taksiran nilai pada x = xi % % ---PENGHITUNGAN INTI: if nargin < 3, xi = X(1); end z0 = 0; n = length(X); S = []; z_lama = z0; % koefisien-koefisien spline linear: for i=1:n-1 z_baru = 2*(Y(i+1)-Y(i))/(X(i+1)-X(i))-z_lama; a = 0.5*(z_baru-z_lama)/(X(i+1)-X(i)); b = z_baru-2*a*X(i+1); c = a*X(i)^2-z_lama*X(i)+Y(i); S = [S; a b c]; z_lama = z_baru;
Bab 6. Interpolasi Spline
42
end % nilai potongan fungsi f(xi): for i=1:n-1 if xi>X(i) & xi
6.3
Spline Kubik S ( x) S n-1 S i S 1 S 0
S i+1
S n-2
S i ( xi +1 )= f ( xi +1 )= S i +1( xi+1 ) S i′ ( xi +1 )= S i′+1 ( xi +1 ) S i′′( xi +1 )= S i′′+1 ( xi +1 )
x0
x1
x2
xi
xi+1
xi+2
xn-2
xn-1
xn
x
Gambar 6.3: Pendekatan dengan polinomial spline kubik.
Diketahui suatu fungsi f ( x) yang dibatasi oleh interval a dan b , dan memiliki sejumlah titik data a = x0 < x1 < x2 < ::: < xn = b. Interpolasi spline kubik S (x) adalah suatu potongan fungsi polinomial berderajat tiga (kubik) yang menghubungkan dua titik yang bersebelahan, lihat Gambar 6.3, dengan ketentuan, untuk i = 0; 1;:::;n 1:
(S0) Potongan fungsi pada subinterval [ xi ; xi+1 ], i = 0; 1;:::;n S i (x) = a i (x
1:
x )3 + b (x x )2 + c (x x ) + d : i
i
i
i
i
i
Bab 6. Interpolasi Spline
43
(S1) Pada setiap titik data x = x i , i = 0; 1;:::;n: S (xi ) = f ( xi ) :
(S2) Nilai-nilai fungsi harus sama pada titik-titik dalam: S i (xi+1 ) = S i+1 (xi+1) ;
i = 0; 1;:::;n
2:
(S3) Turunan-turunan pertama pada titik dalam harus sama: S i0 (xi+1 ) = S i0+1 (xi+1 ) ;
i = 0; 1;:::;n
2:
(S4) Turunan-turunan kedua pada titik dalam harus sama: S i00 (xi+1 ) = S i00+1 (xi+1 ) ;
i = 0; 1;:::;n
2:
(S5) Salah satu syarat batas di antara dua syarat batas x0 dan xn berikut ini harus dipenuhi:
S 00 (x0) = S 00 (x ) = 0 (disebut batas alamiah/ natural boundary ) S 0 (x0) = f 0 (x0) dan S 0 (x ) = f 0 (x ) (disebut batas apitan/ clamped boundn
n
n
ary )
Berikut ini akan dijelaskan bagaimana mencari polinomial spline kubik S . Langkah 1 (Kendala S0): Polinomial spline kubik S untuk suatu fungsi f dide…nisikan oleh (6.6) S i (x) = a i (x xi )3 + bi (x xi )2 + ci (x xi ) + di
dengan i = 0; 1;:::;n
1.
Langkah 2 (Kendala S1): Ketika x = x i maka dipunyai S i (xi ) = ai (xi
x )3 + b (x x )2 + c (x x ) + d i
i
i
i
i
= di = f ( xi ) :
i
i
i
(6.7)
Ini berarti bahwa d i selalu menjadi pasangan titik data dari x i . Dengan pola ini maka pasangan titik data xi+1 adalah ai+1, konsekuensinya S (xi+1) = d i+1. Langkah 3 (Kendala S2): Ketika x = x i+1 disubstitusikan ke persamaan ( S 3) maka diperoleh S i (xi+1 ) = S i+1 (xi+1)
x )3 + b (x +1 x )2 + c (x +1 x ) + d = d +1 dengan i = 0; 1;:::;n 2. Sekarang dimisalkan h = x +1 x , sehingga persamaan ai (xi+1
i
i
i
i
i
i
i
i
i
i
i
i
di atas bisa dituliskan kembali menjadi
di+1 = a i h3i + bi h2i + ci hi + di :
(6.8)
Bab 6. Interpolasi Spline
44
Langkah 4 (Kendala S3): Turunan pertama dari (6.6) yaitu S i0 (x) = 3 ai (x
x )2 + 2b (x x ) + c : i
i
i
i
Ketika x = x i dipunyai S i0 (xi ) = 3 ai (xi
x )2 + 2b (x x ) + c = c ; i
i
i
i
i
i
dan ketika x = x i+1 dipunyai S i0 (xi+1 ) = S i0+1 (xi+1 )
3ai (xi+1
x )2 + 2b (x +1 x ) + c i
i
i
i
= ci+1
i
3ai h2i + 2bi hi + ci = ci+1 :
(6.9)
Langkah 5 (Kendala S4): Turunan kedua dari (6.6) yaitu S i00 (x) = 6 ai (x S 00 (x)
Dengan ketentuan tambahan
2
x ) + 2b : i
i
(6.10)
, persamaan di atas dimodi…kasi menjadi
S i00 (x) = 3 ai (x
x )+b : i
i
Ketika x = x i dipunyai S i00 (xi ) = 3 ai (xi
x ) + b = b ; i
i
i
dan ketika x = x i+1 dipunyai S i00 (xi+1 ) = S i00+1 (xi+1)
3ai (xi+1
x )+b i
i
3ai hi + bi
Dari sini bisa dinyatakan ai =
= bi+1 = bi+1:
1 (bi+1 3hi
b ):
i
(6.11)
Karena itu, persamaan (6.8) dapat dituliskan kembali menjadi di+1
1 (bi+1 3hi
=
h2i
=
3
b ) h3 + b h2 + c h + d i
i
i i
i i
(2bi + bi+1) + ci hi + di ;
i
(6.12)
(6.13)
sedangkan persamaan (6.9) menjadi ci+1
= hi (bi+1 bi ) + 2 bi hi + ci = hi (bi + bi+1) + ci
atau dinyatakan menjadi ci = h i1 (bi1 + bi ) + ci1 :
Bab 6. Interpolasi Spline
45
Dari persamaan (6.12) diperoleh ci =
1 hi
(di+1
d ) h3 (2b + b +1) i
i
i
i
(6.14)
dan untuk ci1 dinyatakan ci1 =
1 hi1
(di
d 1) h 31 (2b 1 + b ) : i
i
i
i
(6.15)
Disubstitusikan ci dan ci1 ke persamaan (6.13) untuk memperoleh 3
hi1 bi1 + 2 (hi1 + hi ) bi + hi bi+1 =
hi
(di+1
d ) h 31 (d d 1) (6.16) i
i
i
i
1 dengan i = 1; 2;:::;n 1. Dalam sistem persamaan ini, nilai hi ni=0 dan nilai n n di i=0 sudah diketahui, sedangkan nilai bi i=0 belum diketahui dan memang nilai inilah yang akan dihitung dari persamaan tersebut.
f g
f g
f g
Langkah 6 (Kendala S5): (1) Batas alamiah. Ketika S 00 (x0 ) = S 00 (xn ) = 0, dari persamaan (6.10) diperoleh S 00 (x0 ) = 6a0 (x0 S 00 (xn ) = 6an (xn
x0) + 2b0 = 0 =) b0 = 0; x ) + 2b = 0 =) b = 0: n
n
n
Sistem persamaan (6.16) dapat dituliskan sebagai perkalian matriks Az = dimana 1
A
=
h0
0 2 (h0 + h1 )
h1
0
h1
2 (h1 + h2 )
2 66 66 4 2 3 66 77 4 5 0
z
=
b0 b1
.. .
;
B =
bn
2 666 64
0
3 h1
3 hn1
(dn
0 0 h2 h 2 0 n
d1) .. .
d 1) n
3 h0
(d1
3 hn2
0 0 0
3 77 7 7 5
2 (hn2 + hn1 ) hn1 0 1
0
(d2
B
d0)
(dn1
0
d 2) n
3 777 75
:
(2) Batas apitan. Dari persamaan (6.14), dengan i = 0, dimana f 0 (x0 ) = S 0 (x0 ) = c 0 diperoleh f 0 (x0 ) =
1 h0
(d1
d0) h30 (2b0 + b1) ;
dan konsekuensinya 2h0 b0 + h0 b1 =
3 h0
(d1
d0) 3c0:
(6.17)
Sementara itu, pada x = xn dan dengan menggunakan persamaan (6.13)
;
Bab 6. Interpolasi Spline
46
diperoleh f 0 (xn ) = c n = c n1 + hn1 (bn1 + bn )
dengan cn1 bisa diperoleh dari persamaan (6.15), dengan i = n cn1 =
1 hn1
(dn
1,
d 1) h 31 (2b 1 + b ) : n
n
n
n
Jadi 1
f 0 (xn ) =
hn1
1
=
hn1
(dn
d 1) h 31 (2b 1 + b ) + h 1 (b 1 + b )
(dn
d 1) + h 31 (b 1 + 2b )
n
n
n
n
n
n
n
n
n
n
n
dan akhirnya diperoleh hn1 bn1 + 2 hn1 bn = 3cn
h 31 (d d 1) n
n
(6.18)
n
Persamaan (6.17) dan (6.18) ditambah persamaan (6.16) membentuk perkalian matriks Az = B dimana
A
=
2h0
h0
h0
2 (h0 + h1 )
0
h1
2 66 664 2 3 66 77 4 5
0
z
=
b0 b1
.. .
bn
;
B =
3
0
0 0 0 h1 2 (h1 + h2 ) h2 0 0 h 2 2 (h 2 + h 1) h 1 2h 1 0 h 1 1 (d1 d0 ) c0 1 (d2 d1 ) 1 (d1 d0 ) n
n
n
n
n
2 666 64
n
3 777 75
h0
h1
1 hn1
(dn
.. .
h0
1
d 1) (d 1 d 2) c 1 (d d 1 ) n
n
hn2 n
hn1
n
n
n
:
Langkah 7: Setelah sistem persamaan diselesaikan untuk b0 , b1 , ..., bn , kita mensubstitusikan nilai-nilai tersebut ke kendala (S1), (6.14), (6.11) untuk mendapatkan koe…sien-koe…sien lain dari spline kubik: (S1)
di = yi ; ci
(6:14)
=
h (2b + b +1) ; a 3
di+1 di hi
i
i
i
i
(6:14)
=
1 (bi+1 3hi
b ) : (6.19) i
Contoh 6.4 Konstruksikan spline kubik untuk 4 titik data berikut: x 0 y 0
1 1
2 4
3 5
terhadap syarat batas: S 0 (x0 ) = S 0 (0) = c 0 = 2 dan S 0 (xn ) = S 0 (3) = c n = 2.
3 77 775
;
Bab 6. Interpolasi Spline
47
Penyelesaian. Lebar subinterval pada sumbu x: h0 = h 1 = h 2 = h 3 = 1
dan beda terbagi pertama, dengan mengingat bahwa di = f ( xi ) = y i , yaitu d1
d0 = 1 ;
d2
h0
d1 = 3;
h1
d3
d2 = 1 :
h2
Persamaan matriks dapat dituliskan sebagai
2 64
2 1 0 0
1 4 1 0
0 1 4 1
0 0 1 2
yang mempunyai penyelesaian b0 =
32 75 64
b0 b1 b2 b3
3 75
=3
2 64
1 3 1 2
2 1 3 1
3 75
=
2 64
3 6 6 3
3 75
;
3; b1 = 3; b2 = 3, dan b3 = 3.
Disubstitusikan penyelesaian tersebut ke persamaan (6.19) untuk memperoleh koe…sienkoe…sien lain dari spline kubik: d0 c0 a0
= 0; d1 = 1; d2 = 4; 1 = 1 (3 + 2 ( 3)) = 2 ; c1 = 3 3 3 ( 3) 3 3 = = 2; a1 = = 3 3
13 (3 + 2 (3)) = 2 ; c2 = 1 13 (3 + 2 (3)) = 2; 2; a2 = 3 3(3) = 2 :
Terakhir, kita dapat menuliskan persamaan spline kubik seperti S 0 (x) = 2x3 S 1 (x) = S 2 (x) =
3x2 + 2x, untuk x 2 [0; 1] ; 2 (x 1)3 + 3 (x 1)2 + 2 (x 1) + 1, untuk x 2 [1; 2] ; 2 (x 2)3 3 (x 2)2 + 2 (x 1) + 4, untuk x 2 [2; 3] .
Fungsi MatLab splinekubik() mengkonstruksi persamaan matriks Ax = b , menyelesaikannya menggunakan iterasi SOR untuk mendapatkan koe…sien-koe…sien spline kubik untuk koordinat-koordinat x , y dari titik-titik data yang dilengkapi syarat batas apitan. function [S,yi] = splinekubik(X,Y,dy0,dyn,xi) % splinekubik Mengkonstruksi polinomial spline kubik untuk titik-titik data (x,y). % % Input: X = [x0,x1,...,xn] % Y = [y0,y1,...,yn] % dy0 = S’(x0) = c0: turunan awal % dyn = S’(xn) = cn: turunan akhir % xi = suatu nilai di antara x0 dan xn % % Output: S = [ai,bi,ci,di], koefisien-koefisien spline kubik
Bab 6. Interpolasi Spline
Gambar 6.4: Spline kuadratik untuk Contoh 6.4.
% Si(x)= ai*(x-xi)^3+bi*(x-xi)^2+ci(x-xi)+di % yi = taksiran nilai pada x = xi % % ---PENGHITUNGAN INTI: if nargin < 5, xi = X(1); end if nargin < 4, dyn = 0; end if nargin < 3, dy0 = 0; end n = length(X); for i = 1:n-1 % lebar subinterval h(i) = X(i+1) - X(i); end % mengkonstruksi matriks A dan B: A = zeros(n,n); B = zeros(n,1); A(1,1) = 2*h(1); A(1,2) = h(1); A(2,1) = h(1); A(n,n) = 2*h(n-1); A(n-1,n) = h(n-1); A(n,n-1) = h(n-1); for i = 2:n-1 for j = 2:n-1 if j==i, A(i,j) = 2*(h(i-1)+h(i)); elseif j==i+1, A(i,j) = h(i); elseif j==i-1, A(i,j) = h(i); else A(i,j) = 0; end end end B(1) = (Y(2)-Y(1))/h(1) - dy0; B(n) = dyn - (Y(n)-Y(n-1))/h(n-1); for i = 2:n-1 B(i) = (Y(i+1)-Y(i))/h(i) - (Y(i)-Y(i-1))/h(i-1); end B = 3*B; % koefisien-koefisien spline kubik: z = iterasiSOR(A,B);
48
Bab 6. Interpolasi Spline bi = z(1:3); di = Y(1:3)’; for i = 1: n-1 ci(i) = (Y(i+1)-Y(i))/h(i) - h(i)/3*(2*z(i)+z(i+1)); ai(i) = (z(i+1)-z(i))/(3*h(i)); end S = [ai’ bi ci’ di]; % orde menurun % nilai potongan fungsi f(xi): for i=1:n-1 if xi>X(i) & xi
49
Bab 7
Regresi Kuadrat Terkecil Tujuan Pembelajaran:
Mengetahui teknik pencocokan kurva yang merepresentasikan trend data. Mengaplikasikan regresi kuadrat terkecil untuk menentukan fungsi hampiran. Jika terdapat banyak galat yang berhubungan dengan data, khususnya data eksperimental, maka interpolasi tidak sesuai dan mungkin memberikan hasil-hasil yang tidak memuaskan jika dipakai untuk meramalkan nilai-nilai antara. Untuk kasus yang demikian, suatu strategi yang sesuai adalah dengan pencocokan kurva. Pencocokan kurva adalah pencarian suatu kurva yang bisa menunjukkan kecenderungan ( trend ) dari himpunan data. Kurva ini tidak harus melalui titik-titik data. Suatu kriteria yang dipakai untuk mengukur kecukupan dari kecocokan yaitu regresi kuadrat terkecil.
7.1
Regresi Linear
Regresi linier digunakan untuk menentukan fungsi linier yang paling sesuai dengan kumpulan titik data (xk ; yk ) : k = 1;:::;n yang diketahui. Pernyataan matematis untuk fungsi linear tersebut yaitu
f
g
y = a 0 + a1 x + e
dengan e dinamakan galat atau sisa. Sisa adalah selisih antara pengamatan dengan garis: e = y a0 + a1 x:
Suatu kriteria untuk pencocokan yang terbaik adalah hampiran kuadrat terkecil yang meminimalkan jumlahan kuadrat dari sisa: n
S r =
n
X X 2
ei =
i=1
i=1
n
(yi;observasi
y
2
i;data )
=
X i=1
(yi
a0 + a1x )2 : i
Kriteria ini menghasilkan suatu garis tunggal untuk himpunan data yang diberikan. Untuk menentukan nilai-nilai a0 dan a1 , diturunkan S r terhadap setiap koe…sien dan
50
Bab 7. Regresi Kuadrat Terkecil
51
selanjutnya disamakan dengan nol: n
@S @a 0
=
@S @a 1
=
X X 2
(yi
a0 a1x ) = 0 ;
(yi
a0 a1x ) x = 0:
i=1 n
2
i=1
i
i
i
Persamaan-persamaan di atas dapat dituliskan kembali menjadi n
n
n
X X X X X X yi
a0
i=1
i=1
n
a0 xi
i=1
= 0;
a1 x2i
= 0;
i=1 n
n
yi xi
a1 xi
i=1
i=1
atau ekivalen dengan na0 + a1
n
X X
=
xi
i=1 n
n
a0
n
X X
X
xi + a1
i=1
yi ;
i=1 n
x2i
=
i=1
yi x i :
i=1
Selanjutnya diselesaikan kedua persamaan untuk memperoleh n
n a1 =
n
n
X X X X X ! xi yi
xi
i=1
i=1
n
i=1
2
n
x2i
n
yi
dan a0 =
xi
i=1
1 n
n
X
yi
i=1
1
a1 n
n
X
xi :
i=1
i=1
Contoh 7.1 Berikut ini akan dicari persamaan kurva linear jika diberikan tujuh data untuk x dan y: 2 3 4 5 6 7 xi 1 yi 0:5 2:5 2:0 4:0 3:5 6:0 5:5 Berdasarkan data diperoleh: n = 7;
7
7
7
7
X
X
X
X
xi yi = 119 :5;
i=1
xi = 28;
i=1
i=1
yi = 24;
x2i = 140 :
i=1
Karena itu a1
=
a0
=
7 119:5 28 24 = 0 :8392857 ; 7 140 (28)2 24 28 0:8392857 = 0:07142857 ; 7 7
Bab 7. Regresi Kuadrat Terkecil
52
sehingga persamaan kurva linear: y = 0:07142857 + 0 :8392857 x:
Fungsi MatLab reglin() menjalankan skema kuadrat terkecil untuk mencari koe…sien dari suatu persamaan kurva linear yang cocok dengan sekumpulan titik data. Fungsi ini juga menampilkan plot titik-titik data dan kurva regresi linear. function Koef = reglin(X,Y) % reglin Mencari koefisien dari persamaan kurva linear untuk sekumpulan % titik data menggunakan kriteria kuadrat terkecil. % % Input: X = data x % Y = data y yang berkorespondensi dengan x % % Output: koefisien dari x^0 dan x^1 secara berturutan % % ---PENGHITUNGAN INTI: n = length(X); sigmaXY = sum(X.*Y); sigmaX = sum(X); sigmaY = sum(Y); sigmaXX = sum(X.^2); koefx1 = (n*sigmaXY-sigmaX*sigmaY)/(n*sigmaXX-sigmaX^2); % koefisien dari x^1 koefx0 = (sigmaY-koefx1*sigmaX)/n; % koefisien dari x^0 % ---OUTPUT: Koef = [koefx0,koefx1]; % orde naik % plot titik-titik data dan kurva linear xreg = X(1):0.01:X(end); yreg = Koef(1)+Koef(2)*xreg; plot(X,Y,’bo’,xreg,yreg,’r-’);
7.2
Regresi Polinomial Derajat Dua
Diberikan kumpulan titik data (xk ; yk ) : k = 1;:::;n . Fungsi pendekatan untuk fungsi polinomial berderajat dua: y = a 0 + a1 x + a2 x2 :
f
g
Jumlahan kuadrat dari sisa yaitu n
S r =
X i=1
yi
a0 a1x a2x2 2 : i
i
Bab 7. Regresi Kuadrat Terkecil
53
Gambar 7.1: Kurva linear dengan metode kuadrat terkecil untuk titik-titik data dalam Contoh 7.1.
Diturunkan S r terhadap semua parameter dan selanjutnya disamakan dengan nol: n
@S r @a 0
=
@S r @a 1
=
@S r @a 1
=
X X X 2
yi
a0 a1x a2x2
yi
a0 a1x a2x2
yi
a0 a1x a2x2
i=1 n
2
i=1 n
2
i=1
i
i
i
i
= 0; xi = 0; x2i = 0
i
i
2
X X X
Persamaan-persamaan disusun kembali menjadi na0 + a1
X X
xi + a1
i=1 n
a0
n
X X X
xi + a2
i=1 n
n
a0
n
X X X
i=1
xi
=
i=1 n
x2i + a2
i=1 n
x2i + a1
n
i=1 n
x3i
=
i=1 n
x3i + a2
i=1
yi ; xi yi ;
i=1 n
x4i
=
i=1
x2i yi ;
i=1
atau dalam bentuk perkalian matriks: n
n
n
2 X X 3 2X 3 66 77 2 3 66 77 66 X X X 77 4 5 66 X 77 66 X X X 77 66 X 77 4 5 4 5 n
i=1 n
n
i=1 n
i=1 n
x2i
xi
i=1 n
x2i
i=1
x2i
xi
x3i
i=1 n
x3i
i=1
yi
a0 a1 a2
i=1 n
xi yi
=
i=1 n
x4i
i=1
x2i yi
i=1
:
Bab 7. Regresi Kuadrat Terkecil
54
Contoh 7.2 Berikut ini akan dicari persamaan kurva polinomial derajat dua untuk data: 1 2 3 4 5 xi 0 yi 2:1 7:7 14 27 41 61 Berdasarkan data diperoleh: n = 6;
6
6
X
X
xi = 15;
i=1
6
X
xi yi = 585:6;
i=1
2
xi = 55 ;
i=1
6
X
4
xi = 979 ;
i=1
6
6
X
X
yi = 152:6;
i=1
6
X
x3i = 225 ;
i=1
x2i yi = 2488 :8:
i=1
Karena itu dipunyai sistem persamaan linear dalam bentuk
2 4
6 15 55 15 55 225 55 225 979
yang mempunyai penyelesaian
32 3 2 54 5 4 a0 a1 a2
=
152:6 585:6 2488:8
3 5
a0 = 2:47857 ; a1 = 2:35929; a2 = 1:86071 ;
sehingga persamaan kurva polinomial derajat dua: y = 2:47857 + 2 :35929 x + 1:86071x2 :
Fungsi MatLab regkuad() mencari koe…sien dari suatu persamaan polinomial derajat dua yang cocok dengan sekumpulan titik data dengan metode kuadrat terkecil. function Koef = regkuad(X,Y) % reglin Mencari koefisien dari persamaan polinomial derajat dua untuk % sekumpulan titik data menggunakan kriteria kuadrat terkecil. % % Input: X = data x % Y = data y yang berkorespondensi dengan x % % Output: koefisien dari x^0, x^1, dan x^2 secara berturutan % % ---PENGHITUNGAN INTI: % penyusun sistem persamaan: A(1,1) = length(X); A(1,2) = sum(X); A(1,3) = sum(X.^2); A(2,1) = A(1,2); A(2,2) = A(1,3); A(2,3) = sum(X.^3); A(3,1) = A(2,2); A(3,2) = A(2,3); A(3,3) = sum(X.^4); b(1) = sum(Y); b(2) = sum(X.*Y); b(3) = sum(X.^2.*Y);
Bab 7. Regresi Kuadrat Terkecil
55
b = b’; % penyelesaian sistem dengan SOR: Koef = SOR(A,b)’; % orde naik % plot titik-titik data dan kurva regresi kuadrat xreg = X(1):0.01:X(end); yreg = Koef(1)+Koef(2)*xreg+Koef(3)*xreg.^2; plot(X,Y,’bo’,xreg,yreg,’r-’);
Gambar 7.2: Kurva polinomial derajat dua dengan metode kuadrat terkecil untuk titiktitik data dalam Contoh 7.2.
7.3
Linearisasi Fungsi Tak Linear
Regresi linear memberikan teknik yang ampuh untuk mencocokkan garis "terbaik" terhadap data. Namun, teknik ini tergantung pada kenyataan bahwa kaitan antara variabel tak bebas dan bebas adalah linear. Dalam analisis regresi, seharusnya langkah pertama adalah penggambaran gra…k untuk memeriksa apakah pada data berlaku suatu hubungan linear. Beberapa data yang tidak linear dapat dilinearkan dengan suatu transformasi data, seperti yang disajikan dalam Tabel 7.1.
Contoh 7.3 Berikut ini akan dicari persamaan kurva eksponensial untuk data: xi yi ln (yi )
1 0:5 0:6931
2 1:7 0:5306
3 3:4 1 :2238
4 5:7 1 :7405
5 8:4 2:1282
Bab 7. Regresi Kuadrat Terkecil
56
Tabel 7.1: Linearisasi dari fungsi tak linear dengan transformasi data.
Fungsi Pencocokan (1) (2) (3) (4) (5) (6) (7) (8) (9)
Fungsi Linearisasi
a y = + b x b y = a+x y = a bx y = b eax y = C b eax y = a xb
y = aX + b
x 1 a Y = ; A = ; B = y b b Y = ln(y) ; A = ln (b) ; B = ln (a) Y = ln(y) ; B = ln (b) Y = ln(C y) ; A = a; B = ln (b) Y = ln(y) ; A = b; X = ln (x) ; B = ln (a) y Y = ln ; A = b; B = ln (a) x C Y = ln 1 ; B = ln (b) y X = ln (x)
1
Y = Ax + B
= Ax + B = ax + B = Ax + B = AX + B
Y Y Y Y
y =
Y = ax + B
Y = Ax + B
C 1 + b eax y = a ln (x) + b
1
X =
y = ax e
bx
Substitusi Variabel/ Perbaikan Parameter
y = aX + b
Jadi di sini digunakan fungsi linearisasi seperti no. 4 dalam Tabel 7.1. Berdasarkan data diperoleh: n = 5;
5
5
5
5
X
X
X
X
xi zi = 21:6425;
i=1
xi = 15;
i=1
zi = 4:93;
i=1
x2i = 55 :
i=1
Karena itu a = B
=
5 21:6425 15 4:93 = 0 :685; 5 55 (15)2 4:93 15 0:685 = 1:069; 5 5
sehingga persamaan kurva linear: Y = 0 :685
1:069x
atau persamaan kurva eksponensial, dengan b = e B = e 1:069 , yaitu y = e 0:685x1:069
Bab 8
Diferensiasi Numerik Tujuan Pembelajaran:
Mencari rumus hampiran beda untuk turunan pertama dan kedua. Mengetahui rumus hampiran turunan yang mempunyai galat lebih baik. Mengaplikasikan rumus hampiran beda untuk suatu fungsi. 8.1
Turunan Pertama
Suatu teori sederhana mengenai hampiran numerik untuk turunan dapat diperoleh melalui ekspansi deret Taylor dari f ( x + h) di sekitar x: 0
f ( x + h) = f ( x) + hf (x) +
h2
2!
00
f (x) +
h3
3!
f 000 (x) +
:
(8.1)
Pengurangan f ( x) dari kedua sisi dan membagi kedua sisi dengan ukuran langkah h menghasilkan f ( x + h) f ( x) h 00 h2 000 0 Df 1 (x; h) = = f (x) + f (x) + f (x) + 2! 3! h 0 = f (x) + O(h);
(8.2)
dimana O(h) menyatakan suku galat pemotongan yang sebanding dengan h untuk 1. Dari sini kita mendapatkan hampiran beda maju ( forward di¤erence aph proximation ) untuk f 0 (x):
j j
Df 1 (x; h) =
f ( x + h) h
f (x)
yang mempunyai galat sebanding dengan ukuran langkah h atau ekivalen dalam orde dari h. Sekarang, untuk menurunkan rumus hampiran yang lain untuk turunan pertama yang mempunyai galat lebih kecil, dihapus suku orde satu terhadap h dari persamaan (8.2) dengan mensubstitusikan 2 h untuk h dalam persamaan: f ( x + 2h) Df 1 (x; 2h) = 2h
f (x) = f 0 (x) + 2 h f 00 (x) + 4 h2 f 000 (x) + 2!
57
3!
Bab 8. Diferensiasi Numerik
58
dan pengurangan hasil ini dari dua kali persamaan (8.2) menghasilkan 2Df 1 (x; h)
f ( x + h) f ( x) f ( x + 2h) f ( x) 2h h f ( x + 2h) + 4 f ( x + h) 3f (x) Df 2 (x; h) = 2h 2 2 h 000 = f 0 (x) f (x) +
D 1 (x; 2h)
= 2
f
3!
= f 0 (x) + O h2
(8.3)
yang dapat dipandang sebagai suatu perbaikan atas persamaan (8.2) karena ini mempunyai galat pemotongan sebanding dengan h2 untuk h 1.
j j
Berikutnya, dengan mensubstitusikan h untuk h dalam persamaan (8.1) akan diperoleh hampiran beda mundur (backward di¤erence approximation ) untuk f 0 (x):
Db1 (x; h) =
f ( x)
f (x h) = D 1 (x; h) f
h
yang galat pemotongan yang sebanding dengan h untuk h 1. Untuk menghasilkan suatu versi perbaikan yang mempunyai galat pemotongan sebanding dengan h2 untuk 1 dapat diproses seperti berikut: h
j j
j j
2Db1 (x; h) Db1 (x; 2h) 3f ( x) = 2 1 0 = f (x) + O h2 :
Db2 (x; h) =
4f (x h) + f (x 2h) 2h
(8.4)
Untuk menurunkan rumus hampiran lain untuk turunan pertama, kita mengambil ekspansi deret Taylor dari f ( x + h) dan f ( x h) sampai orde kelima:
f (x + h) = f ( x) + hf 0 (x) +
1 + h5 f (5) (x) + 5! f (x
h)
1 2 00 1 1 h f (x) + h3 f 000 (x) + h4 f (4) (x) 2! 3! 4!
;
(8.5)
hf 0 (x) + 2!1 h2f 00 (x) 3!1 h3f 000 (x) + 4!1 h4f (4) (x) 5!1 h5f (5) (x) +
= f ( x)
(8.6)
dan membagi selisih antara kedua persamaan dengan 2 h untuk mendapatkan hampiran beda pusat (central di¤erence approximation ) untuk f 0 (x): Dc2 (x; h) =
f ( x + h)
f (x h) = f 0 (x) + 1 h2f 000 (x) + 1 h4f (5) (x) +
2h 0
3!
5!
2
= f (x) + O h
yang mempunyai galat sebanding dengan h 2 , serupa dengan persamaan (8.3) dan (8.4). Ini juga dapat diproses untuk menghasilkan suatu versi perbaikan yang mempunyai
Bab 8. Diferensiasi Numerik
59
galat pemotongan sebanding dengan h4 : 22 Dc2 (x; h)
D 2 (x; 2h)
f (x h) f (x + 2h) f (x 2h) 2h 2 2h 4 12 h (5) f (x) 3f 0 (x) 5! 22 D 2 (x; h) D 2 (x; 2h) 22 1 8f ( x + h) 8f ( x h) f ( x + 2h) + f ( x 2h)
= 4
c
= Dc4 (x; h) =
=
f ( x + h)
c
c
12h
4
0
= f (x) + O h
Contoh 8.1 Hitung nilai turunan f ( x) = x 2 pada x = 2 dengan h = 0:1. Penyelesaian. Hampiran beda maju:
f (2) = 4:41 4 = 4 :1; 0:1 0:1 f (2:2) + 4f (2:1) 3f (2) = 4:84 + 17:64 12 = 4 : D 2 (2; 0:1) =
O (h) : Df 1 (2; 0:1) = O h2
:
f (2 :1)
f
0:2
0:2
Hampiran beda mundur:
f (1:9) = 4 3:61 = 3:9; 0:1 0:1 3f (2) 4f (1 :9) + f (1 :8) 12 14:44 + 3 :24 D 2 (2; 0:1) = = = 4:
O (h) : Db1 (2; 0:1) = O h2
:
f (2)
b
0:2
0:2
Hampiran beda pusat:
f (1:9) = 4:41 3:61 = 4; 0:2 0:2 8f (2 :1) 8f (1 :9) f (2 :2) + f (1 :8) D 4 (2; 0:1) = = 4:
O h2
: Dc2 (2; 0:1) =
O h4
:
c
f (2 :1)
1:2
Berikut ini adalah kode MatLab untuk mencari hampiran beda pusat orde dua untuk fungsi f ( x) = sin (x) di x = 31 untuk h yang bevariasi. f = inline(’sin(x)’); f1 = inline(’cos(x)’); h = 0.1; x = pi/3; fprintf(’ h Galat Mutlak\n’); fprintf(’=========================\n’); for i = 1:6 Dc2 = (f(x+h) - f(x-h))/(2*h); fprintf(’%7.1e %8.1e\n’,h,abs(f1(x)-Dc2))
Bab 8. Diferensiasi Numerik
60
h = h/10; end
Kode MatLab tersebut menghasilkan tabel seperti berikut: h Galat Mutlak ========================= 1.0e-001 8.3e-004 1.0e-002 8.3e-006 1.0e-003 8.3e-008 1.0e-004 8.3e-010 1.0e-005 7.8e-012 1.0e-006 4.1e-011
8.2
Ekstrapolasi Richardson
Prosedur yang sudah diturunkan untuk hampiran beda pada bagian sebelumnya dapat dirumuskan menjadi suatu rumus umum, dinamakan ekstrapolasi Richardson, untuk memperbaiki hampiran beda dari turunan-turunan seperti berikut ini: Df;n +1 (x; h) = Db;n+1 (x; h) = Dc;2(n+1) (x; h) =
2n Df;n (x; h) 2n 2n Db;n (x; h) 2n n 4 Dc;2n (x; h) 4n
f;n (x; 2h)
D 1 D 1 D 1
b;n (x; 2h)
;
[n: orde galat]
;
c;2n (x; 2h)
:
Seperti dalam beda bagi, penghitungan hampiran-hampiran beda juga dapat disusun dalam bentuk tabel piramida. Sebagai contoh, konstruksi ekstrapolasi Richardson untuk beda maju sampai hampiran orde 4 diberikan seperti dalam Tabel ??. Sementara itu, untuk beda pusat sampai hampiran orde 8 diberikan seperti dalam Tabel 8.2. Tabel 8.1: Tabel ekstrapolasi Richardson untuk beda maju sampai hampiran orde 4.
1: 2: 4: 7:
O (h) Df 1 (x; h) Df 1 (x; 2h) Df 1 (x; 4h) Df 1 (x; 8h)
O h2 3: Df 2 (x; h) 5: Df 2 (x; 2h) 8: Df 2 (x; 4h)
O h3 6: Df 3 (x; h) 9: Df 3 (x; 2h)
O h4 10: Df 4 (x; h)
Fungsi MatLab RichDiff1() mengkonstruksi tabel ekstrapolasi dan mendapatkan hampiran beda pusat untuk turunan pertama sampai orde n dimana n adalah bilangan genap positif.
Bab 8. Diferensiasi Numerik
61
Tabel 8.2: Tabel ekstrapolasi Richardson untuk beda pusat sampai hampiran orde 8.
1: 2: 4: 7:
O h2 Dc2 (x; h) Dc2 (x; 2h) Dc2 (x; 4h) Dc2 (x; 8h)
O h4 3: Dc4 (x; h) 5: Dc4 (x; 2h) 8: Dc4 (x; 4h)
O h6 6: Dc6 (x; h) 9: Dc6 (x; 2h)
O h8 10: Dc8 (x; h)
function [Dcf1,Dc] = RichDiff1(f,x0,h,n) % RichDiff1 Mengkonstruksi ekstrapolasi Richardson untuk hampiran beda pusat dari % turunan pertama fungsi f(x). % Input: f = fungsi f(x) % x0 = nilai dimana turunan f’(x0) akan dicari % h = ukuran langkah % n = orde galat pemotongan (harus bilangan genap positif) % Output: Dcf1 = nilai f’(x0) N = n/2; Dc = zeros(N,N); % matriks turunan-turunan Richardson % kolom pertama dalam matriks Dc: Dc2(x,2^k*h) for i = 1:N x = x0 + 2^(i-1)*h; f1 = f(x); x = x0 - 2^(i-1)*h; f2 = f(x); Dc(i,1) = (f1-f2)/(2^i*h); % hampiran turunan-turunan pertama untuk beda pusat end % matriks Dc mulai kolom kedua for j = 2:N for i = 1:N+1-j Dc(i,j) = (4^(j-1)*Dc(i,j-1)-Dc(i+1,j-1))/(4^(j-1)-1); end end % Turunan pertama untuk beda pusat sampai orde n: Dcf1 = Dc(1,end);
8.3
Turunan Kedua
Untuk memperoleh suatu rumus hampiran untuk turunan kedua, kita mengambil ekspansi deret Taylor dari f ( x + h) dan f (x h) sampai orde kelima seperti dalam persamaan (8.5) dan (8.6). Dijumlahkan kedua persamaan (untuk menghapus suku f 0 (x)) dan selanjutnya dilakukan pengurangan 2 f ( x) dari kedua sisi dan dibagi kedua sisi dengan h2 untuk memperoleh hampiran beda pusat untuk turunan kedua:
(2)
Dc2 (x; h) =
f (x + h)
= f 00 (x) +
2f (x) + f (x h) h2
1 2 (4) 2 h f (x) + h4 f (6) (x) + 12 6!
Bab 8. Diferensiasi Numerik
62
yang mempunyai suatu galat pemotongan orde h2 . Ekstrapolasi Richardson dapat digunakan untuk memanipulasi persamaan di atas sehingga suku h2 terhapus. Hasil yang diperoleh adalah suatu versi perbaikan: (2)
22 Dc2 (x; h) 22
D(2)2 (x; 2h) = 1 c
f (x + 2h) + 16f (x + h) 30f (x) + 16f (x h) f (x 2h) 12h2
= f 00 (x)
4
h (5) f (x) + : 90
atau dituliskan kembali menjadi (2)
Dc4 (x; h) =
f (x + 2h) + 16f (x + h) 30f (x) + 16f (x h) f (x 2h) 12h2
= f 00 (x) + O h4
yang mempunyai galat pemotongan orde h4 . Kode MatLab untuk mengkonstruksi tabel ekstrapolasi dan mendapatkan hampiran beda pusat untuk turunan kedua adalah sama seperti untuk turunan pertama tetapi penghitungan hampiran turunan-turunan pertama orde 2 pada kolom pertama diganti dengan rumus turunan kedua: Dc(i,1) = (f1-2*f(x0)+f2)/(2^(i-1)*h)^2;
Lebih lanjut, untuk mendapatkan rumus hampiran dari turunan tingkat tinggi dapat dilihat di [5].
Bab 9
Persamaan Diferensial Biasa Tujuan Pembelajaran:
Mengetahui metode-metode penyelesaian numeris untuk persamaan diferensial biasa. Mengetahui metode yang memberikan hasil paling dekat dengan hasil sesungguhnya. Mengaplikasikan metode untuk menyelesaikan persamaan diferensial biasa. Dalam bab ini, kita akan membicarakan beberapa metode penyelesaian numeris dari persamaan diferensial biasa (disingkat PDB) dimana variabel tak bebas y tergantung pada variabel bebas tunggal t.
9.1
Metode Euler
Diperhatikan persamaan diferensial tingkat satu: dy (t) = f ( t; y (t)) , dt
a
t b,
y (a) = :
(9.1)
Tahap awal penyelesaian pendekatan numerik adalah dengan menentukan titik-titik dalam jarak yang sama di dalam interval [ a; b], yaitu dengan menerapkan ti = a + i h,
i = 0; 1;:::;n
dimana h menyatakan jarak antar titik yang dirumuskan oleh h =
b
a n
yang juga biasa dikenal sebagai lebar langkah ( step size ). Berikutnya, turunan dalam persamaan diferensial diganti dengan suatu turunan numerik (diperkenalkan dalam bab sebelumnya). Untuk kesepakatan diperkenalkan singkatan yi = y (ti ). Metode Euler menghampiri turunan pertama di t = ti dalam persamaan (9.1) dengan persamaan yi0 =
dyi dt
y yt +1 +1 t i
i
i
i
=
yi+1 yi : h
Karena itu, pada saat t = t i , persamaan (9.1) dapat dituliskan sebagai
f (t ; y ) :
yi+1 yi h
i
63
i
Bab 9. Persamaan Diferensial Biasa
64
Jadi, metode Euler mendapatkan barisan numerik yi
n i
f g =0 yang dinyatakan sebagai
y0 yi+1
=
yi + hf ( ti ; yi ) ,
i = 0; 1; 2;:::;n
1:
Pengertian geometris untuk metode Euler diberikan dalam Gambar 9.1. y prediksi galat eksak y(t i) + h·f (t i, y(t i)) y(t i+1) y(t i) t i
t i+1
h
t
Gambar 9.1: Ilustrasi dari penurunan metode Euler.
Contoh 9.1 Selesaikan persamaan diferensial y0 (t) = y (t)
t2 + 1,
0
t 2,
y (0) = 0:5
menggunakan metode Euler dimana n = 10. Penyelesaian. Dicari jarak antar titik dalam interval [0 ; 2] yaitu h =
2
0 = 0 :2,
10
sehingga dipunyai titik-titik diskrit yang dirumuskan oleh ti = 0 + i (0:2) = (0:2) i,
i = 0; 1;:::; 10;
yaitu t0 t6
= 0:0; t1 = 0:2; t2 = 0:4; t3 = 0:6; t4 = 0:8; t5 = 1:0; = 1:2; t7 = 1:4; t8 = 1:6; t9 = 1:8; t10 = 2:0:
Karena diketahui bahwa f (t; y (t)) = y (t)
t2 + 1
dan
y (0) = 0:5,
maka persamaan Euler dapat dinyatakan sebagai y0 yi+1
= 0:5;
yi + 0 :2 yi
= 1:2yi
t2 + 1 i
0:2t2 + 0:2, i
i = 0; 1; 2; :::; 9:
(9.2)
Bab 9. Persamaan Diferensial Biasa
65
Penyelesaian pada saat i = 0: = y (0:2) 1:2y0 0:2t20 + 0 :2 = 1:2 (0:5) 0:2 02 + 0:2 = 0:8;
y1
pada saat i = 1: y2
1:2y1 0:2t21 + 0 :2 1:2 (0:8) 0:2 0:22 + 0:2 = 1:152;
= y (0:4) =
pada saat i = 2: y3
1:2y2 0:2t22 + 0 :2 1:2 (1:152) 0:2 0:42 + 0:2 = 1 :5504;
= y (0:6) =
dan seterusnya sampai pada i = 9: y10
= y (2) 1:2y9 0:2t29 + 0 :2 = 1:2 (1:5504) 0:2 0:82 + 0:2 = 4:8658;
Berikut ini dibuat program MatLab pdb_Euler yang menggunakan metode Euler untuk menyelesaikan persamaan diferensial (9.1). function yt = pdb_Euler(f,t,y0,n) % pdb_Euler Menyelesaikan PDB tingkat satu menggunakan metode Euler, dimana % y’(t) = f(t,y), t = [a,b], y(a) = y0. % Input: f = fungsi turunan pertama dari y, f=inline(’(t,y)’,’t’,’y’) % t = vektor yang memuat ujung-ujung interval dari t % y0 = syarat awal % n = lebar langkah/ banyak subinterval dalam interval t % Output: yt = [t y(t)] h = (t(2)-t(1))/n; tvek(1) = t(1); yvek(1) = y0; for i = 2:n+1 tvek(i) = tvek(i-1)+h; yvek(i) = yvek(i-1)+h*f(tvek(i-1),yvek(i-1)); % persamaan (9.2) end yt(:,1) = tvek’; yt(:,2) = yvek’;
Dicatat bahwa metode Euler mempunyai galat yang sebanding dengan h. Terdapat dua modi…kasi sederhana pada metode Euler yang memberikan keakuratan lebih baik. Dua modi…kasi tersebut yaitu metode Heun dan metode titik tengah. Kedua metode ini mempunyai galat yang sebanding dengan h2 . Metode Heun dikenal juga sebagai metode prediktor-korektor.
Bab 9. Persamaan Diferensial Biasa
9.2
66
Metode Heun
Metode Heun memperbaiki taksiran turunan pertama dengan mengambil rata-rata dari kedua turunan pada titik-titik ujung subinterval. Turunan di titik awal subinterval [ti ; ti+1 ] yaitu yi0 = f (ti ; yi ) : Sementara itu, taksiran untuk yi+1 dihitung menggunakan metode Euler:
y + hf (t ; y )
yi+1
i
i
i
(9.3)
yang selanjutnya digunakan untuk menaksir turunan di titik akhir subinterval: yi0+1 = f ( ti+1; yi+1 )
f (t + h; y + hf (t ; y )) : i
i
i
i
Karena itu, diperoleh rata-rata turunan pertama di t = t i yaitu yi0
f (t ; y ) + f (t +2h; y + hf (t ; y )) : i
i
i
i
i
i
(9.4)
Jadi, metode Heun diperoleh dengan mengganti f ( ti ; yi ) di persamaan (9.2) dengan ruas kanan dari persamaan (9.4): y0
=
y + h2 [ f (t ; y ) + f (t + h; y + hf (t ; y ))] dengan i = 0; 1; 2;:::;n 1. Di sini, persamaan (9.3) dinamakan prediktor dan peryi+1
i
i
i
i
i
i
i
samaan (9.4) dinamakan korektor. Pengertian geometris untuk metode Heun diberikan dalam Gambar 9.2. y
y
turunan= f (t i+1, y(t i+1))
turunan = 0.5( f (t i, y(t i)) + f (t i+1, y(t i+1))) turunan= f (t i, y(t i))
t i h
t i+1
t i
t
h
(a)
t i+1
t
(b)
Gambar 9.2: Ilustrasi dari penurunan metode Heun.
Strategi dari metode Heun dinyatakan dalam Algoritma 7 dan diimplementasikan dalam fungsi MatLab pdb_Heun(). Contoh 9.2 Selesaikan persamaan diferensial y0 (t) = y (t)
t2 + 1,
0
t 2,
y (0) = 0:5
Bab 9. Persamaan Diferensial Biasa
67
Algorithm 7 Algoritma Metode Heun Masukan: fungsi dari turunan pertama untuk y (t): f (t; y ) interval dari t: [t0 ; tn ] syarat awal: y0 = banyak subinterval: n Penghitungan Inti: b
a
dihitung lebar langkah: h = n untuk i = 0 : n 1, dihitung ti+1 = t 0 + ( i + 1) h; H1: yi0 = f (ti ; yi ) ; H2: yi+1 = y i + hyi0 ; H3: yi0+1 = f (ti+1 ; yi+1 ) ;
h
H4: yi+1 = y i + y0 + yi0+1 ; 2 i Hasil akhir: y1 ; y2 ;:::;yn
menggunakan metode Heun dengan lebar langkah h = 0:2. Penyelesaian. Dicatat bahwa t0
= 0:0; t1 = 0:2; t2 = 0:4; t3 = 0:6; t4 = 0:8; t5 = 1:0;
t6
= 1:2; t7 = 1:4; t8 = 1:6; t9 = 1:8; t10 = 2:0:
Strategi dari metode Heun dinyatakan oleh y0
= 0:5;
H1
:
yi0 = y i
H2
:
yi+1
H3
:
H4
:
t2 + 1 i
y + 0 :2y0 y0+1 y +1 t2+1 + 1 y +1 y + 0 :1 y0 + y0+1 i
i
i
i i
i
i
i
i
dengan i = 0; 1; 2; :::; 9. Penyelesaian pada saat i = 0: H1 : y00 = y 0 (0) = y 0
H4 :
02 + 1 = 0 :5 + 1 = 1:5
y0 + 0 :2y00 = 0:5 + 0:2 (1:5) = 0:8 y10 = y0 (0:2) y1 0:22 + 1 = 0 :8 0:04 + 1 = 1 :76 y1 = y (0:2) y0 + 0 :1 y00 + y10 = 0:5 + 0 :1 (3:26) = 0:8260;
H2 : y1 = y (0:2) H3 :
pada saat i = 1: H1 : y10 = y 0 (0:2) = y 1
0:22 + 1 = 0 :8260 0:04 + 1 = 1 :786
H2 : y2 = y (0:4) y1 + 0 :2y10 = 0:8260 + 0 :2 (1:786) = 1 :1832 H3 : y0 2 = y0 (0:4) y2 0:42 + 1 = 1 :1832 0:16 + 1 = 2 :0232
H4 : y2 = y (0:4)
y1 + 0:1
y10 + y20 = 0:8260 + 0 :1 (3:8092) = 1 :2069 :
Bab 9. Persamaan Diferensial Biasa
68
pada saat i = 2: H1 : y20 = y 0 (0:4) = y 2 0:42 + 1 = 1 :2069 0:16 + 1 = 2 :0469 H2 : y3 = y (0:6) y2 + 0 :2y20 = 1:2069 + 0 :2 (2:0469) = 1 :6163
H3 : y0 3 = y0 (0:6) H4 : y3 = y (0:6)
y3 0:62 + 1 = 1 :6163 0:36 + 1 = 2 :2563
y2 + 0:1
y20 + y30 = 1:2069 + 0 :1 (4:3032) = 1 :6372 ;
dan seterusnya sampai pada i = 9 yang diserahkan kepada pembaca sebagai latihan. H
function yt = pdb_Heun(f,t,y0,n) % pdb_Heun Menyelesaikan PDB tingkat satu menggunakan metode Heun % (prediktor-korektor), dimana % y’(t) = f(t,y), t = [a,b], y(a) = y0. % Input: f = fungsi turunan pertama dari y, f=inline(’(t,y)’,’t’,’y’) % t = vektor yang memuat ujung-ujung interval dari t % y0 = syarat awal % n = lebar langkah/ banyak subinterval dalam interval t % Output: yt = [t y(t)] h = (t(2)-t(1))/n; tvek(1) = t(1); yvek(1) = y0; for i = 2:n+1 tvek(i) = tvek(i-1)+h; fi_1 = f(tvek(i-1),yvek(i-1)); % H1 yi = yvek(i-1)+h*f(tvek(i-1),yvek(i-1)); %H2 fi = f(tvek(i),yi); % H3 yvek(i) = yvek(i-1)+h/2*(fi_1+fi); %H4 end yt(:,1) = tvek’; yt(:,2) = yvek’;
9.3
Metode Titik Tengah y
y turunan= f (t i+0.5, y(t i+0.5)) turunan= f (t i+0.5, y(t i+0.5))
t i
t i+0.5 0.5h
(a)
t
t i h
t i+1
t
(b)
Gambar 9.3: Ilustrasi penurunan metode titik tengah.
Metode ini menggunakan metode Euler untuk memprediksi suatu nilai y di titik tengah
Bab 9. Persamaan Diferensial Biasa
69
dari subinterval [ ti ; ti+1 ]: yi+ 1 = y (ti ) +
h
f ( ti ; yi ) : 2 Nilai prediksi tersebut digunakan untuk menghitung turunan di titik tengah: 2
yi0+ 1 2
= f ti+ 1 ; yi+ 1 : 2
2
Selanjutnya turunan tersebut digunakan untuk mengekstrapolasi secara linear dari ti ke ti+1 : yi+1 y (ti ) + h yi0+ 1 = y i + hf ti+ 1 ; yi+ 1 :
2
2
2
Pengertian geometris untuk metode Heun diberikan dalam Gambar 9.3. Strategi dari metode titik tengah dinyatakan dalam Algoritma 8 dan diimplementasikan dalam fungsi MatLab pdb_tengah() . Algorithm 8 Algoritma Metode Titik Tengah Masukan: fungsi dari turunan pertama untuk y (t): f (t; y ) interval dari t: [t0 ; tn ] syarat awal: y0 = banyak subinterval: n Penghitungan Inti: b
a
dihitung lebar langkah: h = n untuk i = 0 : n 1, dihitung ti+ 1 = t 0 + i + 21 h; 2 ti+1 = t 0 + ( i + 1) h;
T1: yi+ 1 = y i + 2
T2:
yi0+ 1 2
h
2
f ( ti ; yi ) ;
= f ti+ 1 ; yi+ 1 ; 2
2
T3: yi+1 yi + hyi0+ 1 ; 2 Hasil akhir: y1 ; y2 ;:::;yn
Contoh 9.3 Selesaikan persamaan diferensial y0 (t) = y (t)
t2 + 1,
0
t 2,
y (0) = 0:5
menggunakan metode titik tengah dengan lebar langkah h = 0:2. Penyelesaian. Dicatat bahwa t0
= 0:0; t1 = 0:2; t2 = 0:4; t3 = 0:6; t4 = 0:8; t5 = 1:0;
t6
= 1:2; t7 = 1:4; t8 = 1:6; t9 = 1:8; t10 = 2:0:
Bab 9. Persamaan Diferensial Biasa
70
Strategi dari metode titik tengah dinyatakan oleh y0
= 0:5;
T1
:
T2
:
2
yi0+ 1 = y i+ 1 2
2
T3
:
t2 + 1
yi+ 1 = y i + 0 :1 yi
t2+ i
i
y + 0 :2y0+
yi+1
i
+1
1 2
i
= 1:1yi
0:1t2 + 0:1 i
1 2
dengan i = 0; 1; 2; :::; 9. Penyelesaian pada saat i = 0: T1 : y 1 = y (0:1) = 1:1y0 2
T2 : y 01 = y 0 (0:1) = y 1 2
2
T3 : y1 = y (0:2)
0:1
02 + 0:1 = 1:1(0:5) + 0:1 = 0:65
0:12 + 1 = 0 :65 0:01 + 1 = 1 :64
y0 + 0:2y0
1 2
= 0:5 + 0 :2 (1:64) = 0:828
pada saat i = 1: T1 : y 3 = y (0:3) = 1:1y1 2
T2 : y 03 = y 0 (0:3) = y 3 2
T3 : y2 = y (0:4)
2
0:1
0:22 + 0:1 = 1:1 (0:828)
0:096 = 1:0068
0:32 + 1 = 1 :0068 0:91 = 1:9168
y1 + 0:2y0
= 0:828 + 0 :2 (1:9168) = 1 :2114;
3 2
pada saat i = 2 : T1 :
y 5 = y (0:5) = 1 :1y2
T2 :
y 05 = y 0 (0:5) = y 5
2
2
T3 :
y3 = y (0:6)
2
0:1
0:42 + 0:1 = 1:1 (1:2114)
0:084 = 1 :4165
0:52 + 1 = 1 :4165 + 0:75 = 2:1665
y2 + 0:2y0
5 2
= 1:2114 + 0 :2 (2:1665) 1:6447;
dan seterusnya sampai pada i = 9 yang diserahkan kepada pembaca sebagai latihan. H
function yt = pdb_tengah(f,t,y0,n) % pdb_tengah Menyelesaikan PDB tingkat satu menggunakan metode titik % tengah, dimana % y’(t) = f(t,y), t = [a,b], y(a) = y0. % Input: f = fungsi turunan pertama dari y, f=inline(’(t,y)’,’t’,’y’) % t = vektor yang memuat ujung-ujung interval dari t % y0 = syarat awal % n = lebar langkah/ banyak subinterval dalam interval t % Output: yt = [t y(t)] h = (t(2)-t(1))/n; tvek(1) = t(1); yvek(1) = y0; for i = 2:n+1 tgh = tvek(i-1)+1/2*h; tvek(i) = tvek(i-1)+h; ymid = yvek(i-1)+h/2*f(tvek(i-1),yvek(i-1)); % T1 fmid = f(tgh,ymid); %T2 yvek(i) = yvek(i-1)+h*fmid; %T3 end
Bab 9. Persamaan Diferensial Biasa
71
yt(:,1) = tvek’; yt(:,2) = yvek’;
9.4
Metode Runge-Kutta
Metode Euler, Heun, dan titik tengah merupakan versi-versi dari suatu klas besar pendekatan satu langkah yang dinamakan metode Runge-Kutta. Metode Runge-Kutta mencapai keakuratan dari suatu pendekatan Taylor tanpa memerlukan turunan-turunan tingkat tinggi. Bentuk umumnya adalah y0 yi+1
= = yi + hF ( ti ; yi ; h)
dengan i = 0; 1; 2;:::. Fungsi F dinamakan fungsi kenaikan yang dapat dituliskan dalam bentuk + am km F = a 1 k1 + a2 k2 +
dimana m dinamakan orde dari metode Runge-Kutta, ai adalah konstanta dan k1
= f ( ti ; yi ) ;
k2
= f ( ti + p1 h; yi + q 11 k1 h) ; = f ( ti + p2 h; yi + q 21 k1 h + q 22 k2 h) ; .. .
k3
km
= f ( ti + pm1 h; yi + q m1;1 k1 h + q m1;2 k2 h +
+ q 1 m
;m1 km1 h) :
Nilai-nilai ai , p i , dan q ij untuk suatu orde m dicari dengan cara menyamakan persamaan tersebut dengan ekspansi deret Taylor orde m. Metode Runge-Kutta orde 1 tidak lain adalah metode Euler. Serupa dengan itu, metode Heun adalah suatu contoh dari metode Runge-Kutta orde 2. Metode Runge-Kutta orde 4 adalah satu dari metode yang banyak digunakan untuk menyelesaikan persamaan diferensial. Metode ini mempunyai suatu galat pemotongan yang sebanding dengan h4 . Algoritmanya digambarkan seperti berikut ini dan diimplementasikan dalam fungsi MatLab pdb_RK4. yi+1 = y i +
h
6
( k1 + 2 k2 + 2 k3 + k4 ) ,
i = 0; 1;:::;n
dimana k1
= f (ti ; yi ) ;
h
k2
= f ti +
k3
= f ti +
k4
= f (ti + h; yi + hk3 ) :
2 h
2
; yi +
h
; yi +
2 h
2
k1 k2
; ;
1
Bab 9. Persamaan Diferensial Biasa
72
Contoh 9.4 Selesaikan persamaan diferensial y0 (t) = y (t)
t2 + 1,
0
t 2,
y (0) = 0:5
menggunakan metode Runge-Kutta orde 4 dengan lebar langkah h = 0:2. Penyelesaian. Dinyatakan bahwa f ( t; y (t)) = y (t)
t2 + 1,
dan t0 t6
= 0:0; t1 = 0:2; t2 = 0:4; t3 = 0:6; t4 = 0:8; t5 = 1:0; = 1:2; t7 = 1:4; t8 = 1:6; t9 = 1:8; t10 = 2:0:
Untuk menghitung y1 = y (0:2), dimana i = 0, tahap-tahap penghitungannya dimulai dari menghitung k1 sampai k5 :
02 + 1 = 1:5;
k1
= f ( t0 ; y0 ) = f (0 ; 0:5) = 0:5
k2
= f ( t0 + 0 :1; y0 + 0 :1k1 ) = f (0 :1; 0:65) = 0:65
k3
=
k4
=
0:12 + 1 = 1 :64; f ( t0 + 0 :1; y0 + 0 :1k2 ) = f (0 :1; 0:664) = 0 :664 0:12 + 1 = 1 :654; f ( t0 + 0 :2; y0 + 0 :2k3 ) = f (0 :2; 0:8308) = 0 :8308 0:22 + 1 = 1 :7908;
sehingga akhirnya diperoleh y1 = y (0:2) = y 0 +
0 :2 (1:5 + 2 1:64 + 2 1:654 + 1 :7908) = 0 :8293: 6
Dengan cara yang sama, y (0:4), :::, y (2) dapat dihitung dan diserahkan kepada pembaca sebagai latihan. H Perbandingan hasil dari metode Runge-Kutta orde 4 dengan metode Euler, metode Heun, dan metode titik tengah untuk contoh di atas diberikan dalam Tabel 9.1. Dicatat bahwa penyelesaian eksak dari persamaan diferensial adalah y (t) = ( t + 1)2
12 e
t
Dari tabel tersebut ditunjukkan bahwa metode Runge-Kutta orde 4 adalah yang terbaik, sedangkan metode Euler adalah yang paling buruk. Selain itu, dari tabel juga dapat dilihat bahwa metode titik tengah lebih baik daripada metode Heun. function yt = pdb_RK4(f,t,y0,n) % pdb_RK4 Menyelesaikan PDB tingkat satu menggunakan metode Runge-Kutta % orde 4, dimana % y’(t) = f(t,y), t = [a,b], y(a) = y0. % Input: f = fungsi turunan pertama dari y, f=inline(’(t,y)’,’t’,’y’) % t = vektor yang memuat ujung-ujung interval dari t % y0 = syarat awal % n = lebar langkah/ banyak subinterval dalam interval t % Output: yt = [t y(t)] h = (t(2)-t(1))/n;
Bab 9. Persamaan Diferensial Biasa
73
Tabel 9.1: Perbandingan hasil dari metode-metode penyelesaian untuk persamaan diferensial y (t) = y (t) t2 + 1, 0 t 2, y (0) = 0:5. 0
yi
ti Euler
0:0 0:2 0:4 0:6 0:8 1:0 1:2 1:4 1:6 1:8 2:0
0:5000 0:8000 1:1520 1:5504 1:9885 2:4582 2:9498 3:4518 3:9501 4:4282 4:8658
Heun 0:5000 0:8260 1:2069 1:6372 2:1102 2:6177 3:1496 3:6937 4:2351 4:7556 5:2331
Titik Tengah 0 :5000 0 :8280 1 :2114 1 :6447 2 :1213 2 :6332 3 :1705 3 :7212 4 :2706 4 :8010 5 :2904
tvek(1) = t(1); yvek(1) = y0; for i=2:n+1 tvek(i) = tvek(i-1)+h; k1 = f(tvek(i-1),yvek(i-1)); k2 = f(tvek(i-1)+0.5*h,yvek(i-1)+0.5*k1*h); k3 = f(tvek(i-1)+0.5*h,yvek(i-1)+0.5*k2*h); k4 = f(tvek(i),yvek(i-1)+k3*h); yvek(i) = yvek(i-1)+(k1+2*k2+2*k3+k4)*h/6; end yt(:,1) = tvek’; yt(:,2) = yvek’;
Runge-Kutta orde 4 0:5000 0:8293 1:2141 1:6489 2:1272 2:6408 3:1799 3:7323 4:2834 4:8151 5:3054
Eksak 0:5000 0:8293 1:2141 1:6489 2:1272 2:6409 3:1799 3:7324 4:2835 4:8152 5:3055
Bab 10
Integrasi Numerik Tujuan Pembelajaran:
Mendapatkan rumus hampiran untuk menghitung integral dari suatu fungsi atas interval terbatas. Mengetahui aturan yang memberikan hasil lebih dekat dengan hasil sebenarnya. Pada bab ini kita tertarik untuk menghitung integral b
I =
Z
f ( x) dx
a
untuk f (x) R (dan tentu saja x R). Bentuk umum integrasi numerik dari suatu fungsi f (x) atas suatu interval [a; b] adalah suatu jumlahan berbobot dari nilai-nilai fungsi di berhingga ( n + 1) titik-titik diskrit, dan dinamakan kuadratur ( quadrature ):
2
2
n
b
Z
f ( x) dx =
a
X
wi f ( xi )
i=0
dengan a = x0 < x1 < < xn = b. Dalam kedua subbab berikut ini akan diturunkan empat bentuk persamaan integral numerik yang dikenal sebagai rumus eksplisit Newton-Cotes. Rumus ini ini didasarkan pada strategi penggantian fungsi f (x) dengan suatu hampiran fungsi f m (x) yang lebih mudah diintegralkan. Fungsi f m (x) tersebut adalah suatu polinomial derajat m.
10.1
Aturan Trapesium
Aturan ini didasarkan pada interpolasi linear dari fungsi f ( x) pada [xi1 ; xi ] untuk x setiap i = 1; 2;:::;n. Jadi, pada setiap subinterval, hampiran untuk x 1 f ( x) dx diberikan oleh luas bidang datar trapesium, lihat Gambar 10.1:
R
i
i
xi
Z
xi1
f (x) dx
12 [f (x ) + f (x 1)] (x x 1) = h2 [ f (x 1) + f (x )] : i
Dicatat bahwa h = x i dipunyai a
i
i
i
i
(10.1)
x 1 = b n a . Diaplikasikan (10.1) untuk i = 1 sampai i = n, i
b
Z
i
f ( x) dx
h
T (n) = 2
74
n
X i=1
[f (xi1 ) + f ( xi )] ;
Bab 10. Integrasi Numerik
75
dan jika jumlahan diekspansikan maka diperoleh h
T ( n) =
2 h
=
2
[ f (x0 ) + 2 f (x1 ) + 2 f ( x2 ) +
2f (x 1) + f (x )] n
n1
"
f ( x0 ) + 2
X
n
#
f ( xi ) + f (xn ) :
i=1
(10.2)
Dicatat bahwa aturan trapesium mempunyai galat pemotongan yang sebanding dengan h3 .
f ( x0)
f ( x1)
f ( x2)
f ( x3)
f ( xn-2) f ( xn-1)
...
f ( xn)
( x)
x0
x1
x2
x3
xn-2 xn-1
xn
x
Gambar 10.1: Ilustrasi dari aturan trapesium.
Contoh 10.1 Hitung luas bidang datar yang dibatasi oleh y = 2x3 dan sumbu x untuk x [0; 1] dengan menggunakan aturan trapesium dan lebar langkah h = 0:1.
2
Penyelesaian. Diperoleh tabel x f ( x)
0 0
0:1 0:002
0:2 0:016
0:3 0:054
0:4 0:128
0:5 0:25
0:6 0:432
0:7 0:686
0:8 1:024
0:9 1:458
1 2
Jadi, 1
T (10) =
Z
2x3 dx
0
0:1 [0 + 2 (0:002 + 0 :016 + 2
+ 1:024 + 1:458) + 2] = 0 :505:
Secara kalkulus, 1
L =
Z 0
1 4 2x dx = x 2 3
1
= 0 :5:
0
Dengan h = 0:1 ternyata terjadi kesalahan e = 0:505
j
0:5j = 0 :005.
Rumus integrasi menggunakan aturan trapesium diimplementasikan dalam fungsi MatLab trapesium() . function Tn = trapesium(f,x,n) % trapesium Menyelesaikan integral f(x) atas interval [a,b] menggunakan % aturan trapesium dengan n subinterval % Input: f = fungsi yang diintegralkan % x = vektor [a b] % n = lebar langkah/ banyak subinterval dalam interval x % Output: Tn = luas bidang datar yang dibatasi oleh f(x) dan sumbu x dalam
Bab 10. Integrasi Numerik
76
% interval [a,b] h = (x(2)-x(1))/n; xvek = x(1):h:x(2); yvek = f(xvek); % nilai f untuk semua titik diskrit Tn = h/2*(yvek(1)+sum(yvek(2:n))+yvek(n+1)); % persamaan(10.2)
10.2
Aturan-aturan Simpson
Taksiran yang lebih akurat dari suatu integral diperoleh jika polinomial derajat tinggi digunakan untuk menghubungkan titik-titik diskrit. Rumus-rumus yang dihasilkan dengan pengambilan integral dari polinomial tersebut dinamakan aturan-aturan Simpson. Pertama kali akan digunakan suatu polinomial derajat dua untuk mencocokkan titiktitik ( xi1 ; f ( xi1 )), (xi ; f ( xi )), dan ( xi+1; f (xi+1)). Dide…nisikan polinomial derajat dua P i1 (x) = a (x xi )2 + b (x xi ) + c: (10.3)
Dicatat bahwa subskrip ( i 1) di sini tidak menyatakan derajat dari polinomial, tetapi menyatakan titik ujung kiri dari interval [ xi1 ; xi+1 ] dimana kita mencocokkan polinomial derajat dua. Untuk kesepakatan diperkenalkan singkatan y i = f (xi ). Oleh karena itu, dari (10.3) kita membentuk tiga persamaan dalam a, b, c:
x )2 + b (x 1 x ) + c a (x x )2 + b (x x ) + c a (x +1 x )2 + b (x +1 x ) + c a (xi1
i
i
i
i
i
i
i
i
= yi1 ; = yi ;
i
= yi+1 :
i
i
Ini adalah suatu sistem persamaan linear, dan kita akan mengasumsikan bahwa h = xi xi1 = x i+1 xi , sehingga
a = b =
2y + y 1 ; 2h2 y +1 y 1 ; yi+1
i
i
c = yi :
i
i
2h
Diintegralkan polinomial derajat dua (10.3) dari x = xi1 sampai x = xi+1 dengan koe…sien-koe…sien tersebut menghasilkan hampiran xi+1
Z
xi1
xi+1
f (x) dx
Z
xi+1
P i1 (x) dx =
xi1
a (x
xi1
=
1 a (x 3
=
2h 3
Z h
yi+1
1 xi )3 + b (x 2
2
x ) i
2y + y 1 + 3y i
2
i
i
x ) i
xi+1
+ cx
=
2
xi1
h
3
+ b (x
i
x )+c i
2 = ah3 + 2ch 3
dx
( yi1 + 4 yi + yi+1 ) :
Seperti dalam bagian sebelumnya, kita ingin mengintegralkan f (x) pada [a; b]. Jadi, seperti sebelumnya, a = x0 dan b = xn . Jika n adalah suatu bilangan genap, maka banyaknya subinterval dalam [a; b] adalah suatu bilangan genap, dan karena itu kita
Bab 10. Integrasi Numerik
77
mempunyai hampiran b
Z Z
I =
a
f ( x) dx
x2
x4
P 0 (x) dx +
x0
=
h
3
Z
xn
P 2 (x) dx +
x2
( y0 + 4 y1 + 2 y2 + 4 y3 + 2 y4 +
+
Z
P n2 (x) dx
xn2
+ 2y 2 + 4 y 1 + y ) n
n
n
Dengan menggunakan notasi jumlahan, dari persamaan terakhir kita memperoleh aturan 31 -Simpson: S 3 (n) =
h
3
1 n 2
0 X @ y0 + 4
1 n 2
1
y2i1 + 2
i=1
X
y2i + yn
i=1
1 A
:
(10.4)
Dicatat bahwa aturan 31 -Simpson mempunyai galat pemotongan yang sebanding dengan h5 . Rumus integrasi menggunakan aturan 31 -Simpson diimplementasikan dalam fungsi MatLab Simpson3(). function S3n = Simpson3(f,x,n) % Simpson3 Menyelesaikan integral f(x) atas interval [a,b] menggunakan % aturan 1/3-Simpson dengan n subinterval % Input: f = fungsi yang diintegralkan % x = vektor [a b] % n = lebar langkah/ banyak subinterval dalam interval x % Output: S3n = luas bidang datar yang dibatasi oleh f(x) dan sumbu x dalam % interval [a,b] h = (x(2)-x(1))/n; xvek = x(1):h:x(2); yvek = f(xvek); % nilai f untuk semua titik diskrit S3n = h/3*(yvek(1)+4*sum(yvek(2:2:n))+2*sum(yvek(3:2:n-1))+yvek(n+1)); % persamaan(10.4)
Selanjutnya, jika digunakan suatu polinomial derajat tiga untuk mencocokkan titiktitik (xi1 ; f ( xi1 )), (xi ; f ( xi )), (xi+1; f (xi+1 )), dan (xi+2; f ( xi+2)), maka dengan cara yang sama seperti sebelumnya akan diperoleh hampiran xi+2
I =
Z
f ( x) dx
xi1
38h (y 1 + 3y + 3y +1 + y +2) : i
i
i
i
Jika n adalah suatu bilangan kelipatan tiga, maka kita mempunyai hampiran b
I =
Z Z a
=
f ( x) dx
x3
x0
x6
P 0 (x) dx +
Z
x3
xn
P 3 (x) dx +
+
Z
P n3 (x) dx
xn3
3h (y0 + 3 y1 + 3 y2 + 2 y3 + 3 y4 + 3 y5 + 2 y6 + 8
+ 2y 3 + 3y 2 + 3 y 1 + y ) : n
n
n
n
Bab 10. Integrasi Numerik
78
Dengan menggunakan notasi jumlahan, dari persamaan terakhir kita memperoleh aturan 83 -Simpson: S 4 (n) =
3h 8
1 n 3
0 X @ y0 + 3
1
(y3i2 + y3i1 ) + 2
i=1
3 8 -Simpson
Dicatat bahwa aturan dengan h5 .
1 n 3
X
y3i + yn
i=1
1 A
:
(10.5)
juga mempunyai galat pemotongan yang sebanding
Rumus integrasi menggunakan aturan 83 -Simpson diimplementasikan dalam fungsi MatLab Simpson4(). function S4n = Simpson4(f,x,n) % Simpson4 Menyelesaikan integral f(x) atas interval [a,b] menggunakan % aturan 3/8-Simpson dengan n subinterval % Input: f = fungsi yang diintegralkan % x = vektor [a b] % n = lebar langkah/ banyak subinterval dalam interval x % Output: S4n = luas bidang datar yang dibatasi oleh f(x) dan sumbu x dalam % interval [a,b] h = (x(2)-x(1))/n; xvek = x(1):h:x(2); yvek = f(xvek); % nilai f untuk semua titik diskrit sum1 = sum(yvek(2:3:n-1)); sum2 = sum(yvek(3:3:n)); sum3 = sum(yvek(4:3:n-2)); S4n = 3*h/8*(yvek(1)+3*(sum1+sum2)+2*sum3+yvek(n+1)); % persamaan(10.5)
Terakhir, jika digunakan suatu polinomial derajat empat untuk mencocokkan titik-titik (xi1 ; f ( xi1 )), ( xi ; f ( xi )), ( xi+1 ; f ( xi+1 )), ( xi+2 ; f ( xi+2 )), dan ( xi+3 ; f ( xi+3)), maka dengan cara yang sama seperti sebelumnya akan diperoleh hampiran xi+3
Z
I =
f (x) dx
xi1
245h (7 y 1 + 32 y + 12y +1 + 32 y +2 + 7 y +3) : i
i
i
i
i
Jika n adalah suatu bilangan kelipatan empat, maka kita mempunyai hampiran b
I =
Z Z a
=
f ( x) dx
x4
x8
P 0 (x) dx +
x0
Z
x4
xn
P 4 (x) dx +
+
Z
P n4 (x) dx
xn4
2h (7 y0 + 32 y1 + 12 y2 + 32 y3 + 14 y4 + 32 y5 + 12 y6 + 32 y7 + 14 y8 + 45 14yn4 + 32 yn3 + 12 yn2 + 32 yn1 + 7 yn ) :
Dengan menggunakan notasi jumlahan, dari persamaan terakhir kita memperoleh atur2 an 45 -Simpson: S 4 (n) =
2h 45
1 n 4
0 X @ 7y0 +
i=1
1 n 4
1
(32y4i3 + 12 y4i2 + 32 y4i1 ) + 14
X i=1
y4i + 7 yn
1 A
: (10.6)
Bab 10. Integrasi Numerik Dicatat bahwa aturan gan h7 .
2 45 -Simpson
79 mempunyai galat pemotongan yang sebanding den-
Rumus integrasi menggunakan aturan MatLab Simpson5() .
2 45 -Simpson
diimplementasikan dalam fungsi
function S5n = Simpson5(f,x,n) % Simpson5 Menyelesaikan integral f(x) atas interval [a,b] menggunakan % aturan 2/45-Simpson dengan n subinterval % Input: f = fungsi yang diintegralkan % x = vektor [a b] % n = lebar langkah/ banyak subinterval dalam interval x % Output: S4n = luas bidang datar yang dibatasi oleh f(x) dan sumbu x dalam % interval [a,b] h = (x(2)-x(1))/n; xvek = x(1):h:x(2); yvek = f(xvek); % nilai f untuk semua titik diskrit sum1 = sum(yvek(2:4:n-2)); sum2 = sum(yvek(3:4:n-1)); sum3 = sum(yvek(4:3:n)); sum4 = sum(yvek(5:3:n-3)); S5n = 2*h/45*(7*yvek(1)+32*sum1+12*sum2+32*sum3+14*sum4+7*yvek(n+1)); % pers.(10.6)
