Komputasi untuk Sains dan Teknik

Komputasi untuk Sains dan Teknik

Supriyanto Suparno ( Website: http://supriyanto.fisika.ui.edu ) ( Email: [email protected] atau [email protected] ) Edisi II Revisi terakhir tgl: 12 Februari 2008

Departemen Fisika-FMIPA, Univeristas Indonesia Dipublikasikan pertama kali pada September 2007

Untuk Nina Marliyani Muflih Syamil dan Hasan Azmi

Ketekunan adalah jalan yang terpercaya untuk mengantarkan kita menuju kesuksesan (Supriyanto, 2007)

Kata Pengantar

Secara garis besar, ilmu fisika dapat dipelajari lewat 3 jalan, yaitu pertama, dengan menggunakan konsep atau teori fisika yang akhirnya melahirkan fisika teori. Kedua, dengan cara eksperimen yang menghasilkan aliran fisika eksperimental, dan ketiga, fisika bisa dipelajari lewat simulasi fisika yang sangat mengandalkan komputer serta algoritma numerik. Tujuan penyusunan buku ini adalah untuk meletakkan pondasi dasar dari bangunan pemahaman akan metode-metode komputasi yang banyak digunakan pada simulasi-simulasi fenomena fisika. Rujukan utama buku ini bersumber pada buku teks standar yang sangat populer di dunia komputasi, yaitu buku yang ditulis oleh Richard L. Burden dan J. Douglas Faires dengan judul Numerical Analysis edisi ke-7, diterbitkan oleh Penerbit Brooks/Cole, Thomson Learning Academic Resource Center. Disamping itu, buku ini dilengkapi oleh sejumlah contoh aplikasi komputasi pada upaya penyelesaian problem-problem fisika. Dalam edisi ke-2 ini, algoritma numerik disalin ke dalam 2 bahasa pemrograman, yaitu Fortran77 dan Matlab. Disamping itu penjelasan lebih terperinci tentang bagaimana menentukan indeks i, j dan k dalam proses looping disajikan pada Bab I, untuk memberi pondasi yang sangat penting bagi berdirinya bangunan pemahaman akan teknik-teknik numerik selanjutnya. Akhirnya saya ingin mengucapkan rasa terima kasih yang tak terhingga kepada Dede Djuhana yang telah berkenan memberikan format LATEX-nya sehingga tampilan tulisan pada buku ini benar-benar layaknya sebuah buku yang siap dicetak. Rasa terima kasih juga ingin saya teruskan kepada Sarah Wardhani yang telah memicu langkah awal penulisan buku ini hingga masuk ke Edisi-2. Tak lupa, saya pun sepatutnya berterima kasih kepada seluruh rekan diskusi yaitu para mahasiswa yang telah mengambil mata kuliah Komputasi Fisika PTA 2006/2007 di Departemen Fisika, FMIPA, Universitas Indonesia. Tiga orang mahasiswi dari Universitas Pakuan yaitu Eni Nurliani, Saidah Al-adawiyah dan Deni Fitri A juga perlu saya tulis disini sebagai ungkapan terima kasih atas pertanyaan-pertanyaan mereka yang turut memperkaya isi buku ini. Walaupun buku ini masih jauh dari sempurna, namun semoga ia dapat menyumbangkan energi bagi terciptanya gelombang kebangkitan ilmu pengetahuan pada bangsa Indonesia yang saat ini sedang terpuruk. Saya wariskan ilmu ini untuk anak bangsa. Saya izinkan anda untuk meng-copy dan menggunakan buku ini selama itu ditujukan untuk belajar dan bukan untuk tujuan komersial. Bagi yang ingin berdiskusi, memberikan masukan, kritikan dan saran, silakan dikirimkan ke email: [email protected]

Depok, 16 September 2007 Supriyanto Suparno v

Daftar Isi

Lembar Persembahan

i

Kata Pengantar

v

Daftar Isi

vii

Daftar Gambar

xi

Daftar Tabel

xiii

1 Matrik dan Komputasi

1

1.1

Pengenalan matrik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Inisialisasi matrik dalam memori komputer . . . . . . . . . . . . . . . . . . . . . .

2

1.3

Macam-macam matrik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3.1

Matrik transpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3.2

Matrik bujursangkar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3.3

Matrik simetrik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3.4

Matrik diagonal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3.5

Matrik identitas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.3.6

Matrik upper-triangular . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.3.7

Matrik lower-triangular . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.3.8

Matrik tridiagonal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.3.9

Matrik diagonal dominan . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.3.10 Matrik positive-definite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.3.11 Vektor-baris dan vektor-kolom . . . . . . . . . . . . . . . . . . . . . . . . .

6

Operasi matematika . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.4.1

Penjumlahan matrik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.4.2

Komputasi penjumlahan matrik . . . . . . . . . . . . . . . . . . . . . . . .

7

1.4.3

Perkalian matrik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

1.4.4

Komputasi perkalian matrik . . . . . . . . . . . . . . . . . . . . . . . . . .

12

1.4.5

Perkalian matrik dan vektor-kolom . . . . . . . . . . . . . . . . . . . . . .

13

1.4.6

Komputasi perkalian matrik dan vektor-kolom . . . . . . . . . . . . . . . .

14

1.5

Penutup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

1.6

Latihan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

1.4

2 Metode Eliminasi Gauss

17

2.1

Sistem persamaan linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

2.2

Triangularisasi dan Substitusi Mundur . . . . . . . . . . . . . . . . . . . . . . . . .

18

vii

viii 2.3

Matrik dan Eliminasi Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

2.4

Algoritma eliminasi Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

2.4.1

Algoritma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

Contoh aplikasi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

2.5.1

Menghitung arus listrik . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

2.6

Menghitung invers matrik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

2.7

Penutup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

2.5

3 Aplikasi Eliminasi Gauss pada Masalah Inversi 3.1

39

Inversi Model Garis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

3.1.1

Script matlab inversi model garis . . . . . . . . . . . . . . . . . . . . . . . .

42

Inversi Model Parabola . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

3.2.1

Script matlab inversi model parabola . . . . . . . . . . . . . . . . . . . . .

48

3.3

Inversi Model Bidang . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

3.4

Contoh aplikasi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

3.4.1

53

3.2

Menghitung gravitasi di planet X . . . . . . . . . . . . . . . . . . . . . . . .

4 Metode LU Decomposition

61

4.1

Faktorisasi matrik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

4.2

Algoritma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65

5 Metode Iterasi

71

5.1

Kelebihan Vektor-kolom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

5.2

Pengertian Norm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

5.2.1

Script perhitungan norm dalam Matlab . . . . . . . . . . . . . . . . . . . .

72

5.2.2

Perhitungan norm-selisih . . . . . . . . . . . . . . . . . . . . . . . . . . . .

73

Iterasi Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

73

5.3.1

Script Matlab untuk menghitung iterasi . . . . . . . . . . . . . . . . . . . .

76

5.3.2

Optimasi script Matlab untuk menghitung iterasi . . . . . . . . . . . . . .

80

5.3.3

Algoritma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

5.3.4

Program dalam Fortran . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

Iterasi Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

84

5.4.1

Script iterasi Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . .

85

5.4.2

Algoritma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87

5.4.3

Script iterasi Gauss-Seidel dalam Fortran . . . . . . . . . . . . . . . . . . .

87

Iterasi dengan Relaksasi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

89

5.5.1

90

5.3

5.4

5.5

Algoritma Iterasi Relaksasi . . . . . . . . . . . . . . . . . . . . . . . . . . .

6 Interpolasi

93

6.1

Interpolasi Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93

6.2

Interpolasi Cubic Spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95

ix 7 Diferensial Numerik

103

7.1

Metode Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

7.2

Metode Runge Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

7.3

Metode Finite Difference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.3.1

Script Finite-Difference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

7.3.2

Aplikasi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

8 Pers. Diferensial Parsial Numerik

123

8.1

Pendahuluan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

8.2

PDP eliptik dalam Finite-Difference . . . . . . . . . . . . . . . . . . . . . . . . . . 124

9 Integral Numerik

131

9.1

Metode Trapezoida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

9.2

Metode Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

9.3

Metode Composite-Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

10 Metode Newton

137

10.1 Penyederhanaan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 11 Metode Monte Carlo

139

11.1 Penyederhanaan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 12 Inversi

143

12.1 Inversi Linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 12.2 Inversi Non-Linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 Daftar Pustaka

149

Indeks

151

Daftar Gambar

3.1

Grafik data pengukuran gerak batu . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

3.2

Grafik hasil inversi parabola . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

6.1

Fungsi f (x) dengan sejumlah titik data . . . . . . . . . . . . . . . . . . . . . . . .

95

6.2

Pendekatan dengan polinomial cubic spline . . . . . . . . . . . . . . . . . . . . . .

95

6.3

Profil suatu object . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

6.4

Sampling titik data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

6.5

Hasil interpolasi cubic spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

6.6

Hasil interpolasi lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

7.1

Metode Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

7.2

Trend error metode euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

7.3

Rangkaian RC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

7.4

Kurva muatan q terhadap waktu t . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

7.5

Kurva suatu fungsi f (x) yang dibagi sama besar berjarak h. Evaluasi kurva yang dilakukan Finite-Difference dimulai dari batas bawah X0 = a hingga batas atas x6 = b . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

8.1

Distribusi temperatur pada lempeng logam . . . . . . . . . . . . . . . . . . . . . . 126

9.1

Metode Trapezoida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

9.2

Metode Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

9.3

Metode Composite Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

10.1 Metode Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 11.1 Lingkaran dan bujursangkar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 11.2 Dart yang menancap pada bidang lingkaran dan bujursangkar . . . . . . . . . . . 140 11.3 Dart yang menancap pada bidang 1/4 lingkaran dan bujursangkar . . . . . . . . 141

xi

Daftar Tabel

3.1

Data temperatur bawah permukaan tanah terhadap kedalaman . . . . . . . . . .

39

3.2

Data temperatur bawah permukaan tanah terhadap kedalaman . . . . . . . . . .

44

3.3

Data ketinggian terhadap waktu dari planet X . . . . . . . . . . . . . . . . . . . .

53

5.1

Hasil akhir elemen-elemen vektor x hingga iterasi ke-10 . . . . . . . . . . . . . . .

79

5.2

Hasil perhitungan norm-selisih (dengan ℓ2 ) hingga iterasi ke-10 . . . . . . . . . .

81

5.3

Hasil Iterasi Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

85

5.4

Hasil perhitungan iterasi Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . .

90

5.5

Hasil perhitungan iterasi Relaksasi dengan ω = 1, 25 . . . . . . . . . . . . . . . . .

90

xiii

Bab 1

Matrik dan Komputasi

✍ Objektif : ⊲ Mengenalkan matrik dan jenis-jenis matrik. ⊲ Mengenalkan operasi penjumlahan dan perkalian matrik. ⊲ Mendeklarasikan elemen-elemen matrik ke dalam memori komputer. ⊲ Membuat script operasi matrik.

1.1

Pengenalan matrik

Notasi suatu matrik berukuran n x m ditulis dengan huruf besar dan dicetak tebal, misalnya An×m . Huruf n menyatakan jumlah baris, dan huruf m jumlah kolom. Suatu matrik tersusun dari elemen-elemen yang dinyatakan dengan huruf kecil diikuti angka-angka indeks, misalnya aij , dimana indeks i menunjukan posisi baris ke-i dan indeks j menentukan posisi kolom ke-j. 



a11

a12

. . . a1m

  a21 A = (aij ) =   ..  .

a22 .. .

 . . . a2m  ..   . 

(1.1)

an1 an2 . . . anm

Contoh 1: Matrik A2×3 A=

"

# 3 8 5 6 4 7

dimana masing-masing elemennya adalah a11 = 3, a12 = 8, a13 = 5, a21 = 6, a22 = 4, dan a23 = 7. Contoh 2: Matrik B3×2 

 1 3   B = 5 9 2 4 1

BAB 1. MATRIK DAN KOMPUTASI

2

dimana masing-masing elemennya adalah b11 = 1, b12 = 3, b21 = 5, b22 = 9, b31 = 2, dan b32 = 4.

1.2

Inisialisasi matrik dalam memori komputer

Dalam bahasa pemrograman Fortran77, cara mengisi memori komputer dengan elemen-elemen matrik A2×3 , sesuai dengan Contoh 1 adalah 1 2 3 4 5 6

A(1,1) A(1,2) A(1,3) A(2,1) A(2,2) A(2,3)

= = = = = =

3 8 5 6 4 7

Sedangkan untuk matrik B3×2 , sesuai Contoh 2 adalah 1 2 3 4 5 6

B(1,1) B(1,2) B(2,1) B(2,2) B(3,1) B(3,2)

= = = = = =

1 3 5 9 2 4

Sementara dalam Matlab, cara mengisi memori komputer dengan elemen-elemen matrik A2×3 , sesuai dengan Contoh 1 adalah 1 2

clear all clc

3 4 5 6 7 8 9 10

A(1,1) A(1,2) A(1,3) A(2,1) A(2,2) A(2,3) A

= = = = = =

3; 8; 5; 6; 4; 7;

Sedangkan untuk matrik B3×2 , sesuai Contoh 2 adalah 1 2

clear all clc

3 4 5 6 7 8 9 10

B(1,1) B(1,2) B(2,1) B(2,2) B(3,1) B(3,2) B

= = = = = =

1; 3; 5; 9; 2; 4;

1.3. MACAM-MACAM MATRIK

1.3

3

Macam-macam matrik

1.3.1 Matrik transpose Operasi transpose terhadap suatu matrik akan menukar elemen-elemen dalam satu kolom menjadi elemen-elemen dalam satu baris; demikian pula sebaliknya. Notasi matrik tranpose adalah AT atau At . Contoh 3: Operasi transpose terhadap matrik A

A=

"

3 8 5 6 4 7

#

  3 6   At = 8 4 5 7

1.3.2 Matrik bujursangkar Matrik bujursangkar adalah matrik yang jumlah baris dan jumlah kolomnya sama. Contoh 4: Matrik bujursangkar berukuran 3x3 atau sering juga disebut matrik bujursangkar orde 3

  1 3 8   A = 5 9 7 2 4 6

1.3.3 Matrik simetrik Matrik simetrik adalah matrik bujursangkar yang elemen-elemen matrik A bernilai sama dengan matrik transpose-nya (At ). Contoh 5: Matrik simetrik   2 −3 7 1   −3 5 6 −2  A= 7 6 9 8   1 −2 8 10



2

−3 7

1



  −3 5 6 −2  A = 7 6 9 8   1 −2 8 10 t

1.3.4 Matrik diagonal Matrik diagonal adalah matrik bujursangkar yang seluruh elemen-nya bernilai 0 (nol), kecuali elemen-elemen diagonalnya. Contoh 6: Matrik diagonal orde 3   11 0 0   A =  0 29 0  0 0 61


4 1.3.5 Matrik identitas

Matrik identitas adalah matrik bujursangkar yang semua elemen-nya bernilai 0 (nol), kecuali elemen-elemen diagonal yang seluruhnya bernilai 1. Contoh 7: Matrik identitas orde 3   1 0 0   I = 0 1 0 0 0 1 1.3.6 Matrik upper-triangular Matrik upper-tringular adalah matrik bujursangkar yang seluruh elemen dibawah elemen diagonal bernilai 0 (nol). Contoh 8: Matrik upper-triangular  3  0 A= 0  0

6 2 1



 4 1 5  0 8 7  0 0 9

1.3.7 Matrik lower-triangular Matrik lower-tringular adalah matrik bujursangkar yang seluruh elemen diatas elemen diagonal bernilai 0 (nol). Contoh 9: Matrik lower-triangular 

12

  32 A= 8 

−5

0 −2 7

10

0

0



 0  11 0  6 9 0

1.3.8 Matrik tridiagonal Matrik tridiagonal adalah matrik bujursangkar yang seluruh elemen bukan 0 (nol) berada disekitar elemen diagonal, sementara elemen lainnya bernilai 0 (nol). Contoh 10: Matrik tridiagonal   3 6 0 0   2 −4 1 0   A= 0 5 8 −7   0 0 3 9

1.3. MACAM-MACAM MATRIK

5

1.3.9 Matrik diagonal dominan Matrik diagonal dominan adalah matrik bujursangkar yang memenuhi |aii | >

n X

j=1,j6=i

|aij |

(1.2)

dimana i=1,2,3,..n. Coba perhatikan matrik-matrik berikut ini   7 2 0   A = 3 5 −1 0 5 −6



6

4

 B= 4

−3

−3

−2



 0 1

0

Pada elemen diagonal aii matrik A, |7| > |2|+|0|, lalu |5| > |3|+|−1|, dan |−6| > |5|+|0|. Maka

matrik A disebut matrik diagonal dominan. Sekarang perhatikan elemen diagonal matrik B, |6| < |4| + | − 3|, | − 2| < |4| + |0|, dan |1| < | − 3| + |0|. Dengan demikian, matrik B bukan matrik diagonal dominan.

1.3.10 Matrik positive-definite Suatu matrik dikatakan positive-definite bila matrik tersebut simetrik dan memenuhi xt Ax > 0

(1.3)

Contoh 11: Diketahui matrik simetrik berikut 

2

 A = −1 0

−1 2

−1

0



 −1 2

untuk menguji apakah matrik A bersifat positive-definite, maka

xt Ax =

h

i



2

−1

0



x1



   x1 x2 x3 −1 2 −1 x2  0 −1 2 x3   2x1 − x2 i h   = x1 x2 x3 −x1 + 2x2 − x3  −x2 + 2x3

= 2x21 − 2x1 x2 + 2x22 − 2x2 x3 + 2x23

= x21 + (x21 − 2x1 x2 + x22 ) + (x22 − 2x2 x3 + x23 ) + x23 = x21 + (x1 − x2 )2 + (x2 − x3 )2 + x23

Dari sini dapat disimpulkan bahwa matrik A bersifat positive-definite, karena memenuhi x21 + (x1 − x2 )2 + (x2 − x3 )2 + x23 > 0


6 kecuali jika x1 =x2 =x3 =0. 1.3.11 Vektor-baris dan vektor-kolom

Notasi vektor biasanya dinyatakan dengan huruf kecil dan dicetak tebal. Suatu matrik dinamakan vektor-baris berukuran m, bila hanya memiliki satu baris dan m kolom, yang dinyatakan sebagai berikut i i h h a = a11 a12 . . . a1m = a1 a2 . . . am

(1.4)

Sedangkan suatu matrik dinamakan vektor-kolom berukuran n, bila hanya memiliki satu kolom dan n baris, yang dinyatakan sebagai berikut 

a11





a1



     a21   a2     a=  ..  =  ..  .   . an1 an

1.4

(1.5)

Operasi matematika

1.4.1 Penjumlahan matrik Operasi penjumlahan pada dua buah matrik hanya bisa dilakukan bila kedua matrik tersebut berukuran sama. Misalnya matrik C2×3 C=

"

9 5 3 7 2 1

#

dijumlahkan dengan matrik A2×3 , lalu hasilnya (misalnya) dinamakan matrik D2×3 D=A+C

D =

"

=

"

7 2 1 # 3+9 8+5 5+3

=

"

13

# 3 8 5 6 4 7

+

"

# 9 5 3

6+7 4+2 7+1 # 12 13 8 6

8

Tanpa mempedulikan nilai elemen-elemen masing-masing matrik, operasi penjumlahan antara matrik A2×3 dan C2×3 , bisa juga dinyatakan dalam indeks masing-masing dari kedua matrik tersebut, yaitu "

d11 d12 d13 d21 d22 d23

#

=

"

a11 + c11 a12 + c12 a13 + c13 a21 + c21 a22 + c22 a23 + c23

#

1.4. OPERASI MATEMATIKA

7

Dijabarkan satu persatu sebagai berikut d11 = a11 + c11 d12 = a12 + c12 d13 = a13 + c13

(1.6)

d21 = a21 + c21 d22 = a22 + c22 d23 = a23 + c23 Dari sini dapat diturunkan sebuah rumus umum penjumlahan dua buah matrik dij = aij + cij

(1.7)

dimana i=1,2 dan j=1,2,3. 1.4.2 Komputasi penjumlahan matrik Berdasarkan contoh operasi penjumlahan di atas, indeks j pada persamaan (1.7) lebih cepat berubah dibanding indeks i sebagaimana ditulis pada persamaan (1.6), d11 = a11 + c11 d12 = a12 + c12 d13 = a13 + c13 Jelas terlihat, ketika indeks i masih bernilai 1, indeks j sudah berubah dari nilai 1 sampai 3. Hal ini membawa konsekuensi pada script pemrograman, dimana looping untuk indeks j harus diletakkan di dalam looping indeks i. Pokoknya yang looping-nya paling cepat harus diletakkan paling dalam; sebaliknya, looping paling luar adalah looping yang indeksnya paling jarang berubah. Dalam matlab, algoritma penjumlahan dua matrik ditulis sebagai berikut: 1 2 3 4 5

for i=1:2 for j=1:3 D(i,j)=A(i,j)+C(i,j); end end

Sedangkan dalam Fortran77, operasi penjumlahan antara matrik ditulis sebagai berikut: A2×3 dan C2×3 adalah 1 2 3 4 5

do i=1,2 do j=1,3 D(i,j)=A(i,j)+C(i,j) end do end do


8

Perhatikan kedua script di atas! Penulisan indeks i harus didahulukan daripada indeks j. Perlu dicatat bahwa ukuran matrik tidak terbatas hanya 2x3. Tentu saja anda bisa mengubah ukurannya sesuai dengan keperluan atau kebutuhan anda. Jika ukuran matrik dinyatakan secara umum sebagai n x m, dimana n adalah jumlah baris dan m adalah jumlah kolom, maka bentuk pernyataan komputasinya dalam matlab menjadi 1 2 3 4 5

for i=1:n for j=1:m D(i,j)=A(i,j)+C(i,j); end end

sedangkan dalam Fortran77 1 2 3 4 5

do i=1,n do j=1,m D(i,j)=A(i,j)+C(i,j) end do end do

Sekarang, mari kita lengkapi dengan contoh sebagai berikut: diketahui matrik A2×3 A=

"

# 3 8 5

C=

"

#

dan matrik C2×3

6 4 7

9 5 3 7 2 1

Program untuk menjumlahkan kedua matrik tersebut dalam matlab adalah: 1 2

clear all clc

3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

A(1,1) = 3; A(1,2) = 8; A(1,3) = 5; A(2,1) = 6; A(2,2) = 4; A(2,3) = 7; C(1,1) = 9; C(1,2) = 5; C(1,3) = 3; C(2,1) = 7; C(2,2) = 2; C(2,3) = 1; n=2 m=3 for i=1:n for j=1:m D(i,j)=A(i,j)+C(i,j); end end


9

sedangkan dalam Fortran77

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

A(1,1) = 3 A(1,2) = 8 A(1,3) = 5 A(2,1) = 6 A(2,2) = 4 A(2,3) = 7 C(1,1) = 9 C(1,2) = 5 C(1,3) = 3 C(2,1) = 7 C(2,2) = 2 C(2,3) = 1 n=2 m=3 do i=1,n do j=1,m D(i,j)=A(i,j)+C(i,j) end do end do

1.4.3 Perkalian matrik Operasi perkalian dua buah matrik hanya bisa dilakukan bila jumlah kolom matrik pertama sama dengan jumlah baris matrik kedua. Jadi kedua matrik tersebut tidak harus berukuran sama seperti pada penjumlahan dua matrik. Misalnya matrik A2×3 dikalikan dengan matrik B3×2 , lalu hasilnya (misalnya) dinamakan matrik E2×2 E2×2 = A2×3 .B3×2

E =

"

=

"

=

"

  # 1 3 3 8 5   5 9 6 4 7 2 4

# 3.1 + 8.5 + 5.2 3.3 + 8.9 + 5.4

6.1 + 4.5 + 7.2 6.3 + 4.9 + 7.4 # 53 101 40

82

Tanpa mempedulikan nilai elemen-elemen masing-masing matrik, operasi perkalian antara matrik A2×3 dan B3×2 , bisa juga dinyatakan dalam indeks masing-masing dari kedua matrik tersebut, yaitu "

e11 e12 e21 e22

#

=

"

a11 .b11 + a12 .b21 + a13 .b31 a11 .b12 + a12 .b22 + a13 .b32 a21 .b11 + a22 .b21 + a23 .b31 a21 .b12 + a22 .b22 + a23 .b32

#


10

Bila dijabarkan, maka elemen-elemen matrik E2×2 adalah e11 = a11 .b11 + a12 .b21 + a13 .b31

(1.8)

e12 = a11 .b12 + a12 .b22 + a13 .b32

(1.9)

e21 = a21 .b11 + a22 .b21 + a23 .b31

(1.10)

e22 = a21 .b12 + a22 .b22 + a23 .b32

(1.11)

Sejenak, mari kita amati perubahan pasangan angka-angka indeks yang mengiringi elemen e, a dan b pada persamaan (1.8) sampai persamaan (1.11). Perhatikan perubahan angka indeks pertama pada elemen e seperti berikut ini e1.. = .. e1.. = .. e2.. = .. e2.. = .. Pola perubahan yang sama akan kita dapati pada angka indeks pertama dari elemen a e1.. = a1.. .b... + a1.. .b... + a1.. .b... e1.. = a1.. .b... + a1.. .b... + a1.. .b... e2.. = a2.. .b... + a2.. .b... + a2.. .b... e2.. = a2.. .b... + a2.. .b... + a2.. .b... Dengan demikian kita bisa mencantumkan huruf i sebagai pengganti angka-angka indeks yang polanya sama ei.. = ai.. .b... + ai.. .b... + ai.. .b... ei.. = ai.. .b... + ai.. .b... + ai.. .b... ei.. = ai.. .b... + ai.. .b... + ai.. .b... ei.. = ai.. .b... + ai.. .b... + ai.. .b... dimana i bergerak mulai dari angka 1 hingga angka 2, atau kita nyatakan i=1,2. Selanjutnya, masih dari persamaan (1.8) sampai persamaan (1.11), marilah kita perhatikan perubahan angka indeks masih pada elemen e dan elemen b, ei1 = ai.. .b..1 + ai.. .b..1 + ai.. .b..1 ei2 = ai.. .b..2 + ai.. .b..2 + ai.. .b..2 ei1 = ai.. .b..1 + ai.. .b..1 + ai.. .b..1 ei2 = ai.. .b..2 + ai.. .b..2 + ai.. .b..2


11

Dengan demikian kita bisa mencantumkan huruf j sebagai pengganti angka-angka indeks yang polanya sama eij = ai.. .b..j + ai.. .b..j + ai.. .b..j eij = ai.. .b..j + ai.. .b..j + ai.. .b..j eij = ai.. .b..j + ai.. .b..j + ai.. .b..j eij = ai.. .b..j + ai.. .b..j + ai.. .b..j dimana j bergerak mulai dari angka 1 hingga angka 2, atau kita nyatakan j=1,2. Selanjutnya, masih dari persamaan (1.8) sampai persamaan (1.11), mari kita perhatikan perubahan angka indeks masih pada elemen a dan elemen b, dimana kita akan dapati pola sebagai berikut eij = ai1 .b1j + ai2 .b2j + ai3 .b3j eij = ai1 .b1j + ai2 .b2j + ai3 .b3j eij = ai1 .b1j + ai2 .b2j + ai3 .b3j eij = ai1 .b1j + ai2 .b2j + ai3 .b3j Dan kita bisa mencantumkan huruf k sebagai pengganti angka-angka indeks yang polanya sama, dimana k bergerak mulai dari angka 1 hingga angka 3, atau kita nyatakan k=1,2,3. eij = aik .bkj + aik .bkj + aik .bkj eij = aik .bkj + aik .bkj + aik .bkj eij = aik .bkj + aik .bkj + aik .bkj eij = aik .bkj + aik .bkj + aik .bkj Kemudian secara sederhana dapat ditulis sebagai berikut eij = aik .bkj + aik .bkj + aik .bkj

(1.12)

Selanjutnya dapat ditulis pula formula berikut eij =

3 X

aik bkj

(1.13)

k=1

dimana i=1,2; j=1,2; dan k=1,2,3. Berdasarkan contoh ini, maka secara umum bila ada matrik An×m yang dikalikan dengan matrik Bm×p , akan didapatkan matrik En×p dimana elemen-elemen matrik E memenuhi eij =

m X k=1

dengan i=1,2,. . . ,n; j=1,2. . . ,p; dan k=1,2. . . ,m.

aik bkj

(1.14)


12 1.4.4 Komputasi perkalian matrik

Komputasi operasi perkalian antara matrik A2×3 dan B3×2 dilakukan melalui 2 tahap; pertama adalah memberikan nilai 0 (nol) pada elemen-elemen matrik E2×2 dengan cara (dalam matlab) 1 2 3 4 5

for i=1:2 for j=1:2 E(i,j)=0.0; end end

dalam Fortran77 1 2 3 4 5

do i=1,2 do j=1,2 E(i,j)=0.0 end do end do

kedua adalah menghitung perkalian matrik dengan cara (dalam matlab) 1 2 3 4 5 6 7

for i=1:2 for j=1:2 for k=1:3 E(i,j)=E(i,j)+A(i,k)*B(k,j); end end end

dalam Fortran77 1 2 3 4 5 6 7

do i=1,2 do j=1,2 do k=1,3 E(i,j)=E(i,j)+A(i,k)*B(k,j) end do end do end do

Sebentar.., sebelum dilanjut tolong perhatikan penempatan indeks i, j dan k pada script di atas. Mengapa indeks i didahulukan daripada indeks j dan k? Ini bukan sesuatu yang kebetulan. Dan ini juga bukan sekedar mengikuti urutan huruf abjad i,j,k. Sekali lagi ingin saya tegaskan bahwa penempatan yang demikian semata-mata mengikuti aturan umum yaitu looping yang indeksnya berubah paling cepat harus diletakkan paling dalam; sebaliknya, looping paling luar adalah looping yang indeksnya paling jarang berubah. Kalau anda perhatikan dengan teliti, pasti anda akan menemukan fakta bahwa indeks k paling cepat berubah. Kemudian disusul oleh indeks j. Lalu yang paling jarang berubah adalah indeks i. Itulah sebabnya, penempatan urutan indeks pada script di atas harus dimulai dari i terlebih dahulu sebagai looping terluar, kemudian indeks j, dan yang terakhir indeks k sebagai looping terdalam.


13

Tentu saja anda bisa mengubah ukurannya sesuai dengan keperluan atau kebutuhan anda. Jika ukuran matrik A dinyatakan secara umum sebagai n x m dan matrik B berukuran m x p, maka bentuk pernyataan komputasinya dalam Matlab menjadi

1 2 3 4 5 6 7 8 9 10 11 12

for i=1:n for j=1:p E(i,j)=0.0; end end for i=1:n for j=1:p for k=1:m E(i,j)=E(i,j)+A(i,k)*B(k,j); end end end

dalam Fortran77

1 2 3 4 5 6 7 8 9 10 11 12

do i=1,n do j=1,p E(i,j)=0.0 end do end do do i=1,n do j=1,p do k=1,m E(i,j)=E(i,j)+A(i,k)*B(k,j) end do end do end do

dimana akan diperoleh hasil berupa matrik E yang berukuran n x p.

1.4.5 Perkalian matrik dan vektor-kolom Operasi perkalian antara matrik dan vektor-kolom sebenarnya sama saja dengan perkalian antara dua matrik. Hanya saja ukuran vektor-kolom boleh dibilang spesial yaitu m x 1, dimana m merupakan jumlah baris sementara jumlah kolomnya hanya satu. Misalnya matrik A, pada contoh 1, dikalikan dengan vektor-kolom x yang berukuran 3 x 1 atau disingkat dengan mengatakan vektor-kolom x berukuran 3, lalu hasilnya (misalnya) dinamakan vektor-kolom y y = Ax


14

y =

"

=

"

=

"

  # 2 3 8 5   3 6 4 7 4

3.2 + 8.3 + 5.4

6.2 + 4.3 + 7.4 # 50

#

52

Sekali lagi, tanpa mempedulikan nilai elemen-elemen masing-masing, operasi perkalian antara matrik A dan vektor-kolom x, bisa juga dinyatakan dalam indeksnya masing-masing, yaitu " # y1 y2

=

"

a11 .x1 + a12 .x2 + a13 .x3 a21 .x1 + a22 .x2 + a23 .x3

#

Bila dijabarkan, maka elemen-elemen vektor-kolom y adalah y1 = a11 .x1 + a12 .x2 + a13 .x3 y2 = a21 .x1 + a22 .x2 + a23 .x3 kemudian secara sederhana dapat diwakili oleh rumus berikut yi =

3 X

aij xj

j=1

dimana i=1,2. Berdasarkan contoh tersebut, secara umum bila ada matrik A berukuran n x m yang dikalikan dengan vektor-kolom x berukuran m, maka akan didapatkan vektor-kolom y berukuran n x 1 dimana elemen-elemen vektor-kolom y memenuhi yi =

m X

aij xj

(1.15)

j=1

dengan i=1,2,. . . ,n. 1.4.6 Komputasi perkalian matrik dan vektor-kolom Sama seperti perkalian dua matrik, komputasi untuk operasi perkalian antara matrik A berukuran n x m dan vektor-kolom x berukuran m dilakukan melalui 2 tahap; pertama adalah memberikan nilai 0 (nol) pada elemen-elemen vektor-kolom y yang berukuran n. Lalu tahap kedua adalah melakukan proses perkalian. Kedua tahapan ini digabung jadi satu dalam program berikut ini 1 2

for i=1:n b(i,1)=0.0;

1.5. PENUTUP 3 4 5 6 7 8

15

end for i=1:n for j=1:m b(i,1)=b(i,1)+A(i,j)*x(j,1); end end

dan dalam Fortran 1 2 3 4 5 6 7 8

do i=1,n b(i,1)=0.0 end do do i=1,n do j=1,m b(i,1)=b(i,1)+A(i,j)*x(j,1) end do end do

1.5

Penutup

Demikianlah catatan singkat dan sederhana mengenai jenis-jenis matrik dasar yang seringkali dijumpai dalam pengolahan data fisika secara numerik. Semuanya akan dijadikan acuan atau referensi pada pembahasan topik-topik numerik yang akan datang.

1.6

Latihan

Diketahui matrik A, matrik B, dan vektor x sebagai berikut 

1

3

−6

−2



  5  9 7 5.6  A= 2 4 8 −1    2.3 1.4 0.8 −2.3



8

1

4

21



  3  10 5 0.1  B= 7 −2 9 −5   2.7 −12 −8.9 5.7



0.4178



  −2.9587  x=  56.3069    8.1

1. Buatlah script untuk menyelesaikan penjumlahan matrik A dan matrik B. 2. Buatlah script untuk menyelesaikan perkalian matrik A dan matrik B. 3. Buatlah script untuk menyelesaikan perkalian matrik A dan vektor x. 4. Buatlah script untuk menyelesaikan perkalian matrik A dan vektor x.

============================================================

16


Bab 2

Metode Eliminasi Gauss

✍ Objektif : ⊲ Mengenalkan sistem persamaan linear. ⊲ Mengenalkan teknik triangularisasi dan substitusi mundur. ⊲ Aplikasi metode Eliminasi Gauss menggunakan matrik. ⊲ Membuat algoritma metode Eliminasi Gauss. ⊲ Menghitung invers matrik menggunakan metode Eliminasi Gauss.

2.1

Sistem persamaan linear

Secara umum, sistem persamaan linear dinyatakan sebagai berikut Pn :

an1 x1 + an2 x2 + ... + ann xn = bn

(2.1)

dimana a dan b merupakan konstanta, x adalah variable, n = 1, 2, 3, .... Contoh pertama Misalnya ada sistem persamaan linear yang terdiri dari empat buah persamaan yaitu P1 , P2 , P3 , dan P4 seperti berikut ini: P1 P2 P3 P4

: : : :

x1 2x1 3x1 −x1

+ + − +

x2 x2 x2 2x2

− − +

x3 x3 3x3

+ + + −

3x4 x4 2x4 x4

= = = =

4 1 -3 4

Problem dari sistem persamaan linear adalah bagaimana mencari nilai pengganti bagi variabel x1 , x2 , x3 , dan x4 sehingga semua persamaan diatas menjadi benar. Langkah awal penyelesaian problem tersebut adalah dengan melakukan penyederhanaan sistem persamaan linear. 17

BAB 2. METODE ELIMINASI GAUSS

18

2.2

Triangularisasi dan Substitusi Mundur

Ada banyak jalan untuk mendapatkan bentuk yang lebih sederhana, namun masalahnya, kita ingin mendapatkan sebuah algoritma program yang nantinya bisa berjalan di komputer, sedemikian rupa sehingga apapun persamaannya, bisa disederhanakan oleh komputer. Kita akan berpatokan pada tiga buah aturan operasi untuk menyederhanakan sistem persamaan linear di atas, yaitu • Persamaan Pi dapat dikalikan dengan sembarang konstanta λ, lalu hasilnya ditempatkan di posisi persamaan Pi . Simbol operasi ini adalah (λPi ) → (Pi ).

• Persamaan Pj dapat dikalikan dengan sembarang konstanta λ kemudian dijumlahkan dengan persamaan Pi , lalu hasilnya ditempatkan di posisi persamaan Pi . Simbol operasi ini adalah (Pi + λPj ) → (Pi ). • Persamaan Pi dan Pj dapat bertukar posisi. Simbol operasi ini adalah (Pi ) ↔ (Pj ). Maka dengan berpegang pada aturan-aturan tersebut, problem sistem persamaan linear di atas akan diselesaikan dengan langkah-langkah berikut ini: 1. Gunakan persamaan P1 untuk menghilangkan variabel x1 dari persamaan P2 , P3 dan P4 dengan cara (P2 − 2P1 ) → (P2 ), (P3 − 3P1 ) → (P3 ) dan (P4 + P1 ) → (P4 ). Hasilnya akan

seperti ini

P1 :

x1 + x2 + 3x4 = 4,

P2 :

−x2 − x3 − 5x4 = −7,

P3 :

−4x2 − x3 − 7x4 = −15,

P4 :

3x2 + 3x3 + 2x4 = 8

2. Gunakan persamaan P2 untuk menghilangkan variabel x2 dari persamaan P3 dan P4 dengan cara (P3 − 4P2 ) → (P3 ) dan (P4 + 3P2 ) → (P4 ). Hasilnya akan seperti ini P1 :

x1 + x2 + 3x4 = 4,

P2 :

−x2 − x3 − 5x4 = −7,

P3 :

3x3 + 13x4 = 13,

P4 :

−13x4 = −13

Kalau x3 masih ada di persamaan P4 , dibutuhkan satu operasi lagi untuk menghilangkannya. Namun hasil operasi pada langkah ke-2 ternyata sudah otomatis menghilangkan x3 . Bentuk akhir dari keempat persamaan di atas, dikenal sebagai bentuk triangular. Sampai dengan langkah ke-2 ini, kita berhasil mendapatkan sistem persamaan linear yang lebih sederhana. Apa yang dimaksud dengan sederhana dalam konteks ini? Suatu sistem persamaan linear dikatakan sederhana bila kita bisa mendapatkan seluruh nilai pengganti variabelnya dengan cara yang lebih mudah atau dengan usaha yang tidak

2.2. TRIANGULARISASI DAN SUBSTITUSI MUNDUR

19

memakan waktu lama dibandingkan sebelum disederhanakan. Sekali kita mendapatkan nilai pengganti bagi variabel x4 , maka x3 , x2 dan x1 akan diperoleh dengan mudah dan cepat, sebagaimana yang dijelaskan pada langkah berikutnya. 3. Selanjutnya kita jalankan proses backward-substitution. Melalui proses ini, yang pertama kali didapat adalah nilai pengganti bagi variabel x4 , kemudian x3 , lalu diikuti x2 , dan akhirnya x1 . P4 : P3 : P2 : P1 :

x4 =

−13 −13

= 1,

1 1 x3 = (13 − 13x4 ) = (13 − 13) = 0, 3 3 x2 = −(−7 + 5x4 + x3 ) = −(−7 + 5 + 0) = 2, x1 = 4 − 3x4 − x2 = 4 − 3 − 2 = −1

Jadi solusinya adalah x1 = −1, x2 = 2, x3 = 0 dan x4 = 1. Coba sekarang anda cek,

apakah semua solusi ini cocok dan tepat bila dimasukan ke sistem persamaan linear yang pertama, yaitu yang belum disederhanakan?

OK, mudah-mudahan ngerti ya... Kalau belum paham, coba diulangi bacanya sekali lagi. Atau, sekarang kita beralih kecontoh yang lain.


20 Contoh kedua

Misalnya ada sistem persamaan linear, terdiri dari empat buah persamaan yaitu P1 , P2 , P3 , dan P4 seperti berikut ini: P1 P2 P3 P4

: : : :

x1 2x1 x1 x1

− − + −

x2 2x2 x2 x2

+ + + +

2x3 3x3 x3 4x3

− −

x4 3x4

+

3x4

= = = =

-8 -20 -2 4

Seperti contoh pertama, solusi sistem persamaan linear di atas akan dicari dengan langkahlangkah berikut ini: 1. Gunakan persamaan P1 untuk menghilangkan x1 dari persamaan P2 , P3 dan P4 dengan cara (P2 − 2P1 ) → (P2 ), (P3 − P1 ) → (P3 ) dan (P4 − P1 ) → (P4 ). Hasilnya akan seperti ini P1 :

x1 − x2 + 2x3 − x4 = −8,

P2 :

−x3 − x4 = −4,

P3 :

2x2 − x3 + x4 = 6,

P4 :

2x3 + 4x4 = 12

Perhatikan persamaan P2 ! Akibat dari langkah yang pertama tadi, x2 hilang dari persamaan P2 . Kondisi ini bisa menggagalkan proses triangularisasi. Untuk itu, posisi P2 mesti ditukar dengan persamaan yang berada dibawahnya, yaitu P3 atau P4 . Supaya proses triangularisasi dilanjutkan kembali, maka yang paling cocok adalah ditukar dengan P3 . 2. Tukar posisi persamaan P2 dengan persamaan P3 , (P2 ↔ P3 ). Hasilnya akan seperti ini P1 :

x1 − x2 + 2x3 − x4 = −8,

P2 :

2x2 − x3 + x4 = 6,

P3 :

−x3 − x4 = −4,

P4 :

2x3 + 4x4 = 12

3. Gunakan persamaan P3 untuk menghilangkan x3 dari persamaan P4 dengan cara (P4 − 2P3 ) → (P4 ). Hasilnya akan seperti ini P1 : P2 : P3 :

x1 − x2 + 2x3 − x4 = −8, 2x2 − x3 + x4 = 6, −x3 − x4 = −4,

P4 : Sampai disini proses triangularisasi telah selesai.

2x4 = 4

2.3. MATRIK DAN ELIMINASI GAUSS

21

4. Selanjutnya adalah proses backward-substitution. Melalui proses ini, yang pertama kali didapat solusinya adalah x4 , kemudian x3 , lalu diikuti x2 , dan akhirnya x1 . P4 : P3 : P2 : P1 :

4 2 −4 + x4 x3 = −1 6 + x3 − x4 x2 = 2 x1 = −8 + x2 − 2x3 + x4 x4 =

= 2, = 2, = 3, = −7

Jadi solusinya adalah x1 = −7, x2 = 3, x3 = 2 dan x4 = 2.

Berdasarkan kedua contoh di atas, untuk mendapatkan solusi sistem persamaan linear, diperlukan operasi triangularisasi dan proses backward-substitution. Kata backward-substitution kalau diterjemahkan kedalam bahasa indonesia, menjadi substitusi-mundur. Gabungan proses triangularisasi dan substitusi-mundur untuk menyelesaikan sistem persamaan linear dikenal sebagai metode eliminasi gauss.

2.3

Matrik dan Eliminasi Gauss

Sejumlah matrik bisa digunakan untuk menyatakan suatu sistem persamaan linear. Sejenak, mari kita kembali lagi melihat sistem persamaan linear secara umum seperti berikut ini: a11 x1 + a12 x2 + . . . + a1n xn = b1 a21 x1 + a22 x2 + . . . + a2n xn = b2 ............... = ... ............... = ... an1 x1 + an2 x2 + . . . + ann xn = bn Sementara, kalau dinyatakan dalam bentuk operasi matrik, maka akan seperti ini: 

a11

a12

. . . a1n

  a21 a22 . . . a2n  . .. ..  . . .  . an1 an2 . . . ann

     

x1 x2 .. . xn





    =    

b1 b2 .. . bn

     

(2.2)

Dalam mencari solusi suatu sistem persamaan linear dengan metode eliminasi gauss, bentuk operasi matrik di atas dimanipulasi menjadi matrik augment, yaitu suatu matrik yang beruku-


22 ran n x (n + 1) seperti berikut ini: 

a11

a12

. . . a1n

| b1

  a21 a22 . . . a2n | b2  . .. .. .  . . . | ..  . an1 an2 . . . ann | bn





a11

a12

. . . a1n

| a1,n+1

    a21 a22 . . . a2n | a2,n+1 = . .. .. ..   . . . | .   . an1 an2 . . . ann | an,n+1

     

(2.3)

Berdasarkan contoh pertama yang ada dihalaman depan catatan ini, saya akan tunjukkan proses triangularisasi dan substitusi-mundur dalam operasi matrik terhadap sistem persamaan linear yang terdiri dari empat persamaan matematika, yaitu (silakan lihat kembali contoh pertama): 

1

1

0

3

  2 1 −1 1   3 −1 −1 2  −1 2 3 −1



x1





4



      x2   1        x  =  −3   3    x4 4

Lalu kita dapat membuat matrik augment sebagai berikut: 

1

1

0

3

|

4

  2 1 −1 1 | 1   3 −1 −1 2 | −3  −1 2 3 −1 | 4

     

Kemudian kita lakukan operasi triangularisai terhadap matrik augment, dimulai dari kolom pertama, yaitu 

1

1

0

3

|

4

|

4



   0 −1 −1 −5 | −7     0 −4 −1 −7 | −15    0 3 3 2 | 8 lalu dilanjutkan ke kolom berikutnya 

1

1

0

3

  0 −1 −1 −5 | −7   0 0 3 13 | 13  0 0 0 −13 | −13

     

Sebelum dilanjutkan ke substitusi-mundur, saya ingin menegaskan peranan angka-angka indeks dari masing-masing elemen matrik augment tersebut. Silakan perhatikan posisi masing-

2.4. ALGORITMA ELIMINASI GAUSS

23

masing elemen berikut ini: 

1

1

0

3

|



4

  0 −1 −1 −5 | −7   0 0 3 13 | 13  0 0 0 −13 | −13



a11 a12 a13 a14 | a15

     →  a21 a22 a23 a24 | a25   a   31 a32 a33 a34 | a35 a41 a42 a43 a44 | a45

     

Dengan memperhatikan angka-angka indeks pada matrik augment di atas, kita akan mencoba membuat rumusan proses substitusi-mundur untuk mendapatkan seluruh nilai pengganti variabel x. Dimulai dari x4 , x4 =

−13 a45 = =1 a44 −13

ini dapat dinyatakan dalam rumus umum, yaitu xn =

an,n+1 ann

lalu dilanjutkan dengan x3 , x2 , dan x1 . a35 − a34 x4 a33 a25 − (a23 x3 + a24 x4 ) x2 = a22 a15 − (a12 x2 + a13 x3 + a14 x4 ) x1 = a11 x3 =

= = =

13 − [(13)(1)] =0 3 (−7) − [(−1)(0) + (−5)(1)] =2 (−1) 4 − [(1)(2) + (0)(0) + (3)(1)] = −1 1

ini juga dapat dinyatakan dalam rumus umum yaitu: xi =

ai,n+1 −

Pn

j=i+1 aij xj

aii

Proses triangularisasi dan substitusi-mundur dibakukan menjadi algoritma metode eliminasi gauss yang dapat diterapkan dalam berbagai bahasa pemrograman komputer, misalnya fortran, C, java, pascal, matlab, dan lain-lain.

2.4

Algoritma eliminasi Gauss

Secara umum, sistem persamaan linear adalah sebagai berikut: a11 x1 + a12 x2 + . . . + a1n xn = b1 a21 x1 + a22 x2 + . . . + a2n xn = b2 .. .. . . . = .. an1 x1 + an2 x2 + . . . + ann xn = bn Algoritma dasar metode eliminasi gauss, adalah sebagai berikut: 1. Ubahlah sistem persamaan linear tersebut menjadi matrik augment, yaitu suatu matrik


24 yang berukuran n x (n + 1) seperti berikut ini: 

a11

a12

. . . a1n

| b1

  a21 a22 . . . a2n | b2  . .. .. .  . . . | ..  . an1 an2 . . . ann | bn





a11

a12

. . . a1n | a1,n+1


     

(2.4)

Jelas terlihat bahwa elemen-elemen yang menempati kolom terakhir matrik augment adalah nilai dari bi ; yaitu ai,n+1 = bi dimana i = 1, 2, ..., n. 2. Periksalah elemen-elemen pivot. Apakah ada yang bernilai nol? Elemen-elemen pivot adalah elemen-elemen yang menempati diagonal suatu matrik, yaitu a11 , a22 , ..., ann atau disingkat aii . Jika aii 6= 0, bisa dilanjutkan ke langkah no.3. Namun, jika ada elemen

diagonal yang bernilai nol, aii = 0, maka baris dimana elemen itu berada harus ditukar posisinya dengan baris yang ada dibawahnya, (Pi ) ↔ (Pj ) dimana j = i + 1, i + 2, ..., n,

sampai elemen diagonal matrik menjadi tidak nol, aii 6= 0. (Kalau kurang jelas, silakan lihat

lagi contoh kedua yang ada dihalaman 3. Sebaiknya, walaupun elemen diagonalnya tidak nol, namun mendekati nol (misalnya 0,03), maka proses pertukaran ini dilakukan juga).

3. Proses triangularisasi. Lakukanlah operasi berikut: Pj −

aji Pi → Pj aii

(2.5)

dimana j = i + 1, i + 2, ..., n. Maka matrik augment akan menjadi: 

a11 a12 a13 . . . a1n

  0   0   .  ..  0

a22 a23 . . . a2n 0 .. .

a33 . . . a3n .. .. .. . . .

0

0

0

4. Hitunglah nilai xn dengan cara: xn =

ann

| a1,n+1



 | a2,n+1   | a3,n+1    ..  . |  | an,n+1

an,n+1 ann

(2.6)

(2.7)

5. Lakukanlah proses substitusi-mundur untuk memperoleh xn−1 , xn−2 , ..., x2 , x1 dengan cara: xi =

ai,n+1 −

Pn

j=i+1 aij xj

aii

(2.8)

dimana i = n − 1, n − 2, ..., 2, 1. Demikianlan algoritma dasar metode eliminasi gauss. Selanjutnya algoritma dasar tersebut perlu dirinci lagi sebelum dapat diterjemahkan kedalam bahasa pemrograman komputer.

2.4. ALGORITMA ELIMINASI GAUSS

25

2.4.1 Algoritma Algoritma metode eliminasi gauss untuk menyelesaikan n x n sistem persamaan linear. P1 :

a11 x1 + a12 x2 + . . . + a1n xn = b1

P2 : .. .

a21 x1 + a22 x2 + . . . + a2n xn = b2 .. .. . . . = ..

Pn :

an1 x1 + an2 x2 + . . . + ann xn = bn

INPUT: sejumlah persamaan linear dimana konstanta-konstanta-nya menjadi elemen-elemen matrik augment A = (aij ), dengan 1 ≤ i ≤ n dan 1 ≤ j ≤ n + 1.

OUTPUT: solusi x1 , x2 , x3 , ..., xn atau pesan kesalahan yang mengatakan bahwa sistem per-

samaan linear tidak memiliki solusi yang unik. • Langkah 1: Inputkan konstanta-konstanta dari sistem persamaan linear kedalam elemenelemen matrik augment, yaitu suatu matrik yang berukuran n x (n + 1) seperti berikut ini: 

a11

a12

. . . a1n

| b1

  a21 a22 . . . a2n | b2  . .. .. .  . . . | ..  . an1 an2 . . . ann | bn





a11

a12

. . . a1n | a1,n+1


     

(2.9)

• Langkah 2: Untuk i = 1, ..., n − 1, lakukan Langkah 3 sampai Langkah 5. • Langkah 3: Definisikan p sebagai integer dimana i ≤ p ≤ n. Lalu pastikan bahwa

api 6= 0. Jika ada elemen diagonal yang bernilai nol (aii = 0), maka program harus

mencari dan memeriksa elemen-elemen yang tidak bernilai nol dalam kolom yang

sama dengan kolom tempat elemen diagonal tersebut berada. Jadi saat proses ini berlangsung, integer i (indeks dari kolom) dibuat konstan, sementara integer p (indeks dari baris) bergerak dari p = i sampai p = n. Bila ternyata setelah mencapai elemen paling bawah dalam kolom tersebut, yaitu saat p = n tetap didapat nilai api = 0, maka sebuah pesan dimunculkan: sistem persamaan linear tidak memiliki

solusi yang unik. Lalu program berakhir: STOP. • Langkah 4: Namun jika sebelum integer p mencapai nilai p = n sudah diperoleh elemen yang tidak nol (api 6= 0), maka bisa dipastikan p 6= i. Jika p 6= i maka

lakukan proses pertukaran (Pp ) ↔ (Pi ).

• Langkah 5: Untuk j = i + 1, .., n, lakukan Langkah 6 dan Langkah 7. • Langkah 6: Tentukan mji ,

mji =

aji aii


26

• Langkah 7: Lakukan proses triangularisasi, (Pj − mji Pi ) → (Pj ) • Langkah 8: Setelah proses triangularisasi dilalui, periksalah ann . Jika ann = 0, kirimkan pesan: sistem persamaan linear tidak memiliki solusi yang unik. Lalu program berakhir:

STOP. • Langkah 9: Jika ann 6= 0, lakukan proses substitusi mundur, dimulai dengan menentukan xn ,

xn =

an,n+1 ann

• Langkah 10: Untuk i = n − 1, ..., 1 tentukan xi , xi =

ai,n+1 −

Pn

j=i+1 aij xj

aii

• Langkah 11: Diperoleh solusi yaitu x1 , x2 , ..., xn . Algoritma telah dijalankan dengan sukses. STOP.

Saya telah membuat program sederhana dalam fortran untuk mewujudkan algoritma eliminasi gauss. Saya berasumsi bahwa anda sudah menguasai dasar-dasar pemrograman dalam fortran. Program ini sudah dicoba di-compile dengan fortran77 under Linux Debian dan visual-fortran under windows-XP. Langkah-langkah yang tercantum pada program ini disesuaikan dengan langkah-langkah yang tertulis di atas. Dalam program ini, ukuran maksimum matrik augment adalah 10 x 11, untuk mencari 10 variabel yang tidak diketahui. Jika anda bermaksud memperbesar atau memperkecil ukuran matrik augment, silakan sesuaikan angka ukuran matrik yang anda inginkan pada statemen pertama dari program ini, yaitu statemen DIMENSION. Inilah programnya, 1 2 3 4 5

C

6 7 8 9 10 11 12 13 14 15 16

60 50

17 18

C

DIMENSION A(10,11), X(10) REAL MJI WRITE (*,*) ’=PROGRAM ELIMINASI GAUSS=’ WRITE (*,*) LANGKAH 1: MEMASUKAN NILAI ELEMEN-ELEMEN MATRIK AUGMENT WRITE (*,’(1X,A)’) ’JUMLAH PERSAMAAN ? ’ READ (*,*) N WRITE (*,*) WRITE (*,*) ’MASUKAN ELEMEN-ELEMEN MATRIK AUGMENT’ M = N + 1 DO 50 I = 1,N DO 60 J = 1,M WRITE (*,’(1X,A,I2,A,I2,A)’) ’A(’,I,’,’,J,’) = ’ READ (*,*) A(I,J) CONTINUE CONTINUE WRITE (*,*) MENAMPILKAN MATRIK AUGMENT

2.4. ALGORITMA ELIMINASI GAUSS 19 20 21 22

110

23 24

C

25 26 27

C

28 29

100

30 31 32 33

200 C

34 35 36 37

C

38 39 40 41 42 43

20

44 45

C

46 47 48

C

49 50

C

51 52 53

40

54 55 56 57

30 10 C

58 59 60 61 62

120 C

63 64

C

65 66 67 68

C

69 70 71 72 73 74 75 76 77

C

WRITE (*,’(1X,A)’) ’MATRIK AUGMENT:’ DO 110 I = 1,N WRITE (*,’(1X,5(F14.8))’) (A(I,J),J=1,M) CONTINUE WRITE (*,*) LANGKAH 2: MEMERIKSA ELEMEN-ELEMEN PIVOT DAN PROSES TUKAR POSISI NN = N-1 DO 10 I=1,NN LANGKAH 3: MENDEFINISIKAN P P = I IF (ABS(A(P,I)).GE.1.0E-20 .OR. P.GT.N) GOTO 200 P = P+1 GOTO 100 IF(P.EQ.N+1)THEN MENAMPILKAN PESAN TIDAK UNIK WRITE(*,5) GOTO 400 END IF LANGKAH 4: PROSES TUKAR POSISI IF(P.NE.I) THEN DO 20 JJ=1,M C = A(I,JJ) A(I,JJ) = A(P,JJ) A(P,JJ) = C CONTINUE END IF LANGKAH 5: PERSIAPAN PROSES TRIANGULARISASI JJ = I+1 DO 30 J=JJ,N LANGKAH 6: TENTUKAN MJI MJI = A(J,I)/A(I,I) LANGKAH 7: MELAKUKAN PROSES TRIANGULARISASI DO 40 K=JJ,M A(J,K) = A(J,K)-MJI*A(I,K) CONTINUE A(J,I) = 0 CONTINUE CONTINUE MENAMPILKAN HASIL TRIANGULARISASI WRITE (*,’(1X,A)’) ’HASIL TRIANGULARISASI:’ DO 120 I = 1,N WRITE (*,’(1X,5(F14.8))’) (A(I,J),J=1,M) CONTINUE LANGKAH 8: MEMERIKSA ELEMEN A(N,N) IF(ABS(A(N,N)).LT.1.0E-20) THEN MENAMPILKAN PESAN TIDAK UNIK WRITE(*,5) GOTO 400 END IF LANGKAH 9: MENGHITUNG X(N) X(N) = A(N,N+1)/A(N,N) LANGKAH 10: PROSES SUBSTITUSI MUNDUR L = N-1 DO 15 K=1,L I = L-K+1 JJ = I+1 SUM = 0.0 DO 16 KK=JJ,N SUM = SUM+A(I,KK)*X(KK)

27


28 78

16

79 80 81

15 C

82 83 84 85 86 87

18 400

CONTINUE X(I) = (A(I,N+1)-SUM)/A(I,I) CONTINUE LANGKAH 11: MENAMPILKAN HASIL PERHITUNGAN WRITE (*,*) WRITE (*,7) DO 18 I = 1,N WRITE (*,’(1X,A,I2,A,F14.8)’) ’X(’,I,’) = ’,X(I) CONTINUE STOP

88 89 90 91

5 7

FORMAT(1X,’SISTEM LINEAR TIDAK MEMILIKI SOLUSI YANG UNIK’) FORMAT(1X,’SOLUSI UNIK’) END

Script eliminasi gauss dalam matlab juga telah dibuat. Namun dalam anda perlu memodifikasi elemen-elemen matrik A agar sesuai dengan data yang hendak anda olah. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

clear all clc A(1,1)=1; A(1,2)=1; A(1,3)=-1; A(1,4)=0; A(2,1)=6; A(2,2)=-4; A(2,3)=0; A(2,4)=24; A(3,1)=6; A(3,2)=0; A(3,3)=2; A(3,4)=10; A n=3 %jumlah persamaan pause

18 19 20

%========== Proses Triangularisasi ========= for j=1:(n-1)

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

%----mulai proses pivot--if (A(j,j)==0) for p=1:n+1 u=A(j,p); v=A(j+1,p); A(j+1,p)=u; A(j,p)=v; end end %----akhir proses pivot--jj=j+1; for i=jj:n m=A(i,j)/A(j,j); for k=1:(n+1) A(i,k)=A(i,k)-(m*A(j,k)); end end end A

2.5. CONTOH APLIKASI 41 42

29

pause %========= Akhir Proses Triangularisasi ===

43 44 45

%------Proses Substitusi mundur------------x(n,1)=A(n,n+1)/A(n,n);

46 47 48 49 50 51 52 53 54

for i=n-1:-1:1 S=0; for j=n:-1:i+1 S=S+A(i,j)*x(j,1); end x(i,1)=(A(i,n+1)-S)/A(i,i); end x

2.5

Contoh aplikasi

2.5.1 Menghitung arus listrik Gunakan metode Eliminasi Gauss untuk menentukan arus i1 , i2 dan i3 yang mengalir pada rangkaian berikut ini

jawab: Berdasarkan Hukum Kirchhoff: I1 + I2 = I3 10 − 6I1 − 2I3 = 0 −14 + 6I1 − 10 − 4I2 = 0 Lalu kita susun ulang ketiga persamaan di atas menjadi seperti ini: I1 + I2 − I3 = 0 6I1 + 2I3 = 10 6I1 − 4I2 = 24


30 Kemudian dinyatakan dalam bentuk matriks: 

1

1

−1



I1





0



    0   I2  =  24  2 I3 10

  6 −4 6 0

Selanjutkan kita susun matriks augmentasi sebagai berikut: 

1

1

  6 −4 6 0

−1 0

2

0



 24  10

Langkah berikutnya adalah menghitung matriks triangularisasi dengan langkah-langkah sebagai berikut: 6 a21 = = 6 a11 1 = 6 − (6).(1) = 0

m= a21 = a21 − m.a11

a22 = a22 − m.a12 = −4 − (6).(1) = −10 a23 = a23 − m.a13 = 0 − (6).(−1) = 6 a24 = a24 − m.a14 = 24 − (6).(0) = 24 6 a31 = = 6 a11 1 = 6 − (6).(1) = 0

m= a31 = a31 − m.a11

a32 = a32 − m.a12 = 0 − (6).(1) = −6 a33 = a33 − m.a13 = 2 − (6).(−1) = 8 a34 = a34 − m.a14 = 10 − (6).(0) = 10 Sampai disini matriks augment mengalami perubahan menjadi 

1

1

  0 −10 0 −6

−1 6

8

0



 24  10

2.6. MENGHITUNG INVERS MATRIK

31

Kelanjutan langkah menuju triangularisasi adalah m=

a32 a22

−6 ).(0) −10 −6 ).(−10) a32 = a32 − m.a22 = −6 − ( −10 −6 a33 = a33 − m.a23 = 8 − ( ).(6) −10 −6 a34 = a34 − m.a24 = 10 − ( ).(24) −10 a31 = a31 − m.a21 = 0 − (

=

−6 −10

= 0 = 0 = 4, 4 = −4, 4

maka matriks triangularisasi berhasil didapat yaitu 

1

1

−1

0



  24   0 −10 6 0 0 4, 4 −4, 4 Sekarang tinggal melakukan proses substitusi mundur a34 a33 a24 − a23 .I3 I2 = a22 a14 − (a13 .I3 + a12 .I2 ) I1 = a11 I3 =

= = =

−4, 4 = −1 4, 4 24 − (6).(−1) = −3 −10 (0 − [(−1).(−1) + (1).(−3)] =2 1

Dengan demikian, besar masing-masing arus pada rangkaian di atas adalah I1 = 2A, I2 = −3A dan I3 = −1A. Tanda minus (-) memiliki arti bahwa arah arus yang sesungguhnya berlawanan

arah dengan asumsi awal yang kita gunakan.

2.6

Menghitung invers matrik

Sekali lagi saya ulangi apa yang pernah kita bahas di awal bab ini yaitu bahwa sistem persamaan linear dapat dinyatakan sebagai berikut: a11 x1 + a12 x2 + . . . + a1n xn = b1 a21 x1 + a22 x2 + . . . + a2n xn = b2 ............... = ... ............... = ... an1 x1 + an2 x2 + . . . + ann xn = bn Sistem persamaan linear tersebut dapat dinyatakan dalam bentuk operasi matrik, Ax = b

(2.10)


32 sehingga bentuknya menjadi seperti ini: 

a11

a12



. . . a1n

  a21 a22 . . . a2n  . .. ..  . . .  . an1 an2 . . . ann

    



x1 x2 .. .



    =    

xn

b1 b2 .. . bn

     

dimana 

a11

a12



. . . a1n

  a21 a22 . . . a2n A= .. ..  .. . .  . an1 an2 . . . ann



  ,  

x1 x2 .. .

  x=  

xn





  ,  

  b=  

b1 b2 .. . bn

     

Dalam kaitannya dengan invers matrik, matrik A disebut matrik non-singular jika matrik A memiliki matrik invers dirinya yaitu A−1 . Atau dengan kata lain, matrik A−1 adalah invers dari matrik A. Jika matrik A tidak memiliki invers, maka matrik A disebut singular. Bila matrik A dikalikan dengan matrik A−1 maka akan menghasilkan matrik identitas I, yaitu suatu matrik yang elemen-elemen diagonalnya bernilai 1.

−1

AA



  =I=  

1 0 ... 0



 0 1 ... 0  .. .. . . ..   . .  . . 0 0 ... 1

(2.11)

Misalnya diketahui, 

 A=

1 2 −1



 0 , 2

2 1 −1 1



 A−1 = 

− 92

4 9 − 13

5 9 − 19 1 3

− 91

2 9 1 3

  

Bila keduanya dikalikan, maka akan menghasilkan matrik identitas, 

 AA−1 = 

1 2 −1



− 29

 0   49 − 13 2

2 1 −1 1

5 9 1 −9 1 3

− 91 2 9 1 3





1 0 0



   = 0 1 0  0 0 1

Lalu bagaimana cara mendapatkan matrik invers, A−1 ? Persamaan (2.11) bisa dijadikan pedoman.. AA−1 = I   

1 2 −1

2 1 −1 1



i11 i12 i13





1 0 0



    0   i21 i22 i23  =  0 1 0  2 i31 i32 i33 0 0 1

2.6. MENGHITUNG INVERS MATRIK

33

dalam hal ini matrik A−1 adalah 

i11 i12 i13



  A−1 =  i21 i22 i23  i31 i32 i33 Elemen-elemen matrik invers, A−1 dapat diperoleh dengan menerapkan metode eliminasi gauss. Diawali dengan membentuk matrik augment:   

1 2 −1 | 1 0 0

2 1



 0 | 0 1 0  2 | 0 0 1

−1 1

Lalu dilanjutkan dengan proses triangularisasi: (P2 −2P1 )→(P2 ) dan (P3 +P1 )→(P3 ), kemudian diikuti oleh (P3 + P2 )→(P3 ): 

1

2 −1 |

  0 −3 0 3

1 0 0





→

 2 | −2 1 0  1 | 1 0 1

1

2 −1 |

  0 −3 0 0

1 0 0



 2 | −2 1 0  3 | −1 1 1

Langkah berikutnya, matrik augment yang telah mengalami triangularisasi tersebut dipecah menjadi tiga buah matrik augment seperti berikut ini: 

1

2 −1 |

  0 −3 0 0

1



 2 | −2  3 | −1



1

2 −1 | 0

  0 −3 0 0



 2 | 1  3 | 1



1

2 −1 | 0

  0 −3 0 0



 2 | 0  3 | 1

Langkah pamungkasnya adalah melakukan proses substitusi mundur pada ketiga matrik augment di atas, sehingga diperoleh: i11 = − 29 i21 =

4 9

i31 = − 31

i12 =

5 9

i13 = − 19

i22 = − 19 i32 =

i23 =

1 3

i33 =

2 9 1 3

Hasil tersebut digabung menjadi sebuah matrik, yaitu matrik A−1 , 

 A−1 = 

− 92

4 9 − 13

5 9 − 19 1 3

− 91 2 9 1 3

  

Keberadaan matrik A−1 bisa digunakan untuk menyelesaikan sistem persamaan linear


34 (mencari nilai x ), dengan cara sebagai berikut Ax = b

A−1 Ax = A−1 b Ix = A−1 b x = A−1 b

(2.12)

Contoh berikut ini akan menjelaskan prosesnya secara lebih rinci. Misalnya diketahui sistem persamaan linear x1 + 2x2 − x3 = 2 2x1 + x2 = 3 −x1 + x2 + 2x3 = 4 Bila dikonversikan kedalam operasi matrik menjadi   

1 2 −1

2 1 −1 1



x1





2



    0   x2  =  3  2 x3 4

Berdasarkan persamaan (2.12), maka elemen-elemen vektor x dapat dicari dengan cara x = A−1 b 

 x=

− 29 4 9 − 13

5 9 − 91 1 3

− 19 2 9 1 3



2





    3  =  4

7 9 13 9 5 3

  

Akhirnya diperoleh solusi x1 = 7/9, x2 = 13/9, dan x3 = 5/3. Penyelesaian sistem persamaan linear menjadi lebih mudah bila matrik A−1 sudah diketahui. Sayangnya, untuk mendapatkan matrik A−1 , diperlukan langkah-langkah, seperti yang sudah dibahas pada contoh pertama di atas, yang berakibat in-efisiensi proses penyelesaian (secara komputasi) bila dibandingkan dengan metode eliminasi gauss untuk memecahkan sistem persamaan linear. Namun bagaimanapun, secara konseptual kita dianjurkan mengetahui cara bagaimana mendapatkan matrik A−1 . Saya telah memodifikasi program eliminasi gauss yang terdahulu, untuk keperluan perhitungan matrik invers. Program ini ditulis dengan bahasa fortran, sudah berhasil dikompilasi dalam Linux Debian (g77) dan Windows XP (Visual Fortran). Inilah programnya, 1 2 3 4 5

DIMENSION A(10,20), D(10,10), X(10) REAL MJI INTEGER TKR, BK, TK, Q WRITE (*,*) ’=PROGRAM INVERS MATRIK DENGAN ELIMINASI GAUSS=’ WRITE (*,*)

2.6. MENGHITUNG INVERS MATRIK 6

C

7 8 9 10 11 12 13 14 15 16 17 18

60 50 C

19 20 21 22 23 24 25 26 27

80 70

28 29

C

30 31 32 33

110

34 35

C

36 37

C

38 39

C

40 41

C

42 43 44

C

45 46

100

47 48 49 50

200 C

51 52 53 54

C

55 56 57 58 59 60 61

20

62 63 64

C

LANGKAH 1: MEMASUKAN NILAI ELEMEN-ELEMEN MATRIK A WRITE (*,’(1X,A)’) ’JUMLAH PERSAMAAN ? ’ READ (*,*) N WRITE (*,*) WRITE (*,*) ’MASUKAN ELEMEN-ELEMEN MATRIK A’ M = N + 1 DO 50 I = 1,N DO 60 J = 1,N WRITE (*,’(1X,A,I2,A,I2,A)’) ’A(’,I,’,’,J,’) = ’ READ (*,*) A(I,J) CONTINUE CONTINUE LANGKAH 2: MENDEFINISIKAN MATRIK IDENTITAS WRITE (*,*) ’MENDEFINISIKAN MATRIK IDENTITAS’ DO 70 I = 1,N DO 80 J = M,N+N A(I,J) = 0 IF (I+N .EQ. J) THEN A(I,J) = 1 END IF CONTINUE CONTINUE WRITE (*,*) MENAMPILKAN MATRIK AUGMENT WRITE (*,’(1X,A)’) ’MATRIK AUGMENT:’ DO 110 I = 1,N WRITE (*,’(1X,5(F14.8))’) (A(I,J),J=1,N+N) CONTINUE WRITE (*,*) MENGHITUNG JUMLAH TUKAR (TKR) POSISI. MULA2 TKR = 0 TKR = 0 MENGHITUNG JUMLAH OPERASI BAGI/KALI (BK). BK = 0 MENGHITUNG JUMLAH OPERASI TAMBAH/KURANG (TK). TK = 0 LANGKAH 3: MEMERIKSA ELEMEN2 PIVOT DAN PROSES TUKAR POSISI NN = N-1 DO 10 I=1,NN LANGKAH 4: MENDEFINISIKAN P P = I IF (ABS(A(P,I)).GE.1.0E-20 .OR. P.GT.N) GOTO 200 P = P+1 GOTO 100 IF(P.EQ.N+1)THEN MENAMPILKAN PESAN SINGULAR WRITE(*,5) GOTO 400 END IF LANGKAH 5: PROSES TUKAR POSISI IF(P.NE.I) THEN DO 20 JJ=1,N+N C = A(I,JJ) A(I,JJ) = A(P,JJ) A(P,JJ) = C TKR = TKR + 1 CONTINUE END IF LANGKAH 6: PERSIAPAN PROSES TRIANGULARISASI JJ = I+1

35

36 65 66

C

67 68 69

C

70 71 72 73 74

40

75 76 77 78

30 10 C

79 80 81 82 83

120 C

84 85

C

86 87 88 89 90 91

C

92 93 94

C

95 96 97 98 99 100 101 102 103 104

16

105 106 107 108 109 110

15 500 C

111 112 113 114 115

220

116 117 118 119 120 121 122 123

400 5 8 9 11

BAB 2. METODE ELIMINASI GAUSS DO 30 J=JJ,N LANGKAH 7: TENTUKAN MJI MJI = A(J,I)/A(I,I) BK = BK + 1 LANGKAH 8: MELAKUKAN PROSES TRIANGULARISASI DO 40 K=JJ,N+N A(J,K) = A(J,K)-MJI*A(I,K) BK = BK + 1 TK = TK + 1 CONTINUE A(J,I) = 0 CONTINUE CONTINUE MENAMPILKAN HASIL TRIANGULARISASI WRITE (*,’(1X,A)’) ’HASIL TRIANGULARISASI:’ DO 120 I = 1,N WRITE (*,’(1X,5(F14.8))’) (A(I,J),J=1,N+N) CONTINUE LANGKAH 9: MEMERIKSA ELEMEN A(N,N) IF(ABS(A(N,N)).LT.1.0E-20) THEN MENAMPILKAN PESAN SINGULAR WRITE(*,5) GOTO 400 END IF DO 500 J = 1,N Q=N+J LANGKAH 10: MENGHITUNG A(N,N) D(J,N) = A(N,Q)/A(N,N) BK = BK + 1 LANGKAH 11: PROSES SUBSTITUSI MUNDUR L = N-1 DO 15 K=1,L I = L-K+1 JJ = I+1 SUM = 0.0 DO 16 KK=JJ,N SUM = SUM+A(I,KK)*D(J,KK) BK = BK + 1 TK = TK + 1 CONTINUE D(J,I) = (A(I,Q)-SUM)/A(I,I) BK = BK + 1 TK = TK + 1 CONTINUE CONTINUE LANGKAH 12: MENAMPILKAN HASIL PERHITUNGAN WRITE (*,*) WRITE (*,’(1X,A)’) ’MATRIK INVERS:’ DO 220 I = 1,N WRITE (*,’(1X,5(F14.8))’) (D(J,I),J=1,N) CONTINUE WRITE(*,8) TKR WRITE(*,9) BK WRITE(*,11) TK STOP FORMAT(1X,’MATRIK A BERSIFAT SINGULAR’) FORMAT(1X,’JUMLAH TUKAR POSISI = ’,3X,I5) FORMAT(1X,’JUMLAH OPERASI BAGI/KALI = ’,3X,I6) FORMAT(1X,’JUMLAH OPERASI JUMLAH/KURANG = ’,3X,I6)

2.7. PENUTUP

37

END

124

2.7

Penutup

Silakan anda coba aplikasikan program di atas dengan berbagai sistem persamaan linear yang pernah dijadikan contoh pada catatan terdahulu. Saya cukupkan sementara sampai disini. Insya Allah akan saya sambung lagi dilain waktu. Kalau ada yang mau didiskusikan, silakan hubungi saya melalui email yang tercantum di halaman paling depan.

Bab 3

Aplikasi Eliminasi Gauss pada Masalah Inversi

✍ Objektif : ⊲ Mengenalkan model garis. ⊲ Mengenalkan model parabola. ⊲ Mengenalkan model bidang.

Pada bab ini, saya mencoba menuliskan aplikasi Metode Eliminasi Gauss sebagai dasardasar teknik inversi yaitu meliputi model garis, model parabola dan model bidang. Uraian aplikasi tersebut diawali dari ketersediaan data observasi, lalu sejumlah parameter model mesti dicari dengan teknik inversi. Mari kita mulai dari model garis.

3.1

Inversi Model Garis

Pengukuran temperatur terhadap kedalaman di bawah permukaan bumi menunjukkan bahwa semakin dalam, temperatur semakin tinggi. Misalnya telah dilakukan sebanyak empat kali (N = 4) pengukuran temperatur (Ti ) pada kedalaman yang berbeda beda (zi ). Tabel pengukuran secara sederhana disajikan seperti ini: Lalu kita berasumsi bahwa variasi temperatur terhadap Tabel 3.1: Data temperatur bawah permukaan tanah terhadap kedalaman Pengukuran ke-i Kedalaman (m) Temperatur (O C) 1 z1 = 5 T1 = 35 2 z2 = 16 T2 = 57 3 z3 = 25 T3 = 75 4 z4 = 100 T4 = 225 kedalaman ditentukan oleh rumus berikut ini: m1 + m2 zi = Ti 39

(3.1)

BAB 3. APLIKASI ELIMINASI GAUSS PADA MASALAH INVERSI

40

dimana m1 dan m2 adalah konstanta-konstanta yang akan dicari. Rumus di atas disebut model. Sedangkan m1 dan m2 disebut model parameter. Jadi pada model di atas terdapat dua buah model parameter, (M = 2). Adapun yang berlaku sebagai data adalah nilai-nilai temperatur T1 , T2 ,..., dan T4 . Berdasarkan model tersebut, kita bisa menyatakan temperatur dan kedalaman masing-masing sebagai berikut: m1 + m2 z1 = T1 m1 + m2 z2 = T2 m1 + m2 z3 = T3 m1 + m2 z4 = T4 Semua persamaan tersebut dapat dinyatakan dalam operasi matrik berikut ini: 

1 z1

  1 z2   1 z 3  1 z4





T1

#  "  m1  T2   =  m  T 2   3 T4

     

(3.2)

Lalu ditulis secara singkat Gm = d

(3.3)

dimana d adalah data yang dinyatakan dalam vektor kolom, m adalah model parameter, juga dinyatakan dalam vektor kolom, dan G disebut matrik kernel. Lantas bagaimana cara mendapatkan nilai m1 dan m2 pada vektor kolom m? Manipulasi berikut ini bisa menjawabnya Gt Gm = Gt d

(3.4)

dimana t disini maksudnya adalah tanda transpos matrik. Selanjutnya, untuk mendapatkan elemen-elemen m, diperlukan langkah-langkah perhitungan berikut ini: 1. Tentukan transpos dari matrik kernel, yaitu Gt 

1 z1



   1 z2    G=  1 z 3   1 z4

⇒

Gt =

"

1

1

1

1

z1 z2 z3 z4

#

2. Tentukan Gt G

Gt G =

"

1

1

1

z1 z2 z3



1 z1

#  1 z2  z4   1 z3 1 z4 1



P #  "  N zi = P P 2  zi zi 

3.1. INVERSI MODEL GARIS

41

dimana N = 4 dan i = 1, 2, 3, 4.

3. Kemudian tentukan pula Gt d

Gt d =

"

1

1



T1



# #  " P  T2  T i  = P  z4  T z i Ti  3  T4

1

1

z1 z2 z3

4. Sekarang persamaan (3.4) dapat dinyatakan sebagai "

# " P # P #" N zi m1 Ti = P P P 2 zi zi m2 z i Ti

(3.5)

5. Aplikasikan metode Eliminasi Gauss dengan Substitusi Mundur. Untuk itu, tentukan matrik augment-nya # P P N zi | Ti P P P 2 z i Ti zi zi |

"

6. Untuk mempermudah perhitungan, kita masukan dulu angka-angka yang tertera pada tabel pengukuran dihalaman depan. "

4

146

|

392

146 10906 | 25462

#

7. Lakukan proses triangularisasi dengan operasi (P2 − (36, 5)P1 ) → P2 . Saya sertakan pula indeks masing-masing elemen pada matrik augment sebagaimana yang telah saya lakukan pada catatan kuliah yang berjudul Metode Eliminasi Gauss. Hasilnya adalah "

4

146

|

392

0 5577 | 11154

#

=

"

a11 a12 | a13

a21 a22 | a23

#

8. Terakhir, tentukan konstanta m1 dan m2 yang merupakan elemen-elemen vektor kolom m, dengan proses substitusi mundur. Pertama tentukan m2 m2 = lalu tentukan m1 m1 =

11154 a23 = =2 a22 5577

392 − (146)(2) a13 − a12 m2 = = 25 a11 4

42


3.1.1 Script matlab inversi model garis Script inversi model garis ini dibangun dari beberapa script yang sudah kita pelajari sebelumnya, yaitu script transpose matriks, perkalian matrik dan script eliminasi gauss. Silakan pelajari maksud tiap-tiap baris pada script ini.

1 2

clc clear all

3 4 5 6 7 8

disp(’Data observasi’) z1=5; z2=16; z3=25; z4=100;

9 10 11 12 13

T(1,1)=35; T(2,1)=57; T(3,1)=75; T(4,1)=225;

14 15 16 17 18 19 20 21 22 23 24 25 26

disp(’Elemen-elemen matriks kernel G’) G(1,1)=1; G(1,2)=z1; G(2,1)=1; G(2,2)=z2; G(3,1)=1; G(3,2)=z3; G(4,1)=1; G(4,2)=z4; G d=T; d

27 28 29

N=4; %jumlah data M=2; %model parameter

30 31 32 33 34 35 36 37

disp(’Mencari G transpos’) for i=1:N for j=1:M GT(j,i)=G(i,j); end end GT

38 39 40 41 42 43 44 45 46 47 48 49 50 51

disp(’Perkalian GT dan G’) for i=1:M for j=1:M GTG(i,j)=0; end end for i=1:M for j=1:M for k=1:N GTG(i,j)=GTG(i,j)+GT(i,k)*G(k,j); end end end

3.1. INVERSI MODEL GARIS 52

GTG

53 54 55 56 57 58 59 60 61 62 63 64 65 66 67

disp(’Perkalian GT dan d’) for i=1:M for j=1:1 GTd(i,j)=0; end end for i=1:M for j=1:1 for k=1:N GTd(i,j)=GTd(i,j)+GT(i,k)*d(k,j); end end end GTd

68 69 70 71 72 73 74 75

A=GTG; %====== Menggabungkan Vektor GTd kedalam matrik A ======== n=M; for i=1:n A(i,n+1)=GTd(i,1); end A

76 77 78 79

%&&&&&& Proses Eliminasi Gauss &&&&&&&&&&&&&&&&&&&&&&&&& %---------Proses Triangularisasi----------for j=1:(n-1)

80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99

%----mulai proses pivot--if (A(j,j)==0) for p=1:n+1 u=A(j,p); v=A(j+1,p); A(j+1,p)=u; A(j,p)=v; end end %----akhir proses pivot--jj=j+1; for i=jj:n m=A(i,j)/A(j,j); for k=1:(n+1) A(i,k)=A(i,k)-(m*A(j,k)); end end end %-------------------------------------------

100 101 102


103 104 105 106 107 108 109 110

for i=n-1:-1:1 S=0; for j=n:-1:i+1 S=S+A(i,j)*x(j,1); end x(i,1)=(A(i,n+1)-S)/A(i,i); end

43


44 111 112 113

%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& disp(’Model parameter yang dicari’) m=x

Sebetulnya, matlab telah menyediakan fungsi-fungsi intrinsik yang bisa digunakan sehingga dapat memperkecil jumlah baris pada script di atas. Dari line 28 sampai line 113 dapat dipangkas menjadi 1

m=inv(G’*G)*G’*d

%Proses inversi linear

Lalu mengapa kita harus bersusah payah membangun script yang begitu panjang bila matlab bisa melakukannya dengan mudah? Karena kita sedang mempelajari teknik-teknik komputasi untuk menyelesaikan problem sains dan teknik. Kita tidak sedang belajar matlab. Jadi teknik-teknik yang dipelajari disini harus bisa diterapkan di selain matlab. Script singkat m = inv(G′ ∗ G) ∗ G′ ∗ d hanya berlaku di matlab, sementara script yang panjangnya 113 line dapat diterjemahkan dengan sangat mudah ke dalam bahasa pemrograman selain matlab.

Demikianlah contoh aplikasi metode Eliminasi Gauss dengan substitusi mundur. Anda bisa mengaplikasikan pada kasus lain, dengan syarat kasus yang anda tangani memiliki bentuk model yang sama dengan yang telah dikerjakan pada catatan ini, yaitu model persamaan garis atau disingkat model garis: y = m1 + m2x. Selanjutnya mari kita pelajari inversi model parabola.

3.2

Inversi Model Parabola

Pengukuran temperatur terhadap kedalaman di bawah permukaan bumi menunjukkan bahwa semakin dalam, temperatur semakin tinggi. Misalnya telah dilakukan sebanyak delapan kali (N = 8) pengukuran temperatur (Ti ) pada kedalaman yang berbeda beda (zi ). Tabel pengukuran secara sederhana disajikan seperti ini: Tabel 3.2: Data temperatur bawah permukaan tanah terhadap kedalaman Pengukuran ke-i Kedalaman (m) Temperatur (O C) 1 z1 = 5 T1 = 21, 75 2 z2 = 8 T2 = 22, 68 3 z3 = 14 T3 = 25, 62 4 z4 = 21 T4 = 30, 87 5 z5 = 30 T5 = 40, 5 6 z6 = 36 T6 = 48, 72 7 z7 = 45 T7 = 63, 75 8 z8 = 60 T8 = 96 Lalu kita berasumsi bahwa variasi temperatur terhadap kedalaman ditentukan oleh rumus berikut ini: m1 + m2 zi + m3 zi2 = Ti

(3.6)

dimana m1 , m2 dan m3 adalah konstanta-konstanta yang akan dicari. Rumus di atas disebut model. Sedangkan m1 , m2 dan m3 disebut model parameter. Jadi pada model di atas terdapat

3.2. INVERSI MODEL PARABOLA

45

tiga buah model parameter, (M = 3). Adapun yang berlaku sebagai data adalah nilai-nilai temperatur T1 , T2 ,..., dan T8 . Berdasarkan model tersebut, kita bisa menyatakan temperatur dan kedalaman masing-masing sebagai berikut: m1 + m2 z1 + m3 z12 = T1 m1 + m2 z2 + m3 z22 = T2 m1 + m2 z3 + m3 z32 = T3 m1 + m2 z4 + m3 z42 = T4 m1 + m2 z5 + m3 z52 = T5 m1 + m2 z6 + m3 z62 = T6 m1 + m2 z7 + m3 z72 = T7 m1 + m2 z8 + m3 z82 = T8 Semua persamaan tersebut dapat dinyatakan dalam operasi matrik berikut ini:                

1 z1 z12





T1



    T2  1 z2 z22        1 z3 z32   T3   m    1 T4  1 z4 z42        m2  =   T 1 z5 z52   5      m 3  T6  1 z6 z62      T  1 z7 z72   7   T8 1 z8 z82

(3.7)


(3.8)

dimana d adalah data yang dinyatakan dalam vektor kolom, m adalah model parameter, juga dinyatakan dalam vektor kolom, dan G disebut matrik kernel. Lantas bagaimana cara mendapatkan nilai m1 , m2 dan m3 pada vektor kolom m? Manipulasi berikut ini bisa menjawabnya Gt Gm = Gt d

(3.9)

dimana t disini maksudnya adalah tanda transpos matrik. Selanjutnya, untuk mendapatkan elemen-elemen m, diperlukan langkah-langkah perhitungan berikut ini:


46

1. Tentukan transpos dari matrik kernel, yaitu Gt 

       G=       

1 z1 z12



 1 z2 z22   1 z3 z32    2 1 z4 z4   1 z5 z52    1 z6 z62   1 z7 z72   1 z8 z82

⇒





1

1

1

1

1

1

1

1

 Gt =  z1

z2

z3

z4

z5

z6

z7

 z8  z82

z12 z22 z32 z42 z52 z62 z72

2. Tentukan Gt G 



1 1 1 1 1 1 1  G G =  z1 z2 z3 z4 z5 z6 z7 z12 z22 z32 z42 z52 z62 z72 t

     1   z8     z82     

1 z1 z12



 1 z2 z22   1 z3 z32  P P 2    N zi z  2 1 z4 z4  P P 2 P i3  = zi z z   1 z5 z52  P 2 P i3 P i4   zi zi zi 1 z6 z62   1 z7 z72   1 z8 z82

dimana N = 8 dan i = 1, 2, 3, ..., 8. 3. Kemudian tentukan pula Gt d 



1 1 1 1 1 1 1  G d =  z1 z2 z3 z4 z5 z6 z7 z12 z22 z32 z42 z52 z62 z72 t

     1   z8     z82     

T1



 T2    T3    P Ti  T4  P  = z i Ti    T5  P 2  z i Ti T6   T7   T8

4. Sekarang persamaan (3.14) dapat dinyatakan sebagai (ini khan least square juga...!?)    P P P 2  N zi zi m1 Ti P 2 P 3   P    P zi zi z i Ti  zi   m2  =   P 2 P 2 P 3 P 4 m3 z i Ti zi zi zi 

(3.10)

5. Aplikasikan metode Eliminasi Gauss dengan Substitusi Mundur. Untuk itu, tentukan

3.2. INVERSI MODEL PARABOLA

47

matrik augment-nya  P P P 2 Ti N zi zi | P 2 P 3 P  P  zi zi z i Ti  zi |  P 2 P 3 P 4 P 2 zi zi zi | z i Ti 

6. Untuk mempermudah perhitungan, kita masukan dulu angka-angka yang tertera pada tabel pengukuran dihalaman depan. 

8

219

8547

|

349, 89



  8547 393423 | 12894, 81   219 8547 393423 19787859 | 594915, 33 7. Lakukan proses triangularisasi dengan operasi (P2 − (219/8)P1 ) → P2 . Hasilnya adalah   

8

219

0 8547

8547

|



349, 89

2551, 88 159448, 88 | 393423

 3316, 57  | 594915, 33

19787859

8. Masih dalam proses triangularisai, operasi berikutnya (P3 − (8547/8)P1 ) → P3 . Hasilnya adalah



8

219

8547

|



349, 89

  159448, 88 | 3316, 57   0 2551, 88 0 159448.88 10656457, 88 | 221101, 6 9. Masih dalam proses triangularisai, operasi berikutnya (P3 − (159448, 88/2551, 88)P2 ) → P3 . Hasilnya adalah



8

219

8547

|

349, 89



(3.11)

   0 2551, 88 159448, 88 | 3316, 57  0 0 693609, 48 | 13872, 19

Seperti catatan yang lalu, saya ingin menyertakan pula notasi masing-masing elemen pada matrik augment sebelum melakukan proses substitusi mundur. 

8

219

8547

|

349, 89





a11 a12 a13 | a14



     0 2551, 88 159448, 88 | 3316, 57  ⇔  a21 a22 a23 | a24  a31 a32 a33 | a34 0 0 693609, 48 | 13872, 19 10. Terakhir, tentukan konstanta m1 , m2 dan m3 yang merupakan elemen-elemen vektor kolom m, dengan proses substitusi mundur. Pertama tentukan m3 m3 =

13872, 19 a34 = = 0, 02 a33 693609, 48


48 lalu m2

m2 =

a24 − a23 m3 3316, 57 − (159448, 88)(0, 02) = = 0, 05 a22 2551, 88

dan m1 m1 =

a14 − (a12 m2 + a13 m3 ) 349, 89 − [(219)(0, 05) + (8547)(0, 02) = = 21 a11 8

3.2.1 Script matlab inversi model parabola Perbedaan utama script ini dengan script inversi model garis terletak pada inisialisasi elemenelemen matrik kernel. Elemen-elemen matrik kernel sangat ditentukan oleh model matematika yang digunakan. Seperti pada script ini, matrik kernelnya diturunkan dari persamaan parabola. 1 2

clc clear all

3 4 5 6 7 8 9 10 11

z1=5; z2=8; z3=14; z4=21; z5=30; z6=36; z7=45; z8=60;

12 13 14 15 16 17 18 19 20

T(1,1)=21.75; T(2,1)=22.68; T(3,1)=25.62; T(4,1)=30.87; T(5,1)=40.5; T(6,1)=48.72; T(7,1)=63.75; T(8,1)=96;

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

G(1,1)=1; G(1,2)=z1; G(1,3)=z1^2; G(2,1)=1; G(2,2)=z2; G(2,3)=z2^2; G(3,1)=1; G(3,2)=z3; G(3,3)=z3^2; G(4,1)=1; G(4,2)=z4; G(4,3)=z4^2; G(5,1)=1; G(5,2)=z5; G(5,3)=z5^2; G(6,1)=1; G(6,2)=z6; G(6,3)=z6^2; G(7,1)=1; G(7,2)=z7;

3.2. INVERSI MODEL PARABOLA 42 43 44 45

G(7,3)=z7^2; G(8,1)=1; G(8,2)=z8; G(8,3)=z8^2;

46 47 48 49

G d=T; d

50 51 52 53

N=8; %jumlah data M=3; %model parameter pause

54 55 56 57 58 59 60 61 62 63

%%%%%===========Proses inversi============== disp(’Mencari G transpos’) for i=1:N for j=1:M GT(j,i)=G(i,j); end end GT

64 65 66 67 68 69 70 71 72 73 74 75 76 77 78

disp(’Perkalian GT dan G’) for i=1:M for j=1:M GTG(i,j)=0; end end for i=1:M for j=1:M for k=1:N GTG(i,j)=GTG(i,j)+GT(i,k)*G(k,j); end end end GTG

79 80 81 82 83 84 85 86 87 88 89 90 91 92 93


94 95 96 97 98 99 100

A=GTG; %====== Menggabungkan Vektor GTd kedalam matrik A ======== n=M; for i=1:n A(i,n+1)=GTd(i,1); end

49


50 101 102

A pause

103 104 105 106 107

disp(’Hasil Eliminasi Gauss’) %&&&&&& Proses Eliminasi Gauss &&&&&&&&&&&&&&&&&&&&&&&&& %---------Proses Triangularisasi----------for j=1:(n-1)

108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127


128 129 130


131 132 133 134 135 136 137 138 139 140 141

for i=n-1:-1:1 S=0; for j=n:-1:i+1 S=S+A(i,j)*x(j,1); end x(i,1)=(A(i,n+1)-S)/A(i,i); end %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& %%%%%%%%%%%=====AKHIR DARI INVERSI MODEL GARIS========== m=x

Demikianlah contoh aplikasi metode Eliminasi Gauss dengan substitusi mundur. Anda bisa mengaplikasikan pada kasus lain, dengan syarat kasus yang anda tangani memiliki bentuk model yang sama dengan yang telah dikerjakan pada catatan ini, yaitu memiliki tiga buah model parameter yang tidak diketahui dalam bentuk persamaan parabola: y = m1 + m2 x + m3 x2 . Pada catatan berikutnya, saya akan membahas model yang mengandung tiga model parameter dalam 2 dimensi.

3.3

Inversi Model Bidang

Dalam catatan ini saya belum sempat mencari contoh pengukuran yang sesuai untuk model 2-dimensi. Maka, saya ingin langsung saja mengajukan sebuah model untuk 2-dimensi berikut

3.3. INVERSI MODEL BIDANG

51

ini: m1 + m2 xi + m3 yi = di

(3.12)

dimana m1 , m2 dan m3 merupakan model parameter yang akan dicari. Adapun yang berlaku sebagai data adalah d1 , d2 , d3 , ..., di . Berdasarkan model tersebut, kita bisa menyatakan temperatur dan kedalaman masing-masing sebagai berikut: m1 + m2 x1 + m3 y1 = d1 m1 + m2 x2 + m3 y2 = d2 m1 + m2 x3 + m3 y3 = d3 .. .. .. .. .. . . . . . m1 + m2 xN + m3 yN = dN Semua persamaan tersebut dapat dinyatakan dalam operasi matrik berikut ini:         





d1



1

x1

y1

1

x2

1 .. .

x3 .. .

     d2  y2     m1   d   y3  = 3    m2    .  ..     . .  m3  .  dN yN

1 xN


(3.13)

dimana d adalah data yang dinyatakan dalam vektor kolom, m adalah model parameter, juga dinyatakan dalam vektor kolom, dan G disebut matrik kernel. Lantas bagaimana cara mendapatkan nilai m1 , m2 dan m3 pada vektor kolom m? Manipulasi berikut ini bisa menjawabnya Gt Gm = Gt d

(3.14)

dimana t disini maksudnya adalah tanda transpos matrik. Selanjutnya, untuk mendapatkan elemen-elemen m, diperlukan langkah-langkah perhitungan berikut ini:

1. Tentukan transpos dari matrik kernel, yaitu Gt 

    G=   



1

x1

y1

1

x2

1 .. .

x3 .. .

 y2   y3   ..  .   yN

1 xN

⇒



1

1

1

···

1



  Gt =  x1 x2 x3 · · · xN  y1 y2 y3 · · · y N


52 2. Tentukan Gt G



1

1

1

···



1

 Gt G =  x1 x2 x3 · · · xN y1

y2

y3 · · · y N

       



1

x1

y1

1

x2

1 .. .

x3 .. .

   P P y2  N xi yi  P 2 P  P  y3  = xi xi xi yi   P P P 2 ..  yi xi yi yi .  yN

1 xN

dimana N = jumlah data. dan i = 1, 2, 3, ..., N .

3. Kemudian tentukan pula Gt d



1

1

1

···

1

 Gt d =  x1 x2 x3 · · · xN y1

y2

y3

· · · yN



       

d1



  P  d2  di    P d3  = xi di   P ..  yi d i .  dN

4. Sekarang, persamaan (3.14) dapat dinyatakan sebagai    P  P P di m1 yi N xi P 2 P    P   P xi di  xi xi xi yi   m2  =   P P P P 2 yi d i m3 yi xi yi yi 

(3.15)

5. Aplikasikan metode Eliminasi Gauss dengan Substitusi Mundur. Untuk itu, tentukan matrik augment-nya 

 P P P N xi yi | di P 2 P P  P  xi xi xi yi | xi di   P P P 2 P yi xi yi yi | yi d i 6. Langkah-langkah selanjutnya akan sama persis dengan catatan sebelumnya (model linear dan model parabola)

Anda bisa mengaplikasikan data pengukuran yang anda miliki, dengan syarat kasus yang anda tangani memiliki bentuk model yang sama dengan yang telah dikerjakan pada catatan ini, yaitu memiliki tiga buah model parameter yang tidak diketahui dalam bentuk persamaan bidang (atau 2-dimensi): d = m1 + m2 x + m3 y. Saya cukupkan sementara sampai disini. Insya Allah akan saya sambung lagi dilain waktu. Kalau ada yang mau didiskusikan, silakan hubungi saya melalui email: supri@f isika.ui.ac.id.

3.4. CONTOH APLIKASI

3.4

53

Contoh aplikasi

3.4.1 Menghitung gravitasi di planet X Seorang astronot tiba di suatu planet yang tidak dikenal. Setibanya disana, ia segera mengeluarkan kamera otomatis, lalu melakukan ekperimen kinematika yaitu dengan melempar batu vertikal ke atas. Hasil foto-foto yang terekam dalam kamera otomatis adalah sebagai berikut Tabel 3.3: Data ketinggian terhadap waktu dari planet X Waktu (dt) Ketinggian (m) Waktu (dt) Ketinggian (m) 0,00 5,00 2,75 7,62 0,25 5,75 3,00 7,25 0,50 6,40 3,25 6,77 0,75 6,94 3,50 6,20 1,00 7,38 3,75 5,52 1,25 7,72 4,00 4,73 1,50 7,96 4,25 3,85 1,75 8,10 4,50 2,86 2,00 8,13 4,75 1,77 2,25 8,07 5,00 0,58 2,50 7,90 Plot data pengukuran waktu vs ketinggian diperlihatkan sebagai berikut 9

8

7

Tinggi (meter)

6

5

4

3

2

1

0

0

0.5

1

1.5

2

2.5 Waktu (detik)

3

3.5

4

4.5

5

Gambar 3.1: Grafik data pengukuran gerak batu Anda diminta untuk membantu pengolahan data di atas. Jika anda menggunakan asumsi

54


model matematik dari Gerak-Lurus-Berubah-Beraturan (GLBB) seperti ini 1 ho + vo t − gt2 = h 2 maka gunakanlah prinsip-prinsip inversi untuk menentukan kecapatan awal, vo dan konstanta gravitasi, g pada planet tersebut.

jawab: Berdasarkan tabel di atas, diketahui terdapat 21 data. Ketinggian pada saat t = 0 adalah ho = 5 m. Untuk mencari vo dan g menggunakan metode inversi, mula-mula kita definisikan terlebih dahulu m1 dan m2 : m1 = vo

1 m2 = − g 2

sehingga persamaan model GLBB menjadi 5 + m1 ti + m2 t2i = hi dimana i menunjukkan data ke-i. Langkah selanjutnya adalah menentukan nilai tiap-tiap elemen matrik kernel, yaitu dengan memasukan semua data kedalam persamaan model GLBB 5 + m1 t1 + m2 t21 = h1 5 + m1 t2 + m2 t22 = h2 5 + m1 t3 + m2 t23 = h3 .. .. . . . = .. 5 + m1 t20 + m2 t220 = h20 Semua persamaan tersebut dapat dinyatakan dalam operasi matrik berikut ini:              

5

t1

5

t2

5

t3

5 .. .

t4 .. .

5 t19 5 t20

t21



 h1      h2 " #    h3  m1   =  .  m  .. 2 ..   .     h19  2 t19  h20 t220 t22 t23 t24

          

Sebelum dilanjut, coba perhatikan dengan teliti operasi matrik di atas. Adakah yang janggal?? Yep.. matrik kernel G berukuran 20x3 sementara vektor m berukuran 2x1, tentu saja operasi perkalian matrik akan gagal. Untuk menghindarinya, kita tambahkan m0 pada vektor m, se-


55

hingga operasi tersebut menjadi              

5

t1

5

t2

5

t3

5 .. .

t4 .. .

5 t19 5 t20

t21



 h1       h2    m0  h3     m1  =   .   .. ..   m2 .     h19  2 t19  h20 t220 t22 t23 t24

          

Namun, perlu dicatat bahwa m0 harus punya syarat, yaitu harus bernilai 1 atau m0 =1. Ini boleh dibilang sebagai sebuah aksioma, atau sesuatu yang tak perlu dibuktikan lagi tapi tak bisa dibantah. Tinggal nanti bisa kita periksa hasil inversinya. Bila m0 bernilai 1, maka proses inversi dianggap sukses. Kemudian operasi matrik tersebut bisa ditulis secara singkat Gm = d Untuk menyelesaikan persamaan matrik ini, diperlukan modifikasi berikut GT Gm = GT d

(3.16)

dimana T disini maksudnya adalah tanda transpos matrik. Selanjutnya, untuk mendapatkan m0 , m1 dan m2 , diperlukan langkah-langkah perhitungan berikut ini: 1. Tentukan transpos dari matrik kernel, yaitu Gt 

      G=      



5

t1

t21

5

t2

5

t3

5 .. .

t4 .. .

 t22   t23    2 t4   ..  .   t219   2 t20

5 t19 5 t20

⇒



5

5

5

5

...

5

5



  GT =  t1 t2 t3 t4 . . . t19 t20  t21 t22 t23 t24 . . . t219 t220

2. Tentukan GT G 



5 5 5 5 ... 5  G G =  t1 t2 t3 t4 . . . t19 t21 t22 t23 t24 . . . t219 T

    5   t20     t220    



5

t1

t21

5

t2

5

t3

5 .. .

t4 .. .

 t22    P P  t23  25N 5 ti 5 t2i    P P 3  P 2 t24  t   =  5 P ti P ti P i4 ..  3 2 ti ti 5 ti .   t219   2 t20

5 t19 5 t20


56

dimana N = 20 dan i = 1, 2, ..., 20. 3. Kemudian tentukan pula GT d 



5 5 5 5 ... 5  G d =  t1 t2 t3 t4 . . . t19 t21 t22 t23 t24 . . . t219 T

h1



 h2    P  h3  5 hi    P  h4   =  P ti hi  ..  t2i hi .   h19   h20

    5   t20     t220    

4. Sekarang persamaan (3.16) dapat dinyatakan sebagai   P  P P  25N 5 ti 5 t2i m0 5 hi P 2 P 3   P   P  ti ti   m1  =  ti hi   5 ti P 2 P 3 P 4 P 2 5 ti ti ti m2 ti hi 



500

262, 5

896, 9



m0





607, 5

(3.17)



      262, 5 179, 4 689, 1   m1  =  273, 7  896, 9 689, 1 2822, 9 m2 796, 3 Hasil operasi matriks ini dapat diselesaikan dengan metode Eliminasi Gauss, yaitu 

m0





0, 9999



     m1  =  3, 2009  −0, 8169 m2 Lihatlah! m0 = 0,9999 atau mendekati 1. Dan ini sesuai dengan aksioma yang telah dinyatakan di awal bahwa memang m0 harus bernilai 1. Jika m0 tidak bernilai 1 berarti teknik inversinya salah total. Hasil inversi juga menunjukkan bahwa kecepatan awal yaitu saat batu dilempar ke atas adalah sebesar m1 = vo = 3,2009 m/dt. Adapun percepatan gravitasi diperoleh dari m2 dimana m2 = − 12 g = -0,8169. Sehingga nilai g adalah sebesar 1,6338 m/dt2 . Gambar 3.2 memperlihatkan grafik kurva hasil inversi. Garis berwarna biru merupakan garis kurva fitting hasil inversi parabola. Sedangkan bulatan berwarna merah adalah data pengukuran ketinggian (m) terhadap waktu (dt). Jelas terlihat bahwa garis kurva berwarna biru benar-benar cocok melewati semua titik data pengukuran. Ini menunjukkan tingkat akurasi yang sangat tinggi. Sehingga nilai kecepatan awal dan gravitasi hasil inversi cukup valid untuk menjelaskan gerak batu di planet X.


57

9

8

7

Ketinggian (m)

6

5

4

3

2

1

0

0

0.5

1

1.5

2

2.5 Waktu (dt)

3

3.5

4

4.5

5

Gambar 3.2: Grafik hasil inversi parabola Berikut adalah script inversi dalam Matlab untuk memecahkan masalah ini 1 2

clear all clc;

3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

N=20; %jumlah data M=3; %model parameter for i=1:N t(i)=i*0.25; end h(1)=5.75; h(2)=6.40; h(3)=6.94; h(4)=7.38; h(5)=7.72; h(6)=7.96; h(7)=8.10; h(8)=8.13; h(9)=8.07; h(10)=7.90; h(11)=7.62; h(12)=7.25; h(13)=6.77; h(14)=6.20; h(15)=5.52; h(16)=4.73; h(17)=3.85; h(18)=2.86; h(19)=1.77; h(20)=0.58;

29 30 31 32 33

for i=1:N G(i,1)=5; G(i,2)=t(i); G(i,3)=t(i)^2;

58 34 35 36 37 38 39


end G for i=1:N d(i,1)=h(i); end d

40 41 42 43 44 45 46 47 48

%%%%%===========Proses inversi============== disp(’Mencari G transpos’) for i=1:N for j=1:M GT(j,i)=G(i,j); end end GT

49 50 51 52 53 54 55 56 57 58 59 60 61 62 63

disp(’Perkalian GT dan G’) for i=1:M for j=1:M GTG(i,j)=0; end end for i=1:M for j=1:M for k=1:N GTG(i,j)=GTG(i,j)+GT(i,k)*G(k,j); end end end GTG

64 65 66 67 68 69 70 71 72 73 74 75 76 77 78


79 80 81 82 83 84 85 86 87

A=GTG; %====== Menggabungkan Vektor GTd kedalam matrik A ======== n=M; for i=1:n A(i,n+1)=GTd(i,1); end A pause

88 89 90 91 92

disp(’Hasil Eliminasi Gauss’) %&&&&&& Proses Eliminasi Gauss &&&&&&&&&&&&&&&&&&&&&&&&& %---------Proses Triangularisasi----------for j=1:(n-1)

3.4. CONTOH APLIKASI 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114

%----mulai proses pivot--if (A(j,j)==0) for p=1:n+1 u=A(j,p); v=A(j+1,p); A(j+1,p)=u; A(j,p)=v; end end %----akhir proses pivot--jj=j+1; for i=jj:n m=A(i,j)/A(j,j); for k=1:(n+1) A(i,k)=A(i,k)-(m*A(j,k)); end end end %------------------------------------------%------Proses Substitusi mundur------------x(n,1)=A(n,n+1)/A(n,n);

115 116 117 118 119 120 121 122 123 124 125

for i=n-1:-1:1 S=0; for j=n:-1:i+1 S=S+A(i,j)*x(j,1); end x(i,1)=(A(i,n+1)-S)/A(i,i); end %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& %%%%%%%%%%%===== AKHIR DARI PROSES INVERSI ========== m=x

126 127 128 129 130 131 132 133 134 135

%-------MENGGAMBAR GRAFIK---------------------plot(t,h,’ro’); xlabel(’Waktu (dt)’);ylabel(’Ketinggian (m)’); hold on; for i=1:20 hi(i)=5+m(2)*t(i)+m(3)*t(i)^2; end plot(t,hi); hold off;

59

Bab 4

Metode LU Decomposition

✍ Objektif : ⊲ Mengenalkan teknik faktorisasi matrik. ⊲ Mengenalkan aplikasi LU Decomposition pada sistem persamaan linear. ⊲ Merumuskan algoritma LU Decomposition.

4.1

Faktorisasi matrik

Pada semua catatan yang terdahulu, telah diulas secara panjang lebar bahwa sistem persamaan linear dapat dicari solusinya secara langsung dengan metode eliminasi gauss. Namun perlu juga diketahui bahwa eliminasi gauss bukan satu-satunya metode dalam mencari solusi sistem persamaan linear, misalnya ada metode matrik inversi seperti yang dijelaskan pada catatan yang paling terakhir. Terlepas dari masalah in-efisiensi penyelesaiannya, yang jelas metode invers matrik bisa digunakan untuk menyelesaikan sistem persamaan linear. Nah, pada catatan kali ini, saya ingin mengetengahkan sebuah metode yang lain untuk menyelesaikan sistem persamaan linear, yaitu metode faktorisasi matrik yang umum dikenal sebagai LU-decomposition. Metode ini sekaligus menjadi pengantar menuju metode Singular Value Decomposition, (SVD), suatu metode yang saat ini paling “handal” dalam menyelesaikan sistem persamaan linear dan merupakan bagian dari metode least square. Seperti biasa, kita berasumsi bahwa sistem persamaan linear dapat dinyatakan dalam operasi matrik Ax = b

(4.1)

Pada metode LU-decomposition, matrik A difaktorkan menjadi matrik L dan matrik U, dimana dimensi atau ukuran matrik L dan U harus sama dengan dimensi matrik A. Atau dengan kata lain, hasil perkalian matrik L dan matrik U adalah matrik A, A = LU 61

(4.2)

BAB 4. METODE LU DECOMPOSITION

62 sehingga persamaan (6.4) menjadi LU x = b

Langkah penyelesaian sistem persamaan linear dengan metode LU-decomposition, diawali dengan menghadirkan vektor y dimana, Ux = y

(4.3)

Langkah tersebut tidak bermaksud untuk menghitung vektor y, melainkan untuk menghitung vektor x. Artinya, sebelum persamaan (4.3) dieksekusi, nilai-nilai yang menempati elemenelemen vektor y harus sudah diketahui. Lalu bagaimana cara memperoleh vektor y? Begini caranya, Ly = b

(4.4)

Kesimpulannya, metode LU-decomposition dilakukan dengan tiga langkah sebagai berikut: • Melakukan faktorisasi matrik A menjadi matrik L dan matrik U → A = LU . • Menghitung vektor y dengan operasi matrik Ly = b. Ini adalah proses forward-substitution atau substitusi-maju.

• Menghitung vektor x dengan operasi matrik U x = y. Ini adalah proses backward-substitution atau substitusi-mundur.

Metode LU-decomposition bisa dibilang merupakan modifikasi dari eliminasi gauss, karena beberapa langkah yang mesti dibuang pada eliminasi gauss, justru harus dipakai oleh LUdecomposition. Untuk lebih jelasnya, perhatikan contoh berikut ini. Diketahui sistem persamaan linear sebagai berikut P1 P2 P3 P4

: : : :

x1 2x1 3x1 −x1

+ + − +

x2 x2 x2 2x2

− − +

x3 x3 3x3

+ + + −

3x4 x4 2x4 x4

= = = =

4 1 -3 4

Sistem tersebut dapat dinyatakan dalam operasi matrik Ax = y,      

1

1

0

3



x1





4

      1    x2  =  1    3 −1 −1 2   x3    −3 −1 2 3 −1 x4 4 2

1 −1

     

(4.5)

Pada metode eliminasi gauss, matrik A dikonversi menjadi matrik triangular melalui urutan operasi-operasi berikut: (P2 − 2P1 ) → (P2 ), (P3 − 3P1 ) → (P3 ), (P4 − (−1)P1 ) → (P4 ), (P3 −

4P2 ) → (P3 ), (P4 − (−3)P2 ) → (P4 ). Disisi lain, vektor b ikut berubah nilainya menyesuaikan

4.1. FAKTORISASI MATRIK

63

proses triangularisasi, 

1

1

0

3

  0 −1 −1 −5   0 0 3 13  0 0 0 −13



x1







4

     x2   −7      x  =  13  3   x4 −13

    

(4.6)

Lain halnya dengan metode LU-decomposition dimana vektor b tidak mengalami perubahan. Yang berubah hanya matrik A saja, yaitu menjadi matrik L dan matrik U, A = LU 

  A=  

1

1

0

3





1

0 0 0



1

1

0

3

     1  1 0 0  = 2   0 −1 −1 −5    3 −1 −1 2   3 4 1 0  0 3 13  0 −1 2 3 −1 −1 −3 0 1 0 0 0 −13 2

1 −1

     

Jadi matrik L dan U masing-masing adalah 

  L=  

1

0 0 0





 1 0 0   3 4 1 0   −1 −3 0 1 2

1

1

0

3

  0 −1 −1 −5 U =  0 0 3 13  0 0 0 −13

     

Coba bandingkan matrik U di atas dengan matrik hasil triangularisasi dari metode eliminasi gauss pada persamaan (4.6), sama persis bukan? Jadi, cara memperoleh matrik U adalah dengan proses triangularisasi! Lantas, bagaimana cara memperoleh matrik L? Begini caranya: (1) elemen-elemen diagonal matrik L diberi nilai 1 (Asal tahu saja, cara ini dikenal dengan metode Doolittle). (2) elemen-elemen matrik L yang berada di atas elemen-elemen diagonal diberi nilai 0. (3) sedangkan, elemen-elemen matrik L yang berada di bawah elemen-elemen diagonal diisi dengan faktor pengali yang digunakan pada proses triangularisasi eliminasi gauss. Misalnya pada operasi (P2 − 2P1 ) → (P2 ), maka faktor pengalinya adalah 2; pada operasi

(P3 − 3P1 ) → (P3 ), maka faktor pengalinya adalah 3, dan seterusnya.

Inilah letak perbedaannya, seluruh faktor pengali tersebut sangat dibutuhkan pada metode LU-decomposition untuk membentuk matrik L. Padahal dalam metode eliminasi gauss, seluruh faktor pengali tersebut tidak dimanfaatkan alias dibuang begitu saja. Disisi lain, vektor b tidak mengalami proses apapun sehingga nilainya tetap. Jadi, proses konversi matrik pada metode LU-decomposition hanya melibatkan matrik A saja! Setelah langkah faktorisasi matrik A dilalui, maka operasi matrik pada persamaan (4.5) menjadi,      

1

0 0 0



1

1

0

3

  1 0 0    0 −1 −1 −5  3 4 1 0  0 3 13  0 −1 −3 0 1 0 0 0 −13 2



x1





4

     x2   1      x  =  −3  3   x4 4

     

(4.7)


64

Langkah berikutnya adalah menentukan vektor y, dimana Ly = b,      

1

0 0 0





y1



4

      1 0 0    y2  =  1    3 4 1 0    y3   −3 −1 −3 0 1 y4 4 2

     

Dengan proses substitusi-maju, elemen-elemen vektor y dapat ditentukan, y1 = 4, 2y1 + y2 = 1, 3y1 + 4y2 + y3 = −3, −y1 − 3y2 + y4 = 4 maka diperoleh y1 = 4, y2 = −7, y3 = 13, y4 = −13. Langkah terakhir adalah proses substitusi-mundur untuk menghitung vektor x, dimana U x = y, 

1

1

0

3

  0 −1 −1 −5   0 0 3 13  0 0 0 −13



x1





4



      x2   −7   =    x   13   3    x4 −13

Melalui proses ini, yang pertama kali didapat solusinya adalah x4 , kemudian x3 , lalu diikuti x2 , dan akhirnya x1 . x4 = 1 1 (13 − 13x4 ) = 0 x3 = 3 x2 = −(−7 + 5x4 + x3 ) = 2 x1 = 4 − 3x4 − x2 = −1 akhirnya diperoleh solusi x1 = −1, x2 = 2, x3 = 0, dan y4 = 1. Demikianlah contoh penyelesaian sistem persamaan linear dengan metode LU-decomposition.

Sekali matrik A difaktorkan, maka vektor b bisa diganti nilainya sesuai dengan sistem persamaan linear yang lain, misalnya seluruh nilai di ruas kanan diganti menjadi

P1 P2 P3 P4

: : : :

x1 2x1 3x1 −x1

+ + − +

x2 x2 x2 2x2

− − +

x3 x3 3x3

+ + + −

3x4 x4 2x4 x4

= = = =

8 7 14 -7

4.2. ALGORITMA

65

Dalam operasi matrik menjadi      

1

1

0

3



x1





8

      1    x2  =  7    3 −1 −1 2    x3   14 −1 2 3 −1 x4 −7 2

1 −1

     

(4.8)

Perhatikan baik-baik! Matrik A sama persis dengan contoh sebelumnya. Perbedaannya hanya pada vektor b. Selanjutnya, dengan metode LU-decomposition, persamaan (4.8) menjadi      

1

0 0 0



1

1

0

3

  1 0 0    0 −1 −1 −5  3 4 1 0  0 3 13  0 −1 −3 0 1 0 0 0 −13 2



x1





8



      x2   7  =     x   14   3    x4 −7

(4.9)

Silakan anda lanjutkan proses perhitungannya dengan mencari vektor y sesuai contoh yang telah diberikan sebelumnya. Pada akhirnya akan diperoleh solusi sebagai berikut: x1 = 3, x2 = −1, x3 = 0, dan y4 = 2.

4.2

Algoritma

Sekarang saatnya saya tunjukkan algoritma metode LU decomposition. Algoritma ini dibuat untuk menyelesaikan sistem persamaan linear, dengan cara menfaktorkan matrik A = (aij ) berukuran n x n menjadi matrik L = (lij ) dan matrik U = (uij ) dengan ukuran yang sama. Algoritma LU-decomposition yang anda lihat sekarang merupakan modifikasi dari algoritma eliminasi gauss. Silakan anda periksa langkah-langkah mana saja yang telah mengalami modifikasi! Tapi asal tahu saja bahwa ini bukan satu-satunya algoritma untuk mendapatkan matrik LU. Sejauh yang saya tahu, ada algoritma lain untuk tujuan yang sama, dimana algoritma tersebut membutuhkan matrik permutasi untuk menggeser elemen pivot yang bernilai nol agar terhindar dari singular. Nah, sedangkan algoritma yang akan anda baca saat ini, sama sekali tidak “berurusan” dengan matrik permutasi. Algoritma ini cuma memanfaatkan “trik” tukar posisi yang sudah pernah dibahas di awal-awal catatan khususnya ketika membahas konsep eliminasi gauss. Satu lagi yang harus saya sampaikan juga adalah bahwa dalam algoritma ini, elemenelemen matrik L dan matrik U digabung jadi satu dan menggantikan seluruh elemen-elemen matrik A. Perhatian! cara ini jangan diartikan sebagai perkalian matrik L dan matrik U menjadi matrik A kembali. Cara ini dimaksudkan untuk menghemat memori komputer. Suatu aspek yang tidak boleh diabaikan oleh para programer. Marilah kita simak algoritmanya bersamasama! INPUT: dimensi n; nilai elemen aij , 1 ≤ i, j ≤ n; nilai elemen bi . OUTPUT: solusi x1 , x2 , x3 , ..., xn atau pesan kesalahan yang mengatakan bahwa faktorisasi

tidak mungkin dilakukan.


66

• Langkah 1: Inputkan konstanta-konstanta dari sistem persamaan linear kedalam elemenelemen matrik A dan vektor b, seperti berikut ini: 

a11

a12

. . . a1n

  a21 a22 . . . a2n A= .. ..  .. . .  . an1 an2 . . . ann





    

  b=  

b1 b2 .. . bn

     

(4.10)

• Langkah 2: Untuk i = 1, ..., n − 1, lakukan Langkah 3 sampai Langkah 5. • Langkah 3: Definisikan p sebagai integer dimana i ≤ p ≤ n. Lalu pastikan bahwa

api 6= 0. Langkah dilakukan bila ditemukan elemen diagonal yang bernilai nol (aii = 0). Ketika ada elemen diagonal yang bernilai nol, maka program harus mencari dan

memeriksa elemen-elemen yang tidak bernilai nol dalam kolom yang sama dengan kolom tempat elemen diagonal tersebut berada. Jadi saat proses ini berlangsung, integer i (indeks dari kolom) dibuat konstan, sementara integer p (indeks dari baris) bergerak dari p = i sampai p = n. Bila ternyata setelah mencapai elemen paling bawah dalam kolom tersebut, yaitu saat p = n tetap didapat nilai api = 0, maka sebuah pesan dimunculkan: sistem persamaan linear tidak memiliki solusi yang

unik. Lalu program berakhir: STOP. • Langkah 4: Namun jika sebelum integer p mencapai nilai p = n sudah diperoleh elemen yang tidak sama dengan nol (api 6= 0), maka bisa dipastikan p 6= i. Jika p 6= i

maka lakukan proses pertukaran (Pp ) ↔ (Pi ).

• Langkah 5: Untuk j = i + 1, .., n, lakukan Langkah 6 dan Langkah 7. • Langkah 6: Tentukan mji ,

mji =

aji aii

• Langkah 7: Lakukan proses triangularisasi, (Pj − mji Pi ) → (Pj ) • Langkah 8: Nilai mji disimpan ke aji , aji = mji • Langkah 9: Nilai b1 dicopy ke y1 , lalu lakukan substitusi-maju. y 1 = b1 Untuk i = 2, ..., n tentukan xi , y i = bi −

i−1 X j=1

aij yj

4.2. ALGORITMA

67

• Langkah 10: Lakukan proses substitusi-mundur, dimulai dengan menentukan xn , xn =

an,n+1 ann

Untuk i = n − 1, ..., 1 tentukan xi , xi =

ai,n+1 −

Pn

j=i+1 aij xj

aii

• Langkah 11: Diperoleh solusi yaitu x1 , x2 , ..., xn . Algoritma telah dijalankan dengan sukses. STOP.

Algoritma di atas telah diimplementasi kedalam program yang ditulis dengan bahasa Fortran. Program tersebut sudah berhasil dikompilasi dengan visual fortran (windows) dan g77 (debian-linux). Inilah programnya:

1 2 3 4 5 6

C

7 8 9 10 11 12 13 14 15

60

16 17 18 19

50

20 21

C

22 23 24 25

110

26 27

C

28 29 30

C

DIMENSION A(10,11), B(10), Y(10), X(10) REAL MJI WRITE(*,*) WRITE(*,*) ’==> FAKTORISASI MATRIK: LU DECOMPOSITION <==’ WRITE (*,*) LANGKAH 1: MEMASUKAN NILAI ELEMEN-ELEMEN MATRIK A DAN VEKTOR B WRITE (*,’(1X,A)’) ’JUMLAH PERSAMAAN ? ’ READ (*,*) N WRITE (*,*) WRITE (*,*) ’MASUKAN ELEMEN-ELEMEN MATRIK A’ DO 50 I = 1,N DO 60 J = 1,N WRITE (*,’(1X,A,I2,A,I2,A)’) ’A(’,I,’,’,J,’) = ’ READ (*,*) A(I,J) CONTINUE WRITE (*,’(1X,A,I2,A)’) ’B(’,I,’) ? ’ READ (*,*) B(I) WRITE (*,*) CONTINUE WRITE (*,*) MENAMPILKAN MATRIK A WRITE (*,’(1X,A)’) ’MATRIK A:’ DO 110 I = 1,N WRITE (*,6) (A(I,J),J=1,N) CONTINUE WRITE (*,*) LANGKAH 2: MEMERIKSA ELEMEN-ELEMEN PIVOT NN = N-1 DO 10 I=1,NN LANGKAH 3: MENDEFINISIKAN P

68 31 32

100

33 34 35 36

200 C

37 38 39 40

C

41 42 43 44 45 46

20

47 48

C

49 50 51

C

52 53

C

54 55 56 57

40 C

58 59 60 61

30 10 C

62 63 64 65

120

66 67

C

68 69 70 71 72 73

16

74 75 76 77

15 C

BAB 4. METODE LU DECOMPOSITION P = I IF (ABS(A(P,I)).GE.1.0E-20 .OR. P.GT.N) GOTO 200 P = P+1 GOTO 100 IF(P.EQ.N+1)THEN MENAMPILKAN PESAN TIDAK DAPAT DIFAKTORKAN WRITE(*,8) GOTO 400 END IF LANGKAH 4: PROSES TUKAR POSISI IF(P.NE.I) THEN DO 20 JJ=1,N C = A(I,JJ) A(I,JJ) = A(P,JJ) A(P,JJ) = C CONTINUE END IF LANGKAH 5: PERSIAPAN PROSES TRIANGULARISASI JJ = I+1 DO 30 J=JJ,N LANGKAH 6: TENTUKAN MJI MJI = A(J,I)/A(I,I) LANGKAH 7: PROSES TRIANGULARISASI DO 40 K=JJ,N A(J,K) = A(J,K)-MJI*A(I,K) CONTINUE LANGKAH 8: MENYIMPAN MJI KE A(J,I) A(J,I) = MJI CONTINUE CONTINUE MENAMPILKAN MATRIK LU WRITE (*,’(1X,A)’) ’MATRIK LU:’ DO 120 I = 1,N WRITE (*,6) (A(I,J),J=1,N) CONTINUE WRITE (*,*) LANGKAH 9: SUBSTITUSI-MAJU Y(1) = B(1) DO 15 I=2,N SUM = 0.0 DO 16 J=1,I-1 SUM = SUM+A(I,J)*Y(J) CONTINUE Y(I) = B(I)-SUM CONTINUE MENAMPILKAN VEKTOR Y WRITE (*,’(1X,A)’) ’VEKTOR Y:’

4.2. ALGORITMA 78 79 80

138

81 82

C

83 84 85 86 87 88 89 90

26

91 92 93

24 C

94 95 96 97

18

98 99 100 101 102 103 104

400 6 8

69

DO 138 I = 1,N WRITE (*,6) Y(I) CONTINUE WRITE (*,*) LANGKAH 10: SUBSTITUSI-MUNDUR X(N) = Y(N)/A(N,N) DO 24 K=1,N-1 I = N-K JJ = I+1 SUM = 0.0 DO 26 KK=JJ,N SUM = SUM+A(I,KK)*X(KK) CONTINUE X(I) = (Y(I)-SUM)/A(I,I) CONTINUE LANGKAH 11: MENAMPILKAN SOLUSI DAN SELESAI WRITE (*,’(1X,A)’) ’SOLUSI:’ DO 18 I = 1,N WRITE (*,’(1X,A,I2,A,F14.8)’) ’X(’,I,’) = ’,X(I) CONTINUE WRITE(*,*) WRITE(*,*) ’SELESAI --> SUKSES’ WRITE(*,*) CONTINUE FORMAT(1X,5(F14.8)) FORMAT(1X,’TIDAK DAPAT DIFAKTORKAN’) END

Demikianlah, sekarang kita punya tiga buah algoritma untuk memecahkan problem sistem persamaan linear, yaitu eliminasi gauss, invers matrik, dan lu-decomposition. Diantara ketiganya, eliminasi gauss adalah algoritma yang paling simpel dan efisien. Dia hanya butuh proses triangularisasi dan substitusi-mundur untuk mendapatkan solusi. Sedangkan dua algoritma yang lainnya membutuhkan proses-proses tambahan untuk mendapatkan solusi yang sama. Saya cukupkan sementara sampai disini. Insya Allah akan saya sambung lagi dilain waktu. Kalau ada yang mau didiskusikan, silakan hubungi saya melalui email.

Bab 5

Metode Iterasi

✍ Objektif : ⊲ Mengenalkan konsep Norm. ⊲ Mengenalkan iterasi Jacobi. ⊲ Mengenalkan iterasi Gauss-Seidel. ⊲ Mengenalkan iterasi Succesive-Over-Relaxation (SOR).

5.1

Kelebihan Vektor-kolom

Sebelum kita membahas metode iterasi untuk menyelesaikan problem sistem persamaan linear, saya ingin menyampaikan satu hal yang sangat sederhana, yaitu tentang cara merepresentasikan elemen-elemen suatu vektor-kolom. Sebagaimana tertulis pada catatan-catatan sebelumnya, biasanya suatu vektor-kolom ditulis sebagai 

x1



   x2   x=  ..   .  xn

(5.1)

Dengan operasi transpose, vektor-kolom tersebut dapat dinyatakan sebagai it h x = x1 x2 . . . xn Contoh:



3



  h it −2  = 3 −2 8 5 x= 8   5 71

(5.2)

BAB 5. METODE ITERASI

72

Cara penulisan seperti ini digunakan untuk menyatakan vektor-kolom pada suatu kalimat di dalam paragraf. Alasannya supaya tidak terlalu menyita banyak ruang penulisan. Sementara, persamaan (5.1), lebih sering digunakan pada penulisan operasi matrik. Satu hal lagi, pada paragraf-paragraf berikutnya, saya persingkat penulisan istilah vektor-kolom menjadi vektor saja.

5.2

Pengertian Norm

Vektor x=(x1 ; x2 ; ...; xn )t memiliki norm ℓ2 dan ℓ∞ yang didefinisikan sebagai n X

x2i }1/2

(5.3)

ℓ∞ = kxk∞ = max |xi |

(5.4)

ℓ2 = kxk2 = {

i=1

dan 1≤i≤n

Contoh: x=(3; −2; 8; 5)t memiliki norm ℓ2 yaitu ℓ2 = kxk2 = dan norm ℓ∞ yaitu

p (3)2 + (−2)2 + (8)2 + (5)2 = 10, 0995

ℓ∞ = kxk∞ = max{(3), (−2), (8), (5)} = 8 Saya menyarankan agar kedua norm ini diingat-ingat dengan baik, karena akan banyak disinggung pada catatan-catatan berikutnya. 5.2.1 Script perhitungan norm dalam Matlab Script berikut ini merujuk pada contoh di atas, dimana vektor x hanya terdiri dari 4 elemen, yaitu x(1, 1),x(2, 1),x(3, 1) dan x(4, 1)

1 2 3 4 5 6 7 8 9 10 11 12 13 14

clear all clc x(1,1)=3; x(2,1)=-2; x(3,1)=8; x(4,1)=5; x %menampilkan vektor x %=========menghitung norm2============= s=0; for i=1:4 s=s+x(i,1)^2; end norm2=sqrt(s) %menampilkan hasil norm2 %======================================

5.3. ITERASI JACOBI

73

Mohon diperhatikan untuk mengganti angka 4 pada statemen for i=1:4 dengan angka yang lain disesuaikan dengan jumlah elemen vektor yang mau dihitung norm2-nya. 5.2.2 Perhitungan norm-selisih Misalnya kita punya vektor bernama xlama. Lalu ada vektor lainnya bernama xbaru. Norm selisih dari xlama dan xbaru dapat dihitung dengan bantuan script berikut ini

1 2 3 4 5 6 7

clear all clc xlama(1,1)=3; xlama(2,1)=-2; xlama(3,1)=8; xlama(4,1)=5; xlama %menampilkan elemen vektor xlama

8 9 10 11 12 13

xbaru(1,1)=9; xbaru(2,1)=4; xbaru(3,1)=6; xbaru(4,1)=1; xbaru %menampilkan elemen vektor xbaru

14 15

n=4;

%jumlah elemen vektor

16 17 18 19 20 21 22 23

%--------menghitung norm2 selisih ------------s=0; for i=1:n s=s+(xbaru(i,1)-xlama(i,1))^2; end norm2=sqrt(s) %----------------------------------------------

Cara perhitungan norm-selisih seperti ini akan diterapkan pada kebanyakan metode iterasi. Jadi tolong diingat baik-baik!!

5.3

Iterasi Jacobi

Sekarang kita mulai pembahasan tentang metode iterasi untuk menyelesaikan problem sistem persamaan linear. Metode ini berbeda dengan metode-metode yang telah dijelaskan sebelumnya, dimana ia dimulai dari penentuan nilai awal (initial value) untuk setiap elemen vektor x. Kemudian berdasarkan nilai awal tersebut, dilakukan langkah perhitungan untuk mendapatkan elemen-elemen vektor x yang baru. x(baru) = T x(lama) + c

(5.5)


74 atau xk = T xk−1 + c

(5.6)

dimana k = 1, 2, 3, ..., n. Mungkin anda kesulitan memahami rumus-rumus di atas, maka untuk lebih jelasnya, marilah kita perhatikan contoh berikut; diketahui sistem persamaan linear Ax = b yaitu 10x1 − x2 + 2x3 = 6 −x1 + 11x2 − x3 + 3x4 = 25 2x1 − x2 + 10x3 − x4 = −11 3x2 − x3 + 8x4 = 15 yang mana solusinya adalah x=(1; 2; −1; 1)t . Silakan simpan dulu solusi ini, anggap saja kita belum tahu. Lalu perhatikan baik-baik bagaimana metode iterasi Jacobi bisa menemukan so-

lusi tersebut dengan caranya yang khas. Langkah pertama dan merupakan langkah terpenting dari metode iterasi Jacobi adalah dengan mengubah cara penulisan sistem persamaan linear di atas menjadi seperti ini 2 6 1 x2 − x3 + 10 10 10 1 1 3 25 = x1 + x3 − x4 + 11 11 11 11 1 1 11 2 = − x1 + x2 + x4 − 10 10 10 10 3 1 15 = − x2 + x3 + 8 8 8

x1 = x2 x3 x4

Kita bisa menyatakan bahwa nilai x1 , x2 , x3 dan x4 yang berada di ruas kiri tanda = (baca: sama dengan) sebagai x(baru) . Sementara nilai x1 , x2 , x3 dan x4 yang berada di ruas kanan tanda = (baca: sama dengan) sebagai x(lama) . Sehingga sistem persamaan tersebut ditulis seperti ini (baru)

x1

(baru)

x2

(baru)

x3

(baru)

x4

2 (lama) 6 1 (lama) x2 − x3 + 10 10 10 1 (lama) 1 (lama) 3 (lama) 25 = x + x3 − x4 + 11 1 11 11 11 2 (lama) 1 1 (lama) 11 = − x1 + x2 + x4 − 10 10 10 10 3 (lama) 1 (lama) 15 = − x2 + x3 + 8 8 8

=

5.3. ITERASI JACOBI

75

yang secara umum dapat diformulasikan sebagaimana persamaan (5.5). Atau dapat pula ditulis seperti ini 2 (k−1) 6 1 (k−1) x − x3 + 10 2 10 10 1 (k−1) 1 (k−1) 3 (k−1) 25 = x1 + x3 − x4 + 11 11 11 11 2 (k−1) 1 (k−1) 1 (k−1) 11 = − x1 + x2 + x4 − 10 10 10 10 3 (k−1) 1 (k−1) 15 = − x2 + x3 + 8 8 8

(k)

x1

=

(k)

x2

(k)

x3

(k)

x4

yang secara umum dapat diformulasikan sebagaimana persamaan (5.6). Pada persamaan (5.6), indeks k menunjukan jumlah berapa kali perhitungan iterasi telah dilakukan. Mari kita fokuskan sejenak pada indeks k ini; Pada k = 1, maka penulisan sistem persamaan linear menjadi seperti ini (1)

x1

(1)

x2

(1)

x3

(1)

x4

2 (0) 6 1 (0) x2 − x3 + 10 10 10 1 (0) 1 (0) 3 (0) 25 = x + x3 − x4 + 11 1 11 11 11 1 (0) 1 (0) 11 2 (0) = − x1 + x2 + x4 − 10 10 10 10 3 (0) 1 (0) 15 = − x2 + x3 + 8 8 8 =

(0)

(0)

(0)

(0)

Jika kita tentukan nilai-nilai awal x(0) sebagai berikut x1 = 0, x2 = 0, x3 = 0 dan x4 = 0. Atau dinyatakan seperti ini x(0) = (0; 0; 0; 0)t . Maka kita akan memperoleh nilai-nilai x(1) , yaitu hasil iterasi pertama, sebagai berikut (1)

x1

(1)

x2

(1)

x3

(1)

x4

6 10 25 = 11 11 = − 10 15 = 8 =

atau x(1) = (0, 6000; 2, 2727; −1, 1000; 1, 8750)t . Setelah diperoleh nilai-nilai x(1) , perhitungan

tersebut diulangi kembali guna mendapatkan hasil iterasi kedua, dimana nilai k = 2. Caranya

adalah dengan memasukan nilai-nilai x(1) = (0, 6000; 2, 2727; −1, 1000; 1, 8750)t ke ruas kanan


76 tanda sama-dengan, (2)

x1

(2)

x2

(2)

x3

(2)

x4

2 (1) 6 1 (1) x − x3 + 10 2 10 10 1 (1) 1 (1) 3 (1) 25 = x1 + x3 − x4 + 11 11 11 11 2 (1) 1 (1) 1 (1) 11 = − x1 + x2 + x4 − 10 10 10 10 3 (1) 1 (1) 15 = − x2 + x3 + 8 8 8

=

maka nilai-nilai x(2) yang kita dapat adalah x(2) = (1, 0473; 1, 7159; −0, 8052; 0, 8852)t . Setelah diperoleh nilai-nilai x(2) , perhitungan tersebut diulangi kembali guna mendapatkan hasil

iterasi ketiga, dimana nilai k = 3. Caranya adalah dengan memasukan nilai-nilai x(2) = (1, 0473; 1, 7159; −0, 8052; 0, 8852)t ke ruas kanan kembali, (3)

x1

(3)

x2

(3)

x3

(3)

x4

1 (2) 2 (2) 6 x2 − x3 + 10 10 10 1 (2) 1 (2) 3 (2) 25 = x + x3 − x4 + 11 1 11 11 11 1 (2) 1 (2) 11 2 (2) = − x1 + x2 + x4 − 10 10 10 10 3 (2) 1 (2) 15 = − x2 + x3 + 8 8 8

=

maka kita akan memperoleh nilai-nilai x(3) = (0, 9326; 2, 0530; −1, 0493; 1, 1309)t . Lalu proses

perhitungan diulangi lagi dengan k = 4. Begitulah seterusnya. Proses ini diulangi lagi berkalikali untuk nilai-nilai k berikutnya. Proses yang berulang ini disebut proses iterasi. Sampai dengan x(3) di atas, kita sudah melakukan tiga kali proses iterasi. Lantas sampai kapan proses iterasi ini terus berlanjut? Jawabnya adalah sampai x(baru) mendekati solusi yang sesungguhnya, yaitu x = (1; 2; −1; 1)t Dengan kata lain, proses iterasi harus dihentikan bila x(baru) sudah mendekati solusi. Lalu kriteria apa yang digunakan sehingga suatu hasil iterasi bisa dikatakan paling dekat dengan solusi yang sebenarnya? OK, simpan dulu pertanyaan ini, sebagai gantinya marilah kita pelajari script Matlab untuk menghitung iterasi. 5.3.1 Script Matlab untuk menghitung iterasi Pertama-tama kita buat script seperti ini

1 2 3 4 5 6

clear all clc %----nilai awal----------xlama(1,1)=0; xlama(2,1)=0; xlama(3,1)=0;

5.3. ITERASI JACOBI 7 8

77

xlama(4,1)=0; xlama

9 10 11 12 13 14 15

%------nilai baru------------xbaru(1,1)=(1/10)*xlama(2,1)-(2/10)*xlama(3,1)+(6/10); xbaru(2,1)=(1/11)*xlama(1,1)+(1/11)*xlama(3,1)-(3/11)*xlama(4,1)+(25/11); xbaru(3,1)=-(2/10)*xlama(1,1)+(1/10)*xlama(2,1)+(1/10)*xlama(4,1)-(11/10); xbaru(4,1)=-(3/8)*xlama(2,1)+(1/8)*xlama(3,1)+(15/8); xbaru

xbaru yang didapat tak lain adalah hasil iterasi pertama, yaitu x(1) = (0, 6000; 2, 2727; −1, 1000; 1, 8750)t . Kemudian, untuk iterasi ke-2, script di atas dimodifikasi menjadi seperti ini

1 2

clear all clc

3 4 5 6 7 8 9 10 11 12 13 14 15

%----nilai awal----------xlama(1,1)=0; xlama(2,1)=0; xlama(3,1)=0; xlama(4,1)=0; xlama %------nilai baru------------xbaru(1,1)=(1/10)*xlama(2,1)-(2/10)*xlama(3,1)+(6/10); xbaru(2,1)=(1/11)*xlama(1,1)+(1/11)*xlama(3,1)-(3/11)*xlama(4,1)+(25/11); xbaru(3,1)=-(2/10)*xlama(1,1)+(1/10)*xlama(2,1)+(1/10)*xlama(4,1)-(11/10); xbaru(4,1)=-(3/8)*xlama(2,1)+(1/8)*xlama(3,1)+(15/8); xbaru

16 17

xlama=xbaru; %xbaru dijadikan xlama untuk iterasi berikutnya

18 19 20 21 22 23

xbaru(1,1)=(1/10)*xlama(2,1)-(2/10)*xlama(3,1)+(6/10); xbaru(2,1)=(1/11)*xlama(1,1)+(1/11)*xlama(3,1)-(3/11)*xlama(4,1)+(25/11); xbaru(3,1)=-(2/10)*xlama(1,1)+(1/10)*xlama(2,1)+(1/10)*xlama(4,1)-(11/10); xbaru(4,1)=-(3/8)*xlama(2,1)+(1/8)*xlama(3,1)+(15/8); xbaru

Sampai disini, xbaru yang didapat adalah hasil iterasi ke-2, yaitu x(2) = (1, 0473; 1, 7159; − 0, 8052; 0, 8852)t . Kemudian, untuk iterasi ke-3, script di atas dimodifikasi menjadi seperti ini 1 2

clear all clc

3 4 5 6

%----nilai awal----------xlama(1,1)=0; xlama(2,1)=0;

78 7 8 9 10 11 12 13 14 15


xlama(3,1)=0; xlama(4,1)=0; xlama %------nilai baru------------xbaru(1,1)=(1/10)*xlama(2,1)-(2/10)*xlama(3,1)+(6/10); xbaru(2,1)=(1/11)*xlama(1,1)+(1/11)*xlama(3,1)-(3/11)*xlama(4,1)+(25/11); xbaru(3,1)=-(2/10)*xlama(1,1)+(1/10)*xlama(2,1)+(1/10)*xlama(4,1)-(11/10); xbaru(4,1)=-(3/8)*xlama(2,1)+(1/8)*xlama(3,1)+(15/8); xbaru

16 17


18 19 20 21 22 23


24 25


26 27 28 29 30 31


Sampai disini, xbaru yang didapat adalah hasil iterasi ke-3, yaitu x(3) = (0, 9326; 2, 0530; − 1, 0493; 1, 1309)t . Kemudian, untuk iterasi ke-4, script di atas dimodifikasi dengan cara yang

sama. Tapi konsekuensinya script tersebut akan bertambah panjang. Dan itu sesuatu yang

tidak baik. Guna menghindari hal itu, script di atas perlu dioptimasi dengan pasangan for-end sebagai berikut

1 2

clear all clc

3 4 5 6 7 8 9

%----nilai awal----------xlama(1,1)=0; xlama(2,1)=0; xlama(3,1)=0; xlama(4,1)=0; xlama

10 11 12 13 14

for i=1:4 %------nilai update------------xbaru(1,1)=(1/10)*xlama(2,1)-(2/10)*xlama(3,1)+(6/10); xbaru(2,1)=(1/11)*xlama(1,1)+(1/11)*xlama(3,1)-(3/11)*xlama(4,1)+(25/11);

5.3. ITERASI JACOBI 15 16 17

79

xbaru(3,1)=-(2/10)*xlama(1,1)+(1/10)*xlama(2,1)+(1/10)*xlama(4,1)-(11/10); xbaru(4,1)=-(3/8)*xlama(2,1)+(1/8)*xlama(3,1)+(15/8); xbaru

18 19 20

xlama=xbaru; %xbaru dijadikan xlama untuk iterasi berikutnya end

Angka 4 pada statemen for i=1:4 dapat diganti sesuai dengan jumlah iterasi maksimal yang kita kehendaki. 1 2

clear all clc

3 4 5 6 7 8 9


10 11

itermaks=4

%jumlah iterasi maksimal

12 13 14 15 16 17 18 19

for i=1:itermaks %------nilai update------------xbaru(1,1)=(1/10)*xlama(2,1)-(2/10)*xlama(3,1)+(6/10); xbaru(2,1)=(1/11)*xlama(1,1)+(1/11)*xlama(3,1)-(3/11)*xlama(4,1)+(25/11); xbaru(3,1)=-(2/10)*xlama(1,1)+(1/10)*xlama(2,1)+(1/10)*xlama(4,1)-(11/10); xbaru(4,1)=-(3/8)*xlama(2,1)+(1/8)*xlama(3,1)+(15/8); xbaru

20 21 22


Untuk mendapatkan hasil iterasi yang ke-10, silakan nyatakan itermaks=10 pada script di atas. Hasil dari keseluruhan iterasi dari iterasi ke-1 hingga iterasi ke-10 disajikan pada tabel berikut

k (k)

x1 (k) x2 (k) x3 (k) x4

Tabel 5.1: Hasil akhir elemen-elemen vektor x hingga iterasi ke-10 0 1 2 3 4 ... 9 10 0,0000 0,0000 0,0000 0,0000

0,6000 2,2727 -1,1000 1,8852

1,0473 1,7159 -0,8052 0,8852

0,9326 2,0530 -1,0493 1,1309

1,0152 1,9537 -0,9681 0,9739

... ... ... ...

0,9997 2,0004 -1,0004 1,0006

1,0001 1,9998 -0,9998 0,9998

Kita bisa saksikan bahwa hasil iterasi ke-1, x(1) = (0, 6000; 2, 2727; −1, 1000; 1, 8852) adalah

hasil yang paling tidak mendekati solusi, x = (1; 2; −1; 1)t . Dibandingkan dengan hasil iterasi

80


ke-2, jelas terlihat bahwa hasil iterasi ke-2 lebih mendekati solusi. Kalau terus diurutkan, maka hasil iterasi ke-10 merupakan hasil yang paling dekat dengan solusi. 5.3.2 Optimasi script Matlab untuk menghitung iterasi Sekarang mari kita hitung norm-selisih dari masing-masing hasil iterasi secara berurutan. Dimulai dari mencari norm-selisih antara hasil iterasi ke-1 dan ke-2. Lalu dilanjutkan dengan hasil iterasi ke-2 dan ke-3, begitu seterusnya hingga antara hasil iterasi yang ke-9 dan ke-10. Dalam prakteknya, kita cukup menambahkan script norm-selisih pada script yang tadi 1 2

clear all clc

3 4 5 6 7 8 9 10 11

%----nilai awal----------xlama(1,1)=0; xlama(2,1)=0; xlama(3,1)=0; xlama(4,1)=0; xlama n=4 %jumlah elemen vektor itermaks=10 %jumlah iterasi maksimal

12 13 14 15 16 17 18 19


20 21 22 23 24 25 26 27

%------norm selisih------------s=0; for i=1:n s=s+(xbaru(i,1)-xlama(i,1))^2; end epsilon=sqrt(s) %-------------------------------

28 29 30


Tabel dibawah ini memperlihatkan hasil norm-selisih hingga iterasi ke-10. Hasil perhitungan norm-selisih tersebut, saya beri nama epsilon, ǫ, dimana semakin kecil nilai epsilon, ǫ, menandakan hasil iterasinya semakin dekat dengan solusi. Hasil norm-selisih yang semakin kecil pada iterasi ke-10 menunjukan bahwa hasil iterasi ke-10 adalah hasil yang paling dekat dengan solusi yang sebenarnya.

5.3. ITERASI JACOBI

81

Tabel 5.2: °Hasil perhitungan ℓ2°) hingga °iterasi ke-10° ° ° (3)norm-selisih ° °(dengan (2) (1) (2) (4) (3) ° ° ° ° ° norm ℓ2 x −x x −x x − x °2 ... °x(10) − x(9) °2 2 2 ǫ 1,2557 0,4967 0,2189 ... 0,0012 Kembali ke pertanyaan penting yang tadi yaitu kriteria apa yang digunakan sehingga suatu hasil iterasi bisa dikatakan paling dekat dengan solusi yang sebenarnya? Jawabnya: tergantung besar kecilnya nilai ǫ. Artinya kalau nilai ǫ ditentukan sebesar 0,2 , maka iterasi akan berhenti pada iterasi ke-4. Atau kalau nilai ǫ ditentukan sebesar 0,001 , maka proses iterasi akan berhenti pada iterasi ke-10. Kesimpulannya, semakin kecil nilai ǫ, semakin panjang proses iterasinya, namun hasil akhirnya semakin akurat. Jadi nilai ǫ berperan penting untuk menghentikan proses iterasi. Dalam hal ini, ǫ lebih umum dikenal dengan istilah stopping-criteria. Di bawah ini adalah script iterasi Jacobi setelah mengalami optimasi beberapa kali, 1 2

clear all clc

3 4 5 6 7 8 9


10 11 12 13

n=4 itermaks=10 sc=0.001

%jumlah elemen vektor %jumlah iterasi maksimal %stopping-criteria

14 15 16 17 18 19 20 21


22 23 24 25 26 27 28

%------norm selisih------------s=0; for i=1:n s=s+(xbaru(i,1)-xlama(i,1))^2; end epsilon=sqrt(s)

29 30 31 32

%------memeriksa stopping criteria, sc-------if epsilon<sc break


82 33

end

34 35 36


Metode yang baru saja kita bahas ini disebut metode Iterasi Jacobi. Metode ini bertujuan mencari nilai-nilai pengganti variabel-variabel x dengan perumusan (k)

xi

=

Pn

j=1

³

(k−1)

−aij xj aii

´

+ bi

dimana i=1,2,3,...,n. 5.3.3 Algoritma • Langkah 1: Tentukan k=1 • Langkah 2: Ketika (k ≤ N ) lakukan Langkah 3-6 – Langkah 3: Untuk i=1,...,n, hitunglah xi =

−

Pn

j=1 (aij XOj )

+ bi

aii

– Langkah 4: Jika kx − XOk < ǫ, maka keluarkan OUTPUT (x1 , ..., xn ) lalu STOP – Langkah 5: Tentukan k=k+1 – Langkah 6: Untuk i=1,...n, tentukan XOi = xi • Langkah 7: OUTPUT (’Iterasi maksimum telah terlampaui’) lalu STOP 5.3.4 Program dalam Fortran 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

62

IMPLICIT NONE DIMENSION A(10,10),B(10),X(10),XO(10) REAL A,B,X,XO,EPS,NORM,S INTEGER N,I,J,K,ITMAX WRITE(*,*) ’==> ITERASI JACOBI UNTUK SISTEM LINEAR <==’ WRITE(*,*) WRITE (*,’(1X,A)’) ’JUMLAH PERSAMAAN ? ’ READ (*,*) N WRITE (*,*) ’MASUKAN ELEMEN-ELEMEN MATRIK A DAN VEKTOR B’ DO 52 I = 1,N DO 62 J = 1,N WRITE (*,’(1X,A,I2,A,I2,A)’) ’A(’,I,’,’,J,’) = ’ READ (*,*) A(I,J) CONTINUE WRITE (*,’(1X,A,I2,A)’) ’B(’,I,’) ? ’

(5.7)

5.3. ITERASI JACOBI 16 17 18

52

19 20 21 22 23 24 25 26 27

72

28 29

C

30 31 32 33

110

34 35

C

36 37 38 39

111

40 41

C

42 43 44 45

C 100 C

46 47 48 49 50 51

20

52 53 54 55

10

56 57 58 59 60 61 62

C

READ (*,*) B(I) WRITE (*,*) CONTINUE WRITE (*,’(1X,A)’) ’JUMLAH ITERASI MAKSIMUM ? ’ READ (*,*) ITMAX WRITE (*,’(1X,A)’) ’NILAI EPSILON ATAU TOLERANSI ? ’ READ (*,*) EPS WRITE (*,*) ’MASUKAN NILAI AWAL UNTUK XO’ DO 72 I = 1,N WRITE (*,’(1X,A,I2,A)’) ’XO(’,I,’) ? ’ READ (*,*) XO(I) CONTINUE WRITE (*,*) MENAMPILKAN MATRIK A WRITE (*,’(1X,A)’) ’MATRIK A:’ DO 110 I = 1,N WRITE (*,6) (A(I,J),J=1,N) CONTINUE WRITE (*,*) MENAMPILKAN VEKTOR B WRITE (*,’(1X,A)’) ’VEKTOR B:’ DO 111 I = 1,N WRITE (*,6) B(I) CONTINUE WRITE (*,*) LANGKAH 1 K = 1 LANGKAH 2 IF(K.GT.ITMAX) GOTO 200 LANGKAH 3 NORM = 0.0 DO 10 I = 1,N S = 0.0 DO 20 J=1,N S = S-A(I,J)*XO(J) CONTINUE S = (S+B(I))/A(I,I) IF (ABS(S).GT.NORM) NORM=ABS(S) X(I) = XO(I)+S CONTINUE WRITE(*,’(1X,A,I3)’) ’ITERASI KE-’, K WRITE(*,’(1X,A,F14.8)’) ’NORM = ’, NORM WRITE(*,’(1X,A,I3,A,F14.8)’) (’X(’,I,’) = ’, X(I),I=1,N) WRITE(*,*) LANGKAH 4 IF(NORM.LE.EPS) THEN WRITE(*,7) K,NORM

83


84 63 64 65

C

66 67

C

68 69 70

30

71 72 73

C 200

74 75

400

GOTO 400 END IF LANGKAH 5 K = K+1 LANGKAH 6 DO 30 I=1,N XO(I) = X(I) CONTINUE GOTO 100 LANGKAH 7 CONTINUE WRITE(*,9) STOP

76 77 78 79

5 6 7

80 81

9

82

5.4

FORMAT(1X,I3) FORMAT(1X,(6(1X,F14.8))) FORMAT(1X,’KONVERGEN PADA ITERASI YANG KE- ’,I3, *’ , NORM= ’,F14.8) FORMAT(1X,’MELEBIHI BATAS MAKSIMUM ITERASI’) END

Iterasi Gauss-Seidel

Metode Iterasi Gauss-Seidel hampir sama dengan metode Iterasi Jacobi. Perbedaannya hanya terletak pada penggunaan nilai elemen vektor xbaru yang langsung digunakan pada persamaan dibawahnya. Untuk lebih jelasnya, perhatikan sistem persamaan linear berikut, yang diturunkan dari contoh terdahulu 2 6 1 lama x2 − xlama + 3 10 10 10 1 baru 1 lama 3 25 = x1 + x3 − xlama + 4 11 11 11 11 1 baru 1 lama 11 2 baru − = − x1 + x2 + x4 10 10 10 10 3 baru 1 baru 15 = − x2 + x3 + 8 8 8

= xbaru 1 xbaru 2 xbaru 3 xbaru 4

Pada baris pertama, xbaru dihitung berdasarkan xlama dan xlama . Kemudian xbaru tersebut 1 2 3 1 langsung dipakai pada baris kedua untuk menghitung xbaru . Selanjutnya xbaru dan xbaru di2 1 2 gunakan pada baris ketiga untuk mendapatkan xbaru . Begitu seterusnya hingga xbaru pun 3 4 diperoleh pada baris keempat. Sistem persamaan tersebut dapat dinyatakan dalam indeks k

5.4. ITERASI GAUSS-SEIDEL

85

seperti dibawah ini dimana k adalah jumlah iterasi. 2 (k−1) 6 1 (k−1) x − x3 + 10 2 10 10 1 (k) 1 (k−1) 3 (k−1) 25 = x1 + x3 − x4 + 11 11 11 11 2 (k) 1 (k) 1 (k−1) 11 = − x1 + x2 + x4 − 10 10 10 10 3 (k) 1 (k) 15 = − x2 + x3 + 8 8 8

(k)

x1

=

(k)

x2

(k)

x3

(k)

x4

(0)

(0)

(0)

Misalnya kita tentukan nilai-nilai awal x(0) sebagai berikut x1 = 0, x2 = 0, x3 = 0 dan

(0)

x4 = 0. Atau dinyatakan seperti ini x(0) = (0; 0; 0; 0)t . Maka pada k = 1 kita akan memperoleh nilai-nilai x(1) sebagai berikut x1

(1)

= 0, 6000

(1) x2 (1) x3 (1) x4

= 2, 3272 = −0, 9873 = 0, 8789

Lalu proses perhitungan diulangi lagi dengan k = 2. Begitu seterusnya proses ini diulangulang lagi untuk nilai-nilai k berikutnya sampai x(k) mendekati solusi yang sesungguhnya, yaitu x = (1; 2; −1; 1)t Marilah kita amati hasil seluruh iterasi. Tabel di bawah ini menampilkan hasil perhitungan hingga iterasi yang ke-5. Kita bisa saksikan bahwa dibandingkan dengan iterasi Jacobi, problem sistem persamaan linear yang sama, bisa diselesaikan oleh metode iterasi Gauss-Seidel hanya dalam 5 kali iterasi. Dari kasus ini, bisa kita simpulkan bahwa iterasi Gauss-Seidel bek-

k (k)

x1 (k) x2 (k) x3 (k) x4

0

Tabel 5.3: Hasil Iterasi Gauss-Seidel 1 2 3 4

0,0000 0,0000 0,0000 0,0000

0,6000 2,3272 -0,9873 0,8789

1,030 2,037 -1,014 0,9844

1,0065 2,0036 -1,0025 0,9983

1,0009 2,0003 -1,0003 0,9999

5 1,0001 2,0000 -1,0000 1,0000

erja lebih efektif dibandingkan iterasi Jacobi. Ya.., memang secara umum demikian, akan tetapi ternyata ditemukan kondisi yang sebaliknya pada kasus-kasus yang lain. 5.4.1 Script iterasi Gauss-Seidel Secara umum, script iterasi Gauss-Seidel yang saya tuliskan disini hampir sama dengan iterasi Jacobi. Perbedaan kecil-nya terletak pada bagian nilai update, dimana elemen xbaru hasil perhitungan dilibatkan langsung untuk menghitung elemen xbaru selanjutnya.


86

1 2

clear all clc

3 4 5 6 7 8 9


10 11 12 13

n=4 itermaks=10 sc=0.001

%jumlah elemen vektor %jumlah iterasi maksimal %stopping-criteria

14 15 16 17 18 19 20 21

for i=1:itermaks %------nilai update------------xbaru(1,1)=(1/10)*xlama(2,1)-(2/10)*xlama(3,1)+(6/10); xbaru(2,1)=(1/11)*xbaru(1,1)+(1/11)*xlama(3,1)-(3/11)*xlama(4,1)+(25/11); xbaru(3,1)=-(2/10)*xbaru(1,1)+(1/10)*xbaru(2,1)+(1/10)*xlama(4,1)-(11/10); xbaru(4,1)=-(3/8)*xbaru(2,1)+(1/8)*xbaru(3,1)+(15/8); xbaru

22 23 24 25 26 27 28

%------norm selisih------------s=0; for i=1:n s=s+(xbaru(i,1)-xlama(i,1))^2; end epsilon=sqrt(s)

29 30 31 32 33

%------memeriksa stopping criteria, sc-------if epsilon<sc break end

34 35 36


Perumusan metode Iterasi Gauss-Seidel dapat dinyatakan sebagai berikut: (k)

xi dimana i=1,2,3,...,n.

=

−

Pi−1 ³ j=1

(k)

aij xj

´

−

Pn

j=i+1

aii

³

(k−1)

aij xj

´

+ bi

(5.8)

5.4. ITERASI GAUSS-SEIDEL

87

5.4.2 Algoritma • Langkah 1: Tentukan k=1 • Langkah 2: Ketika (k ≤ N ) lakukan Langkah 3-6 – Langkah 3: Untuk i=1,...,n, hitunglah xi =

−

Pi−1

j=1 aij xj

−

Pn

j=i+1 aij XOj

+ bi

aii

– Langkah 4: Jika kx − XOk < ǫ, maka keluarkan OUTPUT (x1 , ..., xn ) lalu STOP – Langkah 5: Tentukan k=k+1 – Langkah 6: Untuk i=1,...n, tentukan XOi = xi • Langkah 7: OUTPUT (’Iterasi maksimum telah terlampaui’) lalu STOP 5.4.3 Script iterasi Gauss-Seidel dalam Fortran 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

62

16 17 18 19

52

20 21 22 23 24 25 26 27 28

72

IMPLICIT NONE DIMENSION A(10,10),B(10),X(10),XO(10) REAL A,B,X,XO,EPS,NORM,S1,S2 INTEGER N,I,J,K,ITMAX WRITE(*,*) WRITE(*,*) ’==> ITERASI GAUSS-SEIDEL UNTUK SISTEM LINEAR <==’ WRITE(*,*) WRITE (*,’(1X,A)’) ’JUMLAH PERSAMAAN ? ’ READ (*,*) N WRITE (*,*) ’MASUKAN ELEMEN-ELEMEN MATRIK A DAN VEKTOR B’ DO 52 I = 1,N DO 62 J = 1,N WRITE (*,’(1X,A,I2,A,I2,A)’) ’A(’,I,’,’,J,’) = ’ READ (*,*) A(I,J) CONTINUE WRITE (*,’(1X,A,I2,A)’) ’B(’,I,’) ? ’ READ (*,*) B(I) WRITE (*,*) CONTINUE WRITE (*,’(1X,A)’) ’JUMLAH ITERASI MAKSIMUM ? ’ READ (*,*) ITMAX WRITE (*,’(1X,A)’) ’NILAI EPSILON ATAU TOLERANSI ? ’ READ (*,*) EPS WRITE (*,*) ’MASUKAN NILAI AWAL UNTUK XO’ DO 72 I = 1,N WRITE (*,’(1X,A,I2,A)’) ’XO(’,I,’) ? ’ READ (*,*) XO(I) CONTINUE

88 29 30

C

31 32 33 34

110

35 36

C

37 38 39 40

111

41 42

C

43 44 45 46

C 100 C

47 48 49 50 51

20

52 53 54 55

23

56 57 58

10 C

59 60 61 62

40

63 64 65 66 67 68

C

69 70 71 72 73

C

74 75

C

BAB 5. METODE ITERASI WRITE (*,*) MENAMPILKAN MATRIK A WRITE (*,’(1X,A)’) ’MATRIK A:’ DO 110 I = 1,N WRITE (*,6) (A(I,J),J=1,N) CONTINUE WRITE (*,*) MENAMPILKAN VEKTOR B WRITE (*,’(1X,A)’) ’VEKTOR B:’ DO 111 I = 1,N WRITE (*,6) B(I) CONTINUE WRITE (*,*) LANGKAH 1 K = 1 LANGKAH 2 IF(K.GT.ITMAX) GOTO 200 LANGKAH 3 DO 10 I = 1,N S1 = 0.0 DO 20 J=I+1,N S1 = S1-A(I,J)*XO(J) CONTINUE S2 = 0.0 DO 23 J=1,I-1 S2 = S2-A(I,J)*X(J) CONTINUE X(I) = (S2+S1+B(I))/A(I,I) CONTINUE SAYA PILIH NORM-2. ANDA BOLEH PAKAI NORM YANG LAIN! NORM = 0.0 DO 40 I=1,N NORM = NORM + (X(I)-XO(I))*(X(I)-XO(I)) CONTINUE NORM = SQRT(NORM) WRITE(*,’(1X,A,I3)’) ’ITERASI KE-’, K WRITE(*,’(1X,A,F14.8)’) ’NORM-2 = ’, NORM WRITE(*,’(1X,A,I3,A,F14.8)’) (’X(’,I,’) = ’, X(I),I=1,N) WRITE(*,*) LANGKAH 4 IF(NORM.LE.EPS) THEN WRITE(*,7) K,NORM GOTO 400 END IF LANGKAH 5 K = K+1 LANGKAH 6

5.5. ITERASI DENGAN RELAKSASI 76 77 78

30

79 80 81

C 200

82 83

400

89

DO 30 I=1,N XO(I) = X(I) CONTINUE GOTO 100 LANGKAH 7 CONTINUE WRITE(*,9) STOP

84 85 86 87

5 6 7

88 89

9

90

5.5

FORMAT(1X,I3) FORMAT(1X,(6(1X,F14.8))) FORMAT(1X,’KONVERGEN PADA ITERASI YANG KE- ’,I3, *’ , NORM= ’,F14.8) FORMAT(1X,’MELEBIHI BATAS MAKSIMUM ITERASI’) END

Iterasi dengan Relaksasi

Metode Iterasi Relaksasi (Relaxation method ) dinyatakan dengan rumus berikut: (k)

xi

(k−1)

= (1 − ω) xi

+



ω  bi − aii

i−1 X j=1

(k)

aij xj −

n X

j=i+1



(k−1) 

aij xj

(5.9)

dimana i=1,2,3,...,n.

Untuk lebih jelasnya, marilah kita perhatikan contoh berikut, diketahui sistem persamaan linear Ax = b yaitu 4x1 + 3x2 + = 24 3x1 + 4x2 − x3 = 30 −x2 + 4x3 = −24 memiliki solusi (3, 4, −5)t . Metode Gauss-Seidel dan Relaksasi dengan ω = 1, 25 akan digu-

nakan untuk menyelesaikan sistem persamaan linear di atas dengan x(0) = (1, 1, 1)t . Untuk setiap nilai k = 1, 2, 3, ..., persamaan Gauss-Seidelnya adalah (k)

x1

(k)

x2

(k) x3

(k−1)

= −0, 75x2

(k)

+6 (k−1)

= −0, 75x1 + 0, 25x3 =

(k) 0, 25x2

−6

+ 7, 5


90

Sedangkan persamaan untuk metode Relaksasi dengan ω = 1, 25 adalah (k)

(k−1)

x1

= −0, 25x1

(k) x2 (k) x3

= =

(k−1)

− 0, 9375x2

+ 7, 5

(k) (k−1) (k−1) −0, 9375x1 − 0, 25x2 + 0, 3125x3 (k) (k−1) 0, 3125x2 − 0, 25x3 − 7, 5

+ 9, 375

Tabel berikut ini menampilkan perhitungan dari masing-masing metode hingga iterasi ke-7.

k

0

(k) x1 (k) x2 (k) x3

1 1 1

k (k)

x1 (k) x2 (k) x3

Tabel 5.4: Hasil perhitungan iterasi Gauss-Seidel 1 2 3 4 5 6 5,2500 3,8125 -5,0468

3,1406 3,8828 -5,0293

3,0879 3,9267 -5,0183

3,0549 3,9542 -5,0114

3,0343 3,9714 -5,0072

3,0215 3,9821 -5,0044

7 3,0134 3,9888 -5,0028

Tabel 5.5: Hasil perhitungan iterasi Relaksasi dengan ω = 1, 25 0 1 2 3 4 5 6 7 1 1 1

6,3125 3,5195 -6,6501

2,6223 3,9585 -4,6004

3,1333 4,0102 -5,0967

2,9570 4,0075 -4,9735

3,0037 4,0029 -5,0057

2,9963 4,0009 -4,9983

3,0000 4,0002 -5,0003

Dari kasus ini, bisa kita simpulkan bahwa iterasi Relaksasi memerlukan proses iterasi yang lebih singkat dibandingkan iterasi Gauss-Seidel. Jadi, pada kasus ini (dan juga secara umum), Relaksasi lebih efektif dibandingkan Gauss-Seidel. Pertanyaannya sekarang, bagaimana menentukan nilai ω optimal? Metode Relaksasi dengan pilihan nilai ω yang berkisar antara 0 dan 1 disebut metode underrelaxation, dimana metode ini berguna agar sistem persamaan linear bisa mencapai kondisi konvergen walaupun sistem tersebut sulit mencapai kondisi konvergen dengan metode GaussSeidel. Sementara bila ω nilainya lebih besar dari angka 1, maka disebut metode successive over-relaxation (SOR), yang mana metode ini berguna untuk mengakselerasi atau mempercepat kondisi konvergen dibandingkan dengan Gauss-Seidel. Metode SOR ini juga sangat berguna untuk menyelesaikan sistem persamaan linear yang muncul dari persamaan diferensial-parsial tertentu. 5.5.1 Algoritma Iterasi Relaksasi • Langkah 1: Tentukan k=1 • Langkah 2: Ketika (k ≤ N ) lakukan Langkah 3-6 – Langkah 3: Untuk i=1,...,n, hitunglah

xi = (1 − ω) XOi +

´ ³ P Pn a XO + b a x − ω − i−1 ij j i ij j j=i+1 j=1 aii

5.5. ITERASI DENGAN RELAKSASI

91

– Langkah 4: Jika kx − XOk < ǫ, maka keluarkan OUTPUT (x1 , ..., xn ) lalu STOP – Langkah 5: Tentukan k=k+1 – Langkah 6: Untuk i=1,...n, tentukan XOi = xi • Langkah 7: OUTPUT (’Iterasi maksimum telah terlampaui’) lalu STOP Demikianlah catatan singkat dari saya tentang metode iterasi untuk menyelesaikan problem sistem persamaan linear. Saya cukupkan sementara sampai disini. Insya Allah akan saya sambung lagi dilain waktu. Kalau ada yang mau didiskusikan, silakan hubungi saya melalui email: [email protected].

Bab 6

Interpolasi

✍ Objektif : ⊲ Mengenalkan Interpolasi Lagrange ⊲ Mengenalkan Interpolasi Spline-cubic

6.1

Interpolasi Lagrange

Interpolasi Lagrange diterapkan untuk mendapatkan fungsi polinomial P (x) berderajat tertentu yang melewati sejumlah titik data. Misalnya, kita ingin mendapatkan fungsi polinomial berderajat satu yang melewati dua buah titik yaitu (x0 , y0 ) dan (x1 , y1 ). Langkah pertama yang kita lakukan adalah mendefinisikan fungsi berikut L0 (x) =

x − x1 x0 − x1

L1 (x) =

x − x0 x1 − x0

dan

kemudian kita definisikan fungsi polinomial sebagai berikut P (x) = L0 (x)y0 + L1 (x)y1 Jika semua persamaan diatas kita gabungkan, maka akan didapat P (x) = L0 (x)y0 + L1 (x)y1 x − x1 x − x0 P (x) = y0 + y1 x0 − x1 x1 − x0 dan ketika x = x0 P (x0 ) =

x0 − x0 x0 − x1 y0 + y1 = y 0 x0 − x1 x1 − x0 93

BAB 6. INTERPOLASI

94 dan pada saat x = x1 P (x1 ) =

x1 − x1 x1 − x0 y0 + y1 = y 1 x0 − x1 x1 − x0

dari contoh ini, kira-kira apa kesimpulan sementara anda? Ya.. kita bisa sepakat bahwa fungsi polinomial P (x) =

x − x0 x − x1 y0 + y1 x0 − x1 x1 − x0

(6.1)

benar-benar melewati titik (x0 , y0 ) dan (x1 , y1 ).

Sekarang mari kita perhatikan lagi contoh lainnya. Misalnya ada tiga titik yaitu (x0 , y0 ), (x1 , y1 ) dan (x2 , y2 ). Tentukanlah fungsi polinomial yang melewati ketiganya! Dengan pola yang sama kita bisa awali langkah pertama yaitu mendefinisikan L0 (x) =

(x − x1 )(x − x2 ) (x0 − x1 )(x0 − x2 )

L1 (x) =

(x − x0 )(x − x2 ) (x1 − x0 )(x1 − x2 )

L2 (x) =

(x − x0 )(x − x1 ) (x2 − x0 )(x2 − x1 )

lalu

dan

kemudian kita definisikan fungsi polinomial sebagai berikut P (x) = L0 (x)y0 + L1 (x)y1 + L2 (x)y2 Jika semua persamaan diatas kita gabungkan, maka akan didapat fungsi polinomial P (x) =

(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 ) y0 + y1 + y2 (x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

Kita uji sebentar. Ketika x = x0 P (x0 ) =

(x0 − x0 )(x0 − x2 ) (x0 − x0 )(x0 − x1 ) (x0 − x1 )(x0 − x2 ) y0 + y1 + y2 = y 0 (x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

pada saat x = x1 P (x1 ) =


pada saat x = x2 P (x2 ) =


6.2. INTERPOLASI CUBIC SPLINE

95

Terbukti bahwa fungsi polonomial P (x) =

(x − x0 )(x − x2 ) (x − x0 )(x − x1 ) (x − x1 )(x − x2 ) y0 + y1 + y2 (x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

(6.2)

melewati ketiga titik tadi.

Kalau kita bandingkan antara persamaan (6.1) dan persamaan (6.2), terlihat bahwa derajat persamaan (6.2) lebih tinggi dibandingkan dengan derajat persamaan (6.1). Hal ini terlihat dari x2 pada persamaan (6.2) sementara pada persamaan (6.1) hanya ada x. persamaan (6.2) disebut funsi polinomial berderajat 2, sedangkan persamaan (6.1) disebut fungsi polinomial berderajat 1.

6.2

Interpolasi Cubic Spline

Gambar 6.1: Fungsi f (x) dengan sejumlah titik data

Gambar 6.2: Pendekatan dengan polinomial cubic spline

BAB 6. INTERPOLASI

96

Diketahui suatu fungsi f (x) (Figure 6.1) yang dibatasi oleh interval a dan b, dan memiliki sejumlah titik data a = x0 < x1 < ... < xn = b. Interpolasi cubic spline S(x) adalah sebuah potongan fungsi polinomial kecil-kecil (Figure 6.2) berderajat tiga (cubic ) yang menghubungkan dua titik data yang bersebelahan dengan ketentuan sebagai berikut: 1. Sj (x) adalah potongan fungsi yang berada pada sub-interval dari xj hingga xj+1 untuk nilai j = 0, 1, ..., n − 1; 2. S(xj ) = f (xj ), artinya pada setiap titik data (xj ), nilai f (xj ) bersesuaian dengan S(xj ) dimana j = 0, 1, ..., n; 3. Sj+1 (xj+1 ) = Sj (xj+1 ). Perhatikan titik xj+1 pada Figure 6.2. Ya.. tentu saja jika fungsi itu kontinyu, maka titik xj+1 menjadi titik sambungan antara Sj dan Sj+1 . ′ (xj+1 ) = Sj′ (xj+1 ), artinya kontinyuitas menuntut turunan pertama dari Sj dan Sj+1 4. Sj+1

pada titik xj+1 harus bersesuaian. ′′ (x ′′ 5. Sj+1 j+1 ) = Sj (xj+1 ), artinya kontinyuitas menuntut turunan kedua dari Sj dan Sj+1

pada titik xj+1 harus bersesuaian juga. 6. Salah satu syarat batas diantara 2 syarat batas x0 dan xn berikut ini mesti terpenuhi: • S ′′ (x0 ) = S ′′ (xn ) = 0 ini disebut natural boundary

• S ′ (x0 ) = f ′ (x0 ) dan S ′ (xn ) = f ′ (xn ) ini disebut clamped boundary

Polinomial cubic spline S (polinomial pangkat 3) untuk suatu fungsi f berdasarkan ketentuan di atas adalah Sj (x) = aj + bj (x − xj ) + cj (x − xj )2 + dj (x − xj )3

(6.3)

dimana j = 0, 1, ..., n − 1. Maka ketika x = xj Sj (xj ) = aj + bj (xj − xj ) + cj (xj − xj )2 + dj (xj − xj )3 Sj (xj ) = aj = f (xj ) Itu artinya, aj selalu jadi pasangan titik data dari xj . Dengan pola ini maka pasangan titik data xj+1 adalah aj+1 , konsekuensinya S(xj+1 ) = aj+1 . Berdasarkan ketentuan (3), yaitu ketika x = xj+1 dimasukan ke persamaan (12.7) aj+1 = Sj+1 (xj+1 ) = Sj (xj+1 ) = aj + bj (xj+1 − xj ) + cj (xj+1 − xj )2 + dj (xj+1 − xj )3 dimana j = 0, 1, ..., n − 2. Sekarang, kita nyatakan hj = xj+1 − xj , sehingga aj+1 = aj + bj hj + cj h2j + dj h3j Kemudian, turunan pertama dari persamaan (12.7) adalah Sj′ (x) = bj + 2cj (x − xj ) + 3dj (x − xj )2

(6.4)


97

ketika x = xj , Sj′ (xj ) = bj + 2cj (xj − xj ) + 3dj (xj − xj )2 = bj dan ketika x = xj+1 , bj+1 = Sj′ (xj+1 ) = bj + 2cj (xj+1 − xj ) + 3dj (xj+1 − xj )2 Ini dapat dinyatakan sebagai bj+1 = bj + 2cj (xj+1 − xj ) + 3dj (xj+1 − xj )2 dan dinyatakan dalam hj bj+1 = bj + 2cj hj + 3dj h2j

(6.5)

Berikutnya, kita hitung turunan kedua dari persamaan (12.7) Sj′′ (x) = 2cj + 6dj (x − xj )

(6.6)

tapi dengan ketentuan tambahan yaitu S ′′ (x)/2, sehingga persamaan ini dimodifikasi menjadi Sj′′ (x) = cj + 3dj (x − xj ) dengan cara yang sama, ketika x = xj Sj′′ (xj ) = cj + 3dj (xj − xj ) = cj dan ketika x = xj+1 cj+1 = Sj′′ (xj+1 ) = cj + 3dj (xj+1 − xj ) cj+1 = cj + 3dj hj dan dj bisa dinyatakan dj =

(6.7)

1 (cj+1 − cj ) 3hj

dari sini, persamaan (6.4) dapat ditulis kembali aj+1 = aj + bj hj + cj h2j + dj h3j = aj + bj hj + cj h2j + = aj + bj hj +

h2j (cj+1 − cj ) 3

h2j (2cj + cj+1 ) 3

(6.8)

BAB 6. INTERPOLASI

98 sementara persamaan (6.5) menjadi bj+1 = bj + 2cj hj + 3dj h2j = bj + 2cj hj + hj (cj+1 − cj ) = bj + hj (cj + cj+1 )

(6.9)

Sampai sini masih bisa diikuti, bukan? Selanjutnya, kita coba mendapatkan bj dari persamaan (6.8) bj = dan untuk bj−1 bj−1 =

hj 1 (aj+1 − aj ) − (2cj + cj+1 ) hj 3

(6.10)

hj−1 1 (aj − aj−1 ) − (2cj−1 + cj ) hj−1 3

(6.11)

Langkah berikutnya adalah mensubtitusikan persamaan (6.10) dan persamaan (6.11) kedalam persamaan (6.9), hj−1 cj−1 + 2(hj−1 + hj )cj + hj cj+1 =

3 3 (aj+1 − aj ) − (aj − aj−1 ) hj hj−1

(6.12)

n dimana j = 1, 2, ..., n − 1. Dalam sistem persamaan ini, nilai {hj }n−1 j=0 dan nilai {aj }j=0 su-

dah diketahui, sementara nilai {cj }nj=0 belum diketahui dan memang nilai inilah yang akan dihitung dari persamaan ini.

Sekarang coba perhatikan ketentuan nomor (6), ketika S ′′ (x0 ) = S ′′ (xn ) = 0, berapakah nilai c0 dan cn ? Nah, kita bisa evaluasi persamaan (6.6) S ′′ (x0 ) = 2c0 + 6d0 (x0 − x0 ) = 0 jelas sekali c0 harus berharga nol. Demikian halnya dengan cn harganya harus nol. Jadi untuk

natural boundary, nilai c0 = cn = 0. Persamaan (6.12) dapat dihitung dengan operasi matrik Ax = b dimana 

1

0

0

  h0 2(h0 + h1 ) h1  0 h1 2(h1 + h2 )  A= . . . ... ...  . . . ... ...  0 ... ...

...

...

0

...

h2

0

...

...

. . . hn−2 ...

  c0    c1   x=  ..  . cn

0

...

0



 0   ... 0    ... ...   2(hn−2 + hn−1 ) hn−1   0 1 ...


99





0

  3 3   h1 (a2 − a1 ) − h0 (a1 − a0 )     .. b=  .     3 (a − a 3 n−1 ) − hn−2 (an−1 − an−2 )  hn−1 n 0

Sekarang kita beralih ke clamped boundary dimana S ′ (a) = f ′ (a) dan S ′ (b) = f ′ (b). Nah, kita bisa evaluasi persamaan (6.10) dengan j = 0, dimana f ′ (a) = S ′ (a) = S ′ (x0 ) = b0 , sehingga f ′ (a) =

h0 1 (a1 − a0 ) − (2c0 + c1 ) h0 3

konsekuensinya, 2h0 c0 + h0 c1 =

3 (a1 − a0 ) − 3f ′ (a) h0

(6.13)

Sementara pada xn = bn dengan persamaan (6.9) f ′ (b) = bn = bn−1 + hn−1 (cn−1 + cn ) sedangkan bn−1 bisa didapat dari persamaan (6.11) dengan j = n − 1 bn−1 =

1 hn−1

(an − an−1 ) −

hn−1 (2cn−1 j + cn ) 3

Jadi f ′ (b) = =

hn−1 1 (an − an−1 ) − (2cn−1 j + cn ) + hn−1 (cn−1 + cn ) hn−1 3 1 hn−1 (an − an−1 + (cn−1 j + 2cn ) hn−1 3

dan akhirnya kita peroleh hn−1 cn−1 + 2hn−1 Cn = 3f ′ (b) −

3 hn−1

(an − an−1 )

(6.14)

Persamaan (6.13) dan persamaan (6.14) ditambah persamaan (6.12 membentuk operasi matrik Ax = b dimana  2h0 h0 0   h0 2(h0 + h1 ) h1   0 h1 2(h1 + h2 )  A= ... ... ...  ... ... ...  0

...

...

...

...

0

...

h2

0

...

...

. . . hn−2 ...

0

...

0



   ... 0    ... ...   2(hn−2 + hn−1 ) hn−1   hn−1 2hn−1 ...

0

BAB 6. INTERPOLASI

100

Gambar 6.3: Profil suatu object



  c0    c1   x=  ..  . cn 3 h0 (a1

− a0 ) − 3f ′ (a)



  3 3   h1 (a2 − a1 ) − h0 (a1 − a0 )   .   .. b=      3 (a − a 3 n−1 ) − hn−2 (an−1 − an−2 )  hn−1 n 3 3f ′ (b) − hn−1 (an − an−1 )


Gambar 6.4: Sampling titik data

Gambar 6.5: Hasil interpolasi cubic spline

101

j 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

xj 0,9 1,3 1,9 2,1 2,6 3,0 3,9 4,4 4,7 5,0 6,0 7,0 8,0 9,2 10,5 11,3 11,6 12,0 12,6 13,0 13,3

aj 1,3 1,5 1,85 2,1 2,6 2,7 2,4 2,15 2,05 2,1 2,25 2,3 2,25 1,95 1,4 0,9 0,7 0,6 0,5 0,4 0,25

bj 5,4 0,42 1,09 1,29 0,59 -0,02 -0,5 -0,48 -0,07 0,26 0,08 0,01 -0,14 -0,34 -0,53 -0,73 -0,49 -0,14 -0,18 -0,39

cj 0,00 -0,30 1,41 -0,37 -1,04 -0,50 -0,03 0,08 1,27 -0,16 -0,03 -0,04 -0,11 -0,05 -0,1 -0,15 0,94 -0,06 0 -0,54

dj -0,25 0,95 -2,96 -0,45 0,45 0,17 0,08 1,31 -1,58 0,04 0,00 -0,02 0,02 -0,01 -0,02 1,21 -0,84 0,04 -0,45 0,60

Gambar 6.6: Hasil interpolasi lagrange

Bab 7

Diferensial Numerik

✍ Objektif : ⊲ Mengenalkan metode Euler. ⊲ Mengenalkan metode Runge Kutta orde 4. ⊲ Mengenalkan metode Finite Difference.

7.1

Metode Euler

Suatu persamaan diferensial dapat dinyatakan sebagai berikut: dy = f (t, y), dt

a ≤ t ≤ b,

y(a) = α

(7.1)

Pada kenyataannya, melalui pendekatan numerik, kita tidak akan memperoleh solusi fungsi yang kontinyu; yang mungkin kita dapat adalah solusi diskrit dalam bentuk mesh points di dalam interval [a,b]. Setelah diperoleh solusi numerik pada suatu point, maka point-point yang lainpun bisa dicari dengan cara interpolasi. Tahap awal solusi pendekatan numerik adalah dengan menentukan point-point dalam jarak yang sama di dalam interval [a,b], yaitu dengan menerapkan ti = a + ih,

i = 0, 1, 2, ..., N

(7.2)

Jarak antar point dirumuskan sebagai h=

b−a N

(7.3)

ini disebut step size. Metode Euler diturunkan dari deret Taylor. Misalnya, fungsi y(t) adalah fungsi yang kontinyu 103

BAB 7. DIFERENSIAL NUMERIK

104

Gambar 7.1: Metode Euler

dan memiliki turunan dalam interval [a,b]. Maka dalam deret Taylor y(ti+1 ) = y(ti ) + (ti+1 − ti )y ′ (ti ) +

(ti+1 − ti )2 ′′ y (ξi ) 2

(7.4)

h2 ′′ y (ξi ) 2

(7.5)

Karena h = (ti+1 − ti ), maka y(ti+1 ) = y(ti ) + hy ′ (ti ) + dan, karena y(t) memenuhi persamaan diferensial (7.1), y(ti+1 ) = y(ti ) + hf (ti , y(ti )) +

h2 ′′ y (ξi ) 2

(7.6)

Metode Euler dibangun dengan pendekatan wi ≈ y(ti ) untuk i = 1, 2, 3, ..., N , dengan mengabaikan suku terakhir yang terdapat pada persamaan (7.6). Jadi metode Euler dinyatakan sebagai w0 = α

(7.7)

wi+1 = wi + hf (ti , wi )

(7.8)

dimana i = 0, 1, 2, .., N − 1

Contoh Diketahui persamaan diferensial y ′ = y − t2 + 1,

0 ≤ t ≤ 2,

y(0) = 0, 5

dimana N = 10. Sehingga h=

2−0 b−a = = 0, 2 N 10

7.1. METODE EULER

105

dan ti = a + ih = 0 + i(0, 2)

→

ti = 0, 2i

serta w0 = 0, 5 Dengan demikian persamaan euler dapat dinyatakan sebagai wi+1 = wi + h(wi − t2i + 1)

= wi + 0, 2(wi − 0, 04i2 + 1) = 1, 2wi − 0, 008i2 + 0, 2

dimana i = 0, 1, ..., 9. Pada saat i = 0 dan dari syarat awal diketahui w0 = 0, 5 w1 = 1, 2w0 − 0, 008(0)2 + 0, 2 = 0, 8000000 Pada saat i = 1 w2 = 1, 2w1 − 0, 008(1)2 + 0, 2 = 1, 1520000 Pada saat i = 2 w3 = 1, 2w2 − 0, 008(2)2 + 0, 2 = 1, 5504000 Demikian seterusnya, hingga pada i = 9 w10 = 1, 2w9 − 0, 008(9)2 + 0, 2 = 4, 8657845 Disisi lain, solusi exact persamaan diferensial adalah y(t) = (t + 1)2 − 0, 5et Tabel dibawah ini memperlihatkan solusi metode euler dan solusi exact serta error atau selisih antara keduanya. Trend error menunjukan bahwa ketika i semakin besar, error juga semakin meningkat. Figure (9.3) memperlihatkan kurva peningkatan error ketika i semakin besar. Untuk mengatasi hal ini, salah satu pemecahannya adalah dengan menerapkan deret Taylor berorde lebih tinggi, atau cara lainnya adalah dengan menggunakan metode Runge-Kutta yang akan dijelaskan pada pertemuan berikutnya.


106 ti 0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 1,8 2,0

wi 0,5000000 0,8000000 1,1520000 1,5504000 1,9884800 2,4581760 2,9498112 3,4517734 3,9501281 4,4281538 4,8657845

yi = y(ti ) 0,5000000 0,8292986 1,2140877 1,6489406 2,1272295 2,6408591 3,1799415 3,7324000 4,2834838 4,8151763 5,3054720

|wi − yi | 0,0000000 0,0292986 0,0620877 0,0985406 0,1387495 0,1826831 0,2301303 0,2806266 0,3333557 0,3870225 0,4396874

Gambar 7.2: Trend error metode euler

7.2

Metode Runge Kutta

Pada saat membahas metode Euler untuk penyelesaian persamaan diferensial, kita telah sampai pada kesimpulan bahwa truncation error metode Euler terus membesar seiring dengan bertambahnya iterasi. Dikaitkan dengan hal tersebut, metode Runge-Kutta Orde Empat menawarkan penyelesaian persamaan diferensial dengan pertumbuhan truncation error yang jauh lebih kecil.

7.2. METODE RUNGE KUTTA

107

Persamaan-persamaan yang menyusun metode Runge-Kutta Orde Empat adalah w0 = α k1 = hf (ti , wi ) h 1 k2 = hf (ti + , wi + k1 ) 2 2 h 1 k3 = hf (ti + , wi + k2 ) 2 2 k4 = hf (ti+1 , wi + k3 ) 1 wi+1 = wi + (k1 + 2k2 + 2k3 + k4 ) 6

(7.9) (7.10) (7.11) (7.12) (7.13)

Contoh Diketahui persamaan diferensial y ′ = y − t2 + 1,

0 ≤ t ≤ 2,

y(0) = 0, 5

dengan mengganti y menjadi w, kita bisa nyatakan f (ti , wi ) sebagai f (ti , wi ) = wi − t2i + 1 Jika N = 10, maka h=

b−a 2−0 = = 0, 2 N 10

dan ti = a + ih = 0 + i(0, 2)

→

ti = 0, 2i

serta w0 = 0, 5 Sekarang mari kita demonstrasikan metode Runge-Kutta Orde Empat ini. Untuk menghitung w1 , tahap-tahap perhitungannya dimulai dari menghitung k1 k1 = hf (t0 , w0 ) = h(w0 − t20 + 1)

= 0, 2((0, 5) − (0, 0)2 + 1) = 0, 3 lalu menghitung k2 k1 h , w0 + ) 2 2 h k1 = h[(w0 + ) − (t0 + )2 + 1)] 2 2 0, 2 2 0, 3 ) − (0, 0 + ) + 1)] = 0, 2[(0, 5 + 2 2 = 0, 328

k2 = hf (t0 +


108 dilanjutkan dengan k3

h k2 , w0 + ) 2 2 k2 h = h[(w0 + ) − (t0 + )2 + 1)] 2 2 0, 328 0, 2 2 = 0, 2[(0, 5 + ) − (0, 0 + ) + 1)] 2 2 = 0, 3308

k3 = hf (t0 +

kemudian k4 k4 = hf (t1 , w0 + k3 ) = h[(w0 + k3 ) − t21 + 1]

= 0, 2[(0, 5 + 0, 3308) − (0, 2)2 + 1] = 0, 35816 akhirnya diperoleh w1 1 (k1 + 2k2 + 2k3 + k4 ) 6 1 = 0, 5 + (0, 3 + 2(0, 328) + 2(0, 3308) + 0, 35816) 6 1 = 0, 5 + (0, 3 + 0, 656 + 0, 6616 + 0, 35816) 6 = 0, 8292933

w1 = w 0 +

Dengan cara yang sama, w2 , w3 , w4 dan seterusnya dapat dihitung. Tabel berikut menunjukkan hasil perhitungannya.

i 0 1 2 3 4 5 6 7 8 9 10

ti 0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 1,8 2,0

wi 0,5000000 0,8292933 1,2140762 1,6489220 2,1272027 2,6408227 3,1799942 3,7323401 4,2834095 4,8150857 5,3054720

yi = y(ti ) 0,5000000 0,8292986 1,2140877 1,6489406 2,1272295 2,6408591 3,1799415 3,7324000 4,2834838 4,8151763 5,3053630

|wi − yi | 0,0000000 0,0000053 0,0000114 0,0000186 0,0000269 0,0000364 0,0000474 0,0000599 0,0000743 0,0000906 0,0001089

Dibandingkan dengan metode Euler, tingkat pertumbuhan truncation error, pada kolom |wi − yi |, jauh lebih rendah sehingga metode Runge-Kutta Orde Empat lebih disukai untuk mem-

bantu menyelesaikan persamaan-diferensial-biasa.


109

Gambar 7.3: Rangkaian RC Contoh tadi tampaknya dapat memberikan gambaran yang jelas bahwa metode Runge-Kutta Orde Empat dapat menyelesaikan persamaan diferensial biasa dengan tingkat akurasi yang lebih tinggi. Namun, kalau anda jeli, ada suatu pertanyaan cukup serius yaitu apakah metode ini dapat digunakan bila pada persamaan diferensialnya tidak ada variabel t ? Misalnya pada kasus berikut ini Sebuah kapasitor yang tidak bermuatan dihubungkan secara seri dengan sebuah resistor dan baterry (Figure 7.3). Diketahui ǫ = 12 volt, C = 5,00 µF dan R = 8,00 ×105 Ω. Saat saklar dihubungkan (t=0), muatan belum ada (q=0).

ǫ q dq = − dt R RC

(7.14)

Solusi exact persamaan (7.14) adalah ³ ´ qexact = q(t) = Cǫ 1 − e−t/RC

(7.15)

Anda bisa lihat semua suku di ruas kanan persamaan (7.14) tidak mengandung variabel t. Padahal persamaan-persamaan penyusun metode Runge-Kutta selalu mencantumkan variabel t. Apakah persamaan (7.14) tidak bisa diselesaikan dengan metode Runge-Kutta? Belum tentu. Sekarang, kita coba selesaikan, pertama kita nyatakan m1 = m2 =

ǫ = 1, 5 × 10−5 R 1 = 0, 25 RC

sehingga persamaan (7.14) dimodifikasi menjadi dq = f (qi ) = m1 − qi m2 dt ti = a + ih Jika t0 = 0, maka a = 0, dan pada saat itu (secara fisis) diketahui q0 = 0, 0. Lalu jika ditetapkan


110

h = 0, 1 maka t1 = 0, 1 dan kita bisa mulai menghitung k1 dengan menggunakan q0 = 0, 0, walaupun t1 tidak dilibatkan dalam perhitungan ini k1 = hf (q0 ) = h(m1 − q0 m2 )

= 0, 1((1, 5 × 10−5 ) − (0, 0)(0, 25))

= 0, 150 × 10−5 lalu menghitung k2 k2 = hf (q0 +

k1 ) 2

= h[(m1 − (q0 +

k1 )m2 )] 2

= 0, 1[(1, 5 × 10−5 − ((0, 0) +

0, 15 × 10−5 )(0, 25)] 2

= 0, 14813 × 10−5 dilanjutkan dengan k3 k3 = hf (q0 +

k2 ) 2

= h[(m1 − (q0 +

k2 )m2 )] 2

= 0, 1[(1, 5 × 10−5 − ((0, 0) +

0, 14813 × 10−5 )(0, 25)] 2

= 0, 14815 × 10−5 kemudian k4 k4 = hf (q0 + k3 ) = h[(m1 − (q0 + k3 )m2 )]

= 0, 1[(1, 5 × 10−5 − ((0, 0) + 0, 14815 × 10−5 )(0, 25)]

= 0, 14630 × 10−5 akhirnya diperoleh q1

1 (k1 + 2k2 + 2k3 + k4 ) 6 1 = 0, 0 + (0, 150 + 2(0, 14813) + 2(0, 14815) + 0, 14630) × 10−5 6 = 0, 14814 × 10−5

q1 = q0 +

Selanjutnya q2 dihitung. Tentu saja pada saat t2 , dimana t2 = 0, 2, namun sekali lagi, t2 tidak


111

terlibat dalam perhitungan ini. Dimulai menghitung k1 kembali k1 = hf (q1 ) = h(m1 − q1 m2 )

= 0, 1((1, 5 × 10−5 ) − (0, 14814 × 10−5 )(0, 25))

= 0, 14630 × 10−5 lalu menghitung k2 k2 = hf (q1 +

k1 ) 2

= h[(m1 − (q1 +

k1 )m2 )] 2

= 0, 1[(1, 5 × 10−5 − ((0, 14814 × 10−5 ) +

0, 14630 × 10−5 )(0, 25)] 2

= 0, 14447 × 10−5 dilanjutkan dengan k3 k3 = hf (q1 +

k2 ) 2

= h[(m1 − (q1 +

k2 )m2 )] 2

= 0, 1[(1, 5 × 10−5 − ((0, 14814 × 10−5 ) +

0, 14447 × 10−5 )(0, 25)] 2

= 0, 14449 × 10−5 kemudian k4 k4 = hf (q1 + k3 ) = h[(m1 − (q1 + k3 )m2 )]

= 0, 1[(1, 5 × 10−5 − ((0, 14814 × 10−5 ) + 0, 14449 × 10−5 )(0, 25)]

= 0, 14268 × 10−5 akhirnya diperoleh q2

1 (k1 + 2k2 + 2k3 + k4 ) 6 1 = 0, 14814 × 10−5 + (0, 14630 + 2(0, 14447) + 2(0, 14449) + 0, 14268) × 10−5 6 −5 = 0, 29262 × 10

q2 = q1 +

Dengan cara yang sama, q3 , q4 , q5 dan seterusnya dapat dihitung. Tabel di atas menunjukkan hasil perhitungannya. Kolom qexact diperoleh dari persamaan (7.15). Luar biasa!! Tak ada error sama sekali. Mungkin, kalau kita buat 7 angka dibelakang koma,

error nya akan terlihat. Tapi kalau anda cukup puas dengan 5 angka dibelakang koma, hasil ini


112 i 0 1 2 3 4 5 6 7 8 9 10

ti 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0

qi 0,00000×10−5 0,14814×10−5 0,29262×10−5 0,43354×10−5 0,57098×10−5 0,70502×10−5 0,83575×10−5 0,96326×10−5 1,0876×10−5 1,2089×10−5 1,3272×10−5

qexact = q(ti ) 0,00000×10−5 0,14814×10−5 0,29262×10−5 0,43354×10−5 0,57098×10−5 0,70502×10−5 0,83575×10−5 0,96326×10−5 1,0876×10−5 1,2089×10−5 1,3272×10−5

|qi − qexact | 0,00000 0,00000 0,00000 0,00000 0,00000 0,00000 0,00000 0,00000 0,00000 0,00000 0,00000

sangat memuaskan. Figure 7.4 memperlihatkan kurva penumpukan muatan q terhadap waktu t. −5

6

x 10

4

2

0

0

2

4

6

8

10

12

14

16

18

Gambar 7.4: Kurva muatan q terhadap waktu t Berikut ini adalah script dalam matlab yang dipakai untuk menghitung q 1 2 3 4 5 6 7 8 9

clear all clc E=12; R=800000; C=5e-6; m1=E/R; m2=1/(R*C); b=20.0; a=0.0;

20

7.3. METODE FINITE DIFFERENCE 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

113

h=0.1; n=(b-a)/h; q0=0.0; t0=0.0; for i=1:n t(i)=a+i*h; end for i=1:n if i==1 k1=h*(m1-(m2*q0)); k2=h*(m1-(m2*(q0+(k1/2)))); k3=h*(m1-(m2*(q0+(k2/2)))); k4=h*(m1-(m2*(q0+k3))); q(i)=q0+(k1+(2*k2)+(2*k3)+k4)/6; else k1=h*(m1-(m2*q(i-1))); k2=h*(m1-(m2*(q(i-1)+(k1/2)))); k3=h*(m1-(m2*(q(i-1)+(k2/2)))); k4=h*(m1-(m2*(q(i-1)+k3))); q(i)=q(i-1)+(k1+(2*k2)+(2*k3)+k4)/6; end end q

Sampai disini mudah-mudahan jelas dan bisa dimengerti. Silakan anda coba untuk kasus yang lain, misalnya proses pembuangan (discharging ) q pada rangkaian yang sama, atau bisa juga anda berlatih dengan rangkaian RL dan RLC. Saya akhiri dulu uraian saya sampai disini.

7.3

Metode Finite Difference

Suatu persamaan diferensial dapat dinyatakan sebagai berikut: dy d2 y (x) = p(x) (x) + q(x)y(x) + r(x), 2 dx dx

a ≤ x ≤ b,

y(a) = α,

y(b) = β

(7.16)

atau juga dapat dituliskan dalam bentuk lain y ′′ = p(x)y ′ + q(x)y + r(x)

(7.17)

Persamaan tersebut dapat diselesaikan dengan melakukan pendekatan numerik terhadap y ′′ dan y ′ . Caranya adalah pertama, kita memilih angka integer sembarang yaitu N dimana N > 0 dan membagi interval [a, b] dengan (N + 1), hasilnya dinamakan h (lihat Gambar 7.5) h=

b−a N +1

(7.18)

Dengan demikian maka titik-titik x yang merupakan sub-interval antara a dan b dapat dinyatakan sebagai xi = a + ih,

i = 0, 1, ..., N + 1

(7.19)

Pencarian solusi persamaan diferensial melalui pendekatan numerik dilakukan dengan memanfaatkan polinomial Taylor untuk mengevaluasi y ′′ dan y ′ pada xi+1 dan xi−1 seperti berikut


114

Gambar 7.5: Kurva suatu fungsi f (x) yang dibagi sama besar berjarak h. Evaluasi kurva yang dilakukan Finite-Difference dimulai dari batas bawah X0 = a hingga batas atas x6 = b

ini y(xi+1 ) = y(xi + h) = y(xi ) + hy ′ (xi ) +

h2 ′′ y (xi ) 2

(7.20)

y(xi−1 ) = y(xi − h) = y(xi ) − hy ′ (xi ) +

h2 ′′ y (xi ) 2

(7.21)

dan

Jika kedua persamaan ini dijumlahkan y(xi+1 ) + y(xi−1 ) = 2y(xi ) + h2 y ′′ (xi ) Dari sini y ′′ dapat ditentukan h2 y ′′ (xi ) = y(xi+1 ) − 2y(xi ) + y(xi−1 ) y ′′ (xi ) =

y(xi+1 ) − 2y(xi ) + y(xi−1 ) h2

(7.22)

Dengan cara yang sama, y ′ (xi ) dapat dicari sebagai berikut y ′ (xi ) =

y(xi+1 ) − y(xi−1 ) 2h

Selanjutnya persamaan (7.22) dan (7.23) disubstitusikan ke persamaan (7.17) maka y(xi+1 ) − 2y(xi ) + y(xi−1 ) y(xi+1 ) − y(xi−1 ) = p(xi ) + q(xi )y(xi ) + r(xi ) 2 h 2h y(xi+1 ) − y(xi−1 ) −y(xi+1 ) + 2y(xi ) − y(xi−1 ) = −p(xi ) − q(xi )y(xi ) − r(xi ) h2 2h −y(xi+1 ) + 2y(xi ) − y(xi−1 ) y(xi+1 ) − y(xi−1 ) + p(xi ) + q(xi )y(xi ) = −r(xi ) 2 h 2h

(7.23)

7.3. METODE FINITE DIFFERENCE

115

Sebelum dilanjut, saya nyatakan bahwa y(xi+1 )=wi+1 dan y(xi )=wi serta y(xi−1 )=wi−1 . Maka persamaan di atas dapat ditulis sebagai berikut µ

−wi+1 + 2wi − wi−1 h2

¶

+ p(xi )

µ

wi+1 − wi−1 2h

¶

+ q(xi )wi = −r(xi )

h p(xi ) (wi+1 − wi−1 ) + h2 q(xi )wi 2 h h −wi+1 + 2wi − wi−1 + p(xi )wi+1 − p(xi )wi−1 + h2 q(xi )wi 2 2 h h 2 −wi−1 − p(xi )wi−1 + 2wi + h q(xi )wi − wi+1 + p(xi )wi+1 2 2 ¶ µ ¶ µ ¡ ¢ h h − 1 + p(xi ) wi−1 + 2 + h2 q(xi ) wi − (1 − p(xi ) wi+1 2 2 (−wi+1 + 2wi − wi−1 ) +

= −h2 r(xi ) = −h2 r(xi ) = −h2 r(xi ) = −h2 r(xi )

(7.24)

dimana i=1,2,3...sampai N, karena yang ingin kita cari adalah w1 , w2 , w3 ,..., wN . Sementara, satu hal yang tak boleh dilupakan yaitu w0 dan wN +1 biasanya selalu sudah diketahui. Pada persamaan (7.16), jelas-jelas sudah diketahui bahwa w0 =α dan wN +1 =β; keduanya dikenal sebagai syarat batas atau istilah asingnya adalah boundary value. Topik yang sedang bahas ini juga sering disebut sebagai Masalah Syarat Batas atau Boundary Value Problem. Sampai disini kita mendapatkan sistem persamaan linear yang selanjutnya dapat dinyatakan sebagai bentuk operasi matrik Aw = b

(7.25)

dimana A adalah matrik tridiagonal dengan orde N × N  A

=

2 + h2 q(x1 ) −1 − h2 p(x2 )  0   0   ...  ... 0

p(x1 ) −1 + h 2 2 + h2 q(x2 ) −1 − h p(x3 ) 2 0 ... ... ...



w1

0 −1 + h p(x2 ) 2 2 + h2 q(x3 ) p(x4 ) −1 − h 2 ... ... ...



   w2     w   3     w w= 4    ..   .    w   N −1  wN 7.3.1 Script Finite-Difference 1 2

... 0 −1 + h p(x3 ) 2 2 + h2 q(x4 ) ... ... ...



... ... 0 −1 + h p(x4 ) 2 ... −1 − h p(xN −1 ) 2 ...

¢ ¡ −h2 r(x1 ) + 1 + h2 p(x1 ) w0



    −h2 r(x2 )     2 −h r(x ) 3     2 r(x )   −h b= 4    ..   .     2 −h r(x ) N −1   ¡ ¢ h 2 −h r(xN ) + 1 − 2 p(xN ) wN +1

clear all clc

3 4

... ... ... 0 ... 2 + h2 q(xN −1 ) −1 − h p(xN ) 2

a=1.0; %ganti angkanya sesuai data yang anda miliki



0 0   0   0   ...  h −1 + 2 p(xN −1 ) 2 2 + h q(xN )


116 5 6 7 8 9

b=2.0; %ganti angkanya sesuai data yang anda miliki n=9; %ganti angkanya sesuai data yang anda miliki h=(b-a)/(n+1); alpha=1; %ganti angkanya sesuai data yang anda miliki beta=2; %ganti angkanya sesuai data yang anda miliki

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

%====== Mencari Elemen Matrik A ======== for i=1:n x=a+i*h; A(i,i)=2+h^2*fungsiQ(x); end for i=1:n-1 x=a+i*h; A(i,i+1)=-1+((h/2)*fungsiP(x)); end for i=2:n x=a+i*h; A(i,i-1)=-1-((h/2)*fungsiP(x)); end A %====== Mencari Elemen Vektor b ======== x=a+h; b(1,1)=-h^2*fungsiR(x)+(1+((h/2)*fungsiP(x)))*alpha; for i=2:8 x=a+i*h; b(i,1)=-h^2*fungsiR(x); end xn=a+n*h b(n,1)=-h^2*fungsiR(xn)+(1-((h/2)*fungsiP(xn)))*beta; b

Pada akhirnya, elemen-elemen matrik A dan vektor b sudah diketahui. Sehingga vektor w dapat dihitung dengan berbagai metode pemecahan sistem persamaan linear, seperti Eliminasi Gauss, Gauss-Jourdan, Iterasi Jacobi dan Iterasi Gauss-Seidel.

Contoh Diketahui persamaan diferensial seperti berikut ini 2 sin(ln x) 2 , y ′′ = − y ′ + 2 y + x x x2 memiliki solusi exact y = c1 x + dimana c2 =

1 ≤ x ≤ 2,

y(1) = 1,

y(2) = 2

3 1 c2 − sin(ln x) − cos(ln x), 2 x 10 10

1 [8 − 12 sin(ln 2) − 4 cos(ln 2)] ≈ −0, 03920701320 70

7.3. METODE FINITE DIFFERENCE dan c1 =

117

11 − c2 ≈ 1, 1392070132. 10

Dengan metode Finite-Difference, solusi pendekatan dapat diperoleh dengan membagi interval 1 ≤ x ≤ 2 menjadi sub-interval, misalnya kita gunakan N = 9, sehingga spasi h diperoleh h=

2−1 b−a = = 0, 1 N +1 9+1

Dari persamaan diferensial tersebut juga didapat p(xi ) = − q(xi ) = r(xi ) =

2 xi

2 x2i sin(ln xi ) x2i

Script matlab telah dibuat untuk menyelesaikan contoh soal ini. Untuk memecahkan persoalan ini, saya membuat 4 buah script, terdiri dari script utama, script fungsiP, script fungsiQ dan script fungsiR. Berikut ini adalah script fungsiP yang disimpan dengan nama file fungsiP.m: 1 2

function y = fungsiP(x) y = -2/x;

lalu inilah script fungsiQ yang disimpan dengan nama file fungsiQ.m: 1 2

function y = fungsiQ(x) y = 2/x^2;

kemudian ini script fungsiR yang disimpan dengan nama file fungsiR.m:: 1 2

function y = fungsiR(x) y = sin(log(x))/x^2;

dan terakhir, inilah script utamanya: 1 2

clear all clc

3 4 5

a=1.0; b=2.0;

6 7 8

alpha=1; beta=2;

9 10

%=======jika diketahui n, maka h dihitung ====

118 11 12


n=9; h=(b-a)/(n+1);

13 14 15 16

%=======jika diketahui h, maka n dihitung ==== %h=0.1; %n=((b-a)/h)-1;

17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

%====== Mencari Elemen Matrik A ======== for i=1:n x=a+i*h; A(i,i)=2+h^2*fungsiQ(x); end for i=1:n-1 x=a+i*h; A(i,i+1)=-1+((h/2)*fungsiP(x)); end for i=2:n x=a+i*h; A(i,i-1)=-1-((h/2)*fungsiP(x)); end A %====== Mencari Elemen Vektor b ======== x=a+h; b(1,1)=-h^2*fungsiR(x)+(1+((h/2)*fungsiP(x)))*alpha; for i=2:8 x=a+i*h; b(i,1)=-h^2*fungsiR(x); end xn=a+n*h b(n,1)=-h^2*fungsiR(xn)+(1-((h/2)*fungsiP(xn)))*beta; b %====== Menggabungkan Vektor b kedalam matrik A ======== for i=1:n A(i,n+1)=b(i,1); end A

47 48 49 50

%&&&&&& Proses Eliminasi Gauss &&&&&&&&&&&&&&&&&&&&&&&&& %---------Proses Triangularisasi----------for j=1:(n-1)

51 52 53 54 55 56 57

%----mulai proses pivot--if (A(j,j)==0) for p=1:n+1 u=A(j,p); v=A(j+1,p); A(j+1,p)=u;

7.3. METODE FINITE DIFFERENCE 58 59 60 61 62 63 64 65 66 67 68 69 70

119

A(j,p)=v; end end %----akhir proses pivot--jj=j+1; for i=jj:n m=A(i,j)/A(j,j); for k=1:(n+1) A(i,k)=A(i,k)-(m*A(j,k)); end end end %-------------------------------------------

71 72 73


74 75 76 77 78 79 80 81 82

for i=n-1:-1:1 S=0; for j=n:-1:i+1 S=S+A(i,j)*x(j,1); end x(i,1)=(A(i,n+1)-S)/A(i,i); end %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

83 84 85

%===== Menampilkan Vektor w ================= w=x

Tabel berikut ini memperlihatkan hasil perhitungan dengan pendekatan metode Finite-Difference wi dan hasil perhitungan dari solusi exact y(xi ), dilengkapi dengan selisih antara keduanya |wi − y(xi )|. Tabel ini memperlihatkan tingkat kesalahan (error) berada pada orde 10−5 . Unxi 1,0 1,1 1,2 1,3 1,4 1,5 1,6 1,7 1,8 1,9 2,0

wi 1,00000000 1,09260052 1,18704313 1,28333687 1,38140205 1,48112026 1,58235990 1,68498902 1,78888175 1,89392110 2,00000000

y(xi ) 1,00000000 1,09262930 1,18708484 1,28338236 1,38144595 1,48115942 1,58239246 1,68501396 1,78889853 1,89392951 2,00000000

|wi − y(xi )| 2,88 × 10−5 4,17 × 10−5 4,55 × 10−5 4,39 × 10−5 3,92 × 10−5 3,26 × 10−5 2,49 × 10−5 1,68 × 10−5 8,41 × 10−6

tuk memperkecil orde kesalahan, kita bisa menggunakan polinomial Taylor berorde tinggi. Akan tetapi proses kalkulasi menjadi semakin banyak dan disisi lain penentuan syarat batas


120

lebih kompleks dibandingkan dengan pemanfaatan polinomial Taylor yang sekarang. Untuk menghindari hal-hal yang rumit itu, salah satu jalan pintas yang cukup efektif adalah dengan menerapkan ekstrapolasi Richardson. Contoh Pemanfaatan ekstrapolasi Richardson pada metode Finite Difference untuk persamaan diferensial seperti berikut ini 2 sin(ln x) 2 , y ′′ = − y ′ + 2 y + x x x2

1 ≤ x ≤ 2,

y(1) = 1,

y(2) = 2,

dengan h = 0, 1, h = 0, 05, h = 0, 025. Ekstrapolasi Richardson terdiri atas 3 tahapan, yaitu ekstrapolasi yang pertama Ext1i =

4wi (h = 0, 05) − wi (h = 0, 1) 3

kemudian ekstrapolasi yang kedua Ext2i =

4wi (h = 0, 025) − wi (h = 0, 05) 3

dan terakhir ekstrapolasi yang ketiga Ext3i =

16Ext2i − Ext1i 15

Tabel berikut ini memperlihatkan hasil perhitungan tahapan-tahapan ekstrapolasi tersebut. Jika seluruh angka di belakang koma diikut-sertakan, maka akan terlihat selisih antara solusi exact dengan solusi pendekatan sebesar 6, 3 × 10−11 . Ini benar-benar improvisasi yang luar

biasa. xi 1,0 1,1 1,2 1,3 1,4 1,5 1,6 1,7 1,8 1,9 2,0

wi (h = 0, 1) 1,00000000 1,09260052 1,18704313 1,28333687 1,38140205 1,48112026 1,58235990 1,68498902 1,78888175 1,89392110 2,00000000

wi (h = 0, 05) 1,00000000 1,09262207 1,18707436 1,28337094 1,38143493 1,48114959 1,58238429 1,68500770 1,78889432 1,89392740 2,00000000

wi (h = 0, 025) 1,00000000 1,09262749 1,18708222 1,28337950 1,38144319 1,48115696 1,58239042 1,68501240 1,78889748 1,89392898 2,00000000

Ext1i 1,00000000 1,09262925 1,18708477 1,28338230 1,38144598 1,48115937 1,58239242 1,68501393 1,78889852 1,89392950 2,00000000

Ext2i 1,00000000 1,09262930 1,18708484 1,28338236 1,38144595 1,48115941 1,58239246 1,68501396 1,78889853 1,89392951 2,00000000

Ext3i 1,00000000 1,09262930 1,18708484 1,28338236 1,38144595 1,48115942 1,58239246 1,68501396 1,78889853 1,89392951 2,00000000

7.3.2 Aplikasi Besar simpangan terhadap waktu (y(t)) suatu sistem osilator mekanik yang padanya diberikan gaya secara periodik (forced-oscilations) memenuhi persamaan diferensial seperti dibawah ini

7.3. METODE FINITE DIFFERENCE

121

berikut syarat-syarat batasnya dy d2 y = + 2y + cos(t), 2 dt dt

0≤t≤

π , 2

y(0) = −0, 3,

π y( ) = −0, 1 2

Dengan metode Finite-Difference, tentukanlah besar masing-masing simpangan di setiap interval h = π/8. Buatlah table untuk membandingkan hasil finite-difference dengan solusi analitik 1 yang memenuhi y(t) = − 10 [sin(t) + 3cos(t)].

jawab: Secara umum, persamaan diferensial dapat dinyatakan sbb: dy d2 y (x) = p(x) (x) + q(x)y(x) + r(x), dx2 dx

a ≤ x ≤ b,

y(a) = α,

y(b) = β

Dengan membandingkan kedua persamaan di atas, kita bisa definisikan p(t) = 1

q(t) = 2

r(t) = cos(t)

a=0

b=

π 2

α = −0, 3

β = −0, 1

Adapun persamaan finite-difference adalah µ ¶ µ ¶ ¡ ¢ h h 2 − 1 + p(xi ) wi−1 + 2 + h q(xi ) wi − (1 − p(xi ) wi+1 = −h2 r(xi ) 2 2 Persamaan diatas dikonversi kedalam operasi matriks Aw = b

(7.26)

dimana A adalah matrik tridiagonal dengan orde N × N  A

=

2 + h2 q(x1 ) −1 − h2 p(x2 )  0   0   ...  ... 0

p(x1 ) −1 + h 2 2 + h2 q(x2 ) −1 − h p(x3 ) 2 0 ... ... ...



w1

0 −1 + h p(x2 ) 2 2 + h2 q(x3 ) −1 − h p(x4 ) 2 ... ... ...



   w2     w  3      w w=  4   ..   .     w  N −1  wN

... 0 −1 + h p(x3 ) 2 2 + h2 q(x4 ) ... ... ...



... ... 0 −1 + h p(x4 ) 2 ... −1 − h p(xN −1 ) 2 ...

... ... ... 0 ... 2 + h2 q(xN −1 ) p(xN ) −1 − h 2

¡ ¢ −h2 r(x1 ) + 1 + h2 p(x1 ) w0



0 0   0   0   ...  −1 + h p(xN −1 ) 2 2 + h2 q(xN )



    −h2 r(x2 )     2 −h r(x ) 3     2 r(x )  −h b= 4     ..   .     2 −h r(x )   N −1 ¢ ¡ −h2 r(xN ) + 1 − h2 p(xN ) wN +1

Jumlah baris matrik ditentukan oleh bilangan n. Namun disoal hanya tersedia informasi nilai


122 h = π/8, sehingga n harus dihitung terlebih dahulu: h=

b−a n+1

n=

π −0 b−a −1= 2 −1=3 h π/8

perhitungan ini dilakukan didalam script matlab. Selanjutnya seluruh elemen matrik A dan vektor b dihitung dengan matlab 

2, 3084

−0, 8037

0



w1





−0, 5014



      −1, 1963 2, 3084 −0, 8037   w2  =  −0, 1090  0 −1, 1963 2, 3084 w3 −0, 1394 Proses diteruskan dengan metode Eliminasi Gauss dan didapat hasil akhir berikut ini w1 = −0.3157

w2 = −0.2829

w3 = −0.2070

Bab 8

Pers. Diferensial Parsial Numerik

✍ Objektif : ⊲ Mengenalkan Persamaan Diferensial Parsial Eliptik ⊲ Mengenalkan Persamaan Diferensial Parsial Hiperbolik ⊲ Mengenalkan Persamaan Diferensial Parsial Parabolik

8.1

Pendahuluan

Dalam bab ini, penulisan ’persamaan diferensial parsial’ akan dipersingkat menjadi PDP. PDP dapat dibagi menjadi 3 jenis, yaitu persamaan diferensial eliptik, parabolik dan hiperbolik. PDP eliptik dinyatakan sebagai berikut ∂2u ∂2u (x, y) + (x, y) = f (x, y) ∂x2 ∂y 2

(8.1)

Di bidang fisika, persamaan (8.1) disebut juga dengan Persamaan Poisson. Jika f (x, y)=0, maka diperoleh persamaan yang lebih sederhana ∂2u ∂2u (x, y) + (x, y) = 0 ∂x2 ∂y 2

(8.2)

yang dikenal sebagai Persamaan Laplace. Contoh masalah PDP eliptik di bidang fisika adalah problem distribusi panas pada kondisi steady-state di dalam obyek 2-dimensi dan 3-dimensi. Adapun, PDP parabolik dinyatakan sebagai berikut ∂2u ∂u (x, t) − α2 2 (x, t) = 0 ∂t ∂x

(8.3)

Fenomena fisis yang bisa dijelaskan oleh persamaan ini adalah masalah aliran panas dalam suatu obyek. Sedangkan PDP hiperbolik dinyatakan sebagai berikut α2

∂2u ∂2u (x, t) = (x, t) ∂2x ∂t2 123

(8.4)

BAB 8. PERS. DIFERENSIAL PARSIAL NUMERIK

124

, biasa digunakan untuk menjelaskan fenomena gelombang. Sekarang, mari kita bahas lebih dalam satu-persatu, difokuskan pada bagaimana cara mengkonversi semua PDP menjadi PDP versi Finite-Difference, (FD).

8.2

PDP eliptik dalam Finite-Difference

Kita mulai dari persamaan aslinya ∂2u ∂2u (x, y) + (x, y) = f (x, y) ∂x2 ∂y 2

(8.5)

dimana R = [(x, y)|a < x < b, c < y < d]. Maksudnya, variasi titik-titik x dibatasi antara a dan b. Demikian pula dengan variasi titik-titik y, dibatasi antara c dan d. Jika h adalah jarak antar titik yang bersebelahan pada titik-titik dalam rentang a dan b, maka titik-titik variasi di antara a dan b dapat diketahui melalui rumus ini xi = a + ih,

dimana i = 1, 2, . . . , n

(8.6)

dimana a adalah titik awal pada sumbu horisontal x. Demikian pula pada sumbu y. Jika k adalah jarak antar titik yang bersebelahan pada titik-titik dalam rentang c dan d, maka titiktitik variasi di antara c dan d dapat diketahui melalui rumus ini yj = c + jk,

dimana j = 1, 2, . . . , m

(8.7)

dimana c adalah titik awal pada sumbu vertikal y. Garis-garis yang sejajar sumbu horisontal, y = yi dan garis-garis yang sejajar sumbu vertikal, x = xi disebut grid lines. Sementara titiktitik perpotongan antara garis-garis horisontal dan vertikal dinamakan mesh points. Turunan kedua sebagaimana yang ada pada persamaan (8.5) dapat dinyatakan dalam rumus centered-difference sebagai berikut u(xi+1 , yj ) − 2u(xi , yj ) + u(xi−1 , yj ) h2 ∂ 4 u ∂2u (x , y ) = − (ξi , yj ) i j ∂x2 h2 12 ∂x4

(8.8)

u(xi , yj+1 ) − 2u(xi , yj ) + u(xi , yj−1 ) k 2 ∂ 4 u ∂2u (xi , yj ) = − (xi , ηj ) 2 ∂y k2 12 ∂y 4

(8.9)

Kemudian, untuk menyederhanakan analisis menggunakan metode Finite-Difference, biasanya suku yang terakhir diabaikan. u(xi+1 , yj ) − 2u(xi , yj ) + u(xi−1 , yj ) ∂2u (xi , yj ) = ∂x2 h2

(8.10)

u(xi , yj+1 ) − 2u(xi , yj ) + u(xi , yj−1 ) ∂2u (xi , yj ) = 2 ∂y k2

(8.11)

Pengabaian suku terakhir otomatis menimbulkan error yang dinamakan truncation error. Jadi, ketika suatu persamaan diferensial diolah secara numerik dengan metode Finite-Difference,

8.2. PDP ELIPTIK DALAM FINITE-DIFFERENCE

125

maka solusinya pasti meleset/keliru "sedikit", karena pada hasilnya pasti mengandung truncation error. Akan tetapi, keberadaan error tersebut dapat ditolerir hingga batas-batas tertentu yang uraiannya akan dikupas pada bagian akhir bab ini. Ok. Mari kita lanjutkan! Sekarang persamaan (8.10) dan (8.11) disubstitusikan ke persamaan (8.5) u(xi+1 , yj ) − 2u(xi , yj ) + u(xi−1 , yj ) u(xi , yj+1 ) − 2u(xi , yj ) + u(xi , yj−1 ) + = f (xi , yj ) (8.12) h2 k2 dimana i = 1, 2, ..., n − 1 dan j = 1, 2, ..., m − 1 dengan syarat batas sebagai berikut u(x0 , yj ) = g(x0 , yj )

u(xn , yj ) = g(xn , yj )

j = 0, 1, ..., m;

u(xi , y0 ) = g(xi , y0 )

u(xi , ym ) = g(xi , ym )

i = 1, 2, ..., n − 1.

Dalam metode Finite-Difference, persamaan (8.12) dinyatakan sebagai berikut wi,j+1 − 2wi,j + wi,j−1 wi+1,j − 2wi,j + wi−1,j + 2 h k2 2 h wi+1,j − 2wi,j + wi−1,j + 2 (wi,j+1 − 2wi,j + wi,j−1 ) k h2 h2 h2 wi+1,j − 2wi,j + wi−1,j + 2 wi,j+1 − 2 2 wi,j + 2 wi,j−1 k k k h2 h2 −2[1 + 2 ]wi,j + (wi+1,j + wi−1,j ) + 2 (wi,j+1 + wi,j−1 ) k k h2 h2 2[1 + 2 ]wi,j − (wi+1,j + wi−1,j ) − 2 (wi,j+1 + wi,j−1 ) k k

= f (xi , yj ) = h2 f (xi , yj ) = h2 f (xi , yj ) = h2 f (xi , yj ) = −h2 f (xi , yj )

(8.13)

dimana i = 1, 2, ..., n − 1 dan j = 1, 2, ..., m − 1, dengan syarat batas sebagai berikut w0,j = g(x0 , yj )

wn,j = g(xn , yj )

j = 0, 1, ..., m;

wi,0 = g(xi , y0 )

wi,m = g(xi , ym )

i = 1, 2, ..., n − 1.

Persamaan (8.13) adalah rumusan akhir metode Finite-Difference untuk PDP Eliptik. Contoh pertama Misalnya kita diminta mensimulasikan distribusi panas pada lempengan logam berukuran 0, 5 m x 0, 5 m. Temperatur pada 2 sisi tepi lempengan logam dijaga pada temperatur 0◦ C, sementara pada 2 sisi tepi lempengan logam yang lain temperaturnya diatur meningkat secara linear dari 0◦ C hingga 100◦ C. Problem ini memenuhi PDP Eliptik: ∂2u ∂2u (x, y) + 2 (x, y) = 0 2 ∂x ∂y dengan syarat-syarat batas u(0, y) = 0,

u(x, 0) = 0,

u(x, 0.5) = 200x,

u(0.5, y) = 200y

(8.14)

Jika n = m = 4, masalah ini dipecahkan dengan grid sebagaimana yang diperlihatkan pada


126 Y

U(0,y)=0

W0,3 W0,2 W0,1

W2,4

W1,4

W3,4

W1,3

W2,3

W3,3

W1,2

W2,2

W3,2

W1,1

W2,1

W3,1

W1,0

W2,0

W3,0

W4,3 W4,2 W4,1

0.5

U(0.5,y)=200y

0.5

U(x,0.5)=200x

X

U(x,0)=0 Gambar 8.1: Distribusi temperatur pada lempeng logam

Gambar 8.1. Adapun persamaan finite-differencenya adalah 4wi,j − wi+1,j − wi−1,j − wi,j−1 − wi,j+1 = 0

(8.15)

dimana i = 1, 2, 3 dan j = 1, 2, 3. Dengan menerapkan persamaan (8.15) pada setiap mesh point yang belum diketahui temperaturnya, diperoleh 4w1,3 − w2,3 − w1,2 = w0,3 + w1,4 4w2,3 − w3,3 − w2,2 − w1,3 = w2,4 4w3,3 − w3,2 − w2,3 = w4,3 + w3,4 4w1,2 − w2,2 − w1,1 − w1,3 = w0,2 4w2,2 − w3,2 − w2,1 − w1,2 − w2,3 = 0 4w3,2 − w3,1 − w2,2 − w3,3 = w4,2 4w1,1 − w2,1 − w1,2 = w0,1 + w1,0 4w2,1 − w3,1 − w1,1 − w2,2 = w2,0 4w3,1 − w2,1 − w3,2 = w3,0 + w4,1 Syarat batas yang telah ditentukan adalah w1,0 = w2,0 = w3,0 = w0,1 = w0,2 = w0,3 = 0, w1,4 = w4,1 = 25,

w2,4 = w4,2 = 50,

and

w3,4 = w4,3 = 75

8.2. PDP ELIPTIK DALAM FINITE-DIFFERENCE

127

Dengan menyertakan syarat batas, sistem persamaan linear dapat ditulis kembali sebagai berikut 4w1,3 − w2,3 − w1,2 = 25 4w2,3 − w3,3 − w2,2 − w1,3 = 50 4w3,3 − w3,2 − w2,3 = 150 4w1,2 − w2,2 − w1,1 − w1,3 = 0 4w2,2 − w3,2 − w2,1 − w1,2 − w2,3 = 0 4w3,2 − w3,1 − w2,2 − w3,3 = 50 4w1,1 − w2,1 − w1,2 = 0 4w2,1 − w3,1 − w1,1 − w2,2 = 0 4w3,1 − w2,1 − w3,2 = 25 Kemudian sistem persamaan linear tersebut dikonversikan ke dalam operasi perkalian matrik                  

4 −1

−1 4

0 −1

−1 0

0

0

0

0

−1

0

0

0

0



 0    0 −1 4 0 0 −1 0 0 0    −1 0 0 4 −1 0 −1 0 0    0 −1 0 −1 4 −1 0 −1 0    0 0 −1 0 −1 4 0 0 −1     0 0 0 −1 0 0 4 −1 0    0 0 0 0 −1 0 −1 4 −1   0 0 0 0 0 −1 0 −1 4

w1,3





  w2,3      w3,3      w1,2      w2,2  =   w3,2      w1,1      w2,1   w3,1

25



 50   150    0   0    50   0    0  25

Mari kita perhatikan sejenak, susunan elemen-elemen angka pada matrik berukuran 9x9. Terlihat jelas bahwa angka diagonal selalu bernilai 4. Ini sama sekali bukan ketidaksengajaan. Melainkan susunan itu sengaja direkayasa sedemikian rupa sehingga elemen-elemen tri-diagonal terisi penuh dan pada diagonal utamanya diletakkan angka yang terbesar. Akhirnya dengan metode Eliminasi Gauss diperoleh hasil w1,3

w2,3

w3,3

w1,2

w2,2

w3,2

w1,1

w2,1

w3,1

18.7500 37.5000 56.2500 12.5000 25.0000 37.5000 6.2500 12.5000 18.7500 Catatan: inilah solusi yang "ditawarkan" oleh Finite-Difference. Kalau diamati dengan teliti, angka-angka distribusi temperatur pada 9 buah mesh points memang logis dan masuk akal. Dalam kondisi real, mungkin kondisi seperti ini hanya bisa terjadi bila lempengan logam tersebut terbuat dari bahan yang homogen isotropik. Berikut adalah script Matlab untuk menghitung nila-nilai w menggunakan metode Eliminasi Gauss. 1 2

clear all clc

128 3 4 5 6 7 8 9 10 11 12


n=9; A=[ 4 -1 0 -1 0 0 0 0 0; -1 4 -1 0 -1 0 0 0 0; 0 -1 4 0 0 -1 0 0 0; -1 0 0 4 -1 0 -1 0 0; 0 -1 0 -1 4 -1 0 -1 0; 0 0 -1 0 -1 4 0 0 -1; 0 0 0 -1 0 0 4 -1 0; 0 0 0 0 -1 0 -1 4 -1; 0 0 0 0 0 -1 0 -1 4];

13 14

b=[25; 50; 150; 0; 0; 50; 0; 0; 25];

15 16 17 18 19 20 21

%&&&&&& Proses Eliminasi Gauss &&&&&&&&&&&&&&&&&&&&&&&&& %====== Menggabungkan Vektor b kedalam matrik A ======== %====== sehingga terbentuk matrik Augmentasi. ======== for i=1:n A(i,n+1)=b(i,1); end

22 23 24

%---------Proses Triangularisasi----------for j=1:(n-1)

25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44


45 46 47


48 49 50 51 52 53 54 55 56

for i=n-1:-1:1 S=0; for j=n:-1:i+1 S=S+A(i,j)*x(j,1); end x(i,1)=(A(i,n+1)-S)/A(i,i); end %&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

57 58 59

%===== Menampilkan Vektor w ================= w=x

Sementara berikut ini adalah script Matlab untuk menghitung nila-nilai w menggunakan metode

8.2. PDP ELIPTIK DALAM FINITE-DIFFERENCE Iterasi Gauss-Seidel. 1 2

clear all clc

3 4 5 6 7 8 9 10 11 12 13

n=9; A=[ 4 -1 0 -1 0 0 0 0 0; -1 4 -1 0 -1 0 0 0 0; 0 -1 4 0 0 -1 0 0 0; -1 0 0 4 -1 0 -1 0 0; 0 -1 0 -1 4 -1 0 -1 0; 0 0 -1 0 -1 4 0 0 -1; 0 0 0 -1 0 0 4 -1 0; 0 0 0 0 -1 0 -1 4 -1; 0 0 0 0 0 -1 0 -1 4];

14 15

b=[25; 50; 150; 0; 0; 50; 0; 0; 25];

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57

%&&&&&&& ITERASI GAUSS-SEIDEL &&&&&&&&&&&&&&&&&& itermax=100; %iterasi maksimum %----nilai awal----------xl=[0; 0; 0; 0; 0; 0; 0; 0; 0]; xb=xl; %----stopping criteria----------sc=0.001; %----memulai iterasi------------for iterasi=1:itermax smtr1=0; for j=2:n smtr1=smtr1+A(1,j)*xl(j,1); end xb(1,1)=(-smtr1+b(1,1))/A(1,1); %---------------------------------------------for i=2:n-1 smtr2=0; for j=i+1:n smtr2=smtr2-A(i,j)*xl(j,1); end smtr3=0; for k=1:i-1 smtr3=smtr3-A(i,k)*xb(k,1); end xb(i,1)=(smtr3+smtr2+b(i,1))/A(i,i); end %---------------------------------------------smtr4=0; for k=1:n-1 smtr4=smtr4-A(n,k)*xb(k,1); end xb(n,1)=(smtr4+b(n,1))/A(n,n); %------perhitungan norm2 ------------s=0; for i=1:n s=s+(xb(i,1)-xl(i,1))^2; end epsilon=sqrt(s); %------------------------------------xl=xb; %------memeriksa stopping criteria--------

129


130

if epsilon<sc w=xb break end %-----------------------------------------

58 59 60 61 62 63

end

Tabel berikut memperlihatkan hasil pemrosesan dengan metode Eliminasi Gauss (disingkat: EG) dan iterasi Gauss-Seidel (disingkat: GS) w1,3

w2,3

w3,3

w1,2

w2,2

w3,2

w1,1

w2,1

w3,1

EG 18.7500 37.5000 56.2500 12.5000 25.0000 37.5000 6.2500 12.5000 18.7500 GS

18.7497 37.4997 56.2498 12.4997 24.9997 37.4998 6.2498 12.4998 18.7499

Memang berbeda. Tapi tak ada perbedaan yang cukup berarti. Keduanya nyaris sama persis. Namun perlu saya tegaskan disini bahwa jika sistem persamaan linear yang diperoleh dari Finite Difference berorde 100 atau kurang dari itu, maka lebih baik memilih metode Eliminasi Gauss sebagai langkah akhir Finite-Difference. Alasannya karena, direct method seperti eliminasi Gauss lebih stabil dibandingkan metode iterasi. Tapi jika ordenya lebih dari 100, disarankan memilih metode iterasi seperti iterasi Gauss-Seidel, atau menggunakan metode SOR yang terbukti lebih efisien dibanding Gauss-Seidel. Jika matrik A bersifat positive definite, metode Court Factorization adalah pilihan yg paling efisien dan menghemat memori komputer.

Bab 9

Integral Numerik

✍ Objektif : ⊲ Mengenalkan metode Trapezoida. ⊲ Mengenalkan metode Simpson. ⊲ Mengenalkan metode Composite-Simpson.

9.1

Metode Trapezoida

Suatu persamaan integral Z

b

f (x)dx

(9.1)

a

disebut numerical quadrature. Pendekatan numerik untuk menyelesaikan integral tersebut adalah

i=0 X

ai f (xi )

(9.2)

n

Metode pendekatan yang paling dasar dalam integral numerik adalah metode Trapezoida Z

b

f (x)dx =

a

h3 h [f (x0 ) + f (x1 )] − f ′′ (ξ) 2 12

(9.3)

dimana x0 = a, x1 = b dan h = b − a. Karena bagian error pada Trapezoida adalah f ′′ , maka

pendekatan Trapezoida bekerja efektif pada fungsi-fungsi yang turunan kedua-nya bernilai nol.

9.2

Metode Simpson

Metode pendekatan yang lebih baik dalam integral numerik adalah metode Simpson Z

b

f (x)dx = a

h5 h [f (x0 ) + 4f (x1 ) + f (x2 )] − f 4 (ξ) 3 90 131

(9.4)

BAB 9. INTEGRAL NUMERIK

132

Gambar 9.1: Metode Trapezoida

Gambar 9.2: Metode Simpson dengan x0 = a, x2 = b, dan x1 = a + h dimana h = (b − a)/2. Contoh

Metode Trapezoida untuk fungsi f pada interval [0,2] adalah Z

0

2

f (x)dx ≈ f (0) + f (2)

dimana x0 = 0, x1 = 2 dan h = 2 − 0 = 2,

sedangkan metode Simpson untuk fungsi f pada interval [0,2] adalah Z

0

2

f (x)dx ≈

1 [f (0) + 4f (1) + f (2)] 3

dengan x0 = 0, x2 = 2, dan x1 = a + h = 1 dimana h = (b − a)/2 = 1.

9.3. METODE COMPOSITE-SIMPSON

133

Tabel berikut ini memperlihatkan evaluasi integral numerik terhadap beberapa fungsi dalam interval [0,2] beserta solusi exact-nya. Jelas terlihat, metode Simpson lebih baik dibanding Trapezoida. x2 2,667 4,000 2,667

f (x) Nilai exact Trapezoida Simpson

x4 6,400 16,000 6,667

1/(x + 1) 1,099 1,333 1,111

√

1 + x2 2,958 3,326 2,964

sin x 1,416 0,909 1,425

ex 6,389 8,389 6,421

Kalau diamati lebih teliti, akan kita dapatkan bahwa interval [0,2] telah dibagi 2 pada metode Simpson, sementara pada metode Trapesoida tidak dibagi sama sekali. Sebenarnya dengan membagi interval lebih kecil lagi, maka error -nya akan semakin kecil. Misalnya, banyaknya pembagian interval dinyatakan dengan n ketika n = 1: Trapesioda Z

x1

h h3 [f (x0 ) + f (x1 )] − f ′′ (ξ) 2 12

(9.5)

h h5 [f (x0 ) + 4f (x1 ) + f (x2 )] − f 4 (ξ) 3 90

(9.6)

f (x)dx = x0

ketika n = 2: Simpson Z

x2

f (x)dx =

x0

ketika n = 3: Simpson tiga-per-delapan Z

x3

f (x)dx =

x0

3h 3h5 4 [f (x0 ) + 3f (x1 ) + 3f (x2 ) + f (x3 )] − f (ξ) 8 80

ketika n = 4: Z x4 8h7 6 2h [7f (x0 ) + 32f (x1 ) + 12f (x2 ) + 32f (x3 ) + 7f (x4 )] − f (ξ) f (x)dx = 45 945 x0

(9.7)

(9.8)

Keempat bentuk persamaan integral numerik di atas dikenal dengan closed Newton-Cotes formulas. Keterbatasan metode Newton-Cotes terlihat dari jumlah pembagian interval. Di atas tadi pembagian interval baru sampai pada n = 4. Bagaimana bila interval evaluasinya dipersempit supaya solusi numeriknya lebih mendekati solusi exact? Atau dengan kata lain n > 4.

9.3

Metode Composite-Simpson

Coba perhatikan contoh berikut ini. Tentukan solusi numerik dari

R4 0

ex dx. Metode Simpson

dengan h = 2 (atau interval evaluasi integral dibagi 2 , n = 2) memberikan hasil Z

4 0

ex dx ≈

¢ 2¡ 0 e + 4e2 + e4 = 56, 76958 3


134

Gambar 9.3: Metode Composite Simpson

Padahal solusi exact dari integral tersebut adalah e4 − e0 = 53, 59815, artinya terdapat er-

ror sebesar 3,17143 yang dinilai masih terlampau besar untuk ditolerir. Bandingkan dengan metode yang sama namun dengan h = 1 (atau interval evaluasi integral dibagi 4 , n = 4) Z

4

x

e dx =

0

Z

2

x

e dx + 0

Z

4

ex dx

2

¢ 1¡ 2 ¢ 1¡ 0 e + 4e + e2 + e + 4e3 + e4 3 3 ¢ 1¡ 0 2 e + 4e + 2e + 4e3 + e4 = 3 = 53, 86385 ≈

Hasil ini memperlihatkan error yang makin kecil, yaitu menjadi 0,26570. Jadi dengan memperkecil h, error menjadi semakin kecil dan itu artinya solusi integral numerik semakin mendekati solusi exact. Sekarang kita coba kecilkan lagi nilai h menjadi h =

1 2

(atau interval evaluasi in-

tegral dibagi 8 , n = 8), Z

4

x

1

Z

2

Z

3

Z

4

ex dx e dx + e dx + 3 2 1 0 ´ 1³ ´ 1³ 0 1/2 3/2 e + 4e + e + e + 4e + e2 + ≈ 6 6 ´ 1³ ´ 1³ 2 e + 4e5/2 + e3 + e3 + 4e7/2 + e4 6 6 ´ 1³ 0 1/2 e + 4e + 2e + 4e3/2 + 2e2 + 4e5/2 + 2e3 + 4e7/2 + e4 = 6 = 53, 61622

e dx = 0

Z

x

e dx +

x

x

dan seperti yang tadi kita simpulkan, error -nya semakin kecil menjadi 0,01807.

9.3. METODE COMPOSITE-SIMPSON

135

Prosedur ini dapat digeneralisir menjadi suatu formula sebagai berikut Z

b

f (x)dx = a

n/2 Z X j=1

=

x2j

f (x)dx

x2j−2

n/2 ½ X h j=1

¾ h5 (4) [f (x2j−2 ) + 4f (x2j−1 ) + f (x2j )] − f (ξj ) 3 90

(9.9)

dimana h = (b − a)/n dan xj = a + jh, untuk j = 1, ..., n/2, dengan x0 = a dan xn = b. Formula ini dapat direduksi menjadi Z

a

b

  n/2 (n/2)−1 n/2 X X h5 X (4) h f (x2j−1 ) + f (xn ) − f (x2j ) + 4 f (x0 ) + 2 f (ξj ) f (x)dx = 3 90 j=1

j=1

Formula ini dikenal sebagai metode Composite Simpson.

j=1

(9.10)


136 Tugas#4

1. Hitunglah integral-integral berikut ini dengan metode Composite Simpson! a. b. c. d. e.

Z

Z

Z

2

x ln xdx,

n=4

1 2 0 3 1 2

Z

2 dx, +4 x dx, 2 x +4

n=6

x2

x3 ex dx,

−2 Z 3π/8

n=8 n=4

tan xdx,

n=8

0

f.

Z

5

√

3

1

x2 − 4

dx,

n=8

2. Tentukan nilai n dan h untuk mengevaluasi Z

2

e2x sin 3xdx

0

dengan metode Composite Simpson, bila error yang ditolerir harus lebih kecil dari 10−4 . 3. A particle of mass m moving through a fluid is subjected to a viscous resistance R, which is a function of the velocity v. The relationship between the resistance R, velocity v, and time t is given by the equation: t=

Z

v(t)

v(t0 )

m du R(u)

√ Suppose that R(v) = −v v, for a particular fluid, where R is in newtons and v in meters/second.

If m = 10 kg and v(0) = 10 m/s, approximate the time required for the particle to slow to v = 5 m/s.

Bab 10

Metode Newton


10.1

Penyederhanaan

Metode Newton sangat populer dan powerfull untuk mencari akar suatu fungsi yang kontinyu. Ada banyak jalan untuk memperkenalkan metode ini. Salah satunya bisa didahului mulai dari deret Taylor atau polinomial Taylor. Suatu fungsi yang kontinyu dapat dinyatakan dalam deret Taylor sebagai berikut f (x) = f (¯ x) + (x − x ¯)f ′ (¯ x) + 0 = f (¯ x) + (p − x ¯)f ′ (¯ x) +

(x − x ¯)2 ′′ f (ξ(x)) 2

(p − x ¯)2 ′′ f (ξ(p)) 2

0 = f (¯ x) + (p − x ¯)f ′ (¯ x) p−x ¯=− p≈x ¯− pn = pn−1 −

f (x) f ′ (¯ x)

f (x) f ′ (¯ x)

f (pn−1 ) f ′ (pn−1 ) 137

,

n≥1

BAB 10. METODE NEWTON

138

Gambar 10.1: Metode Newton Demikianlah catatan singkat dari saya tentang metode Finite-Difference dengan ekstrapolasi Richardson untuk problem linear. Saya cukupkan sementara sampai disini. Insya Allah akan saya sambung dengan problem non-linear dan konsep ekstrapolasi. Kalau ada yang mau didiskusikan, silakan hubungi saya melalui email: [email protected].

Bab 11

Metode Monte Carlo


11.1

Penyederhanaan

Kita awali pembahasan metode Monte Carlo dengan mengetengahkan contoh yang sangat terkenal yaitu menghitung luas suatu lingkaran. Fugure 1 memperlihatkan lingkaran dengan radius r = 1 berada di dalam kotak bujursangkar. Luas lingkaran adalah πr2 = π(1)2 = π sementara luas bujursangkar adalah (2)2 = 4. Rasio antara luas lingkaran dan luas bola adalah ρ=

π luas lingkaran = = 0, 7853981633974483 luas bujursangkar 4

(11.1)

Jadi, dengan mengetahui nilai ρ, maka kita bisa menghitung luas lingkaran dengan cara luas lingkaran = ρ × luas bujursangkar

(11.2)

Bayangkan anda punya satu set permainan dart. Anda lemparkan sejumlah dart ke arah lingkaran tadi. Misalnya, total dart yang menancap di papan dart ada 1024 buah. Sebanyak 812 dart berada di dalam lingkaran, dan yang lainnya di luar lingkaran. Rasio antara keduanya ρ=

total

812 dart di dalam lingkaran = = 0, 79296875 dart di dalam bujursangkar 1024

139

(11.3)

BAB 11. METODE MONTE CARLO

140

Gambar 11.1: Lingkaran dan bujursangkar

Gambar 11.2: Dart yang menancap pada bidang lingkaran dan bujursangkar

Dengan pendekatan ke persamaan (11.2) maka luas lingkaran adalah luas lingkaran = ρ × luas bujursangkar = 0, 79296875 × 4 = 3, 171875 Apakah angka ini make sense ? Mungkin anda masih ragu. Sekarang mari kita coba hitung nilai π dengan mengacu pada rumus di atas. Kita sepakati saja bahwa dart yang berada di dalam lingkaran mesti memenuhi x2i + yi2 ≤ 1. Dalam perhitungan, semua dart diganti dengan bi-

langan acak (random number ). Dari 1000 dart, yang masuk lingkaran ada 787 buah, sehingga,

mengacu persamaan (11.3) ρ=

787 = 0, 787 1000

11.1. PENYEDERHANAAN

141

Gambar 11.3: Dart yang menancap pada bidang 1/4 lingkaran dan bujursangkar maka berdasarkan persamaan (11.1) π = ρ × 4 = 0, 787 × 4 = 3, 148 Lumayan akurat bukan? Semakin banyak jumlah dart, semakin akurat nilai π yang anda peroleh. Sekarang mari kita kembangkan metode Monte Carlo ini untuk menghitung luas suatu area yang terletak di bawah garis kurva suatu fungsi f (x). Atau sebut saja menghitung integral suatu fungsi f (x) yang dievaluasi antara batas a dan b. Luas kotak R yang melingkupi luas bidang integral A adalah R = {(x, y) : a ≤ x ≤ b dan

0 ≤ y ≤ d}

(11.4)

dimana d = maksimum f (x)

,

a≤x≤b

(11.5)

Bab 12

Inversi


12.1

Inversi Linear

Diketahui data eksperimen tersaji dalam tabel berikut ini xi 1 2 3 4 5

yi 1,3 3,5 4,2 5,0 7,0

xi 6 7 8 9 10

yi 8,8 10,1 12,5 13,0 15,6

Lalu data tersebut di-plot dalam sumbu x dan y. Sekilas, kita bisa melihat bahwa data yang telah di-plot tersebut dapat didekati dengan sebuah persamaan garis, yaitu a1 xi + a0 . Artinya, kita melakukan pendekatan secara linear, dimana fungsi pendekatan-nya adalah P (xi ) = a1 xi + a0

(12.1)

Problemnya adalah berapakah nilai konstanta a1 dan a0 yang sedemikian rupa, sehingga posisi garis tersebut paling mendekati atau bahkan melalui titik-titik data yang telah di-plot di atas? 143

BAB 12. INVERSI

144 16 14 12 10

Y 8 6 4 2 0 1

2

3

4

5

6

7

8

9

10

X

Dengan kata lain, sebisa mungkin yi sama dengan P (xi ) atau dapat diformulasikan sebagai m X

yi − P (xi ) = 0

(12.2)

yi − (a1 xi + a0 ) = 0

(12.3)

i=1

m X i=1

dimana jumlah data, m = 10. Suku yang berada disebelah kiri dinamakan fungsi error (error function), yaitu E(a0 , a1 ) =

m X i=1

yi − (a1 xi + a0 )

(12.4)

Semua data yang diperoleh melalui eksperimen, fungsi error-nya tidak pernah bernilai nol. Jadi, tidak pernah didapatkan garis yang berhimpit dengan semua titik data ekperimen. Namun demikian, kita masih bisa berharap agar fungsi error menghasilkan suatu nilai, dimana nilai tersebut adalah nilai yang paling minimum atau paling mendekati nol. Harapan tersebut diwujudkan oleh metode least square dengan sedikit modifikasi pada fungsi error-nya sehingga menjadi E(a0 , a1 ) =

m X i=1

[yi − (a1 xi + a0 )]2

(12.5)

Agar fungsi error bisa mencapai nilai minimum, maka syarat yang harus dipenuhi adalah: ∂E(a0 , a1 ) =0 ∂ai

(12.6)

dimana i = 0 dan 1, karena dalam kasus ini memang cuma ada a0 dan a1 . Maka mesti ada dua

12.1. INVERSI LINEAR

145

buah turunan yaitu: m ∂E(a0 , a1 ) ∂ X = [yi − (a1 xi + a0 )]2 = 0 ∂a0 ∂a0 i=1

m X (yi − a1 xi − a0 )(−1) = 0 2 i=1

a0 .m + a1

m X

xi =

m X

yi

(12.7)

xi yi

(12.8)

i=1

i=1

dan m ∂ X ∂E(a0 , a1 ) = [yi − (a1 xi + a0 )]2 = 0 ∂a1 ∂a1 i=1

2

m X i=1

(yi − a1 xi − a0 )(−xi ) = 0 a0

m X i=1

xi + a1

m X

x2i

=

i=1

m X i=1

Akhirnya persamaan (12.7) dan (12.8) dapat dicari solusinya berikut ini: a0 =

Pm

dan a1 =

Pm Pm Pm i=1 yi − i=1 xi yi i=1 xi ¡ Pm 2 ¢ Pm 2 m i=1 xi − ( i=1 xi )

2 i=1 xi

m

Pm Pm Pm yi i=1 xi i=1 xi yi − ¡ Pm 2 ¢ Pm i=12 m i=1 xi − ( i=1 xi )

(12.9)

(12.10)

Coba anda bandingkan kedua hasil di atas dengan rumus least square yang terdapat pada buku Praktikum Fisika Dasar keluaran Departemen Fisika-UI. Mudah-mudahan sama persis. OK, berdasarkan data ekperimen yang ditampilkan pada tabel diawal catatan ini, maka didapat: a0 =

385(81) − 55(572, 4) = −0, 360 10(385) − (55)2

(12.11)

10(572, 4) − 55(81) = 1, 538 10(385) − (55)2

(12.12)

dan a1 =

Jadi, fungsi pendekatan-nya, P (xi ), adalah P (xi ) = 1, 538xi − 0, 360

(12.13)

Solusi least square dengan pendekatan persamaan garis seperti ini juga dikenal dengan nama lain yaitu regresi linear. Sedangkan nilai a0 dan a1 disebut koefisien regresi. Gambar di bawah ini menampilkan solusi regresi linear tersebut berikut semua titik datanya Tentu saja anda sudah bisa menduga bahwa selain regresi linear, mungkin saja terdapat regresi

BAB 12. INVERSI

146 16

P(x) = 1.538*x − 0.36

14 12 10 8 6 4 2 0 −2 0

2

4

6

8

10

parabola atau quadratik dimana fungsi pendekatannya berupa persamaan parabola, yaitu: P (xi ) = a2 x2i + a1 xi + a0

(12.14)

dimana koefisien regresinya ada tiga yaitu a0 , a1 dan a2 . Kalau anda menduga demikian, maka dugaan anda benar! Bahkan sebenarnya tidak terbatas sampai disitu. Secara umum, fungsi pendekatan, P (xi ), bisa dinyatakan dalam aljabar polinomial berikut ini: P (xi ) = an xni + an−1 xn−1 + ... + a2 x2i + a1 xi + a0 i

(12.15)

Namun untuk saat ini, saya tidak ingin memperluas pembahasan hingga regresi parabola, dan polinomial. Saya masih ingin melibatkan peranan metode eliminasi gauss dalam menyelesaikan problem least square seperti yang selalu saya singgung pada catatan-catatan kuliah saya yang terdahulu. Nah, kalau metode eliminasi gauss hendak digunakan untuk mencari solusi regresi linear, kita bisa mulai dari persamaan (12.7) dan (12.8), yaitu: a0 .m + a1 a0

m X

xi + a1

i=1

m X

i=1 m X

xi = x2i =

m X

i=1 m X

yi xi yi

i=1

i=1

Keduanya bisa dinyatakan dalam operasi matrik: "

m Pm

i=1 xi

Pm

i=1 xi Pm 2 i=1 xi

#"

a0 a1

#

=

" P m

yi Pmi=1 i=1 xi yi

#

(12.16)

Kalau anda mengikuti catatan-catatan terdahulu, pasti anda tidak asing lagi dengan dengan semua elemen-elemen matrik di atas. Semua sudah saya ulas pada catatan yang berjudul Aplikasi Elimininasi Gauss: Model Garis. Silakan anda lanjutkan perhitungan matrik tersebut hingga diperoleh koefisien regresi a0 dan a1 . Selamat mencoba!

12.2. INVERSI NON-LINEAR

12.2

147

Inversi Non-Linear

Persamaan least squares linear adalah sebagai berikut: [Gt G]δm = Gt δd

(12.17)

Persamaan least squares non-linear dapat dinyatakan sebagai berikut: [Gt G + λI]δm = Gt δd

(12.18)

dimana G adalah matrik kernel, namun dia juga biasa dikenal dengan sebutan matrik Jacobian, sementara λ adalah faktor pengali Lagrange, dan I adalah matrik identitas yang ordenya disesuaikan dengan Gt G. Adapun definisi δm dan δd akan dijelaskan pada bagian akhir catatan ini. Langkah-langkah untuk menyelesaikan problem least squares non-linear adalah: 1. Menentukan model, misal f (x) = xm 2. Menghitung jacobian, G. Caranya adalah menghitung turunan pertama dari model terhadap model-parameter, m. Sesuai permisalan pada point 1, didapat A=

∂f (m) = xm ln(x) ∂m

(12.19)

3. Membuat perhitungan simulasi, misalnya ditentukan m = 2. Nilai m adalah nilai yang hendak dicari. Dalam simulasi, nilai m dianggap sudah diketahui bahkan ditentukan. Lalu hitunglah f (x) = xm dengan x bergerak dari x = 1, 2, 3.., 10. Jadi, nanti akan didapat 10 buah f (x). Mau lebih dari 10 juga boleh, terserah saja. Hasil hitungannya dikasih nama d, jadi d = f (x). Karena dalam simulasi ini x-nya bergerak hanya sampai 10, maka hasilnya mesti ada 10 d, yaitu d1 , d2 , .., d10 . 4. Buatlah perhitungan untuk m sembarang, misal mula-mula dipilih m = 5. Ini adalah nilai awal dari m yang akan diiterasikan sedemikian rupa hingga nantinya m akan menuju 2 sesuai dengan nilai m pada simulasi (point 3). Bagusnya dibedakan penulisannya, atau tulis saja m0 = 5, dimana m0 maksudnya adalah m mula-mula. Lalu hitung lagi nilai 0

f (x) = xm . Sekarang dinamakan dc = f (x). Jangan lupa bahwa saat perhitungan, nilai x bergerak dari 1 sampai 10. Jadi, nanti didapat 10 dc . 5. Hitunglah δd, dimana δd = dc − d. Sebelumnya sudah dinyatakan bahwa dc ada 10 buah, demikian juga d ada 10 buah, maka δd harus ada 10 buah juga.

6. Selanjutnya hitung ||δd|| yang rumusnya seperti ini ||δd|| =

1 1 Σ(dc − d)2 = Σδd2 N N

(12.20)

dimana N = 10 karena δd-nya ada 10. Rumus ini tidak mutlak harus demikian, anda bisa juga menggunakan norm 2, ℓ2 .

BAB 12. INVERSI

148

7. Tentukan nilai epsilon, ǫ, misal ǫ = 0.000001. Lalu lakukan evaluasi sederhana. Cek, apakah ||δd|| < ǫ ? Pasti awalnya ||δd|| > ǫ, kenapa? Karena m 6= m0 . Kalau begini

situasinya, δd yang ada 10 biji itu dimasukan kedalam proses berikutnya. 8. Hitunglah operasi matriks berikut ini untuk mendapatkan δm [Gt G + λI]δm = Gt δd

(12.21)

dengan λ-nya dikasih nilai sembarang antara 0 dan 1, misalnya λ = 0.005. Perhitungan ini bisa diselesaikan dengan metode eliminasi gauss. 9. Ganti nilai m0 menjadi m1 sesuai dengan rumus m1 = m0 + δm

(12.22)

Nah, m1 ini dimasukan ke proses yang dijelaskan pada point 4 kemudian proses diulangi hingga point 9, begitu seterusnya. Dari sinilah dimulai proses iterasi. Iterasi akan berhenti bila ||δd|| < ǫ. Pada saat itu, nilai mk akan mendekati m = 2 sesuai dengan m simulasi.

Selamat mencoba! Saya juga telah menulis beberapa persamaan non-linear sebagai bahan latihan. Lihat saja di Latihan 1. Tapi tolong diperiksa lagi, apakah jacobiannya sudah benar atau ada kekeliruan. Selanjutnya, kalau ada pertanyaan atau komentar, silakan kirim ke [email protected]

Daftar Pustaka

[1] Burden, R.L. and Faires, J.D., (2001), Numerical Analysis, Seventh Edition, Brooks/Cole, Thomson Learning Academic Resource Center. [2] Haliday and Resnick, (2001), Fundamental of Physics, Brooks/Cole, Thomson Learning Academic Resource Center.

149

Indek

Positive-definite, 5 Transpose, 3 Tridiagonal, 4 Vektor-baris, 6 Vektor-kolom, 6

151

Komputasi untuk Sains dan Teknik

Recommend Documents