INVERS Z-MATRIKS. Pangeran Sianipar

Bulletin of Mathematics Vol. 01, No. 01 (2009), pp. 1–14.

INVERS Z-MATRIKS

Pangeran Sianipar Abstract. Discussions about convergence of exponents of a matrix generally is not free from the problem of spectral radius. In nonnegative matrices, the characteristics of inverses play an important role in the development of knowledge that had relations with convergence. Some of Z-Matrices that had the Principal minor had the big contribution in the determination of the basic characteristics from this convergence. The analysis in the problem invers Z-Matrices that had the Principal minor very useful to study the characteristics of the transition function in the discussions about the Markov Chain.

1. PENDAHULUAN Persoalan Z-Matriks yang mempunyai principal minor biasanya dikaji dalam permasalahan matriks tidak negatif, dimana selalu kita jumpai misalnya dalam problema-problema Rantai Markov dalam Statistika dan Operasi Riset, Model Input-Output dalam Ekonomi Perusahaan dan Convergensi dalam Teknik Numerik dalam menyelesaikan sistem linier maupun tidak linier serta banyak lagi pemakaian yang lainnya [1]. Penentuan invers suatu matriks dalam suatu kasus berperan penting dalam menyelesaikan problema seperti di atas. Sedangkan dalam permasalahan yang dihadapi seperti di atas untuk menentukan suatu harga pendekatan yang diinginkan tidak terlepas dari permasalahan konvergensi. Received 12-03-2008, Accepted 15-09-2008. 2000 Mathematics Subject Classification: 15A23, 15A48 Key words and Phrases: Z-Matriks.

1

Sianipar – Invers Z-Matriks

2

Z-Matriks yang mempunyai Principal minor didefinisikan oleh Franklin A. Graybill tahun 1983 adalah suatu Z-matriks (aij ≤ 0 untuk i 6= j ) dan semua prinsipal minor adalah positip dikenal dengan M -Matriks adalah semi convergen. Jadi dari keadaan tersebut timbul suatu hasrat untuk meneliti sifat-sifat yang dipunyai oleh invers-invers dari Z-Matriks yang mempunyai prinsipal minor . Dalam masalah ini haruslah diingat bahwa semua Z-Matriks yang mempunyai prinsipal minor positip, tetapi tidak semua ZMatriks memberi sifat ini. Oleh karena itu, permasalahan disini adalah penentuan karakter dari Z-Matriks yang mempunyai prinsipal minor sehingga akan didapatkan inversnya. 2. LANDASAN TEORI Definisi 2.1 Misalkan Am×n adalah sebuah matriks. Jika dihapus sebanyak r baris dan s kolom, maka matriks tersebut dinamakan submatriks A. Jika Amxn sebuah matriks dengan baris dan kolom ke-k dari matriks A dihapus, maka matriks tersebut dinamakan principal matriks A orde k. sedangkan determinannya tersebut principal minor A orde k. jika Am×n sebuah matriks dengan baris dan kolom n − k terakhir dihapus, maka A disebut leading principal matriks orde k dan determinannya disebut leading principal minor order k dari A. Definisi 2.2 Misalkan A sebuah matriks kuadrat sedemikian hingga aij ≤ 0 untuk semua i 6= j dan a > 0 untuk semua i = j , maka A disebut Z-Matriks. Teorema 2.1 Jika A adalah sebuah Z-Matriks dan D sebuah matriks diagonal dengan semua entri diagonal utamanya bernilai positif, maka AD adalah sebuah Z-Matriks. Teorema 2.2 Jika A adalah sebuah Z-Matriks dan terdapat sebuah matriks segitiga bawah L dan sebuah matriks segitiga atas U dengan setiap entri diagonal utamanya bernilai positip sedemikian hingga A = LU . Maka tidak ada entri dari L−1 dan U −1 yang negatif dan entri diagonal dari L−1 dan U −1 semuanya positif. Definisi 2.3 Suatu matriks A = [ai,j ] ∈ Rn×n disebut matriks ultrametric umum jika 1. A mempunyai entri taknegatif. 2. ai,j ≥ min{ai,k ; ak,j } untuk semua i, j, k ∈ N.


3

3. Tiap triple {q, s, t} ⊆ N3 dapat berorder triple {i, j, k} sedemikian hingga 3.1. ai,j = ai,k dan ak,j = ak,i ai,j = ai,k dan ak,j = ak,i ai,j = ak,i 3.2. maks{ai,j ; aj,i } ≥ maks{aj,i ai,k ; ak,i } 4. ai,i ≥ maks{ai,j ; aj,i } untuk semua i, j ∈ N Suatu matriks A disebut secara umum kuat ultrametrik matriks jika ketaksamaan (4) dipenuhi untuk semua i 6= j di dalam N, bila jika n = 1 dapat disajikan bahwa ai,j > 0. Definisi 2.4 Suatu matriks A = [ai,j ] ∈ Rn×n adalah ultrametrik matrik umum kuat jika ada bilangan real r sedemikian sehingga A + rξn ξnT adalah ultrametrik umum kuat. Jika pertaksamaan (4) dari definisi 2.5. dipenuhi untuk semua i ∈ N , maka A adalah ultrametrik kuat umum yang digeser. Jika A = [ai,j ] adalah ultrametrik kuat umum yang digeser yang adalah bukan ultrametrik matrik (−1) umum sehimgga r(A) < 0, maka A disebut dari tipe Up,n−p jika 1. ai,j ≤ 0 untuk i, j ∈ N dengan i 6= j. 2. p adalah bilangan positip terbesar sedemikian sehingga ada p × p principal submatriks dari A yang adalah nonsingular M -Matriks. Teorema 2.3 Misalkan A = [ai,j ] ∈ Rn×n matriks nonsingular ultrametrik (−1) umum yang digeser dari type Up,n−p . Jika 0 ≤ p ≤ n maka det(A) < 0. Selanjutnya tiap principal minor dari order m dengan p < m < n adalah takpositif dan ada principal minor positip dari order p. Jika p = n (−1) (yakni A adalah tipe Up,0 ) atau jika A adalah matriks taksingular ultrametrik umum maka det(A) > 0. Teorema 2.4 Misalkan A = [ai,j ] ∈ Rn×n , n > 1 matriks ultrametrik umum yang digeser. Maka ada matriks permutasi P ∈ Rn×n dan bilangan bulat positif r dengan 1 ≤ r < n sedemikan sehingga A A 1,1 1,2 P AP T = A2,1 A2,2 C 0 T = + r(A)ξn ξnT ξr ξn−r = M + r(A)ξn ξnT δ(A)ξn−r ξrT D


4

dimana A1,1 ∈ Rn×n dan A2,2 ∈ Rn×n adalah matriks ultrametrik umum yang digeser. Selanjutnya matriks A1,1 dan A2,2 dapat direduksi dengan jalan yang sama seperti pada A. Blok dari diagonal di luar diagonal utama T adalah A1,1 = r(A)ξr ξn−r dan A2,2 = w(A)ξn−r ξrT . Matriks M ∈ Rn×n , r×r n−r×n−r C ∈ R dan D ∈ R adalah matriks ultrametrik umum dan δ(A)δ(A) = w(A) − r(A) ≥ 0. Selanjutnya W (A1,1 ) ≥ w(A) dan w(A2,2 ) ≥ w(A). Jika A adalah matriks ultrametrik umum kuat yang digeser, maka M adalah matrik ultrametrik umum yang digeser. Akibat 2.1 Jika Z-Matriks A ∈ Rn×n adalah matriks ultrametrik umum yang digeser , maka A adalah nonsingular M -Matriks jika dan hanya jika det(A) > 0. Teorema 2.5 Misalkan A = [ai,j ] ∈ Rn×n nonsingular matriks ultrametriks umum digeser. Jika A adalah matrik ultrametrik umum, maka A−1 ∈ Ln . (−1) Jika A adalah tipe Un−m,m , untuk m memenuhi 1 ≤ m ≤ n, maka A−1 ∈ Ln , −1 A dapat ditulis dalam bentuk A−1 = tI − B dengan B ≥ 0 dan ρm−1 (B) ≤ t ≤ ρm (B). ρm−1 (B) = t jika dan hanya jika ada principal submatriks order n − m + 1 dalam A yang singular. Teorema 2.6 Jika nonsingular Z-matriks A ∈ Rn×n adalah matriks ultrametrik umum yang digeser, maka A−1 ∈ Lm−1 jika dan hanya jika A adalah (−1) tipe Un−m,m untuk m dengan 1 ≤ m ≤ n. Teorema 2.7 Misalkan Z-matriks A ∈ Rn×n dan andaikan ada matriks segitiga bawah R dan matriks segitiga atas S, dengan elemen diagonal positip sedemikian sehingga A = RS. Maka 1. A mempunyai leading principal minors positif 2. R1 dan S −1 ada 3. Elemen di luar diagonal dari R dan S nonnegative 4. Tidak ada elemen R−1 dan S −1 negatif dan semua elemen diagonal R−1 dan S −1 adalah positif. Teorema 2.8 Misalkan A suatu Z-matriks sedemikangga akar real karakteristik dari A adalah positif. Misalkan B adalah suatu Z-matriks sedemikan sehingga A ≤ B. Maka


5

1. A−1 dan B −1 ada 2. 0 ≤ B −1 ≤ A−1 3. Setiap akar real karakteristik dari B adalah positif 4. det(B) ≥ det(A) > 0

3. KELAS Z-MATRIKS Definisi 3.5 Misalkan Ls (untuk s = 0, . . . , n) menyatakan kelas matriks yang memuat matriks n × n yang mempunyai bentuk

A = tI − B dengan B ≥ 0 dan ρs (B) ≤ t < ρs+1 (B).

(1)

˜ :B ˜ adalah sub matriks prinsipal dari B} dan Di sini ρs (B) := maks{ρ(B) kita set ρ0 (B) := −∞ dan ρn+1 (B) := ∞. Skalar t dan matriks B dalam (1) tidak tunggal, tetapi setiap Z-matriks masuk secara eksak ke dalam satu himpunan Ls . Teorema 3.9 Misalkan A ∈ Rn×n suatu Z-matriks. Maka, A ∈ Ls untuk satu s, s = 0, . . . , n − 1, jika dan hanya jika semua prinsipal submatriks berorder s adalah M -matriks, tetapi ada satu prinsipal submatriks berorder s + 1 yang bukan suatu M -matriks. Akibat 3.2 Jika A ∈ Ls dan D suatu matriks diagonal positif n × n, maka DA dan AD masuk ke Ls untuk s = 0, . . . , n. Teorema 3.10 Misalkan A ∈ Ls , B ∈ Lt dengan A ≤ B. Maka, i) s ≤ t, dan ii) apabila suatu matriks C dipenuhi A ≤ C ≤ B, maka C ∈ Lq dengan s ≤ q ≤ t. Dalam yang berikut diperhatikan tanda dari determinan dari A ∈ Ls . Telah diketahui bahwa det A ≥ 0 jika A ∈ Ln , yaitu A merupakan suatu M -matriks.


6

Teorema 3.11 Misalkan A ∈ Rn×n suatu matriks tak negatif dan misalkan λ suatu nilai eigen real dari A dengan λ 6= ρ(A). Maka λ ≤ ρbn/2c (A). Teorema 3.12 Misalkan A ∈ Ls . Maka (i) det(A) ≥ 0 jika s = n (ii) det(A) ≤ 0 jika bn/2c ≤ s < n Bukti. Mudah diketahui bahwa (i) berlaku (Berman, 1979). Misalkan A ∈ Ls dengan A = tI − B, B ≥ 0, dan suatu nilai tetap t ∈ R. Sekarang diperhatikan polinomial karakteristik dari B, χ(z) = det(zI − B). Nilai nol riil dari χ(z) mengindikasikan perubahan dari tanda dari determinan dari Z-matriks yang bersesuaian dengan matriks B. Tetapi nilai nol dari χ(z) adalah nilai eigen dari B. Sehingga, dengan Teorema 3.5, (ii) mengikut secara langsung. Determinan dari matriks yang berbeda di dalam Ls dengan s < bn/2c tidak perlu mempunyai tanda yang sama. Setiap Z-matriks A, kecuali M matriks, mempunyai sedikitnya satu nilai eigen yang negatif. Tetapi, jika A ∈ Ls dengan bn/2c ≤ s < n, maka teorema berikut menyatakan bahwa A mempunyai secara tepat satu nilai eigen yang negatif. Teorema 3.13 Misal A ∈ Ls dengan bn/2c ≤ s < n. Maka A mempunyai tepat satu nilai eigen negatif. Bukti. Anggap A = tI − B, dengan ρs (B) ≤ t < ρs+1 (B), B ≥ 0. Jika µ ∈ σ(B) dimana σ(B) menyatakan spektrum dari B, maka t − µ ∈ σ(A). Tetapi dengan Teorema 3.5, sembarang nilai eigen riil µ dari B dengan µ 6= ρ(B) memenuhi µ ≤ ρbn/2c (B). Sehingga, t − ρ(B) merupakan satusatunya nilai eigen negatif untuk A. Nilai eigen riil terkecil dari A dinyatakan dengan n(A). Seperti diterangkan di atas, n(A) < 0 jika A ∈ Ls untuk semua s = 0, . . . , n − 1. Fiedler dan Markham (1992) menyebut suatu matriks A dalam Ls murni jika setiap prinsipal minor berorder s + 1 negatif. Tambahan, suatu matriks dalam Ls murni kuat jika setiap prinsipal minor berorder lebih besar atau sama dengan s + 1 adalah negatif. Komplemen Schur dari Z-matriks. Jika A dipartisi sebagai


A=

A1,1 A1,2 A2,1 A2,2

7

(2)

dengan A11 dan A22 matriks bujur sangkar, dan jika A11 tak singular, maka komplemen matriks dari A11 dalam A adalah A/A11 := A22 − A21 A−1 11 A12

(3)

Perhatikan komplemen Schur dari matriks A ∈ Ls dimana prinsipal submatriks A11 berorder kurang atau sama dengan s. Maka komplemen Schur A/A11 selalu merupakan Z-matriks. Umumnya, komplemen Schur dari Z-matriks tak perlu merupakan Z-matriks. Yang berikut diberikan karakterisasi matriks dalam Ls menggunakan inversnya. Teorema 3.14 Misal A ∈ Rn×n Z-matriks tak singular. Maka A ∈ Ls jika dan hanya jika kasus alternatif berikut (a) atau (b) berlaku: (a)

(i) det(A) < 0, (ii) semua prinsipal minor dari A−1 berorder lebih besar atau sama de-ngan n − 1 adalah tak positif. (iii) ada suatu prinsipal minor dari A−1 berorder n−s−1 yang positif;

(b)

(i) det(A) > 0, (ii) semua prinsipal minor dari A−1 berorder lebih besar atau sama de-ngan n − 1 adalah tak negatif. (iii) ada suatu prinsipal minor dari A−1 berorder n−s−1 yang negatif;

Ide Teorema 3.8 telah digunakan dalam (Nabben & Varga, 1995) untuk mencari kelas khusus dari invers Ls -matriks untuk semua s. Suatu Z-matriks n × n dengan determinan positif merupakan Z-matriks dari salah satu kelas Ls dengan 0 ≤ s < bn/2c, karena matriks dari kelas Ls lain mempunyai determinan tak positif. Karakterisasi yang diberikan dalam Teorema 3.8 didasarkan pada tanda dari determinan dari Z-matriks. Suatu matriks tak singular C disebut invers Z-matriks jika C −1 merupakan Z-matriks. Persisnya, suatu matriks tak singular C ∈ Rn×n invers Ls -matriks jika C −1 ∈ Ls untuk suatu s dengan s ∈ {0, . . . , n}. Jadi, invers M -matriks adalah invers Ln -matriks, invers N0 -matriks merupakan invers Ln−1 -matriks, dan invers F0 -matriks merupakan invers Ln−2 -matriks.


8

4. KELAS INVERS Z-MATRIKS Dalam bagian ini Ls menyatakan kelas dari invers Ls -matriks untuk s = 0, . . . , n, sementara Z menyatakann seluruh kelas dari invers Z-matriks. Teorema 4.15 Misal C ∈ Ls . Maka det(C) ≥ 0 jika s = n; dan det(C) ≤ 0 jika bn/2c ≤ s < n. Teorema 4.16 Misal C ∈ Ls dengan bn/2c ≤ s < n. Maka C mempunyai secara eksak satu nilai eigen negatif. Teorema 4.17 Misal C ∈ Rn×n suatu invers Z-matriks. Maka C ∈ Ls jika dan hanya jika satu dari kasus alternatif (a) atau (b) berlaku: (a)

(i) det (C) < 0, (ii) semua prinsipal minor dari C berorder lebih besar atau sama dengan n − s adalah tak positif. (iii) ada suatu prinsipal minor dari C berorder n − s − 1 yang positif;

(b)

(i) det (C) > 0, (ii) semua prinsipal minor dari C berorder lebih besar atau sama dengan n − s adalah tak negatif. (iii) ada suatu prinsipal minor dari C berorder n − s − 1 yang negatif;

Karena tanda dari determinan suatu matriks adalah invarian di bawah pergandaan dengan matriks diagonal positif, diperoleh invarian dari kelas CLs di bawah operasi ini. Akibat 4.3 Jika C ∈ Ls , s ∈ {0, . . . , n} dan jika D adalah matriks diagonal positif n × n, maka DC dan CD masuk ke Ls . Kelas invers M -matriks adalah invarian di bawah penjumlahan dari matriks diagonal positif. Tetapi, ini tidak benar untuk Z. Perhatikan sebagai contoh invers N0 -matriks. Invers N0 -matriks merupakan matriks tak positif, sehingga merupakan Z-matriks. Jika ditambahkan sI, s ∈ R cukup besar, invers N0 -matriks menjadi M -matriks, yang bukan merupakan invers Zmatriks. Setiap prinsipal submatriks dari suatu Z-matriks merupakan Z-matriks lagi. Telah dibuktikan oleh (Markham, 1972), setiap prinsipal submatriks

9


dari suatu M -matriks merupakan M -matriks lagi. Secara umum hal ini tidak benar untuk semua invers Z-matriks. Tetapi jika A berada dalam Ls dan s ≥ bn/2c, maka ada prinsipal submatriks yang merupakan invers Z-matriks lagi. Teorema 4.18 Misalkan matriks n×n C berada dalam Ls dengan bn/2c ≤ s < n, yaitu C −1 ∈ Ls,n . Misalkan C˜ merupakan prinsipal submatriks tak singular dari C yang berorder k, dengan k ≥ n − s. Maka C˜ merupakan invers Z-matriks dalam Lj−n+k,k dengan n − k ≤ j ≤ s. Jika C˜ memuat paling sedikit satu prinsipal minor positif berorder n − s − 1, maka matriks n×n C˜ memenuhi C˜ ∈ Lj−n+k,k , dimana didefinisikan prinsipal minor 0×0 adalah positif. Bukti. Tanpa mengurangi keumuman, diasumsikan bahwa C˜ ∈ Rn×n prinsipal submatriks unggul dari C, yang kita nyatakan dengan C11 . Maka dipartisi C dan C −1 sebagai C=

C1,1 C1,2 C2,1 C2,2

,

C

−1

:= A =

A1,1 A1,2 A2,1 A2,2

(4)

Karena C −1 ∈ Ls,n dan order dari A22 adalah kurang dari atau sama dengan s, kita punya det (A22 ) ≥ 0. Selain itu, kita punya det C < 0, karena bn/2c ≤ s < n. Tetapi karena detC11 = (detA22 )(detA−1 ), maka diperoleh A22 > 0 dan det C11 < 0. Bentuk eksplisit dari A−1 diberikan oleh

A

−1

=

−1 −1 −1 −1 −1 A−1 1,1 + A1,1 A1,2 (A/A1,1 ) A2,1 A1,1 −A1,1 A1,2 (A/A1,1 ) −1 −1 −1 (A/A1,1 ) −(A/A1,1 ) A2,1 A1,1

(5)

Sehingga C11 = (A/A22 )−1 . Tetapi A/A22 adalah suatu Z-matriks. Oleh karena itu C11 merupakan invers Z-matriks. Lebih lanjut, karena C11 merupakan prinsipal submatriks dari C, maka semua prinsipal minor dari C11 yang berorder lebih besar atau sama dengan n − s adalah tak positif. Teorema 4.3 menjamin eksistensi dari prinsipal minor positif berorder n−s− 1 dari C. Jika C11 memuat prinsipal submatriks, maka dikenakan Teorema −1 4.3 dan mengganti n dengan k. Jadi matriks k ×kC11 memenuhi C11 ∈ Lt,m dengan t = s − n + k dan m = k. Jika C11 tidak memuat minor yang demikian, maka C11 ∈ Lj−n+k,k untuk n − k ≤ j ≤ s.


10

Teorema 4.19 Misalkan C suatu invers Z-matriks yang dipartisi sebagai C=

C1,1 C1,2 C2,1 C2,2

dengan C11 prinsipal submatriks tak singular dari C yang berorder sembarang. Maka C/C11 juga merupakan Z-matriks. Bukti. Misalkan A := C −1 dipartisi secara konformal dengan C sedemikian sehingga A=

A1,1 A1,2 A2,1 A2,2

Dengan menggunakan bentuk eksplisit dari C −1 seperti dalam (3.2), diperoleh A22 = (C/C11 )−1 . Sehingga, A22 adalah tak singular. Tetapi, A22 merupakan prinsipal submatriks dari Z-matriks. Oleh karena itu, A22 sendiri adalah Z-matriks. Jadi, C/C11 = A−1 22 .

5. SYARAT CUKUP INVERS Z-MATRIKS Pada Bagian 2 telah dijelaskan bahwa suatu Z-matriks yang mempunyai prinsipal minor positif adalah M -matriks. Oleh karena itu dalam mengkaji invers Z-matriks yang mempunyai prinsipal minor positif, cukup diselidiki invers dari M -matriks yang tak singular. Sifat-sifat terinci dari M -matriks dibicarakan dalam (Berman & Plemmons, 1979). Nabben (1995) telah memecahkan masalah penentuan matriks yang mempunyai invers yang merupakan Z-matriks yang disebut invers Z-matriks. Didefinisikan suatu sistem kelas matriks yang untuknya invers dari setiap matriks dari setiap kelas masuk ke satu kelas dari klasifikasi Z-matriks yang didefinisikan oleh Fiedler dan Markham. Teorema 5.20 Misal A ∈ Rn×n suatu matriks ultrametrik terumum kuat. Maka A adalah tak singulir, dan inversnya A−1 ∈ Rn×n merupakan M matriks dominan secara diagonal baris kuat dan kolom kuat, dengan sifat bahwa ω(A)ξnT A−1 ξn < 1

(6)


11

Jika A suatu matriks ultrametrik terumum tak singulir, maka inversnya A−1 ∈ Rn×n merupakan M -matriks dominant secara diagonal baris dan kolom, dengan sifat bahwa ω(A)ξnT A−1 ξn ≤ 1

(7)

Untuk sembarang Z-matriks A, det A > 0 merupakan syarat perlu untuk A merupakan M -matriks tak singulir. Teorema 5.21 Misal A ∈ Rn×n suatu matriks ultrametrik terumum ter(−1) geser tak singulir berjenis Up,n−p . Jika 0 ≤ p ≤ n, maka det A < 0. Lebih lanjut, setiap prinsipal minor berorder m dengan p < m < n adalah tak positif, dan ada suatu prinsipal minor positif berorder p. Jika p = n atau jika A suatu matriks ultrametrik terumum tak singulir, maka det A > 0. Bukti. Jika A suatu matriks ultrametrik terumum tak singulir, maka A−1 merupakan M -matriks dominan secara diagonal baris dan kolom, sedemikian bahwa det A−1 > 0. Ketika 1 = det(AA−1 ) = detAdetA−1 , maka det A > 0. (−1) Secara sama, jika A berjenis Up,0 , maka det A > 0. Selanjutnya anggap bahwa A ∈ Rn×n suatu matriks ultrametrik teru(−1) mum tergeser tak singulir berjenis Up,n−p dengan 0 ≤ p < n. Jika n = 1 sedemikian bahwa p = 0, diakibatkan bahwa A = [α] ∈ R1×1 dengan α < 0, yang mengakibatkan det A < 0. Untuk 1 ≤ m < n, anggap hipotesis induktif bahwa sembarang B ∈ Rm×m yang merupakan matriks ultrametrik teru(−1) mum tergeser berjenis Uq,m−q untuk suatu q dengan 0 ≤ q < m memenuhi det B < 0. Untuk n > 1, dengan matriks permutasi, A dapat dinyatakan sebagai A=

A1,1 A1,2 A2,1 A2,2

(8)

dengan A11 ∈ Rr,r dan A22 ∈ Rn−r,n−r , 1 ≤ r < n, merupakan matriks T ultrametrik terumum tergeser. Kemudian A12 = τ (A)ξr ξn−r dan A21 = T ω(A)ξn−r ξr . Karena A dianggap tak singulir, maka juga A + I untuk semua > 0 yang cukup kecil. Sehingga, penggeseran membolehkan untuk menganggap bahwa kedua A11 dan A22 dalam (8) adalah tak singulir. Dengan menggunakan komplemen Schur dari A terhadap A11 diperoleh det A = det A11 det(A/A11 )


12

T = det A11 det[A22 − τ (A)ω(A)(ξrT A−1 11 ξr )ξn−r ξn−r ]. (−1)

Jika A11 berjenis Ur,0 , yaitu A11 M -matriks tak singulir, maka det A11 > 0. Anggap bahwa komplemen Schur T A/A11 = A22 − τ (A)ω(A)(ξrT A−1 11 ξr )ξn−r ξn−r

suatu M -matriks tak singulir. Jika A˜ sembarang prinsipal submatriks unggul dari A yang ordernya lebih besar daripada order A11 , maka komple˜ 11 merupakan prinsipal submatriks dari A/A11 . Sehingga, men Schur A/A ˜ 11 ) > 0, dan karena det A = det A11 det(A/A ˜ 11 ), maka det A > 0. det(A/A Kemudian, ketika semua prinsipal minor unggul dari A adalah positif, maka A merupakan M -matriks tak singulir, yang kontradiksi dengan asumsi pada A. Jadi, A/A11 bukan merupakan M -matriks tak singulir, tetapi lagi, ketika det A = det A11 det(A/A11 ) dengan det A 6= 0 dan dengan det A11 > 0, maka A/A11 tak pernah tak singulir. Sehingga dengan hipotesis induktif, diakibatkan det(A/A11 ) < 0, yang mengakibatkan det A < 0. (−1)

Jika A11 tidak berjenis Ur,0 , yaitu A11 bukan M -matriks tak singulir, maka dengan hipotesis induktif det A11 < 0. Lebih lanjut, A/A11 merupakan matriks ultrametrik terumum, sedemikian bahwa det(A/A11 ) > 0. Jadi det A < 0. Karena p adalah order dari prinsipal submatriks terbesar yang merupakan M -matriks tak singulir, maka pernyataan yang tersisa diakibatkan secara langsung. Untuk sembarang Z-matriks A, det A > 0 merupakan syarat perlu untuk A menjadi M -matriks tak singulir. Selanjutnya akan ditunjukkan bahwa dengan suatu struktur tertentu dari elemen-elemen di luar diagonal, det A > 0 juga cukup.

KESIMPULAN Dari hasil pembicaraan dalam bagian sebelumnya dapat ditarik kesimpulansekimpulan berikut. Klasifikasi dari Z-matriks secara langsung menentukan klasifikasi dari invers Z-matriks. Semua matriks yang masuk ke dalam kelas Z-matriks yang sama membentuk sub kelas dari kelas semua invers Zmatriks. Jika suatu Z-matriks tak singulir A ∈ Rn,n merupakan matriks ultrametrik terumum tergeser, maka A−1 ∈ Lm−1 jika dan hanya jika A (−1) berjenis Un−m,m untuk suatu m dengan 1 < m ≤ n.


13

DAFTAR PUSTAKA 1. Ben-Israel, A. Generalized Inverses of Matrices and their Application. 2. 3.

4. 5.

6.

7. 8.

9.

10. 11. 12.

13. 14. 15.

Springer-Verlag, New York, 1980. Berman, A. 1978. The Spectral Radius of a Nonnegative Matrix, Canad. Math. Bull. 21(1):113-114. Berman, A. and Jain, S. K.. 1988. Nonnegative Generalized Invers of Powers of Nonnegative Matrices, Linear Algebra and Its Appl. 107 : 169 - 179. Berman, A. & Plemmons, R. 1979. Nonnegative Matrices in the Mathematical Sciences. New York: Academic Press. Fallat, S. M. & Tsatsomeros, M. J. 2002. On The Cayley Transform of Positivity Classes of Matrices. The Electronic Journal of Linear Algebra. Volume 9, pp. 190-196 Fiedler, M. and Markham, T. L. 1990. Some Connections Between the Drazin Inverse, P-Matrices and the Closure of Invers M-Matrices, Linear Algebra and Its Appl. 132 : 163 - 172. Fiedler, M & Markham, T. L. 1992. A classification of matrices of class Z. Linear Algebra Appl. 173: 115-124. Fiedler, M. and Ptak, V. 1962. On Matrices with Non-positive OffDiagonal Elements and Positive Principal Minors, Czech. Math. J. 12 : 382 - 400. Fiedler, M. and Ptak, V. 1966. Some Result on Matrices of Class K and their Application to the Convergence Rate of Iteration Procedures, Czech. Math. J. 16 : 260 - 273. Hershkowitz, D. and Schneider, H. 1988. On the Generalized Nullspace of M-Matrices and Z-Matrices, Linear Algebra and Its Appl. 106 : 5 - 23. Johnson, G.A. 1985. Inverse N0 -matrices. Linear Algebra Appl. 64:215222. Johnson C. R. 1996. The Combinatorially Simmetric P-Matrix Completion Problem. The Electronic Journal of Linear Algebra. Volume 1, pp. 59-63 Li, B. 1997. Generalizations of Diagonal Dominance in Matrix Theory. Ph.D Thesis. University of Regina. Markham, T. L. 1972. Nonnegative matrices whose inverses are Mmatrices. Proc. Amer. Math. Soc. 36:326-330. McDonald, J. J., Olesky, D. D., Schneider, H., Tsatsomeros, M. J. & Van Den Driessche, P. 1998. Z-Pencils. The Electronic Journal of Linear Algebra. Volume 4, pp. 32-38


14

16. Meyer. Jr. C. D. and Plemmons, R. J. 1977. Convergent Powers of a

17. 18. 19. 20. 21. 22. 23. 24.

25. 26. 27.

Matrix with Application to Iterative Methods for Singular Linear System, SIAM J. Numerical Analysis. 14 : 699 - 705. Meyer. Jr. C. D. and Stadelmaier, M. W. 1978. Singular M-Matrices and Inverse Positivity, Linear Algebra and Its Appl. 22 : 139 - 156. Nabben, R. 1997. Z-matrices and inverse Z-matrices. Linear Algebra and Its Applications, 256: 31-48. Nabben, R. & Varga, R. S. 1995a. On Classes of Inverse Z-Matrices. Linear Algebra and Its Applications, 223/224: 521-552. Nabben, R. & Varga, R. S. 1995b. Generalized ultrametric matrices–a class of inverse M-matrices. Linear Algebra Appl. 220:365-390. Oldenburger, R. 1940. Infinite Powers of Matrices and Characteristic Root, Duke Math. J. 6:357-361. Plemmons, R. J. 1976. M-Matrices Leading to Semiconvergent Splittings, Linear Algebra and Its Appl. 15 : 243 - 252. Plemmons, R. J. 1977. M-Matrices Characterizations. I - Nonsingular M-Matrices, Linear Algebra and Its Appl. 18 : 175 - 188. Rothblum, U. G. 1980. Convergence Properties of Powers of Matrices with Applications to Iterative Methods for Solving Linear Algebra, Lecture Notes in Economics and Mathematical System, Vol. 174 : 231 - 247, Springer-Verlag, Berlin. Rudin, W. 1956. Principles of Mathematical Analysis, McGraw-Hill International Editions, 1976. Schroder, J. 1978. M-Matrices and Generalization Using an Operator Theory Approach, SIAM Review Vol. 20(2) : 213 - 244. Wielandt, H. 1950. Unzerlegbare, nicht Negative Matrizen, Mathematische Zeitchrift, Bd. 52.

Sianipar: Jurusan Matematika, FMIPA, Universitas Sumatera Utara, Jl. Biote-

knologi No. 1 Kampus USU, Medan 20155, Indonesia.

E-mail: [email protected]

Mardiningsih – Digrap Jarang

15


DIGRAPH JARANG DENGAN EKSPONEN BESAR

Mardiningsih dan Astri Syafrianti Abstract. A strong connected digraph D is primitive if there exist a positive integer number k such that for every couple of vertices u and v in D there exist a walk from u to v that its length is k. The smallest positive integer k is called the exponent D and notated as exp(D). A strong connected digraph on n vertices is called to have a big exponent if w2n + 2 ≤ exp(D) ≤ wn , which wn = (n − 1)2 + 1. The sparse digraph on n vertices is strong connected digraph with (n + 1) arcs. Given a digraph with a big exponent which have two cycle in different length that is k-cycle and j-cycle where k > j. Consider D is a digraph with a big exponent which only have one k-cycle. For k = n and k = n − 1, it will be shown that D is a sparse digraph if and only if exp(D) = n + j(k − 2). Given exp(D) = γn then presented an algorithm to find a sparse digraph on n vertex with a big exponents.

1. PENDAHULUAN Suatu digraph terhubung kuat D adalah primitif jika terdapat bilangan bulat positif k sehingga untuk setiap pasangan verteks u dan v di D terdapat walk dari u ke v yang panjangnya k. Bilangan bulat positif terkecil k disebut eksponen D, dinotasikan dengan γ(D). Diketahui bahwa digraph D dikatakan primitif jika dan hanya jika D terhubung kuat dan pembagi persekutuan terbesar dari panjang-panjang cycle di D adalah 1[1]. Kondisi ini memenuhi bahwa suatu digraph primitif D atas n verteks mempunyai Received 12-03-2008, Accepted 27-09-2008. 2000 Mathematics Subject Classification: 05C20, 05C40 Key words and Phrases: digraph primitif, eksponen besar.


16

sedikitnya dua cycle, maka digraph D mempunyai sedikitnya n + 1 arc. Digraph jarang atas n verteks diartikan sebagai digraph terhubung kuat yang terdiri dari n + 1 arc. Wielandt [7] memperlihatkan bahwa untuk suatu digraph primitif atas n ≥ 3 verteks γ(D) ≤ (n − 1)2 + 1. Batas atas diperoleh jika dan hanya jika D isomorfik dengan digraph Wielandt Wn , yaitu suatu digraph atas n verteks yang terdiri dari cycle Hamilton dan suatu arc sehingga terbentuk sebuah (n + 1)-cycle. Maka eksponen terbesar digraph primitif atas n verteks adalah (n − 1)2 + 1, nilai ini dinotasikan dengan wn . Suatu digraph primitif D atas n verteks dikatakan mempunyai eksponen besar jika γ(D) ≥ w2n + 2. Terdapat beberapa penelitian tentang digraph dengan eksponen besar. Kirkland et.al. [4] membahas syarat perlu dan cukup untuk digraph dengan eksponen besar yang berhubungan dengan panjang-panjang cycle di D. Mereka memperlihatkan bahwa terdapat suatu digraph primitif atas n verteks mempunyai hanya dua panjang cycle yaitu dengan panjang k dan j dan eksponen exp(D) ≥ w2n + 2 jika dan hanya j(k − 2) ≥ w2n + 2. MacGillivray et.al. [5] mendiskusikan digraph primitif dengan eksponen w n 2 + 2 yang terdiri dari n + 1 atau n + 2 arc. Penelitian ini memfokuskan pada pembahasan syarat perlu dan cukup untuk digraph dengan eksponen besar yang mempunyai exp(D) = n + j(k − 2). Penelitian ini memfokuskan pembahasan terutama untuk golongan digraph yang terdiri hanya satu kcycle dan (n + 1) arc. 2. PENGERTIAN DASAR Andaikan D adalah suatu digraph terhubung kuat atas n verteks terdiri dari cycle-cycle dengan panjang k dan j. Untuk n ≥ k > j ≥ 2, dengan k + j ≥ n + 2 dan gcd(k, j) =, W (n, k, j) dinotasikan digraph primitif atas n ≥ 3 verteks terdiri tepat satu k-cycle dan satu j-cycle. Andaikan k dan j bilangan bulat positif sedemikian hingga gcd(k, j) = 1. Frobenius-Schur indeks terhadap k dan j, dinotasikan oleh φ(k, j), adalah bilangan bulat positif terkecil N sedemikian sehingga untuk setiap bilangan bulat n0 ≥ N dapat ditulis sebagai ak + bj untuk beberapa bilangan bulat nonnegatif a dan b. Diketahui bahwa φ(k, j) = (k − 1)(j − 1), [1]. Andaikan D digraph primitif dan andaikan u dan v adalah dua verteks di D. Didefinisikan suatu bilangan ruv menjadi 0 jika u = v, dan menjadi panjang terpendek dari uv-walk yang melewati cycle untuk masing-masing panjang di D. Bilangan r = max{ruv } atas semua pasangan verteks u dan v di D adalah circumdiameter dari D.


17

Figure 1: Digraph W (n, k, j) Bagian ini dimulai dari beberapa hasil terdahulu tentang digraph dengan eksponen besar dan kemudian membahas suatu formula untuk eksponen dari digraph primitif yang terdiri dari dua cycle. Teorema 1 Andaikan D digraph primitif atas n ≥ 3 verteks dengan eksponen besar. Maka D mempunyai cycle tepat dengan dua panjang berbeda yaitu yang panjangnya j dan k dengan j < k ≤ n. Hasil berikut, perhatikan Teorema 1 dari [3], memberikan suatu batas bawah untuk panjang dari cycle terkecil. Teorema 2 Andaikan D digraph primitif atas n verteks dengan eksponen besar dan dua cycle dengan dua panjang berbeda yaitu j and k, j < k ≤ n. Maka j ≥ n−1 2 Hasil berikut, lihat [5, Proposisi 2.3], memperlihatkan bahwa suatu digraph jarang dengan eksponen besar harus terdiri dari tepat dua cycle yang prima relatif. Proposisi 1 Andaikan D digraph primitif atas n ≥ 4 verteks yang mempunyai eksponen besar dan n+1 arc. Maka Dterdiri dari dua panjang cycle n−1 yang prima relatif yaitu j dan k dengan 2 j < k ≤ n dan k + j ≥ n + 2. Kirkland et.al [4] telah mengkarakterisasikan semua digraph dengan eksponen besar terdiri dari n-cycle dan (n + 1)-cycles. Teorema 2.1 dan


18

2.2 di atas menyajikan hasilnya untuk digraph yang terdiri dari n-cycle dan (n + 1)-cycles. Teorema 3 Andaikan D digraph primitif atas n ≥ 3 verteks terdiri dari satu n-cycle. (a) Anggap bahwa j > n/2. Digraph D mempunyai eksponen besar dan cycle yang panjangnya n dan j jika dan hanya jika D isomorfik ke suatu subdigraph dari digraph yang terdiri dari cycle 1 → n → n−1 → . . . → 2 → 1, dan arc i → I + j − 1 untuk 1 ≤ I ≤ n − j + 1. (b) Anggap bahwa n ganjil dan j = (n + 1)/2. Digraph D mempunyai eksponen besar dan cycle yang panjangnya n dan j jika dan hanya jika D isomorfik ke suatu subdigraph dari digraph yang terdiri dari cycle 1 → n → n − 1 → ... → 2 → 1, dan arc i → i + j − 1 untuk 1 ≤ I ≤ j. Teorema 4 Andaikan D digraph primitif atas n ≥ 6 verteks terdiri dari (n + 1)-cycle dan j-cycle untuk (n + 1) > j. (a) Anggap bahwa j > n/2. Digraph D mempunyai eksponen besar dan cycle yang panjangnya (n + 1) dan j jika dan hanya jika D isomorfik ke suatu subdigraph dari digraph yang terdiri dari (n + 1)-cycle 1 → n − 1 → n − 2 → ... → 2 → 1, arc i → i + j − 1 untuk 1 ≤ I ≤ n − j dan salah satu dari yang berikut: (i) arc 2 → n → n − 1, atau (i + 1) → n → i − 1 untuk 2 ≤ I ≤ n − 2, (ii) arc 1 → n , n → n − 2, dan n → j − 1. (b) Anggap bahwa n genap dan j = n/2. Digraph D mempunyai eksponen besar dan cycle yang panjangnya (n + 1) dan j jika dan hanya jika D isomorfik ke suatu subdigraph dari digraph yyang terdiri (n + 1)-cycle 1 → n−1 → n−2 → ... → 2 → 1, arc i → i+j−1 untuk 1 ≤ i ≤ n/2−3 dan salah satu dari yang berikut: (i) arc 2 → n → n − 1, atau (i + 1) → n → i − 1 untuk 2 ≤ I ≤ n − 2, (ii) arc 1 → n, n → n − 2, dan n → j − 1 Lemma 1 Andaikan D digraph primitif atas n verteks terdiri dari dua cycle yang panjangnya k dan j, dengan j < k ≤ n. Maka exp(D) = n+j(k+2).


19

BUKTI. Pertama akan diperlihatkan bahwa exp(D) ≤ n + j(k + 2), dengan memperlihatkan bahwa untuk setiap pasangan verteks u dan v di D terdapat suatu walk dari u ke v yang panjangnya n + j(k + 2). Andaikan puv adalah suatu path dari u ke v dan p adalah path terpanjang yang terletak pada kedua cycle di D dan andaikan v0 adalah verteks inisial dari p. (Suatu verteks v dapat dianggap sebagai suatu path trivial yang panjangnya 0.) Andaikan `(puv ) dan `, berturut-turut, panjang dari path puv dan p. Perhatikan bahwa `(puv ) ≤ k + j − 2 − ` dan ` = k + j − n − 1. Definisikan e = kj +k −j −1−` = n+j(k −2). Karena `(puv ) ≤ k +j −2−`, e−`(puv ) ≥ kj −(k +j)+1. Tetapi k dan j adalah prima relatif, maka terdapat bilangan bulat nonnegatif a and b sedemikian hingga ak + bj = e − `(puv ), lihat [1, halaman 83-85]. Maka e = ak + bj + `(puv ), artinya walk wuv dimulai dari u, kemudian bergerak a kali mengelilingi k-cycle dan b kali mengelilingi j-cycle dan kembali ke u, dan akhirnya bergerak menuju verteks v adalah suatu walk dari u ke v yang panjangnya e = n + j(k − 2). Untuk setiap verteks u di D, panjang walk tertutup dari u ke dirinya sendiri dapat dinyatakan dalam bentuk ak + bj untuk beberapa bilangan bulat nonnegatif a dan b. Karena φ(k, j) = (k − 1)(j − 1) adalah bilangan bulat paling kecil sedemikian hingga setiap bilangan bulat q ≥ φ(k, j) dapat dinyatakan dalam bentuk ak + bj, maka exp(D) ≥ φ(k, j). Sekarang anggap pasangan verteks ` + 2 dan n. Walk tertutup terpendek dari ` + 2 ke k yang panjang paling kecilnya φ(k, j) mempunyai panjang 2k − `2 = n + k − j − 1. Maka exp(D) ≥ φ(k, j) + n + k − j − 1 = n + j(k − 2). 3. HASIL UTAMA Pada bagian ini membahas tentang digraph D atas n verteks dengan eksponen besar dan mempunyai dua cycle yang panjangnya k dan j, k > j, untuk D yang mempunyai tepat satu k-cycle. Ditentukan syarat perlu dan cukup untuk digraph dengan eksponen besar yang mempunyai eksponen n = n + j(k + 2). Kasus pertama yang akan dibahas ketika k = n dan k = n − 1. Lemma 2 Andaikan D digraph primitif atas n verteks dengan eksponen besar dan terdiri dari cycle-cycle yang panjangnya n dan j. Maka D mem punyai n + 1 arc jika dan hanya jika j(n + 2) = γn − n dengan j ≥ n−1 2 BUKTI. Karena D mempunyai n + 1 arc, Proposisi 1 menjamin bahwa digraph D terdiri dari tepat dua cycle yang panjangnya n dan j. Oleh


20

Lemma 1, γn = n + j(n + 2). Maka j(n + 2) = γn − n. Karena D adalah suatu digraph dengan eksponen besar, maka diperlukan n−1 ≤ j < n. w2 n Sebaliknya asumsikan bahwa n = n + j(n + 2) ≥ 2 + 2. Anggap dua kasus yang mana j ≥ n/2 dan j = (n + 1)/2 dan n ganjil. Jika j ≥ n/2, maka digraph D isomorfik ke suatu subdigraph primitif pada digraph primitif yang terdiri dari cycle 1 → n → n − 1 → ... → 2 → 1 dan arc a → a + j − 1 for 1 ≤ a ≤ n − j + 1, lihat [4, Teorema 2.2]. Klaim bahwa D terdiri dari tepat satu cycle dengan panjang j. Sebaliknya asumsikan bahwa D mempunyai lebih dari satu j-cycle, dan tanpa kehilangan sifat keumuman asumsikan cycle 1 → j → j − 1 → ... → 2 → 1 adalah salah satu j-cycles di D. Maka circumdiameter D adalah 2n + j − t, yang mana 2 ≤ t ≤ n − j + 1. Maka exp(D) = χ(n, j) + r ≤ φ(n, j) + 2n − j − t untuk beberapa bilangan bulat 2 ≤ t ≤ n − j + 1. Artinya exp(D) < n + j(n + 2), kontradiksi dengan fakta bahwa exp(D) = n+j(n+2) = χ(n, j)+2n+j −1. Maka D harus terdiri dari tepat satu j-cycle. Jika n ganjil dan j = (n + 1)/2, maka digraph D isomorfik dengan subdigraph dari digraph yang terdiri dari cycle 1 → n → n−1 → ... → 2 → 1 dan arc a → a + j − 1 untuk 1 ≤ a ≤ n − j + 1, lihat [4, Teorema 2.3]. Jika D mempunyai lebih dari satu j-cycle, argumen yang sama untuk kasus sebelumnya memperlihatkan bahwa exp(D) ≤ χ(n, j) + 2n − j − t untuk beberapa bilangan bulat 2 ≤ t ≤ j. Kontradiksi dengan fakta exp(D) = n + j(n + 2) = χ(n, j) + 2n − j − 1. Maka D harus terdiri dari satu j-cycle. Kedua kasus memenuhi bahwa D terdiri dari satu n-cycle dan satu j-cycle. Oleh karena itu dapat disimpulkan bahwa D has n + 1 arc. Anggap kasus berikut ketika digraph D terdiri dari tepat satu k-cycle yang mana k = n − 1. Lemma 3 Andaikan D digraph primitif atas n verteks dengan eksponen besar. Anggap D terdiri dari hanya satu (n + 1)-cycle dan beberapa j-cycles. Digraph D mempunyai n + 1 arc jika dan hanya jika n = n + j(n + 3) dengan j ≥ n−1 . 2 BUKTI. Jika D mempunyai n + 1 arc, maka D terdiri dari satu j-cycle. Maka oleh Lemma 2 eksponen D adalah n = n + j(k + 2) = n + j(n + 3), pada kasus ini k = n − 1. w Sekarang anggap eksponen D adalah n = n + j(k + 2) = n + j(n + 3) ≥ n 2 + 2. Akan diperlihatkan D mempunyai n + 1 arc dengan menganggap kasus ketika j > n/2 dan ketika j = n/2 dan n genap. Karena D mempunyai hanya satu (n+1)-cycle, jika j > n/2 maka D isomorfik ke suatu subdigraph


21

dari digraph yang terdiri dari (n−1)-cycle 1 → n−1 → n−2 → ... → 2 → 1, path 1 → n → j − 1, dan arc i → I + j − 1, 1 ≤ I ≤ n − j, lihat [4, Teorema 3.3]. Jika D mempunyai lebih dari satu j-cycle, maka circumdiameter D adalah r = 2n − j − 1 − t untuk beberapa bilangan bulat 1 ≤ t ≤ j − 1. Maka exp(D) ≤ χ(n − 1, j) + 2n − j − 2 < n + j(n − 3). Suatu kontradiksi, sehingga D harus mempunyai hanya satu j-cycle. Karena D mempunyai hanya satu (n + 1)-cycle, jika j = n/2 maka D isomorfik ke suatu subdigraph dari digraph yang terdiri dari (n + 1)-cycle 1 → n − 1 → n − 2 → ... → 2 → 1, path 1 → n → j − 1, dan arc i → i + j − 1, 1 ≤ i ≤ j − 3, lihat[4, Teorema 3.4]. Jika D mempunyai lebih dari satu jcycle, maka circumdiameter D adalah 2n+j +1−t untuk beberapa bilangan integer 1 ≤ t ≤ j −3. Artinya exp(D) ≤ χ(n−1, j)+2n−j −2 < n+j(n−3). Suatu kontradiksi, maka D harus mempunyai hanya satu j-cycle. Sekarang dapat disimpulkan bahwa D mempunyai satu cycle dengan panjang j. Karena D adalah terhubung kuat dan terdiri dari hanya dua cycle, maka D mempunyai n + 1 arc. Perhatikan bahwa Lemma 3 tidak benar ketika D mempunyai lebih dari satu (n + 1)- cycle. Berikut diberikan contoh penyangkalan, digraph D atas n = 8 verteks yang mempunyai dua 7-cycles dan satu 5-cycle.

Figure 2: Digraph atas 8 verteks dengan dua 7-cycles dan satu 5-cycle Eksponen dari D adalah exp(D) = 33. Perhatikan bahwa exp(D) = n + j(k + 2), tetapi D mempunyai n + 2 arc. Pada dasarnya Lemma 2 menyatakan jika terdapat bilangan bulat positif yang prima relatif yaitu k and j dengan j(k + 2) = γn − n, maka dapat dibangun suatu digraph primitif terhubung kuat atas n = k + j + z verteks yang terdiri dari suatu k-cycle dan suatu j-cycle sedemikian hingga eksponen D adalah exp(D) = n+j(k +2), yang mana z adalah banyaknya verteks diantara k-cycle dan j-cycle. Diberikan n dan γn , algoritma berikut mencari


22

bilangan bulat prima relatif k dan j sedemikian hingga n = n + j(k + 2). Maka, algoritma menentukan suatu digraph atas n verteks dengan n + 1 arc yang terdiri dari satu k-cycle dan satu j-cycle untuk exp(D) = n + j(k + 2).

4. ALGORITMA ALGORITHM 3.1 1 t ← γn − n 2 for j = n−1 to n − 1 2 3 if j membagi t 4 d ← t/j 5 c ← n + jd 6 if c = n 7 k ←d+2 8 if j < k dan gcd(k, j) = 1 9 z ←k+j−n 10 print (n, k, j) Hasil dari eksekusi Algorithm 3.1 memperlihatkan bahwa untuk beberapa n terdapat lebih dari satu digraph primitif atas n vertices yang tidak isomorfik dengan n + 1 arc dan mempunyai eksponen γn . Untuk n = (3n2 − n + 2)/4 ditemukan dua digraph atas 13 verteks dan 14 arc mempunyai eksponen sama dengan γn , digraph tersebut adalah W(13, 13, 10) dan W(13, 12, 11). Meskipun, terdapat tiga digraph atas n = 1577 verteks yang tidak isomorfik dengan n + 1 arc dan mempunyai eksponen sama dengan γn , digraph tersebut adalah W(1577, 1183), W(1577, w 1577, n 1523, 1225) dan W(1577, 1367, 1365). Untuk n = 2 + 2, bilangan bulat ganjil terkecil n yang terdapat lebih dari satu digraph yang tidak isomorfik dengan order n mempunyai n+1 arc dan eksponen γn adalah n = 933. Hasil ini sesuai dengan hasil yang disajikan oleh MacGillivray, et all [5].

DAFTAR PUSTAKA 1. R.A. Brualdi, H.J Ryser, Combinatorial Matrix Theory, Cambridge Uni-

versity Press, 1991. 2. A.L. Dulmage, N.S. Mendelsohn, Gaps in the exponent set of primitive

matrices, Illinois J. Math. 8 (1964) 642-656.


23

3. S. Kirkland, A note on the eigenvalues of a primitive matrix with large 4. 5.

6. 7.

exponent, Linear Algebra Appl. 253 (1997) 103-112. S. Kirkland, D.D. Olesky, P. van den Driessche, Digraphs with large exponent, Electro. J. Linear Algebra 7 (2000) 30-40. G. MacGillivray, S. Nasserasr, D.D. Olesky, P. van den Driessche, Primitive digraph with smallest large exponent, Linear Algebra Appl. 428 (2008) 1740-1752. M. Lewin, Y. Vitek, A system gaps in the exponent set of primitive matrices, Illinois J. Math 25 (1981) 87-98. H. Wielandt, Unzerlegbare nicht negative Matrizen, Math. Zeit. 52 (1950) 642-645.

Mardiningsih: Jurusan Matematika, FMIPA, Universitas Sumatera Utara, Jl.

Bioteknologi No. 1 Kampus USU, Medan 20155, Indonesia.


Syafrianti: Jurusan Matematika, FMIPA, Universitas Sumatera Utara, Jl. Biote-




24

Rosmaini – Menentukan titik awal

25


MENENTUKAN TITIK AWAL PENYELESAIAN QUADRATIC ASISGNMENT PROBLEM

Elly Rosmaini Abstract Quadratic Assignment Problem (QAP) is a combinatorial problem in determine the placement of facilities in certain locations such that the minimize objective function is expressed by the distance between the location and flow between facility. First of all do a simple heuristic approach to obtain starting point assignments. Tabu Search Heuristic methods used to obtain starting point for other assignments that will be used to approach the starting point Nonlinear Programming. Further problems in Nonlinear Programming process by using other heuristic methods for solving deserve a round of the QAP. Convex nature of the problem is not caused need a good starting point in order to obtain the optimal solution globally. This procedure is applied to the problems with dimensions Backboard Wiring 20 x 20.

1. PENDAHULUAN Telah lama dikenal dalam ekonomi bahwa adanya bahan baku yang tak terbagi (indivisible) dapat menimbulkan problema yang serius bila alokasi yang efisien dari bahan baku tadi hendak dicapai. Ketakterbagian (indivisibility), merupakan suatu kasus khusus dari problema alokasi bahan baku di mana terdapatnya non-konveksitas dari problema. Received 12-03-2008, Accepted 27-09-2008. 2000 Mathematics Subject Classification: 90C55, 90C59 Key words and Phrases: Quadratic Assignment Problem, Tabu Search.


26

Salah satu problema yang berkaitan dengan situasi ketakterbagian (bulat) dari variabel keputusan dan adanya non-konveksitas dari fungsi objektif adalah problema Quadratic Assignment Problem (QAP) yang dikembangkan oleh Lawler (1963). Model problema ini memegang peranan penting dalam banyak pemakaian, misalnya teori lokasi, penjadwalan, backboard wiring, control panel, lay out, tata kota, dan lain-lain. Secara matematika model problema dapat ditulis sebagai berikut : Pandang problema melokasikan n fasilitas dalam n lokasi yang diberikan, di mana setiap lokasi dapat ditugaskan hanya terhadap satu fasilitas, dan tiap fasilitas dapat ditugaskan hanya pada satu lokasi. Jika fIk merupakan alur antara fasilitas i dan fasilitas k dan Cjq adalah ongkos transport per unit antara lokasi j dan q maka problema menjadi : n

n

X 1 XXX Minimumkanφ = = 1n = 1n fik Cjq Xij Xkq 2 q

(9)

i=1 k=1 j

dengan kendala n X

Xij = 1, i = 1, ..., n

(10)

Xij = 1, j = 1, ..., n

(11)

j=1 n X i=1

0 6 Xij 6 1, Xij bilangan cacah QAP telah terbukti jauh lebih sulit untuk dapat diselesaikan dibanding dengan problema yang linier. Problema taklinier assignment dapat diatasi dengan mengendorkan (relaxing) program bilangan cacahnya menjadi program linier dan kemudian diselesaikan. Hasilnya selalu bilangan cacah. QAP tidak hanya non linier, tetapi juga tidak uni modal . QAP dikenal sebagai problema sulit Non Polinomial (NP-Hard Problem, istilah komputasi ). Hingga kini belum ada suatu algoritma yang efisien yang dapat menyelesaikan QAP dengan dimensi n > 15 secara optimal. Umumnya dalam pemakaian dimensi yang lebih besar akan muncul karena perangai fungsi tujuan adalah nonkonveks mengakibatkan diperlukannya suatu metode heuristik yang baik untuk memperoleh titik awal penyelesaian yang tepat sehingga penyelesaian optimum global dapat diperoleh. Perkembangan baru-baru ini dalam program non linier berskala besar membuka jalan untuk memperlakukan QAP sebagai suatu program kwadrat berskala besar dari pada metode kombinatorial.


27

Metode Heuristik akan digabungkan dengan program non linier tersebut untuk dapat menyelesaikan QAP dengan dimensi n > 15 Metode penyelesaian eksak yang biasa dipakai untuk menylesaikan problema program bilangan cacah linier yaitu metode Branch and Bound dan metode Cutting Plane telah pernah diajukan peneliti (lihat Bazzara dan Elshafei [1979]; Gendron and Craine [1994] ; Anstreicher dan Brixius [1999]; Brixius dan Anstreicher [2000]; Anstreicher et al [2002] ), namun karena QAP adalah program sulit non polinomial maka tidaklah mengherankan bahwa metode yang diajukan tersebut hanya berfaedah untuk problema dengan dimensi kecil. Juga beberapa metode Heuristik telah pernah disarankan, misalnya algoritma Genetika, Tabu Search, Ant Colonies telah pernah diajukan oleh para peneliti untuk menyelesaikan problema ini (lihat Bruijs [1984]; Mawengkang dan Murtagh [1985]; Bunhard et al [1997]; Gambardella et al [1999]; Taillard [1991]; Thoreman dan Bolte [1999] dan Hahn [2000].Namun sebegitu jauh belum dapat memberikan penyelesaian yang optimal untuk problema QAP berskala besar terutama untuk n > 15. Metode Heuristik dipakai untuk menghasilkan suatu prosedur sederhana penentuan titik awal . Titik awal dikembangkan dengan metode Tabu search. Tabu Search merupakan improvement dari improving search yang digunakan untuk menghindari infinite loop dan short-term cycling akibat pemilihan improving moves. Bentuk kwadrat dari fungsi objektif lebih cenrung untuk tidak konveks sehingga penyelesaian optimal kontiniu yang diperoleh melalui program non linier, tidak merupakan optimal global. Jadi salah satu usaha untuk memperoleh penyelesaian yang baik banyak dipengaruhi titik awal. Kemudian prosedur heuristik lainnya akan dipakai untuk menetukan penyelesaian bulat fisibel dari penyelesaian yang diperoleh dari program non linier. Sebagai studi kasus penerapan metode yang diajukan akan diseleaikan jaringan instalasi kabel komunikasi (backboard wiring) dengan jumlah variabel 30 × 30, dengan menggunakan model QAP, maka panjang total kawat akan diminimumkan. Jarak dari tiap tiap lokasi yang akan dihubungkan dinyatakan dalam satu matriks D dan matriks ini merupakan mtriks simetri dimana bilangan pada diagonal utamanya adalah nol. Jarak antara lokasi i dan lokasi j dinyatakan dengan dij . QAP dikenal sebagai problem sulit non-polinomial (NP-Hard problem, istilah komputasi). Hingga kini belum ada suatu algoritma yang efisien yang dapat menyelesaikan QAP dengan dimensi n > 15 secara optimal. Umumnya dalam pemakaian, dimensi yang lebih besar akan muncul. Karena itu diperlukan metode Heuristik yang baik. Perkembangan baru-baru ini


28

dalam program non linier berskala besar telah membuka jalan untuk memperlakukan QAP sebagai suatu program kuadrat berskala besar dari pada memakai metode kombinatorial.Metode Heuristik akan digabung dengan Program non linier tersebut untuk dapat menyelesaikan QAP dengan dimensi n > 15. Metode penyelesaian eksak yang biasa dipakai untuk menyelesaikan problem penggunaan bilangan cacah linier yaitu Branch and Bound (B&B) dan metode Cutting Plane telah pernah diajukan oleh para peneliti lihat Bazarra dan Elshafei [1979]; Gendron dan Craininc [1994]; Anstreicher dan Baixius [1999]; Baixius dan Anstreicher [2000]; Anstreicher et al [2002]), namun karena QAP adalah problem sulit non polinomial maka tidaklah mengherankan bahwa metode-metode yang diajukan tersebut hanya berfaedah untuk problema dengan dimensi kecil. Suatu faktor yang penting dalam kinerja algoritma B & B adalah metode pemilihan batas bawah. Beberapa metode Heuristik, misalnya algoritma Genetika, Tabu Search, Ant Colonies telah pernah diajukan oleh para peneliti untuk menyelesaikan problema ini lihat, Bruijs [1984]; Mawengkang dan Murtagh [1985]; Burkard et al. [1997]; Gambardella et al. [1999]; Taillard [1991]; Thoremann dan Bolte [1999]; dan Hahn [2000]. Namun sebegitu jauh belum dapat memberikan penyelesaian yang optimal untuk problem QAP berskala besar. Tujuan penelitian ini menghasilkan suatu penyelesaian Quadratic Assignment Problem dengan menggabungkan metode heuristik dan program non linier sehingga ditemukan penugasan yang meminimumkan biaya dari suatu penugasan. Masalah QAP banyak memegang peranan penting dalam berbagai pemakaian, sehingga dengan adanya suatu algoritma yang efisien maka model masalah tersebut dapat diterapkan pada permasalahan yang muncul dalam dunia nyata, asalkan berkaitan dengan model QAP. 2. TABU SEARCH DAN HEURISTIK TITIK AWAL Penelitian ini merupakan penelitian kepustakaan untuk keperluan eksperimen. Langkah awal untuk menyelesaikan masalah QAP adalah dengan menggunakan metode heuristik untuk mendapatkan titik awal penyelesaian yang fisibel. Filosopi yang umum untuk membuat suatu penugasan awal fisibel terhadap lokasi adalah menempatkan fasilitas tersibuk dalam lokasi paling utama. Suatu prosedur yang menggambarkan prisip tersebut dapat diajukan sbb: Strategi pegembangan titik titik awal yang disajikan adalah menggunakan metode heuristik Tabu Search yaitu :


1. Hitung F1 =

n P

29

fik , untuk fasilitas i, i = 1, 2, . . . , n.

k=1

2. Membuat daftar fasilitas yang diurutkan secara menurun (descending) dari F1 . 3. Hitung Cj =

n P

Cj1 untuk setiap lokasi j, j = 1, 2, . . . , n.

j=1

4. Membuat daftar lokasi yang berturutan secara menaik (ascending) dari Cj . Prosedur sebagai berikut : 1. Inisialisasi. Pilih penyeleaian titik awal x(0) dan satu limit iterasi tmax . Set solusi Sementara x ← x(0) dan indeks t ← 0 Tidak terdapat tabu moves. 2. Kriteria stop Jika tidak terdapat non-tabu move ∆x dimove set µ yang menghasilkan feasible neighbor dari solusi x(t) , atau jika t = tmax , maka stop. Solusi sementara x adalah optimum lokal. 3. Move Pilih non-tabu feasible move ∆x ∈ µ sebagai ∆x(t+1) 4. Step Perbaharui x(t+1) ← x(t) + ∆x(t+1) 5. Solusi sementara Jika nilai fungsi objektif dari x(t+1) lebih baik dari x, gantikan x ← x(t+1) 6. Tabu list Hapus tabu move yang berada cukup lama di tabu list dan tambah tabu moves yang mengakibatkan cycling dari x(t+1) ke x(t) . 7. Increment Set t ← t + 1 dan kembali ke langkah 1. Setelah titik awal diperoleh kemudian diselesaikan relaksasi dari QAP dengan titik awal tadi menjadi titik awal untuk menyeleaikan nonlinier programing. Metode heuristik lainnya akan dipergunakan untuk memperoleh penyelesaian layak cacah dari problem awal QAP Masalah kombinatorial yang sering terjadi pada aplikasi seperti lokasi fasilitas, rencana layout dan backboard wiring. Dengan memperhatikan problema dari melokasikan n fasilitas pada n lokasi. Misalkan fik merupakan

30


alur (flow) antara fasilitas i dan fasilitas k dan cjq adalah biaya (jarak) pemindahan per unit antara lokasi j dan q sehingga masalahnya menjadi : Minimumkan φ =

n X n X n X n X

fik cjq xij xkq

(12)

i=1 k=1 j=1 q=1

Dengan kendala

n X

xij = 1i = 1, ..., n

(13)

xij = 1j = 1, ..., n

(14)

j=1 n X i=1

0 6 xij 6 1, xij : bilangan cacah Dengan batasan bahwa setiap lokasi hanya dapat ditempati oleh satu fasilitas saja dan setiap fasilitas hanya dapat menempati pada satu lokasi saja. Penentuan Titik Awal Penyelesaian Bentuk kuadrat dari fungsi objektif (4.1) lebih cenderung untuk tidak konveks sehingga penyelesaian optimal kontinu yang diperoleh melalui metode program non linier tidak merupakan optimal global. Jadi salah satu usaha untuk memperoleh penyelesaian yang baik banyak dipengaruhi oleh titik awal. Penentuan titik awal QAP dilakukan secara heuristik, secara umum untuk memperoleh sebuah penugasan awal fasilitas pada lokasi adalah dengan menempatkan fasilitas tersibuk pada lokasi yang utama, setelah tabel matriks aliran fasilitas (f ) dan matrik jarak (c) dibentuk. Maka titik awal penyelesaian berdasarkan algoritma berikut Algoritma 1 1. Hitung Fi =

n P

fik untuk setiap fasilitas i, i = 1, . . . , n

k=1

2. Membentuk daftar fasilitas yang diurutkan secara menurun dari Fi 3. Hitung Cj =

n P

untuk setiap lokasi j, j = 1, . . . , n

t=1

4. Membentuk daftar lokasi yang diurutkan secara menaik dari Cj


31

5. Tugaskan fasilitas pertama dari urutan fasilitas pada lokasi pertama dari daftar lokasi, kemudian fasilitas berikutnya sampai semua fasilitas ditugaskan. Algoritma ini cukup sederhana dan mudah untuk diimplementasikan dalam mendapatkan penugasan awal. Titik awal penugasan yang dihasilkan berdasarkan algoritma 1 adalah bukan merupakan hasil penyelesaian yang mendekati optimal atau mungkin juga bukan merupakan hasil penyelesaian yang baik oleh karena itu penyelesaian titik awal masih perlu dikembangkan hingga mendapat suatu hasil penyelesaian yang layak. Strategi pengembangan titik awal yang disajikan berikut menggunakan metode heuristik tabu search. Tabu search memerlukan titik awal, titik awal diperoleh dari algoritma 1. Tabu search harus mendefinisikan move set yaitu µ = { single comlement move } dan mentabukan single complement move yang sama. Prosedur Tabu Search sebagai berikut : Algoritma 2 1. Inisialisasi Titik awal penyeleaian diperoleh dari algoritma 1 sebut x(0) dan satu limit iterasi tmax . Set solusi sementara x ← x(0) dan indeks t ← 0. Tidak terdapat tabu moves. 2. Kriteria stop Jika tidak terdapat non-tabu move ∆x dimove set µ yang menghasilkan feasible neighbor dari solusi x(t) , atau jika t = tmax , maka stop. Solusi sementara x adalah optimum lokal. 3. Move Pilih non-tabu feasible move ∆x ∈ µ sebagai ∆x( t + 1) 4. Step Perbaharui. x(t+1) ← x(t) + ∆x(t+1) 5. Solusi sementara Jika nilai fungsi objektif dari x(t+1) lebih baik dari x, gantikan x ← x(t+1) 6. Tabu list Hapus tabu move yang berada cukup lama di tabu list dan tambahkan tabu moves yang mengakibatkan cycling dari x(t+1) ke x(t) 7. Increment. Set t ← t + 1 dan kembali ke langkah 1.


32

Setelah titik awal diperoleh dari tabu search kemudian diselesaikan relaksasi dari QAP dengan titik awal tadi menjadi titik awal untuk penyelesaian Nonlinier Programing (NLP). Metode heuristik lainnya akan dipergunakan untuk memperoleh penyelesaian layak cacah dari problema QAP. Penyelesaian optimum dari model nonlinier programming masih belum membentuk penyelesaian bulat terhadap peubah. Untuk itu perlu diperoleh suatu strategi yang dapat menentukan peubah mana yang dapat bernilai satu dan yang mana bernilai nol. Strategi untuk memperoleh penyelesaian demikian tertera pada langkahlangkah heuristik dibawah ini. Algoritma 3 1. Perhatikan semua peubah yang belum memperoleh nilai cacah, misal∗. nya Xij ∗ , i = 1, . . . , n; j = 1, . . . , n 2. Hitung σij = 1Xij

3. Tentukan himpunan i dengan nilai σij terkecil, misalnya i = i. 4. Buat peubah Xij dengan σij terkecil menjadi 1 (satu), maka σij = 0, misalnya untuk j = j. 5. Buat semua peubah lainnya xij , j 6= j, dalam himpunan terakhir menjadi 0 (nol) dan σij yang berkaitan menjadi 1 (satu) (semua lokasi lainnnya dikeluarkan untuk fasilitas ini) 3. HASIL DAN PEMBAHASAN Masalah Qudratic Assignment Problem (QAP) ini dikutip dari Christoper E.Nugent (1968) dan N. Christofides (1989). Untuk menyelesaikan masalah ini digunakan prosedur penentuan titik awal menggunakan metode heuristik (algoritma 1) sebagai titik awal untuk menyelesaikan metode Tabu Search (algoritma 2), dengan menggunakan bahasa pemograman C++ yang hasilnya merupakan titik awal dari penyelesaian Nonlinier Programming dengan menggunakan bahasa pemograman Fortran kemudian diteruskan dengan metode heuristik lainnya (algoritma 3). Problema yang diselesaikan adalah permasalahan permasalahan Christofides berdimensi N= 20. Problema QAP berdimensi N=20 merupakan masalah bilangan cacah yang cukup besar dimana jumlah variabel binernya adalah 400. Berikut ini disajikan tabel matriks alur F (fik ) yang menyatakan aliran dari


33

pasangan fasilitas i dan k serta matrik ongkos transport C(cjq ) yang menyatakan jarak antara lokasi j dan q Penentuan titik awal menggunakan algoritma 1 1. Menghitung F i untuk setiap fasilitas i, i = 1, . . . , n. dari matriks flow (F ) tabel 9 2. Membentuk daftar fasilitas yang diurutkan secara menurun dari Fi adalah sebagai berkut: Tabel Daftar lokasi menurun untuk N = 20 P No Fi = fik Urutan Menurun 1 101 10 2 226 1 3 82 15 4 204 2 5 148 6 6 145 7 7 19 19 8 24 18 9 113 9 10 172 3 11 91 13 12 128 8 13 6 20 14 168 4 15 157 5 16 92 12 17 51 16 18 93 11 19 84 14 20 40 17 3. Menghitung Cj untuk setiap lokasi j, j = 1, . . . , n dari matriks jarak (C) 4. Membentuk daftar lokasi yang diurutkan secara menaik dari Cj adalah sebagai berikut Tabel Daftar Fasilitas menurun untuk N = 20

34


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

P Cj = cj 101 91 106 82 110 97 105 104 104 80 93 92 87 98 94 87 97 74 107 89

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

5. Tugaskan fasilitas pertama dari tabel 11 dengan lokasi pertama dari tabel 12 diperoleh titik awal i j i j

: : : :

1 15 14 16

2 18 15 13

3 9 16 6

4 10 17 8

5 20 18 17

6 2 19 1

7 19 20 6

8 3

9 11

10 4

11 14

12 12

13 5

Berdasarkan penugasan fasilitas ke lokasi seperti diatas maka fungsi objektif Z yang diperoleh adalah Z = 4819. Penentuan Titik Awal Menggunakan Tabu Search 1. Titik awal penyelesaian yang diperoleh dari algoritma 1 : x(0) :

15 6

18 8

9 17

10 1

20 6

2

19

3

2. Move set µ: { pertukaran dua posisi } 3. Tabu move

11

4

14

12

5

16

13


35

Pertambahan edge yaitu tabu aktif selama satu iterasi dan pengurangan edge yaitu tabu aktif selama dua iterasi. Dengan melakukan langkah algoritma 2 menggunakan bahasa pemograman C++ dengan menggunakan titik awal x(0) dengan iterasi 81840 diperoleh penugasan fasilitas adalah sebagai berikut: i: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 j: 4 10 6 18 8 12 15 19 20 11 16 7 14 9 2 i: 16 17 18 19 20 j: 5 3 1 13 17 yang merupakan titik awal dari Nonlinier Programing. Berdasarkan penugasan fasilitas ke lokasi diperoleh fungsi objektif Z = 2232 Penugasan fasilitas ke lokasi diatas dikonversikan dalam bentuk matrik penugasan berukuran 20 × 20 dimana xij bernilai 0 apa bila tidak ada penugasan dari fasillitas i ke j dan bernilai 1 apabila terdapat penugasan fasilitas i ke lokasi j kemudian matrik penugasan tersebut diselesaikan dengan menggunakan Nonlinier Programming dengan bahasa pemorograman Fortran dan diteruskan dengan algoritma 3 untuk memperoleh nilai fungsi objektif yang lebih baik Z.optimum =1096 dimana penugasan fasilitasnya adalah sebagai berikut: i: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 j : 3 20 7 18 9 12 19 4 10 11 1 6 15 8 2 i : 16 17 18 19 20 j : 5 14 16 13 17

KESIMPULAN Problema Quadratic Assignment Problem (QAP) merupakan permasalahan sulit nonpolinomial. Penyelesaian awal dari problema dilakukan dengan heuristik yaitu metode Tabu Search selanjutnya digunakan Nonlinier Programming Dengan menentukan titik awal menggunakan heuristik Tabu Search, prosedur penyelesaian masalah dapat diperoleh hasil yang optimal . Sifat tidak konveks dari QAP mengakibatkan setiap penyelesaian tergantung pada titik awal penugasan. Untuk proses pengembangan selanjutnya. pengalaman komputasi juga menunjukkan untuk setiap titik awal yang berbeda maka penyelesaian akhir juga berbeda. Sehingga kesimpulan akhir


36

adalah titik awal penugasan merupakan masalah penting dalam penentuan penyelesaian optimal. DAFTAR PUSTAKA 1. Anstreicher K.M and N.W. Brixius (1999), A new Bound for the Quadratic

2.

3.

4.

5. 6. 7. 8.

9.

10.

11.

12.

Assignment Problem Based on Convex Quadratic Programming, Technical Report, Dept of Magement Sciences, University of Iowa. Anstreicher K.M and H.Wolkowigs (1988),. On Lagrangian Relaxation of Quadratic Matrix Constrains, Research Report CORR 98 - 24, Combinatorics and Optimization, Universityof Waterloo,Waterloo,Ontario,Canada. Bazara, M.S and Elshafei, A.N., (1979), An Exact Branch and Bound Procedure for the Quadratic Assignment Problem, NRLQ, Vol.26, pp109120. Brixius N.W and K.M. Anstreicher, (2000), Solving Quadratic Assignment Problem Using Convexs Quadratic Programming Relaxation. Technical Report, Dept of Management Science, Universiy of Iowa. Bruijs, P.A (1984), On the Quality of Heuristic Solution to a 19 x 19 Quadratic Assignment Problem, EJOR, Vol.17, pp.21 30. Burkhad R.E, S.E,.Karisch and F.Rendl (1997), QAPLIB, Quadratic Assignment Problem Library”, Journal of Global Optimization. Christofides, N., Mingozzi, A and Toth, P., (1980)., Constributions to The Quadratic Assignment Problem, EJOR, Vol.4, pp.243 247. Gambardella L.M, E.D.Tailard dan M.Dorigo (1999), ”Colonies for the Quadratic Assignment Problem”, Journal of the Operational Research Society, 50:167-176. Hahn P.M (2000), Progress in Solving The Nugent Insta Of Quadratic Assignment Problem. Working Paper, System Enginering, University of Pennsylvania. Mawengkang, H and Murtagh, B.A (1985), Solving Nonlinier Integer Programs with Larger Scale Optimization Software, Annalas of Operation Research, Vol.5, pp. 425 437. Murtagh B.A and Souders, M.A (1982), A Projected Lagrangian Algoritm and Its Implementation fors Sparse Nonlinier Constrain, Mathemathical Progamming Study, Vol.16, pp.84 117. Resende MG.C, P.M.Pardalos and Y.Li (1996), Algoritm 754: Fortran Soubroutines for Approximate Solution of Dense Quadratic Assigment Problem Using GRASP. ACM Transaction on Mathematical Software, 22:104-118.


37

13. Tailard E.D (1991), Robust Taboo Search for the Quadrati Assgnmet Prob-

lem, Paralel Computing, 17:443-455. 14. U.V. Thonemann and A.Bolte (1994), An Improved Stimulatin Annealing Algoritm for the Quadratic Assignment Problem, Working Paper, School of Business, Departement of Production And Operations Research, University of Paderborn, Germany.

Rosmaini: Jurusan Matematika, FMIPA, Universitas Sumatera Utara, Jl. Biote-




38

Irvan – Developing sensitivity analysis method

39


DEVELOPING SENSITIVITY ANALYSIS METHOD IN MINIMUM COST FLOW PROBLEM

Irvan and Herman Mawengkang Abstract. The minimum cost flow problem is a linear optimization model which can be solved by Simplex method. This paper proposed a sensitivity analysis method for minimum cost flow problem such that maintains the structure of the optimal solution. We improved a sensitivity analysis method which is formerly applicable to a tree solution. the proposed method preserves solution of upper bounds at upper bounds and those of lower bounds at lower bounds.

1. INTRODUCTION Let G = (N, E) be a directed network defined by a set N = (n1 , n2 , · · · , nn ) nodes and a set E of a directed arcs. Each arc (i, j) ∈ E has an associated cost cij that denotes the cost per unit flow on that arc. We assume that the flow cost varies linearly with the amount of flow. We also associate with each arc (i, j) ∈ E a capacity aij that denotes the maximum amount that can flow on the arc and lower bound lij that denotes the minimum amount that must flow m the arc. We associate with each node i ∈ N an integer value bi representing its supply/demand. If bi > 0, node i is a supply node; if bi < 0 node i is a demand node with a demand of −bi ; and if bi = 0 Received 12-04-2008, Accepted 04-09-2008. 2000 Mathematics Subject Classification: 49Q12, 78M50 Key words and Phrases: cost flow, sensitivity, simplex method


40

node i is a transhipment node. The decision variables m the minimum cost flow problem are arc flow which is represented by xij on (i, j) ∈ E. The model of the minimum cost flow problem (MCFP) can be expressed as follow. X cij xij ........................... (P ) min (i,j)∈E

subject to X

xij −

j:(i,j)∈E

X

xji = bi , ∀i ∈ N

j:(i,j)∈E

ìj ≤ xij ≤ aij , ∀(i, j) ∈ E where

n X

bi = 0

i=1

We would like to consider sensitivity analysis with respect to changes in the cost coefficients in problem (P). The purpose of sensitivity analysis of the problem mentioned above is to determine changes in the optional solution of the problem resulting from changes in the data (supply/demand vets, capacity or cost of any arc) [4,5,7,8,9,10,11]. Traditionally, researchers have inducted this sensitivity analysis using primal is dual simplex method [7]. There is, however a conceptual draw bad to this approach. The simplex based approach maintains a basis free at every iteration and conduct sensitivity analysis by determining changes in the basis free precipitated by changes in the data. The basis in the simplex method in often degenerate and consequently changes in the basis tree are not necessarily translated into the changes in the optional solution. Therefore, the simplex based approach does not give information about the changes in the solution as the data changes; instead it tells us about the changes in the basis free [1,2,3,7]. To avoid there short comings, we propose another approach in the sensitivity analysis. Let x∗ij be an optimal solution of the minimum cost flow problem this optimal solution partitions the arc set E in to three subsets as follows: B = {(i, j) ∈ E; ìj < x∗ij < aij } L = {(i, j) ∈ E; x∗ij = ìj } U = {(i, j) ∈ E; x∗ij = aij } We call the triple (B, L, U ) as the optional solution structure. Definition 1.1 Given an optimum solution x∗ij of the minimum cost flow problem. SA is to determined changer in the data that the optional solution structure remains unchanged


41

2. COST SENSITIVITY ANALYSIS Here we consider SA with respect to changes in the cost coeffiecents. Suppose first that an optimal solution x∗ij of the MCFP is given. If the cost cij of an arc (i, j) changes to cîj , then its corresponding optional solution x∗ij or the optional solution structure (B, L, U ) may also change. Hence, cost sensitivity analysis should determine changes in the cost of any arc such that the given optional solution stretare (B, L, U ) is unchanged. First, given a feasible solution xij of the MCFP, xij is an optimum solution of the MCFP if and only if for some set of node potentials w, the reduced costs and flow values satisfy the following complementary slackness optimality condition for every are (i, j) ∈ E [1,6,7] If c¯ij = cij − wi + wj > 0, the xij = ìj ; If c¯ij = cij − wi + wj < 0, the xij = uij ; If c¯ij = cij − wi + wj = 0, the xij ≤ ìj ≤ uij . Thus, given an optimum solution x∗ij of the MCFP, SA with respect to changes in the cost coefficient in equivalent to determining cost range cij satisfying the following conditions : If x∗ij = ìj , the cîj − wi∗ − wj∗ ≥ 0 If x∗ij = uij , the cîj − wi∗ − wj∗ ≤ 0 If ìj ≤ x∗ij ≤ uij , the cîj − wi−1 = 0 Now, we consider a concept for the method of calculating SA for cost coefficient. Given an optimal solution x∗ij of the MCFP, suppose that the cost cij of an arc (i,j) is changed to cij . If there is a path P (i, j) that can reroute the flow x∗ij from node i to node j with less cost than cij without violating any of the optimality conditions, we should reroute the flow x∗ij along P (i, j). This rerouting changes the optimal solution structure (B, L, U ) due to a change in arc flow. Thus range of cost coefficients that maintains the optimal structure (B, L, U ) is the range right before the point where the flow of arc is changed. Now suppose that the cost cij of an arc (i,j) increases by cîj = cij + θ. Then we want to know the range of θ, that is θe ≤ θ ≤ θu , where θu (θe ) denote the amount of maximum increasing (decreasing) flow that preserves the optimal solution structure (B, L, U ). We first consider the following network G = (N, E) with multiple path P (i, j) from node i to node j. Suppose that the cost cij of an arc (i, j) is changed to cîj . This changes the reduced cost optimality conditions to restore the optimality condition


42

of the arc, we must change the flow on arc (i, j). We can both increase and decreare the flow in an arc (i, j) while knowing the bounds an arc flows. However, in an arc (i, h) at its lower bound (i.e., xih = ìh ) we can only increase the flow, and for flow on an arc (i, m) at its upper bound (i.e., xim = uim ) we can only decrease the flow. Using there concepts, we can find two paths P 1 (i, j) = {(i, h), (h, j)} and P 2 (i, j) = {(k, i), (k, j)} that can send the flow from node i to node j. In our example, per unit cost to send one unit of flow along the path P 1 (i, j) is cih + chj . If cih + chjk < c¯ij , then the flow is sent along P 1 (i, j). In this case the optimal structure (B, L, U ) may change. But, if the per unit cost to send one unit of flow along P 1 (i, j) is greates than the per unit cost c¯ij , the flow is not sent along P 1 (i, j). In this case since the flow along P 1 (i, j) is not created, we can maintain the given optimal solution structure (B, L, U ) in the interval of θ ≤ −cij + eih + ehj . Similarly, θu1 the per unit cost to send one unit of flow along P 2 (i, j), is −cki + chj . If −cki + chj > cîj the flow 2 and the structure (B, L, U ) can be preserved in will not be sent along Pi,j the interval of θ ≤ −cij + eih + ehj . Here, let θu1 and θu2 be the upper bound of θ with respect to the path P 1 (i, j) and P 2 (i, j) can preserve the given structure (B, L, U ). Then θu1 = −eij + cih + chj and θu2 = −eij − ehi + ehj . Consequently, θu that maintains the optimal structure (B, L, U ) is restricted in value of θu = min(θu1 , θu2 ). In order to do there calculations easily consider a definition of transformed network. Definition 2.2 The transformed network G0 = (N, E 0 ) corresponding to a network G = (N, E) of the MCFP is defined as follows. We replace each arc (i, j) ∈ E with flow lij < x∗ij < uij by two arcs (i, j) and (j, i) with the cost c¯ij = cij and c¯ji = −cji respectively. We also replace each arc (i, j) ∈ E with the flow x∗ij = uij by arc (j, i) with the cost c¯ji = −cji . Finally, we replace each arc (i, j) ∈ E with the flow x∗ij = ìj by arc (i, j) with the cost c¯ij = cij . Using the above transformed network G0 = (N, E 0 ) , we can calculate θu by summing up the length of the directed path P (i, j) from node i to node j and −cij in G0 = (N, E 0 ). If there exists multiple directed paths P (i, j) form node i to node j, calculate θu by summing the minimum length of the directed path from node i to node j, calculate θu by summing the minimum length of the directed path among various directed path and −cij in G0 = (N, E 0 ). Since among the various directed path the length of the directed path associated with the shortest path from node i to node j is


43

the minimum, we need to calculate the length of the shortest path for the calculation of θu . Lemma 2.1 The transformed network G0 = (N, E 0 ) contains a negative cycles. Proof Assume the transformed network G0 = (N, E 0 ) contains a negative cycle. This means that the network G = (N, E) does not satisfy the optimality condition, and we can improve the current optimal solution with respect to this negative cycle. Because a feasible solution which is less than the objective function value of the given optimal solution x∗ij exists, it contradicts the assumption that the optimum solution in given. Therefore, the transformed network G0 = (N, E 0 ) contains no negative cycles. Lemma 2.2 If the directed path P (i, j) from node i to node j in G0 = (N, E 0 ), the θu =∼. Proof If the directed path P (i, j) from node i to node j in G0 = (N, E 0 ) does not exist, there is no path P (i, j) that can send the flow from node i to node j in the network G = (N, E). Therefore, even if the cost cij of an arc (i, j) increases infinitely, it is optimal to send the current flow x∗ij along the arc (i, j). So if the directed path P (i, j) does not exist in G0 = (N, E 0 ), the θu =∼. Ultimately, if the shortest path exist in G0 = (N, E 0 ) we obtain a constant as the upper bound value θu . But, if the shortest path does not exist in G0 = (N, E 0 ) we obtain an infinite value as the upper bound. Theorem 2.1 If G0 = (N, E 0 ) contains a shortest path from node i to node j and ` is the length of the shortest path, then θu = ` − cij . If G0 = (N, E 0 ) contaion no shortest path from node i to node j, then θu =∼ Proof By lemma 2.1, G0 = (N, E 0 ) contains no negative cycle. Therefore if the shortest path exists in G0 = (N, E 0 ), we can solve it. If the shortest path does not exist in G0 = (N, E 0 ), the directed paths P (i, j) does not exist and by lemma 2.2. θu =∼. If the are multiple directed path P (i, j) from node i to node j, the θu is obtained by summing the minimum length of the directed path among them, and −cij in G0 = (N, E 0 ). let Ph (i, j) be the length of the both directed path among the multiple directed paths. Then l = min{Ph (i, j)} and θu = ` − cij . Next we consider cij = cij − θ. Because per unit cost to send are unit of flow along the arc (i, j) is decreased by cîj ,


44

more flow along the arc (i, j) sent instead of sending flow along the path from node i to node j. Here, the meaning of sending more flow along arc (i, j) is equivalent to sending the flow along the path from node j to node i. Therefore in this case after comparing the per unit cost to send one unit of flow along the path from node j to node i with the per unit cost to send one unit of flow along the arc (i, j) from node j to node i, a flow along the path associated with the cheaper cost of the two cases is created. In this case, θl is the sum of the length of the directed path P (j, i) from node j to node i and Cij in G0 = (N, E 0 ). If there are multiple directed paths P (j, i) from node j to node i, θl is calculated by summing the minimum length of the directed path and cij G0 = (N, E 0 ). Here, let Ph (j, i) be the length of the both directed path among the multiple directed paths. Then the length l, the shortest path from node j to node i is ` = min{Ph (j, i)}. Therefore, if the shortest path from node j to node i exists in G0 = (N, E 0 ) we obtain θ` = ` + cij . But if the shortest path does not exist in G0 = (N, E 0 ), θ` =∼ by lemma 2.2. In case the number of arcs, with a flow of the ìj < xij < uij , is greater than or equal to n − 1, we can easily compute the upper bound or the lower bound. Given a spanning tree optimal solution where the number of arcs with a flow of ìj < xij < uij is exactly equal to n − 1, we can obtain the node potentials w associated with the spanning tree optimal solution because arcs with a flow of ìj < xij < uij consist of the optimum basis of the network G = (N, E). By using node potentials w we can compute the length of the shortest path more easily because the node potentials w is the length of the tree path with respect to the optimum basis O. Theorem 2.2 If O is an optimal basis and w is the node potentials with respect to O, then the length of the tree path P (r, s) from node r to node s is wr − ws . Proof We compute P (r, s) using the fact that c¯ij = cij − wi + wj = 0, for every arc (i, j) ∈ O X X c¯ij = (cij − wi + ws ) (i,j)∈P (r,s)

(i,j)∈P (r,s)

X

=

(i,j)∈P (r,s)

X

=

(i,j)∈P (r,s)

=0

cij −

X (i,j)∈P (r,s)

cij − wr + ws

(wi − ws )


45

because all w corresponding to the nodes in the path, P other than the terminal nodes r and s, cancel each ofher in the expression (i,j)∈P (r,s) (wi −wj ) since P the expression (i,j)∈P (r,s) cij is associated with the length of path P (r, s), the length of the tree path P (r, s) is wr − ws . 3. DISCUSSION Generally, given the optimal basis associated with spanning tree O, because the length of the shortest path from node i to node j is equivalent to the length of the tree path P (i, j) from node i to node j, we can compute the sensitivity analysis more easily by calculating the length of the tree path using the node potentials. By this result, SA2 and SA1 are the same in spanning tree solution [2]. Given an optimal solution that the number of arcs with a flow of ìj < xij < uij is greater than n − 1, there arcs with the intermediate flow always form cycle. Let c be a cycle that consist of arcs with intermediate flow, ìj < xij < uij . Then the results of the sensitivity analysis with respect to an arc (i, j) that belong to cycle c is given in theorem 3. Theorem 3.3 Given a non-tree optimal solution where the number of arcs with a flow of ìj < xij < uij is greater than n − 1, the range θ of SA2 with respect to an arc (i, j) that belongs to a cycle c is zero. Proof Given a non-tree optimal solution arcs with intermediate flow always form cycles. PLet C be a cycle that consists of arcs with the intermediate flow. Then (i,j)∈c cij = 0 if the cost of cij of an arc (i, j) is changed to c¯ij = cij + θ, then the value of the cycle C is X X (cij + θ) + clh = clh + θ = θ (l,h)∈c

(l,h)∈c

If θ > 0(orθ < 0) then it is better to send the flow in the opposite direction (or same direction) of cycle to improve the objective function value. In this case, since the flow along cycle C is created, the given optimal solution structure (B, L, U ) may change. Therefore in order to maintain the given optimal solution structure (B, L, U ), the range of θ has to be zero. CONCLUSION In this paper, we have categorized the sensitivity analysis of the minimum cost flow problem into two types. We define SA1 to be the acquirement of a region where the given optimum basis is unchanged. This is the


46

well known method applicable to a tree solution. However SA1 can not be applied to a non-tree solution or a degenerate tree solution. So we proposed the SA2 that finds the region where the upper bound valued arcs in the optimum solution maintain upper bound valued, low bound valued arcs maintain lower bound valued, and intermediate valued arcs maintain intermediate valued. This SA2 provider the generalized concept of the sensitivity analysis for the optimum solution of the minimum cost flow problem.

REFERENCES 1. Ravindra K. Ahuji, T.L. Magranti, and J.B. Orlin (1993). 2.

3. 4. 5. 6. 7. 8.

9. 10.

11.

”Network

Flows”, Prentice - Hall.Inc. Hoyeon Chung (1994). ”A Study on tje Sensitivity Analysis for a Non-tree Solution and the Performance Improvement of the minimum Cost Flow Inable”, Dissertation of the Seul Nation University. W.H. Cunningham (1979). ”Theoritical Properties of the Network Simplex Method”, Math. of. open. Res. Vol. 4, No. 2. pp196 - 20P. Tomas Gal (1979). ”Postoptimal Analysis, Parametric Programing and Related Topics”, McGraw Hill. Dan Gusfield (1983). ”A Note on Arc Tolerances in Sparse Shortest Path and Network Flow Problems”, Networks, Vol. 13, pp. 191 - 196. K.G. Marty (1992). ”Network Programming”, Prentice - Hall Inc. G.L. Nemhauser and Kan (1989). ”Hand Books in Operations Research and Management ”, Sci. Vol. 1, North - Holland, pp. 321 - 323. N. Ravi and R.E. Wendell (1988). ”The Tolerance Approach to Sensitivity Analysis in Network Linear Programming”, Networks, Vol. 18, pp. 159 -171. D.R. Shier and Ch. Witzgall (1980). ”Arc Tolerance in Shortest Path and Network Flow Problems”, Networks, Vol. 10, pp. 227 - 291. Robert E. Tarjan (1982). ”Sensitivity Analysis of Minimum Spanning Tree and Shortest Path Trees”, Information Processing Letters, Vol. 14, No. 1, pp. 30 - 33. R.E. Wendell (1985). ”The Tolerance Approach to Sensitivity Analysis in Linear Programming”, Management Sci, Vol. 31, No. 5, pp. 564 - 578.


47

Irvan: Program Studi Pendidikan Matematika, FKIP, Universitas Muhammadiyah

Sumatera Utara Medan 20238, Indonesia.


Mawengkang: Jurusan Matematika, FMIPA, Universitas Sumatera Utara, Jl.

Bioteknologi No. 1 Kampus USU, Medan 20155, Indonesia.



48

Marwan – Teori linear pembangkit gelombang

49

Vol. 01, No. 01 (2009), pp. 49–57.

TEORI LINEAR PEMBANGKIT GELOMBANG MULTI ARAH TIPE SNAKE UNTUK FLAP TUNGGAL

Marwan Abstract This paper discusses a mathematical model of three-dimensional linear wave generation in multi-direction type of snake for a single flap. Wave generation of this type is performed by moving some of the flap which is not uniform. The referred model is the Laplace equation and the boundary condition, that is lateral boundary condition at the generation of waves, boundary condition at the bottom and the dynamic boundary condition and kinematic free surface by using linear wave theory approach. For simplicity, assume that the water in the form of an ideal fluid; homogeneous, inviscid, irrational, incrompressible. By using the linear theory of water waves, the relationship between the deviation of the wave generator with generated wave parameters, such as wave amplitude and the number will be derived. In addition, the influence of the position flap was also investigated.

1. PENDAHULUAN Kondisi geografis Indonesia, yang hampir delapan puluh persen wilayahnya terdiri atas lautan, menyebabkan topik-topik yang berkaitan dengan gelombang air merupakan permasalahan strategis yang sangat relevan untuk disimak, ditelaah dan dikembangkan. Hal ini disebabkan karena permasalahan gelombang air, secara tidak langsung, sangat besar pengaruhnya bagi kehidupan perekonomian di Indonesia. Pengembangan lalu lintas laut, penambangan minyak dan gas lepas pantai, penangkapan ikan, masalah abrasi, sedimentasi dan pengikisan garis pantai adalah contoh-contoh kasus Received 10-04-2008, Accepted 09-09-2008. 2000 Mathematics Subject Classification: 76B15, 76B70 Kata kunci: fluida ideal, gelombang multi arah, teori linear.


50

yang tidak terlepas dari studi gelombang air. Pada kasus yang terkait dengan lalu lintas laut, sering kali dikontruksi suatu model fisis terskala yang akan diuji di kolam pengujian pada laboratorium hidrodinamika. Sebagai contoh, sebelum suatu kapal dibangun, diperlukan rancangan model fisis yang berukuran kecil, yaitu berupa miniature kapal untuk diuji kelayakannya melalui pengujian dengan gelombang laut tiruan yang dibangkitkan di laboratorium. Pada pengujian model fisis ini, perlu dibangkitkan gelombang yang skalanya disesuaikan dengan karakteristik gelombang di mana kapal ini akan dioperasikan. Keakuratan pengujian ini sangat bergantung pada ketetapan pembangkitan gelombang yang dilakukan di laboratorium hidrodinamika. Pembangkit gelombang pada laboratorium hidrodinamika ada dua jenis, yang pertama pembangkit gelombang multi arah yang dilakukan di wave basin. Salah satu contoh pembangkitan gelombang multi arah ini adalah tipe snake. Tipe ini terdiri dari beberapa flap, yang masingmasing membuat gerakan yang berbeda-beda. Contoh kolam pengujian (wave basin) pada laboratorium hidrodinamika disajikan pada Gambar 1.

Gambar 1: Kolam pengujian untuk jenis gelombang multi arah pada laboratorium hidrodinamika (www.hazama.co.jp). Yang kedua adalah pembangkit gelombang satu arah yang dilakukan pada towing tank. Pembangkitan gelombang jenis ini dilakukan dengan cara menggerakkan beberapa flap secara seragam. Untuk kasus pembangkitan gelombang satu arah tipe flap yang diletakkan pada dasar kolam telah diteliti Dean dan Dalrymple (1991), Karjanto (2001) dan Iskandar (2004). Sementara itu, untuk kasus di mana flap diletakkan pada posisi tertentu


51

telah pula diselidiki oleh Kusumawinahyu, dkk (2005) dan Marwan (2006). Selain itu, Kusumawinahyu telah pula menyelidiki untuk kasus flap ganda [6]. Tulisan ini memaparkan model matematis pembangkitan gelombang multi arah tipe snake. Secara khusus, di sini akan dipresentasikan ketertarikan antara simpangan pembangkit gelombang dengan amplitude, bilangan, arah rambat dan posisi di mana flap diletakkan. Perlu ditambahkan bahwa dalam kasus ini, flap diletakkan pada posisi dengan jarak tertentu di antara permukaan rata-rata air (SWL) dan dasar kolam.

2. MODEL MATEMATIKA Asumsikan air adalah fluida ideal: tak kental, tak termampatkan dan tak berotasi secara lokal [5]. Dengan demikian, perambatan gelombang air dalam arah horizontal dapat dimodelkan oleh persamaan Laplace, yang dituliskan dalam bentuk ∂2Φ ∂2Φ ∂2Φ + + = 0, ∂x2 ∂y 2 ∂z 2

(15)

dengan −h ≤ z ≤ η(x, y, t). Di sini P hi dan η berturut-turut menotasikan potensial air dan elevasi di permukaan. Syarat batas di dasar kolam, di permukaan bebas dan di pembangkit dituliskan sebagai berikut ini.

• Di dasar kolam ∂Φ = 0, z = −h ∂z

(16)

• di permukaan bebas z = η(x, y, t) ada dua syarat, yaitu syarat batas kinematik ∂η ∂η ∂Φ ∂η ∂Φ ∂Φ + + = ∂t ∂x ∂x ∂y ∂y ∂z

(17)

dan syarat batas dinamik gη +

∂η 1 ∂Φ 2 ∂Φ 2 ∂Φ 2 + [( ) +( ) +( ) ]=0 ∂t 2 ∂x ∂y ∂z

(18)


52

• di pembangkit gelombang x = 0, biasanya disebut juga sebagai syarat batas lateral. Posisi flap S(z) dituliskan sebagai x=

S(z) sin(ky y − ωt) 2

(19)

atau dapat juga dituliskan dalam bentuk F (x, y, z, t) = x −

S(z) sin(ky y − ωt) = 0, 2

(20)

dengan ω adalah frekuensi gerakan flap.

Gambar 2: Pembangkit gelombang flap tunggal dengan persamaan pengatur yang disertai dengan syarat batasnya. Selanjutnya dengan menghitung turunan total dari (20) diperoleh u−

D Dt

=

∂y ∂ ∂ ∂x ∂ ∂z ∂ ∂t + ∂t ∂x + ∂t ∂y + ∂t ∂z

S(z) 1 dS(z) ω cos(ky y − ωt) = w sin(ky y − ωt), 2 2 dz

(21)

∂Φ ∂Φ dengan (u, v, w) = ( ∂Φ ∂x , ∂y , ∂z ) adalah kecepatan. Persamaan (21) dapat juga dituliskan sebagai ∂Φ 1 ∂Φ dS(z) = sin(ky y − ωt) + ωS(z) cos(ky y − ωt) . (22) ∂x 2 ∂z dz


53

3. PELINEARAN Beradasarkan uraian di atas terlihat bahwa syarat batas pada permukaan dan di pembangkit gelombang adalah tak linear. Secara umum, cukup sulit untuk menentukan solusi untuk menentukan eksak dari permasalahan syarat batas di atas. Oleh karena itu, solusi diperoleh dengan menggunakan pendekatan linear dari syarat batas. Berikut ini dipaparkan pelinearan terhadap syarat batas di pembangkit gelombang. Untuk jarak yang cukup kecil mengakibatkan S(z) bergerak dengan kecepatan yang cukup kecil. Dengan demikian, suku pertama bagian kanan dari persamaan (22) dapat diabaikan. Akibatnya diperoleh ∂Φ S(z) S(z) −ω cos(ky y − ωt), x = sin(ky y − ωt). ∂x 2 2

(23)

Ekspansi Taylor bagian kiri persamaan (23) di sekitar x = 0 memberikan ∂Φ S(z) − ω cos(k y − ωt) y S(z) ∂x 2 x= 2 sin(ky y−ωt) ∂Φ S(z) = −ω cos(ky y − ωt) ∂x 2 x=0 d ∂Φ S(z) S(z) sin(ky y − ωt) −ω cos(ky y − ωt) + 2 dx ∂x 2 x=0 + ... (24) Terlihat bahwa hanya suku pertama dari ekspansi Taylor tersebut yang linear dalam ∂Φ ∂x dan S(z). Dengan demikian diperoleh syarat batas lateral linear, yang dapat dituliskan sebagai u(0, y, z, t) = ω

S(z) cos(ky y − ωt). 2

Dengan menggunakan cara yang sama diperoleh syarat batas kinematik linear dan dinamik linear yang berturut-turut dapat dituliskan sebagai ∂η ∂Φ = ∂t ∂z dan η=−

1 ∂Φ , g ∂t


54

pada z = 0. Akibatnya permasalahan syarat batas di atas dapat pula dituliskan sebagai berikut ∂2Φ ∂2Φ ∂2Φ + + = 0, ∂x2 ∂y 2 ∂z 2

(25)

dengan −h ≤ z ≤ η(x, y, t) ∂Φ = 0, z = −h ∂z ∂η ∂Φ = ,z = 0 ∂t ∂z 1 ∂Φ η=− ,z = 0 g ∂t u(0, y, z, t) = ω

S(z) cos(ky y − ωt). 2

(26) (27) (28) (29)

4. SOLUSI PERSAMAAN PENGATUR Secara umum, solusi persamaan Laplace dengan syarat batas di dasar dan permukaan bebas dapat dituliskan sebagai [1] Φ(x, y, z, t) = Ap cosh(kp (h + z)) cos(kx x − ky y − ωt) +Ce−ks x cos(ks (h + z)) cos(ky y − ωt).

(30)

Suku pertama dari Φ(x, y, z, t) berkorespondensi dengan gelombang berjalan, sedangkan suku yang kedua berasosiasi dengan gelombang sendiri. Bilangan gelombang dari gelombang berjalan kp dan bilangan gelombang dari gelombang berdiri ks berelasi dengan frekuensi ω oleh ω 2 = gkp tanh kp h

(31)

ω 2 = −gks tanh ks h.

(32)

dan Grafik dari solusi kedua persamaan ini disajikan dalam Gambar 3, untuk 2 kasus ωg h = 1. Jelas bahwa solusi persamaan (32) tak berhingga banyaknya. Akibatnya solusi untuk permasalahan syarat batas di atas dipilih berbentuk


55

Gambar 3: Pembangkit gelombang flap tunggal dengan persamaan pengatur yang disertai dengan syarat batasnya.

Φ(x, y, z) = Ap cosh(kp (h + z)) cos(kx x − ky y − ωt) −ks (n)x + Σ∞ cos(ks (n)(h + z)) cos(ky y − ωt), (33) n=1 Cn e

dengan Ap dan Cn akan ditentukan. Substitusi kecepatan potensial ke dalam syarat batas di pembangkit gelombang (linear) memberikan q ωS(z) ∞ = Ap kx cosh(kp (h + z)) − Σn=1 Cn ks2 (n) + ky2 cos(ks (n)(h + z)). 2 (34) Diketahui bahwa berdasarkan Sturn-Liouville himpunan (cosh(kp (h + z)), cos(ks (n)(h + z)), n = 1, 2, ) adalah himpunan yang ortogonal, sehingga Z 0 cosh(kp (h + z)) cos(ks (n)(h + z)) = 0 −h

dan Z

0

cosh(kp (m)(h + z)) cos(ks (n)(h + z)) = 0, −h

m 6= n. Kemudian nilai Ap ditentukan dengan cara mengalikan kedua ruas persamaan (34) dengan cosh(kp (h + z)) dan dintegralkan atas doaminnya R 0 ωS(z) cosh(kp (h + z))dz . (35) Ap = −h R 0 2 kx −h cosh2 (kp (h + z))dz


56

Selanjutnya nilai Cn ditentukan dengan cara mengalikan kedua ruas persamaan (34) dengan cos(ks (n)(h + z)) dan dintegralkan atas doaminnya, sehingga R 0 ωS(z) cos(ks (n)(h + z))dz Cn = q −h 2 R . (36) 0 ks2 (n) + ky2 −h cos2 (ks (n)(h + z))dz Terlihat bahwa nilai Ap dan Cn bergantung pada S(z). Dalam kasus ini S(z) = 1 + dz s, −d ≤ z ≤ 0 dan S(z) = 0, 0 ≤ z ≤ d. Dengan demikian diperoleh 4s[kp d sinh(kp h) − cosh(kp h) + 1] Ap = (37) kx kp d[sinh(2kp h) + 2kp h] dan

−4s[ks (n)d sin(ks (n)h + cos(ks (n)h)) + 1] Cn = p . ks2 (n) + λ2 ks (n)d[sin(2ks (n)h) + 2ks (n)h]

(38)

KESIMPULAN Telah diturunkan formula proses pembangkitan gelombang multi arah tipe snake di laboratorium hidrodinamika. Snake di sini adalah gabungan dari beberapa flap. Pendekatan teori gelombang linear untuk pembangkit gelombang tipe flap dengan menggunakan persamaan lengkap tak linear untuk gelombang air telah dipaparkan. Di sini didiskusikan teori linear untuk kasus flap tunggal. Dengan menyelesaikan persamaan pengatur menggunakan syarat batasnya, diperoleh elevasi gelombang permukaan yang memuat gelombang berjalan dan gelombang berdiri. Gelombang berdiri bersifat meluruh seiring dengan perambatannya menjauhi pembangkit gelombang. Diperoleh bahwa relasi yang mengaitkan amplitude gelombang terbangkit dan simpangan pembangkit gelombang bergantung pada ke dalaman, posisi flap, frekuensi gelombang dan arah rambat gelombang. UCAPAN TERIMA KASIH Penulis mengucapkan terima kasih kepada Fajariani, S.Si, Nanda Maulina, S.Si dan Laila Qadriah, S.Si atas segala bantuannya sehingga penelitian ini dapat terlaksana dengan baik. Terima kasih juga disampaikan kepada Ketua Jurusan Matematika FMIPA Unsyiah atas dukungan fasilitas pada pelaksanaan penelitian ini.


57

REFERENSI 1. R. G. Dean and R. A. Dalrymple, Water Wave Mechanics for Engineers and Scientists, Advanced series on Ocean Engineering, World Scientific, Singapore (1991). 2. C. J. Galvin, “Wave Height Prediction for Wave Generators in Shallow Water”, Tech Memo 4, U.S. Army, Coastal Engineering Research Center (1964). 3. N. Karjanto, Teori Pembangkit Gelombang Dua Dimensi Tipe Flap, undergraduate final project, Department of Mathematics Institut Teknologi Bandung (2001). 4. Iskandar, Studi Pembangkitan Gelombang Satu Arah Tipe Flap, Report of Student Research Grant, TPSDP Unsyiah (2004). 5. G. B. Whitham, Linear and Non-Linear Waves, John Wiley and Sons, New York (1974). 6. W.M. Kusumawinahyu, N. Karjanto, G. Klopman, “Linear theory for single and double flap wavemakers”, J. Indones. Math. Soc., Vol. 12 No. 1, (2005). 7. W.M. Marwan, “Linear wave maker theory for single flap”, Jurnal Natural, FMIPA Unsyiah, Vol. 6 No. 2 (2006).

Marwan: Jurusan Matematika, FMIPA, Universitas Syiah Kuala, Jl. Tgk. Daud

Beureuh No. 4 Darussalam-Banda Aceh 23111, Indonesia.



58

Taufiq – Prediksi pasang surut laut

59

Vol. 01, No. 01 (2009), pp. 59–71.

PREDIKSI PASANG SURUT LAUT DI SELAT MALAKA DENGAN MENGGUNAKAN HAMSOM MODEL

Taufiq Iskandar Abstrak. Simulasi Model Pasang Surut Laut (Tides) Tiga Dimensi di Selat Malaka dengan menggunakan HAMSOM model memperlihatkan bahwa jenis pasang surut laut di Selat Malaka adalah semi diurnal. Pasang surut laut di daerah tenggara lebih tinggi dibandingkan dengan pasang surut di barat laut, dimana puncak ketinggian pasang surut laut terjadi ketika bulan purnama dan bulan baru. Arus laut di selat malaka secara umum bergerak dari tengggara menuju barat laut.

1. PENDAHULUAN Selat Malaka merupakan selat yang terletak di antara semenanjung Malaysia dan pulau Sumatra. Dari segi ekonomi dan strategis, selat Malaka merupakan salah satu jalur pelayaran terpenting di dunia, sama pentingnya seperti terusan Suez atau terusan Panama. Selat Malaka membentuk jalur pelayaran terusan antara samudra Hindia, laut Cina Selatan dan samudra Pasifik serta menghubungkan tiga dari negara-negara dengan jumlah penduduk terbesar di dunia: India, Indonesia dan Republik Rakyat Cina. Pada tahun 1993 dan 1995 lebih 100.000 kapal melintasi selat Malaka setiap tahunnya, serta mengangkut 3.23 juta barel minyak per hari. Kecelakaan kapal juga sering terjadi di selat Malaka yang menandakan bahwa selat Malaka Received 30-03-2008, Accepted 10-11-2008. 2000 Mathematics Subject Classification: 86A05, 76Q05 Kata kunci: Selat malaka, pasang surut laut, HAMSOM model


60

merupakan salah satu dari kemacetan lalu lintas terpenting di dunia (ThiaEng, Gorre, Ross, and Regina, 2000). Secara umum penelitian di selat Malaka secara ilmiah (Scientific) sangat jarang dilakukan sehingga peneliti ingin mengkaji secara mendalam mengenai pasang surut di selat tersebut. Adapun penelitian terdahulu yang telah dilakukan di selat Malaka mengenai sedimen yang menjelaskan karakteristik sedimen di lingkungan laut (Keller and Richards, 1967). Sedangkan penelitian mengenai pasang surut laut dengan menggunakan model tiga dimensi juga sudah pernah dilakukan oleh Rizal (1993) yaitu dengan menggunakanan model hidro dinamik tiga dimensi. Begitu juga dengan energi pasang surut dan momemtum di selat Malaka yang merupakan bagian dari perairan Indonesia telah diteliti oleh Mihardja (1991) dengan menggunakan model semi implisit dua dimensi. Pada penelitian kali ini penulis ingin meneliti masalah pasang surut laut di selat Malaka dengan menggunakan model Hamburg Shelf Ocean Model (HAMSOM). Dimana model ini sudah pernah diaplikasikan pada beberapa tempat seperti di Bohai sea, laut Jawa dan laut Utara. 2. MODEL HAMSOM Model HAMSOM pertama kali disusun tahun 1983 oleh Prof. J Backhaus (Backhaus, 1983) dan (Backhaus, 1985) dan tahun-tahun selanjutnya dikembangkan baik secara teoritis persamaan dan perhitungan masingmasing sukunya, maupun penyelesaian numeriknya, hingga penambahan input dan initialisasi modelnya sesuai dengan kajian yang diinginkan (Huang, 1995 ; Pohlmann, 1996; Putri, 2005). Model HAMSOM ini telah sukses dikembangkan di Bohai sea oleh Huang (1995), begitu juga di laut Jawa, selat Sunda dan bagian timur samudra Hindia oleh Putri (2005). Persamaan dasar yang digunakan pada penelitian ini adalah persamaan gerak dinamika oseanografi (Pond and Pickard 1983). Persamaan dinamika oseanografi dalam arah x dan y dapat ditulis sebagai berikut ∂u ∂u ∂u ∂u 1 ∂p ∂ ∂u +u +v +w − fv = − + AH ∂t ∂x ∂y ∂z ρ ∂x ∂x ∂x ∂ ∂u ∂ ∂u + AH + Av ∂y ∂y ∂z ∂z


61

∂v ∂v ∂v ∂v 1 ∂p ∂ ∂v AH +u +v +w + fu = − + ∂t ∂x ∂y ∂z ρ ∂y ∂x ∂x ∂ ∂v ∂ ∂v + AH + Av (39) ∂y ∂y ∂z ∂z Dalam arah z persamaan hidrostatika disederhanakan dalam bentuk ∂p = −ρg ∂z Persamaan differensial parsial diintegrasikan dalam arah vertikal untuk setiap lapisan sehingga diperoleh solusi persamaan medan transpor rata-rata lapisan (U, V ). Penyelesaian persamaan gerak ratarata lapisan dalam literatur dapat diperoleh pada Pohlmann (1991). Selanjutnya, persamaan ini diselesaikan secara numerik menggunakan metode beda hingga dan memakai Arakawa C-grid (Arakawa and Lamb, 1977). 3. METODOLOGI Pada penelitian ini dimasukkan lima komponen utama pasang surut laut yaitu komponen pasang surut M2 , S2 , N2 , K1 dan O1 , terdiri atas amplitudo dan fase. Pada simulasi ini, selain data syarat batas pasang surut berupa nilai amplitudo dan fase dari (Zahel, Gavio and Seiler, 2000). Juga dimasukkan data topografi yang bersumber dari http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO5/TOPO/. Disamping itu juga dalam penelitian ini digunakan domain model wilayah terletak pada lintang astronomis 95o 30103o 30 BT dan 1o 30 LU 5o 30 LU, model didiskritisasi 5 dalam arah x dan y. Dalam arah vertikal, model terbagi menjadi 17 lapisan yaitu, 0-10 m, 10-20 m, 20-30 m, 30-50 m, 50-75 m, 75100 m, 100-125 m, 125-150 m, 150-200 m, 200-250 m, 250-300 m, 300-400 m, 400-500 m, 500-600 m, 600-700 m, 700-800 m, 800-900 m, juga digunakan nilai AH = 500m2 /s, Av = 0, 001m2 /s, e = 1e − 5, ∆t = 300, 0s, α = 0, 5.

4. DISKUSI Untuk mendapatkan gambar pasang surut laut kita harus menyimpan data ς(x, y, t), dimana (x, y) adalah posisi dan t adalah waktu. Untuk penelitian ini dipilih 3 lokasi yaitu Lhokseumawe dengan koordinat


62

(97o 9BT, 5o 12LU ), Belawan dengan koordinat (98o 42BT, 3o 48LU ), dan Bagan siapi-api dengan koordinat (100o 36BT, 2o 12LU ), dan t dipilih pada bulan maret dan oktober 2008. Pada tanggal 9 maret 2008 bertepatan dengan tanggal 1 rabiul awal 1429 H yang merupakan bulan baru, tanggal 16 maret 2008 bertepatan dengan tanggal 8 rabiul awal 1429 H yang merupakan quarter pertama bulan, tanggal 23 maret 2008 bertepatan dengan tanggal 15 rabiul awal 1429 H yang merupakan bulan purnama dan 31 maret 2008 bertepatan dengan tanggal 23 rabiul awal 1429 H yang merupakan quater ketiga. Dari gambar 1 samapai dengan gambar 3 terlihat bahwa ketinggian pasang surut mencapai puncak pada bulan baru. Selanjutnya ketinggian pasang surut juga sangat tinggi pada saat bulan purnama. Sedangkan pada quarter pertama dan quarter ketiga pasang surut mengalami penurunan dibanding pada saat bulan purnama dan bulan baru.

Gambar 1: Pasang surut bulat maret 2008 di Lhokseumawe.

Gambar 2: Pasang surut bulat maret 2008 di Belawan. Pada bulan oktober 2008 tanggal 1 oktober 2008 bertepatan dengan tanggal 1 syawal 1429 H dimana posisi bulan ialah bulan baru, tanggal 8


63

Gambar 3: Pasang surut bulat maret 2008 di Bagan siapi-api. oktober 2008 bertepatan dengan tanggal 8 syawal 1429 H dimana posisi bulan quarter pertama, tanggal 15 oktober 2008 bertepatan dengan tanggal 15 syawal 1429 H dimana posisi bulan ialah bulan purnama, dan tanggal 23 oktober 2008 bertepatan dengan tanggal 23 syawal 1429 H dimana posisi bulan quarter ketiga. Dari gambar 4 samapai dengan gambar 6 terlihat bahwa ketinggian pasang surut mencapai puncak pada bulan purnam. Selanjutnya ketinggian pasang surut juga sangat tinggi pada saat bulan baru. Sedangkan pada quarter pertama dan quarter ketiga pasang surut mengalami penurunan dibanding pada saat bulan purnama dan bulan baru.

Gambar 4: Pasang surut bulat oktober 2008 di Lhokseumawe. Dari gambar 1 sampai dengan gambar 6 juga terlihat bahwa ketinggian pasang surut di Bagan siapi-api lebih tinggi dibandingkan dengan ketinggian pasang surut di Belawan. Begitu juga ketinggian pasang surut di Belawan lebih tinggi dibandingkan dengan ketinggian pasang surut di Lhokseumawe. Hal ini diakibatkan di dekat Bagan siapi-api mempunyai selat sempit sehingga gradien ketebalan lapisan pertama di Bagan siapi-api lebih


64

Gambar 5: Pasang surut bulat oktober 2008 di Belawan.

Gambar 6: Pasang surut bulat oktober 2008 di Bagan siapi-api. besar dibandingkan dengan gradien ketebalan lapisan pertama di Belawan. Hal yang sama juga terjadi antara di Belawan dengan Lhokseumawe, dimana gradien ketebalan lapisan pertama di Belawan lebih besar dibandingkan dengan gradien ketebalan lapisan pertama di Lhokseumawe. Hal ini berakibat laju vertikal lapisan pertama di Bagan siapi-api lebih besar dibandingkan dengan laju vertikal lapisan pertama di Belawan. Begitu juga laju vertikal lapisan pertama di Belawan lebih besar dibandingkan dengan laju vertikal lapisan pertama di Lhokseumawe. Hal ini mengakibatkan ketinggian pasang surut di Bagan siapi-api kebih tinggi dibandingkan dengan ketinggian pasang surut di Belawan. Begitu juga ketinggian pasang surut di Belawan lebih tinggi dibandingkan dengan ketinggian pasang surut di Lhokseumawe. Disamping itu juga terlihat dari gambar 1 sampai 6 bahwa di daerah selat Malaka mempunyai tipe pasang surut semidiurnal, dimana dalam satu hari, mempunyai dua kali pasang naik dan dua kali pasang surut. Hal ini sesuai dengan penelitian yang pernah dilakukan oleh Wyrtki seperti terlihat


65

pada gambar 7, dimana tipe pasang surut di selat Malaka ialah semidiurnal

Gambar 7: Tipe pasang surut laut di Asia Tenggara berdasarkan Wyrki (1961) Untuk menggambarkan arus maka harus didapatkan nilai u(x, y, z, t) dan v(x, y, z, t) pada semua nilai x dan y. Selanjutnya untuk waktu t kita sesuaikan dengan bulan dan tahun yang akan kita gambarkan. Karena ∆t = 300, 0s, maka kita harus mencari nilai rata-rata dari u dan v pada bulan tersebut. Misalkan u ¯ dan v¯ adalah rata-rata kecepatan pada bulan tersebut. Sehingga kita bisa menghitung rata-rata besar kecepatan dan arah kecepatan pada bulan tersebut dengan menggunakan persamaan p K = u2 + v 2 (40) −1 v (41) ϕ = tan u Dimana V adalah besar kecepatan dan ϕ adalah arah kecepatan. Selanjutnya untuk arus permukaan kita ambil nilai z = 1, sehingga kita mempunyai besar kecepatan dan arah dari kecepatan arus pada permukaan. Gambar dari arus di lapisan permukaan di selat Malaka hasil dari simulasi yang penulis lakukan pada bulan maret 2008 dan bulan oktober 2008 dapat dilihat pada gambar 8 dan gambar 9. Untuk arus di selat Malaka baik pada southeast monsoon maupun northwest monsoon secara umum bergerak ke arah barat laut, hal ini diakibatkan perbedaan ketinggian pasang surut laut di daerah tenggara mempunyai ketinggian pasang surut yang lebih tinggi di bandingkan dengan di daerah barat laut. Perbedaannya ketika south east monsoon arus laut di selat Malaka berasal dari laut Cina Selatan, sedangkan ketika north west


66

Gambar 8: Arus di lapisan permukaan pada bulan maret 2008

Gambar 9: Arus di lapisan permukaan pada bulan oktober 2008 monsoon arus di selat Malakan berasal dari laut Jawa. Hal ini dapat dilihat dari hasil simulasi pada gambar 8 dan gambar 9 dan beberapa penelitian sebelumnnya yang pernah dilakukan di selat Malaka.

Gambar 10: Arus laut di selat Malaka pada saat northeast monsoon dan southwest monsoon berdasarkan kantor navigasi Amerika Serikat (1944) dalam Keller dan Richard (1967). Dari gambar 10 terlihat bahwa arus secara umum bergerak ke arah barat laut, baik pada saat southwest monsoon maupun northeast monsoon. Terlihat juga bahwa arus di selat Malaka ketika northeast monsoon berasal dari laut Cina Selatan dan ketika southwest monsoon berasal dari laut Jawa Hal yang sama juga terlihat dari hasil penelitian yang dilakukan oleh Wyrtki (1961), gambar 11 dan gambar 12 menunjukkan sirkulasi arus laut pada bulan juni dan oktober. Untuk di selat Malaka terlihat bahwa arus laut secara umum bergerak menuju ke arah barat laut. Demikian juga halnya dengan penelitian yang dilakukan oleh Hennesey


67

Gambar 11: Sirkulasi arus laut pada bulan Juni berdasarkan Wyrtki (1961).

Gambar 12: Sirkulasi arus laut pada bulan Oktober berdasarkan Wyrtki (1961).

Gambar 13: Sirkulasi arus laut pada bulan Februari berdasarkan Hennesey (1971).


68

Gambar 14: Sirkulasi arus laut pada bulan Agustus berdasarkan Hennesey (1971). (1971) seperti terlihat pada gambar 13 dan 14, menunjukkan arus di permukaan sekat Malaka pada saat northeast monsoon dan south west monsoon secara umum bergerak ke arah barat laut. Terlihat juga pada saat northeast monsson arus di selat Malaka berasal dari laut Cina Selatan, sedangkan pada saat southwest monsoon berasal dari laut Jawa. Begitu juga dengan simulasi yang dilakukan oleh Thomas Pohlmann (1985) menunjukkan arus di permukaan sekat Malaka pada saat northeast monsoon dan south west monsoon secara umum bergerak ke arah barat laut. Terlihat juga pada saat northeast monsson arus di selat Malaka berasal dari laut Cina Selatan, sedangkan pada saat southwest monsoon berasal dari laut Jawa. Simulasi yang dilakukan oleh Thomas Pholmann tersebut dapat dilihat pada gambar 15 dan gambar 16.

KESIMPULAN Berdasarkan hasil penelitian yang telah dilakukan maka dapat disimpulkan bahwa 1. Puncak ketinggian pasang surut terjadi ketika bulan baru atau bulan purnama. 2. Pasang surut laut di daerah tenggara selat Malaka mempunyai ketinggian pasang surut laut yang lebih tinggi di bandingkan dengan pasang surut laut di daerah barat laut selat Malaka.


69

Gambar 15: Simulasi arus laut pada saat northeast monsoon berdasarkan Pohlmann (1985).

Gambar 16: Simulasi arus laut pada saat southwest monsoon berdasarkan Pohlmann (1985). 3. Tipe pasang surut yang mendominasi di selat Malaka ialah tipe semi diurnal, dimana dalam satu hari terjadi dua kali pasang naik dan dua kali pasang surut. 4. Secara umum arus di selat Malaka bergerak dari arah tenggara selat Malaka menuju ke arah barat laut selat Malaka.

REFERENSI 1. Arakawa, A. and V.R. Lamb (1977), Computational design of the basic


2.

3.

4. 5. 6.

7.

8.

9.

10.

11.

12.

13.

70

dynamic processes of the UCLA general circulation model. Methods in Computational Physics 17, 173-263, 1977 Backhaus, J.O. (1983), A Semi-implicit Scheme for the Shallow Water for Application to the Shelf Sea Modeling. Continental Shelf Research, 2, 243-254. Backhaus, J.O. (1985), A Three-Dimensional Model for the Simulation of Shelf Sea Dynamics. Deutsche Hydrographische Zeitschrifft, 38(4), 165 187. Keller, G.H. and A.F. Richards (1967), Sediments of the Malacca Strait, Southeast Asia, Journal Sedimentary Petrology, 102-127 Hennesey, S.J. (1971), Malacca Strait and West Coast of Sumatra Pilot, Hydrographer of the Navy, England. Huang, D. (1995) Modeling Studies of Barotropic and Baroclinic Dynamics in the Bohai Sea. Ph.D Thesis Institut fuer Meereskunde, Hamburg University. Mihardja, D.K. (1991), Energy and Momentum Budget of the Tides in Indonesian Waters, Ph.D Thesis Institut fuer Meereskunde, Hamburg University. Pohlmann, T. (1985), Simulation von Bewegungsvorgngen im Sdchinesischen Meer, Dimplomarbeit Institut fuer Meereskunde, Hamburg University. Pohlmann, T. (1991), Untersuchung Hydro- und thermodynamischer Prozesse in der Nordsee mit einemdreidimensionalem numerischen in der Nordsee mit einemdreidimensionalem numerischen Modell. Ph.D Thesis Institut fuer Meereskunde, Hamburg University. Pohlmann, T. (1996), Predicting the Thermocline in a Circulation Model of North Sea Part I : Model Description, Calibration and Verivication. Continental Shelf Research, 16(2), 131-146. Pond, S. and G.L. Pickard. (1983), Introductory Dynamical Oceanography, Departement of Oceanography, University of Columbia, Vancouver, Canada, second edition. Putri, M.R. (2005), Study of Ocean Climate Variability (1959-2002) in the Eastern Indian Ocean, Java Sea and Sunda Strait using the HAMburg Shelf Ocean Model, Ph.D Thesis Institut fuer Meereskunde, Hamburg University. Rizal,S. (1994), Numerical Study on the Malacca Strait (Southeast Asia) with a Three-Dimensional Hydro dynamical Model, Ph.D Thesis Institut fuer Meereskunde, Hamburg University.


71

14. Thia-Eng, C., I.R.L. Gorre, S.A. Ross and S. Regina. (2000), The Malacca

Straits. Marine Pollution Bulletin, 41, 160-177. 15. Wyrtki, K. (1961), Scientific Results of Marine Investigations of the South China Sea and the Gulf of Thailand 1959-1961, Naga Report Volume 2, The University of California, Scripps Institutions of Oceanography, La Jolla, California. 16. Zahel, W., J.H. Gavio and U. Seiler (2000), Angular Momentum and Energy Budget of a Global Ocean Tide Model with Data Assimilation, GEOS, Ensenada, 20(4), 400-413.

Taufiq Iskandar: Jurusan Matematika, FMIPA, Universitas Syiah Kuala, Jl. Tgk.

Daud Beureuh No. 4 Darussalam-Banda Aceh 23111, Indonesia.



72

Suwilo – Vertex Exponent of Two-Cycles

73


VERTEX EXPONENT OF TWO-CYCLES

Saib Suwilo and Masna Meisaroh Abstract. A digraph D is primitive provided there is a positive integer ` such that for each pair of vertices u and v in D there exists a walk from u to v of length `. The smallest of such integer ` is the exponent of D. The exponent of a vertex v in D is defined to be the smallest positive integer ` such that there is a walk of length ` from v to any vertices in D. In this paper we discuss the vertex exponent of primitive digraph consisting of two cycles of length k and j. We show that for any vertex v in D the exp(v) = exp(v1 ) + d(v, v1 ) where v1 is the vertex in D with outdegree 2 and exp(v1 ) = n − k + 1 + j(k − 2).

1. INTRODUCTION We begin with digraph concepts from [2] that lead to discussion of vertex exponents. Let D be a digraph on n vertices. The outdegree of a vertex v is the number of arcs starting at v. A walk from u to v of length m is a sequence of m arcs of the form (u = v0 , v1 ), (v1 , v2 ), . . . , (vm−1 , vm = v). The notation v0 → v1 → v2 → vm−1 → vm is used for a walk from v0 to vm . In general, the notation uv-walk is used for a walk from u to v. A uv-path is a uv-walk with all vertices distinct except possibly u = v. A uv-path is open whenever u 6= v and is closed otherwise. A cycle is a closed path. The distance from vertex u to vertex v in D, denoted d(u, v), is the length of the shortest uv-path. A digraph D is strongly connected provided that for Received 12-10-2008, Accepted 14-11-2008. 2000 Mathematics Subject Classification: 05C20, 05C40 Key words and Phrases: primitive, vertex exponent, two cycles


74

each pair of vertices u and v in D there exists a uv-walk. By a two-cycles we mean a strongly connected digraphs consisting of two cycles. A strongly connected digraph D is primitive if there exists a positive integer ` such that for each pair of vertices u and v there is a walk from u to v of length `. It is well known that a strongly connected digraph is primitive if and only if the greatest common divisor of cycle lengths in D is 1 [2]. The smallest integer ` such that for each pair of vertices D there exists a walk of length ` connecting them is the exponent of D. The exponent of a digraph D is denoted by exp(D). The exponent of a vertex v in D is the smallest positive integer ` such that there is a walk of length ` from v to every vertex in D. The exponent of a vertex v is denoted by exp(v). We note that for any vertex v in D, exp(v) ≤ exp(D). Moreover, one can show that exp(D) = maxv∈V (D) {exp(v)}, see [2] for a proof. Research on exponent of digraphs is initiated by Wielandt in 1950, see [8]. Wielandt shows that for a primitive digraph D on n vertices exp(D) ≤ n2 − 2n + 2. Brualdi and Liu [1] introduce the notion of vertex exponent of digraphs and discuss an application in memoryless communication network. Liu [5], Shao and Li [6] discuss vertex exponents of symmetric primitive digraphs. Liu especially discusses a formula for vertex exponent of lollipop with cycle length 3 only. In this paper we discuss vertex exponent of primitive directed two-cycles. In Section 2, we discuss some results regarding primitive two-cycles. In Section 3, we give explicit formula for exponent of any vertex v in a two-cycles in term of the vertex exponent of a distinguished vertex that has outdegree 2. 2. PRELIMINARIES We consider primitive two-cycles D on n vertices and cycle lengths k and j as in Figure 1. We note that n = k + j − t where t ≤ j < k is the number of vertices belong to both cycles of length k and j. Let k and j be relatively prime integers. We define φ(k, j) to be the smallest positive integer such that for each integer m ≥ φ(k, j) there are nonnegative integers a and b with t = ak + bj. Lemma 4.3 Let k and j be two relatively prime positive integers. Then φ(k, j) = (k − 1)(j − 1) = kj − k − j + 1. Proof. We first show that φ(k, j) ≥ kj − k − j + 1. Suppose on the contrary that there are nonnegative integers a and b such that kj − k − j = ak + bj. Hence k(j − 1) = ak + (b + 1)j. Since k and j are relatively prime, we have


75

Figure 17: A two-cycles with cycle lengths k and j. k must divide b + 1. Similarly, j must divide a + 1. This implies a ≥ j − 1 and b ≥ k − 1. Now we have kj − k − j = ak + bj ≥ (j − 1)k + (k − 1)j = kj − k − j + kj, a contradiction. Therefore, it must be the case that φ(k, j) ≥ (k − 1)(j − 1). We show that every positive integer m > kj can be expressed as rk+sj where r and s are positive integers. This implies that every positive integer m > kj − k − j can be expressed as a nonnegative integer linear combination of k and j. Since m > kj, there exists a positive integer r with 1 ≤ r ≤ j such that m ≡ rk mod j. We define s to be s = (m − rk)/j. Then s is a positive integer and we have that m = rk + sj. Let D be a primitive digraph containing cycles of two different lengths k and j, and let u and v be any two vertices in D. We define a number ruv to be 0 if u = v and if u is on a cycle of length k and a cycle of length j, otherwise ruv is the length of the shortest uv-walk that passes through cycle of each length in D. The number r = max{ruv } over all pairs of vertices u and v in D is the circumdiameter of D. An ordered pair of vertices (u, v) in a primitive digraph D has the unique path property if every uv-walk that has length at least ruv consists of some walk w of length ruv augmented by a number of cycles which has a vertex in common with the walk w. The following lemma is a special case of a result in Dulmage and Mendelsohn [4], see Corollary 2 of Dulmage and Mendelsohn. Lemma 4.4 Let D be a primitive two-cycles on n vertices with lengths k and j. Let x and y be vertices in D such that the ordered pair (x,y) has the unique path property and rxy = r, then exp(D) = φ(k, j) + r. Proof. We first show that exp(D) ≤ φ(k, j) + r. Let u and v be a pair of vertices in D. For any set of non-negative integers a and b there exists a


76

walk from vertex u to vertex v of length ak + bj + ruv . This implies for any pair of vertices u and v in D, there is a walk of length φ(k, j) + ruv + M for every positive integer M . Taking M = r − ruv there exists a walk from u to v of length φ(k, j) + r. Thus, exp(D) ≤ φ(k, j) + r. Notices that by definition exp(D) is the smallest positif integer such that for each pair of vertices u and v in D there is a walk from u to v of length k, for all integers k ≥ exp(D). Considering closed walk from a vertex u to itself, then exp(D) ≥ φ(k, j). Let x and y be vertices in D such that the pair (x, y) has a unique path property. Then there is no walk from x to y of length φ(k, j) − 1 + r, since the existence of such walk would imply that there is nonnegative integers a1 and a2 such that φ(k, j) = a1 k + a2 j. Hence exp(D) ≥ φ(k, j) + r. As a direct consequence of Lemma 4.4 we have a formula for exponent of primi tive two-cyles on n vertices with cycle lengths k and j (see also [3, 7].) Theorem 4.4 Let D be a primitive two-cycles on n vertices and with cycle lengths k and j. Then exp(D) = n + j(k − 2). Proof. Notice that the pair of vertices (t + 1, k) has a unique path property with rt+1,k = r = n + k − j − 1. Hence by Lemma 4.4 we have exp(D) = φ(k, j) + n + k − j − 1 = n + j(k − 2). 3. THE VERTEX EXPONENT Let D be a primitive digraph consisting of two cycles C1 and C2 of relatively prime length k and j with k > j. Suppose t is the number of vertices belong to both cycles C1 and C2 . Then t = k + j − n. Let v1 = t be the vertex in D with outdegree 2. Lemma 4.5 Let D be a primitive two cycles of length k and j. Let v1 be the vertex in D with outdegree 2. Then exp(v1 ) = n − k + 1 + j(k − 2). Proof. By definition we have that exp(v1 ) ≥ φ(k, j). Let vn be the vertex in D such that d(v1 , vn ) = k − t and let wv1 vn be the shortest walk from v1 to vn with `(wv1 vn ) ≥ φ(k, j). Hence `(v1 vn ) = φ(k, j) + k − t. This implies exp(v1 ) ≥ φ(k, j) + k − t. We show that exp(v1 ) ≤ φ(k, j)+k −t by showing that for any vertices v in D there is a walk from v1 to v of length φ(k, j) + k − t. We first note


77

that by definition of φ(k, j) that for each nonnegative integer k − t there are nonnegative integers a and b such that φ(k, j) + k − t = ak + bj. This implies there is a closed walk from v1 to v1 of length exactly φ(k, j) + k − t. Namely, the walk that starts at v1 , moves a times around C1 and backs to v1 and then moves b times around C2 and ends at v1 . Let v be any vertex in D where v 6= v1 . Then d(v1 , v) ≤ k − t. This implies φ(k, j) + (k − t) − d(v1 , v) ≥ φ(k, j), and hence there are nonnegative integers r and s such that φ(k, j) + (k − t) − d(v1 , v) = rk + sj. Notice that the walk that starts at v1 , moves r times around C1 , then moves s times around C2 and finally moves to v along the v1 v-path is a walk from v1 to v of length φ(k, j) + k − t. Therefore, exp(v1 ) ≤ φ(k, j) + k − t. We now conclude that exp(v1 ) = φ(k, j) + k − t = (k − 1)(j − 1) + k − (k + j − n) = n − k + 1 + j(k − 2). The following theorem gives a formula for vertex exponent of any vertex v in a two-cycles in term of the vertex exponent of v1 . Theorem 4.5 Let D be a primitive two cycles and let v1 be the vertex with outdegree 2 in D. Then For any vertex v in D we have exp(v) = n − k + 1 + j(k − 2) + d(v, v1 ). Proof. We first show that exp(v) ≤ n − k + 1 + j(k − 2). For any vertex w, the walk that starts at v follows the vv1 -path of length d(v, v1 ) to v1 , and then follows the walk from v1 to w of length n − k + 1 + j(k − 2) created in Lemma 4.5 is a vw-walk of length n − k + 1 + j(k − 2) + d(v, w). We note that exp(v) ≥ φ(k, j). Notice that the vvn -path can be decomposed into vv1 -path of length d(v, v1 ) and v1 vn -path of length k − t. This implies the shortest vvn -walk with length at least φ(k, j) is of length φ(k, j) + k − t + d(v, v1 ). Hence exp(v) ≥ n − k + 1 + j(k − 2) + d(v, v1 ). We note that in a two-cycles D of length k and j, we have d(v1 , v) ≤ k − 1. This implies for any vertex v in D, exp(v) ≤ n + j(k − 2) = exp(D).

REFERENCES 1. R.A. Brualdi and B. Liu. Generalized exponents of primitive directed

graphs. J. Graph Theory 14(4)(1990),483-499 . 2. R.A. Brualdi and H.J. Ryser. Combinatorial matrix theory. Cambridge University Press, Cambridge (1991).


78

3. L.F. Dame. The exponent and circumdiameter of primitive directed graphs.

4. 5. 6. 7. 8.

M.Sc. Thesis. Departement Mathematics and Statistics, University of Victoria, Victoria BC, Canada (2004). A.L. Dulmage, N.S. Mendelsohn, Gaps in the exponent set of primitive matrices, Illinois J. Math. 8 (1964) 642–656. B. Liu. Generalized exponents of primitive simple graphs. Acta Math. Appl. Sin, 9:1 (1993), 36-43. J. Shao and B. Li. The set of generalized exponents of primitive simple graphs. Linear Algebra Appl. 258, (1997) 95-127. S. Suwilo and N. Herawati. Eksponen dari dua cycle. Sain dan Teknologi. 8:3 (2002), 123–132. H. Wielandt. Unzerlegbare nicht negative Matrizen. Math. Zeit 52 (1950), 642–648.

Saib Suwilo

Department of Mathematics, Faculty of Mathematics and Natural Sciences, University of Sumatera Utara, Indonesia E-mail: [email protected] Masna Meisaroh

Department of Mathematics, Faculty of Mathematics and Natural Sciences, University of Sumatera Utara, Indonesia E-mail: [email protected]

Tulus – Metoda Elemen Batas

79

Vol. 01, No. 01 (2009), pp. 79–99.

METODE ELEMEN BATAS UNTUK MASALAH POTENSIAL DENGAN HAMPIRAN ELEMEN KUBIK

Tulus Abstract. Boundary element method is a method for approximation solution of boundary value problems. One of problems that can be good approximated by this method is potential problem in which potential and flux are considered. In this paper, cubic element approximation for potential problems. The results are presented in the form of program scheme. The method used to determine the values at boundaries is collocation points method. The integral is calculated using quadrature Gauss formula. The results are compared with exact values.

1. PENDAHULUAN Metoda elemen batas adalah suatu metoda untuk penyelesaian hampiran dari masalah nilai batas dua atau tiga dimensional [5] atau juga untuk dimensi yang lebih tinggi. Metoda ini mengandalkan pembagian masalah tersebut ke dalam bentuk persamaan integral batas. Batas dipilah menjadi potongan-potongan (kita sebut ’elemen’) yang masing-masingnya didekati oleh suatu suku banyak derajat rendah. Penyelesaian hampiran untuk masalah semula dapat diperoleh dengan mengambil suatu integrasi hampiran yang sesuai. Dalam hal konvergensi, Hayakawa dan Iso [4] telah menunjukkan konvergensi seragam orde tinggi dari skema elemen batas ini. Dalam tulisan ini dikaji perihal elemen kubik secara numerik untuk masalah potensial. Tujuan akhirnya adalah pembuatan skema program Received 25-03-2008, Accepted 12-12-2008. 2000 Mathematics Subject Classification: 65Q38, 65R20 Kata kunci: hampiran kubik, metode elemen batas.

80


komputer untuk penyelesaian numerik masalah potensial menggunakan metoda elemen batas dengan elemen kubik. Bahasa pemrograman yang dipakai adalah FORTRAN. Dengan penurunan yang sama seperti dalam penurunan elemen kuadratik, hasil dari elemen kubik akan dibandingkan dengan hasil elemen kuadratik yang telah ada. Pembandingan dilakukan terhadap beberapa bentuk batas dari masalah potensial. Dalam tulisan ini pembandingan masih dibatasi terhadap masalah batas yang telah diketahui penyelesaian eksaknya, yaitu masalah hantaran panas pada suatu lempeng berbentuk empat persegi panjang yang terdapat pada [3] dan masalah hantaran panas pada suatu penampang silinder yang terdapat pada [2].

2. PERSAMAAN INTEGRAL DASAR Penurunan dari persamaan integral dasar ada dijelaskan pada beberapa kepustakaan, seperti halnya dalam [1] dan [6]. Di sini kita akan menurunkan persamaan menggunakan teorema Green. Misalkan Ω suatu daerah buka di R2 yang terhubung dan terbatas dengan batas Γ, dengan Γ = Γ1 + Γ2 seperti gambar berikut. Γ2

Ω Γ1

Gambar 1 Daerah Ω dengan batas Γ Pandang fungsi u ∈ C 2 yang didefinisikan pada Ω + Γ dan memenuhi syarat ∇2 u = 0 di Ω (42) serta kondisi pada batas u = u di Γ1

dan q =

∂u = q di Γ2 ∂n

(43)

81


dimana ∇2 adalah operator Laplace dan n adalah normal terhadap batas. Asumsikan sekarang bahwa u merupakan nilai hampiran. Akibatnya akan timbul galat. Galat yang ditimbulkan dalam persamaan di atas, bilamana nilai eksak dari u dan q merupakan penyelesaian hampiran, dapat diminimalkan dengan ortogonalisasi terhadap fungsi bobot u∗ , dengan tu∗ runan pada batas q ∗ = ∂u ∂n . Dengan mengalikan persamaan (42) dengan u∗ dan kemudian mengintegralkan pada Ω diperoleh Z (∇2 u)u∗ dΩ = 0 (44) Ω

Dengan menggunakan teorema Green dua kali dan persyaratan bahwa u∗ ∈ C 2 dan ∇2 u∗ = 0 kita punya Z Z Z ∗ 2 ∗ (45) uq dΓ − qu∗ dΓ (∇ u )u dΩ = Γ

Γ

Ω

Persamaan ini merupakan titik awal dari pemakaian metoda elemen batas. Pada tulisan ini u∗ akan diambil yang mrupakan penyelesaian fundamental. 2.1. Penyelesaian Fundamental Penyelesaian fundamental u∗ memenuhi persamaan Laplace dan merepresentasikan lapangan yang dibangun oleh muatan satuan terkonsentrasi yang beraksi pada suatu titik i. Penyelesaian tersebut dapat ditulis dengan ∇2 u ∗ + ∆ i = 0

(46)

dimana ∆i menyatakan fungsi Dirac delta yang menuju tak hingga pada titik x = xi dan sama dengan nol untuk lainnya. Integral dari ∆i sama dengan satu. Integral dari fungsi Dirac delta yang dikalikan dengan sebarang fungsi lain sama dengan fungsi tersebut pada titik xi , sehingga Z Z (∇2 u∗ )u dΩ = (−∆i )u dΩ = −ui (47) Ω

Ω

Dari persamaan (45) diperoleh Z Z i ∗ u + uq dΓ = qu∗ dΓ Γ

(48)

Γ

Dari persamaan ini dapat ditentukan nilai ui untuk suatu xi di dalam Ω bilamana nilai-nilai u dan q diketahui pada batas.

82


Untuk daerah isotropik dua dimensional nilai u∗ adalah sebagai berikut. 1 1 u∗ = ln (49) 2π r dimana r adalah jarak dari titik xi dari ke sebarang titik yang sedang dipertimbangkan. Persamaan (49) memenuhi persamaan Laplace dua dimensional. Sekarang akan dicari nilai u dan q pada batas dengan membawa xi menuju batas. Dengan memperhatikan kondisi dari persamaan (48) cukup dicari nilai u di Γ2 dan q di Γ1 . Karena persamaan (48) berlaku di daerah Ω, maka Ω diperbesar menjadi Ω + Ω . Akibatnya diperoleh batas Γ− + Γ . Karena sekarang xi telah berada di dalam Ω + Ω , maka persamaan (48) dapat digunakan. Ambil → 0. Maka Ω + Ω → Ω, Γ− + Γ → Γ dan Γ− → Γ. Persamaan (48) sekarang dapat ditulis menjadi Z Z Z Z i ∗ ∗ ∗ u + uq dΓ + uq dΓ = qu dΓ + qu∗ dΓ (50) Γ−

Γ

Γ−

Γ

Dari persamaan (50) dengan mengambil → 0 diperoleh Z Z ci ui + uq ∗ dΓ = qu∗ dΓ Γ

(51)

Γ

3. INTEGRASI NUMERIK DALAM ELEMEN KUBIK Ekspresi (51) dapat didiskritkan untuk mendapatkan sistem persamaan yang darinya nilai batas dapat dicari. Setelah pendiskritan batas ke dalam N elemen, persamaan (51) dapat ditulis menjadi ci ui +

N Z X j=1

Γj

uq ∗ dΓ =

N Z X j=1

qu∗ dΓ

(52)

Γj

Titik-titik dimana nilai-nilai tak diketahui diperhatikan disebut ’simpul’. 3.1. Variasi potensial dan fungsi interpolasi dalam elemen kubik Ada beberapa bentuk elemen yang bergantung pada cara penempatan dari titik-titik simpul di atas. Bentuk-bentuk tersebut antara lain metoda elemen

83


batas yang disebut elemen konstan, elemen linier, elemen kuadratik, elemen kubik, dan elemen dengan derajat lebih tinggi. Pada elemen kubik setiap elemen didekati oleh suatu kurva kubik, sehingga diperlukan empat titik simpul, dua pada ujung-ujung elemen dan dua di dalam elemen. Untuk setiap elemen ke-j Γj keempat titik ini masing-masing ditandai dengan (1), (2), (3), dan (4). Dan untuk problema dimensi dua setiap titik simpul ke-k, k = 1, 2, 3, 4 ditandai dengan (xk , y k ). Perhatikan gambar 2 di bawah. Pada elemen kubik nilai u dan q bervariasi secara kubik atas masingmasing elemen Γj . Nilai u (dan analog untuk q) pada sebarang titik pada elemen Γj dapat didefinisikan dalam suku-suku dari titik simpul mereka dan empat fungsi interpolasi kubik, sebut φ1 , φ2 , φ3 , dan φ4 , yang diberikan dalam suku-suku dari koordinat homogen ξ sebagai berikut. Elemen Simpul ) j u u= ?B 6 uB BN

u

u u j

j u

u

Ω u

u u j

u

u

u

j u

u

Ω

u

(1)

u

(4) u

(3)

(2)

Gambar 2 Elemen Kubik

84


u(ξ) = φ1 u1 + φ2 u2 + φ3 u3 + φ4 u4 =

 1    u    u2  φ4 u3     4   u

φ1 φ2 φ3

(53)

dimana ξ kordinat tak berdimensi yang bervariasi dari −1 sampai 1. Keempat fungsi interpolasi tersebut adalah φ1 =

1 2 16 (9ξ

φ2 =

9 16 (3ξ

− 1)(1 − ξ),

φ3 =

9 16 (3ξ

φ4 =

1 2 16 (9ξ

+ 1)(1 − ξ 2 ), (54)

−

1)(ξ 2

− 1),

− 1)(1 + ξ),

3.2. Matriks Integral Sekarang kita pandang integral atas elemen ke-j. Dengan memperhatikan persamaan (53), maka integral pada sisi kiri persamaan (52) dapat ditulis menjadi  1  u     2   u uq ∗ dΓ = [ φ1 φ2 φ3 φ4 q ∗ dΓ  u3   Γj Γj   4  u  1   u    2  o n u ij ij ij ij = h1 h 2 h3 h4 u3      4  u

Z

Z

dimana untuk setiap elemen ke-j kita punya keempat bentuk Z ij hk = φk q ∗ dΓ, k = 1, 2, 3, 4

(55)

(56)

Γj

Dengan cara sama, integral pada sisi kanan persamaan (52) dapat ditulis


85

menjadi  1  q     2   q qu∗ dΓ = [ φ1 φ2 φ3 φ4 u∗ dΓ 3  Γj Γj  q     q4  1    q2    n o q ij ij ij ij = g1 g2 g3 g4 q3     4   q

Z

Z

dimana gkij

Z =

φk u∗ dΓ,

k = 1, 2, 3, 4

(57)

(58)

Γj

Dengan mensubstitusi persamaan (??) (dan hal yang sama untuk q) untuk semua elemen j ke persamaan (52) kita punya persamaan untuk simpul xi sebagai berikut.

˜ i1 H ˜ i2 ci ui + [ H

=

Gi1 Gi2

 1  u        u2  i4N ˜ ... H ..  .      4N   u  1  q       q2   . . . Gi4N ..  .      4N   q

(59)

Karena pada setiap simpul nilai potensial sebelum dan sesudah simpul adalah sama, ini berarti nilai potensial pada simpul terakhir suatu elemen sama dengan nilai potensial pada simpul pertama pada elemen berikutnya. Dengan demikian kita punya u1 = u4j dan u4j = u4j+1 , j = 1, . . . , N − 1

(60)

Hal yang berbeda terjadi pada flux, yaitu nilai flux sebelum dan sesudah suatu simpul bisa berbeda. Dari kenyataan ini persamaan (59) dapat ditulis

86


menjadi

˜ i1 H ˜ i2 ci ui + [ H

=

Gi1 Gi2

 1  u       u2   i3N ˜ ... H ..  .      3N   u  1  q     2    q  i4N ... G ..  .      4N   q

(61)

i(j+1)

ˆ i(3j+1) sama dengan hij dari elemen ke-j tambah h dimana H 4 1 elemen ke-j + 1. Persamaan (61) dapat ditulis menjadi i

ci u +

3N X

ˆ ij uj = H

4N X

j=1

Gij q j

dari

(62)

j=1

Karena nilai ui nol untuk i 6= j dan 1 untuk lainnya, maka kita dapat memisalkan ˆ ij , H jika i 6= j H ij = ˆ ij H + ci , jika i = j

(63)

sehingga persamaan (62) dapat ditulis menjadi 3N X j=1

ij j

H u =

4N X

Gij q j

(64)

j=1

Dalam bentuk matriks persamaan ini dapat ditulis menjadi HU=GQ.

(65)

dengan H matriks 3N × 3N , U vektor dengan ukuran 3N , G matriks 3N × 4N , dan Q vektor dengan ukuran 4N . Persamaan di atas merupakan sistem persamaan simultan. Penghitungan dari komponen-komponen matrik G dan H di atas dibahas pada sub bab di bawah. Pada bagian ini kita menganggap, bahwa semua komponen-komponen tersebut telah diketahui.

87


3.3. Perlakuan terhadap sistem persamaan akibat syarat batas Setelah seluruh komponen matriks G dan H diketahui masih perlu dilakukan pertukaran antar komponen dari kedua matriks tersebut, sehingga diperoleh bentuk persamaan AX = B (66) dengan A matriks 3N × 3N , B vektor dengan komponen yang diketahui, dan X vektor yang komponennya belum diketahui nilainya. Ini dikarenakan pada persamaan (65) komponen dari vektor U dan Q ada yang belum dan ada yang sudah diketahui. Dalam pengurutan kembali di atas perlu diperhatikan empat kasus berkenaan dengan syarat batas untuk potensial dan flux pada simpul ujung elemen. Empat kasus tersebut adalah Kasus Nilai diketahui Nilai tak diketahui (a) flux(-), flux(+) potensial (b) flux(-), potensial flux(+) (c) flux(+), potensial flux(-) (d) potensial flux(-), flux(+) Keterangan: (-) = sebelum simpul, (+) setelah simpul Untuk kasus (a) kita tidak perlu mempertukarkan komponsen dari G dan H. Untuk kasus (b) dan (c), kolom matriks H yang sesuai dengan nilai yang diketahui hanya dipertukarkan dengan kolom matriks G sekawannya. Namun jika yang tak diketahui ada dua nilai flux, seperti pada kasus (d), maka kita perlu mengganti kolom matriks H yang sesuai dengan titik simpul tersebut dengan jumlah dua kolom matriks G yang juga bersesuaian dengan titik simpul tersebut. Kemudian pada matriks G kolom yang diganti oleh kolom matriks H tersebut hanya satu kolom, yaitu kolom yang pertama dari keduanya, sedangkan yang bersesuaian dengan nilai flux setelah titik simpul diganti dengan nol. 3.4. Penghitungan Matriks G dan H Untuk menentukan komponen-komponen matrik pada persamaan (65) kita perlu mendefinisikan sub matrik GW dan HW yang masing-masing mempunyai indeks 1, 2, 3 dan 4. Ini dikarenakan pengintegralan dilakukan terhadap batas dengan daerah integrasi elemen Γj . Untuk setiap elemen j kita punya Z Z φk q ∗ dΓ =

HW (k) = Γj

1

−1

φk q ∗ J dξ

(67)


Z GW (k) = Γj

φk u∗ dΓ =

Z

1

φk u∗ J dξ

88

(68)

−1

dengan k = 1, 2, 3, 4. Dalam penghitungan pengintegralan terdapat dua kasus utama, yaitu kasus nonsingular dan kasus singular. Kasus nonsingular terjadi bilamana titik sanding tidak berada pada keempat titik simpul pada elemen ke-j. Sedangkan kasus singular terjadi bilamana titik sanding berada pada salah satu titik simpul pada elemen ke-j. Kasus singular hanya mempengaruhi submatriks GW , karena kernel dari integralnya memuat faktor logaritma. Kasus Nonsingular Pandang submatriks yang isinya integral pada persamaan (67) dan (68). Apabila titik xi tidak berada pada elemen ke-j, maka nilai r tidak akan pernah nol. Dengan demikian kedua integral yang terdapat pada sisi kanan persamaan (67) dan (68) di atas dapat dipakai. Kasus ini yang disebut kasus non singular. Tetapi apabila titik xi merupakan salah satu titik yang terdapat pada elemen ke-j akan timbul nilai r = 0, sehingga kasus singular akan timbul. Kasus singular ini akan dibicarakan pada bagian berikutnya. Pada persamaan (67) harus juga dicari turunan ln 1r terhadap vektor normal. Untuk keperluan penghitungan secara numerik kita perlukan bentuk vektor normal dari vektor yang dibentuk oleh titik xi dengan titik yang sedang ditinjau dalam pengintegralan, sebut xp , yaitu ~n = (xi − xp )n~1 + (y i − y p )n~2

(69)

dimana n~1 dan n~2 dihitung dari komponen-komponen nilai Jacobi. Untuk elemen kubik, Jacobi dihitung dengan mengambil turunan dari ekspresi untuk koordinat x dan y, yaitu x = φ1 x1 + φ2 x2 + φ3 x3 + φ4 x4 y = φ1 y 1 + φ2 y 2 + φ3 y 3 + φ4 y 4

(70)

Dengan mensubstitusi bentuk φi , i = 1, 2, 3, 4 yang terdapat pada persamaan (54) ke persamaan (70), maka ekspresi yang terdapat pada (70)


89

dapat ditulis menjadi x =

9 27 27 9 − x1 + x2 + x3 + x4 ξ 3 16 16 16 16 9 9 9 9 + x1 − x2 − x3 + x4 ξ 2 16 16 16 16 1 27 2 27 3 1 1 + x − x − x − x4 ξ 16 16 16 16 9 2 9 3 1 4 1 1 − x + x + x − x 16 16 16 16

(71)

Penulisan untuk y analog. Kemudian dengan menurunkan persamaan (71) terhadap ξ, maka diperoleh bentuk Jacobi 1 hn 1 J = 27 −x + 3x2 + 3x3 + x4 ξ 2 + 18 x1 − x2 − x3 + x4 ξ 16 o2 n + x1 − 27x2 − 27x3 − x4 + 27 −y 1 + 3y 2 + 3y 3 + y 4 ξ 2 o2 i1/2 +18 y 1 − y 2 − y 3 + y 4 ξ + y 1 − 27y 2 − 27y 3 − y4 (72) Kasus Singular Jika titik sanding berada pada salah satu dari simpul yang terdapat pada elemen, maka kasus singular akan terjadi. Dengan demikian penghitungan integral dengan menggunakan metoda kuadratur Gauss baku tidak dapat dipakai. Untuk kasus ini dipakai metoda Gauss kuadratur khusus. Untuk menghitung integral pada bentuk (68) diperlukan transformasi variabel demikian sehingga titik sanding menjadi titik nol dan salah satu titik ujung elemen menjadi titik satu. Pada persoalan elemen kubik dipandang empat kasus yang berikut. Kasus 1: Titik sanding berada pada simpul (1) Pada kasus ini dilakukan transformasi koordinat, sehingga titik simpul (1) pada elemen Γj menjadi titik nol dan titik simpul (4) menjadi titik satu. Dengan perubahan koordinat dari variabel ξ ke variabel η, η = ξ+1 2 maka dengan memperhatikan bahwa Z Z 1 1 GW = φu∗ dΓ = φ ln J dξ (73) r Γj −1 diperoleh Z GW = 2 0

1

1 φ ln J dη r

(74)

90


Untuk mengganti r menjadi suatu fungsi dari η pandang sistem koordinat lokal berikut. Misalkan x2 = x2 − x1 ,

x3 = x3 − x1 ,

x4 = x4 − x1 ,

(75)

yaitu beda komponen x dari titik-titik simpul pada Γj terhadap titik simpul (1) pada elemen yang sama. Hal yang sama untuk komponen y. Pandang fungsi pangkat tiga r1 (η) = A1 η 3 + B1 η 2 + C1 η + D1

(76)

dengan r1 (0) = 0,

1 r1 ( ) = x2 , 3

2 r1 ( ) = x3 , 3

r1 (1) = x4 ,

(77)

Untuk jelasnya perhatikan gambar berikut. x4 6

x3

-s 6 (4)

-s 6 (3)

x2 -s 6 (2)

y4

y3

y2 s

-

(1) s

−1

ξ

s

s

− 13

1 3

s

s

0

1

η

s

s

s

1 3

2 3

1

Gambar 3 Titik Collocation pada simpul (1) Dengan menyelesaikan sistem persamaan yang dihasilkan dari sistem persamaan (77) diperoleh D1 = 0, dan 9 27 27 9 45 9 A1 = x4 − x3 + x2 , B1 = − x4 + 18x3 − x2 , C1 = x4 − x3 + 9x2 , 2 2 2 2 2 2 (78) Dengan demikian kita punya jarak rn n o o r(η) = η

A1 η 2 + B1 η + C1

2

+ A2 η 2 + B2 η + C2

2

(79)

91


dimana faktor berindeks 2 adalah ekspresi untuk y yang bersesuaian dengan x. Dengan mensubtitusi persamaan (79) ke (74), serta mengingat aturan dalam logaritma dan integral, maka diperoleh Z 1 1 GW = 2 φ ln J dη (80) η 0 rn Z o n o 1

−2

A1 η 2 + B1 η + C1

φ ln

2

+ A2 η 2 + B2 η + C2

2

J dη

0

Karena integral kedua tidak memuat faktor ln(1/r), maka penghitungan integral kedua ini dapat dilakukan dengan menggunakan metoda Gauss kuadratur baku seperti yang terdapat pada bagian sebelumnya dengan transformasi koordinat di atas. Dengan demikian, dari (80) kita punya Z 1 Z 1 1 hn ξ + 1 2 ξ + 1 o2 GW = 2 φ ln J dη − φ ln A1 + B1 + C1 η 2 2 0 −1 n ξ + 1 2 ξ + 1 o2 i 1 2 + A2 + B2 + C2 J dξ (81) 2 2 Z 1 1 J dη = 2 φ ln η 0 Z 1 1 rn o2 n o2 ˜1 ξ + C˜1 + A˜2 ξ 2 + B ˜2 ξ + C˜2 − φ ln A˜1 ξ 2 + B J dξ 8 −1 dengan A˜1 = 9x4 − 27x3 + 27x2 , A˜2 = 9y4 − 27y3 + 27y2 , ˜1 = 18x3 − 36x2 , ˜2 = 18y3 − 36y2 , B B C˜1 = −x4 + 9x3 + 9x2 , C˜2 = −y4 + 9y3 + 9y2 .

(82)

Bentuk φk , k = 1, 2, 3, 4 juga akan berubah sesuai dengan transformasi yang dipakai. Kasus 2: Titik sanding berada pada simpul (2) Berbeda dengan kasus pertama, pada kasus ini dilakukan perubahan koordinat, sehingga titik simpul (2) pada elemen Γj menjadi titik nol. Dengan demikian perlu dilakukan pemisahan daerah integrasi, yaitu Z (4) Z (2) Z (4) ∗ ∗ GW = φu dΓ = φu dΓ + φu∗ dΓ (83) (1)

(1)

(2)

Ada dua transformasi yang diperlukan sehingga dua integral terakhir pada persamaan (83) dapat dihitung dengan menggunakan metoda Gauss kuadratur

92


khusus, yaitu untuk masing-masing subelemen η=−

3ξ + 1 2

dan

η0 =

3ξ + 1 4

(84)

Jadi bentuk persamaan (83) dapat ditulis menjadi Z Z 1 1 2 1 4 1 φ ln φ ln GW = J dη + J dη 0 3 0 r 3 0 r

(85)

Sistem koordinat lokal untuk kasus ini adalah x1 = x1 − x2 ,

x3 = x3 − x2 ,

x4 = x4 − x2 ,

(86)

Untuk y dilakukan hal yang sama. x4 6

2 x

x3

-s 6 (4)

-s 6 (3)

y4

y3

s

-

(2)

y2 s ?

?

(1) s

ξ

−1 s

1

η

s

s

2 3

1 3

/

s

0

s g

s

− 13

1 3

s

1 s

A A 0 η AU

0

s

s

s

1 3

2 3

1

Gambar 4 Titik Collocation pada simpul (2) Proses pengerjaan untuk masing-masing subelemen dilakukan analog dengan kasus 1. Dengan memandang fungsi pangkat tiga masing-masing untuk subelemen pertama r11 (η) = A11 η 3 + B11 η 2 + C11 η + D11

(87)

dengan r11 (0) = 0,

r11 (1) = x1 ,

r11 (−1) = x3 ,

r11 (−2) = x4 .

(88)

93


dan untuk subelemen kedua r21 (η) = A21 η 3 + B21 η 2 + C21 η + D21

(89)

dengan 1 r21 (− ) = x1 , 2

r21 (0) = 0,

1 r21 ( ) = x3 , 2

r21 (1) = x4 .

dan kemudian menyelesaikannya, kita punya bentuk berikut. Z 1 2 1 φ ln GW = J dη 3 0 η Z −1/3 hn o2 ˜11 ξ + C˜11 φ ln A˜11 ξ 2 + B +

(90)

(91)

−1

Z 1 o2 i 1 n 1 4 2 2 ˜12 ξ + C˜12 + A˜12 ξ + B φ ln 0 J dη 0 J dξ + 3 0 η Z 1 hn + φ ln A˜21 ξ 2 −1/3

˜21 ξ + C˜21 +B

o2

n o2 i 1 2 ˜22 ξ + C˜22 + A˜22 ξ 2 + B J dξ

dengan A˜11 = 18 (9x3 + 3x1 − 3x4 ), A˜21 = − 41 (9x3 + 3x1 − 3x4 ), ˜21 = 1 (2x1 + x4 ), ˜11 = − 1 (2x1 + x4 ), B B 4 2 C˜11 = 18 (−9x3 + x1 + x4 ), C˜21 = − 41 (−9x3 + x1 + x4 )

(92)

Sekali lagi, indeks i2, i = 1, 2 adalah bentuk y yang sesuai dengan bbentuk x. Bentuk ke dua di atas mempunyai faktor −2 dari bentuk pertama, sehingga integral ke dua dan ke tiga pada persamaan (91) dapat digabungkan dengan tambahan suku, seperti yang terdapat pada suku terakhir pada persamaan berikut. Z Z 1 1 2 1 4 1 GW = φ ln J dη + φ ln 0 J dη 0 (93) 3 0 η 3 0 η Z 1 hn o2 ˜11 ξ + C˜11 + φ ln A˜11 ξ 2 + B −1

n o2 i 1 2 ˜12 ξ + C˜12 + A˜12 ξ 2 + B J dξ Z

1

+

φ ln(1/2) J dξ −1/3

94


Integral terakhir pada persamaan (93) dapat dihitung dengan metoda Gauss kuadratur dengan transformasi koordinat ξ 0 = 3ξ−1 2 Kasus 3: Titik sanding berada pada simpul (3) Kasus ini analog dengan kasus 2. Perubahan koordinat dilakukan demikian sehingga titik simpul (3) pada elemen Γj menjadi titik nol, dan kita punya Z

(4)

GW =

φu∗ dΓ =

(1)

Z

(3)

φu∗ dΓ +

(1)

Z

(4)

φu∗ dΓ

(94)

(3)

Transformasi yang digunakan masing-masing subelemen adalah η=

1 − 3ξ 4

dan

η0 =

3ξ − 1 2

(95)

Dengan pengerjaan yang sama dengan kasus 2 kita dapat menuliskan persamaan (94) menjadi

GW

Z 1 1 2 1 φ ln J dη + φ ln 0 J dη 0 = η 3 0 η 0 Z 1 hn o2 ˜11 ξ + C˜11 φ ln A˜11 ξ 2 + B + 4 3

1

Z

(96)

−1

n o2 i 1 2 ˜12 ξ + C˜12 + A˜12 ξ 2 + B J dξ Z

1

φ ln(2) J dξ

+ 1/3

dengan A˜11 = − 14 (9x2 + 3x4 − 3x1 ), A˜12 = − 14 (9y2 + 3y4 − 3y1 ), ˜11 = − 1 (2x1 + x4 ), ˜12 = − 1 (2y1 + y4 ), B B 2 2 1 2 1 4 ˜ ˜ C11 = 4 (9x − x − x ), C12 = 41 (9y2 − y1 − y4 )

(97)

Integral terakhir pada persamaan (96) dapat dihitung dengan metoda Gauss kuadratur dengan transformasi koordinat ξ 0 = 3ξ + 2 Kasus 4: Titik sanding berada pada simpul (4) Kasus ini analog dengan kasus 1, namun transformasi yang digunakan adalah η = 1−ξ 2

95


4. APLIKASI MASALAH POTENSIAL Untuk melihat keakuratan dari program di sini akan ditinjau hasil penghitungan program untuk permasalahan potensial pada daerah berikut. q=0 6y

? 6

u = 300-

T2

6

?

flux nol Temp.

0 6

x

T1

q=0

Temperatur

u=0

flux nol

-

Gambar 5 Definisi permasalahan Hasil tersebut akan dibandingkan dengan metoda elemen batas yang menggunakan elemen kuadratik terhadap penyelesaian analitik. 4.1. Aliran panas pada lempeng bujur sangkar Pandang pada suatu daerah tertutup berbentuk bujur sangkar dengan panjang sisi 6 satuan panjang, seperti terlihat pada gambar 5 kiri. Pada sisi tegak diberikan temperatur, yaitu sisi kiri 300 satuan panas dan sisi kanan 0 satuan panas. Pada kedua sisi datar diketahui nilai flux nol. Penyelesaian analitik dari masalah ini mengikuti aturan T = −50x + 300, Untuk membandingkan dengan hasil yang diperoleh dengan hasil dari elemen kuadratik perlu ditinjau hasil pada pendiskritan dengan jumlah simpul yang sama. Ini dapat dilakukan pada 8 elemen kubik dengan 12 elemen kuadratik, dan pada 16 elemen kubik dengan 24 elemen kuadratik, yaitu dengan dengan jumlah simpul 24 dan 48. Dengan elemen-elemen ini diperoleh perbandingan hasil seperti pada tabel di bawah ini.


X 0.0 1.0 2.0 3.0 4.0 5.0 6.0

Y

Eksak 300 250 200 150 100 50 0

2.0 2.0 3.0 4.0 4.0

2.0 4.0 3.0 2.0 4.0

200 200 150 100 100

X 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0

Y

Eksak 300 275 250 225 200 175 150 125 100 75 50 25 0

2.0 2.0 3.0 4.0 4.0

2.0 4.0 3.0 2.0 4.0

200 200 150 100 100

96

Potensial 8 Kubik 12 Kuadratik .300000000E+03 .300000000E+03 .250000000E+03 .250000000E+03 .200000000E+03 .200000000E+03 .150000000E+03 .150000000E+03 .999999998E+02 .999999998E+02 .499999998E+02 .499999998E+02 .000000000E+00 .000000000E+00 Titik di dalam .200000000E+03 .200000003E+03 .200000000E+03 .200000003E+03 .150000000E+03 .150000002E+03 .999999996E+02 .100000001E+03 .999999996E+02 .100000001E+03 Potensial 16 Kubik 24 Kuad. .300000000E+03 .300000000E+03 .275000000E+03 .275000000E+03 .250000000E+03 .250000000E+03 .225000000E+03 .225000000E+03 .200000000E+03 .200000000E+03 .175000000E+03 .175000000E+03 .150000000E+03 .150000000E+03 .125000000E+03 .125000000E+03 .100000000E+03 .999999999E+02 .750000000E+02 .749999999E+02 .500000000E+02 .499999998E+02 .250000001E+02 .249999998E+02 .000000000E+00 .000000000E+00 Titik di dalam .199999999E+03 .200000003E+03 .199999999E+03 .200000003E+03 .149999999E+03 .150000002E+03 .999999996E+02 .100000001E+03 .999999996E+02 .100000001E+03

Dari tabel terlihat bahwa untuk jumlah simpul 24 hasil potensial yang diperoleh kedua metoda untuk titik-titik pada batas adalah sama. Untuk titik-titik dalam dengan elemen kubik hanya terjadi galat 0.0000004% pada titik-titik dengan x = 4.0. Sedangkan dengan elemen kuadratik terjadi galat pada titik-titik dengan x = 2, x = 3 sebesar 0.0000002% dan pada titiktitik dengan x = 1 sebesar 0.0000001%. Dengan jumlah simpul 48 hasil yang

97


diperoleh dari elemen kubik relatif lebih baik dari elemen kuadratik. Letak terjadinya galat serupa dengan dengan jumlah 24 simpul. Jadi kesalahan pada elemen kubik untuk titik-titik dalam hanya terjadi pada satu sisi saja, sedangkan pada elemen kuadratik terjadi secara merata. Ini berarti elemen kubik untuk kasus ini relatif lebih baik dari pada elemen kuadratik. Hasil flux untuk kedua kasus di atas adalah sebagai berikut. Galat yang terjadi pada elemen kubik dengan 8 elemen untuk titik-titik pada elemen batas adalah antara 0.0000004% sampai 0.0000054%, sementara dengan 12 elemen kuadratik terjadi galat dibawah 0.0000022%. Kemudian untuk elemen kubik dengan 16 elemen terjadi galat antara 0.0000004% sampai 0.0000022%, sedangkan dengan 24 elemen kuadratik terjadi galat antara 0.0000002% sampai 0.0000020%. Dengan demikian untuk kedua kasus ini hampiran untuk nilai flux oleh elemen kuadratik lebih baik dari pada elemen kubik. Tabel hasil lengkap dari program yang dibuat terdapat pada lampiran. 4.2. Aliran panas dalam silinder Sekarang pandang contoh yang berkenaan dengan silinder seperti pada gambar 5 kanan dengan jari-jari r1 pada sisi dalam dan r2 pada sisi luar dengan temperatur yang diberikan masing-masing T1 dan T2 . Bila nilai numerik yang diberikan r1 = 1.0, r2 = 2.0, T1 = 10.0 dan T2 = 6.0 dengan konduktivitas 3.5, maka diperoleh penyelesaian analitik T = −5.77078016 ln(r) + 10.0

(98)

Untuk membandingkan hasil diambil pendiskritan masing dengan 10 dan 18 elemen kubik serta 15 dan 27 elemen kuadratik. Dengan demikian titik-titik simpul yang diperhatikan adalah sama. Hasil dari pendiskritan ini seperti pada tabel berikut. X 1.000 1.167 1.333 1.500 1.667 1.833 2.000

Y 0.0 0.0 0.0 0.0 0.0 0.0 0.0

Eksak 10.00000000 9.11043015 8.33985015 7.66014999 7.05213751 6.50212336 6.00000000

Temperatur 10 Kubik 15 Kuad. 10.00000000 10.00000000 9.11099303 9.10997766 8.34045962 8.33937059 7.66075605 7.65968983 7.05268478 7.05173871 6.50252604 6.50184447 6.00000000 6.00000000

Galat 10 Kubik diberikan 0.0062 0.0073 0.0079 0.0078 0.0062 diberikan

(%) 15 Kuad. diberikan 0.0050 0.0058 0.0060 0.0057 0.0043 diberikan

Pada kasus ini elemen kuadratik relatif lebih baik dari pada elemen kubik. Ini disebabkan batas tali busur lebih baik didekati dengan bentuk

98


kuadratik. Ini berarti elemen kubik tidak baik sebagai hampiran untuk mencari penyelesaian numerik dari persoalan aliran panas dalam silinder. X 1.000 1.083 1.167 1.250 1.333 1.417 1.500 1.583 1.667 1.750 1.833 1.917 2.000

Y 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

Eksak 10.00000000 9.53809131 9.11043015 8.71228762 8.33985015 7.98999850 7.66014999 7.34814007 7.05213751 6.77058031 6.50212364 6.24560208 6.00000000

Temperatur 18 Kubik 27 Kuad. 10.00000000 10.00000000 9.53812373 9.53805940 9.11046732 9.11039383 8.71232869 8.71224914 8.33989394 8.33981084 7.99004377 7.98995935 7.66019554 7.66011181 7.34818467 7.34810366 7.05217988 7.05210376 6.77061893 6.77055030 6.50215645 6.50209886 6.24562581 6.24558512 6.00000000 6.00000000

Galat 18 Kubik diberikan 0.00034 0.00041 0.00047 0.00053 0.00057 0.00059 0.00061 0.00060 0.00057 0.00050 0.00038 diberikan

(%) 27 Kuad. diberikan 0.00033 0.00040 0.00044 0.00047 0.00049 0.00050 0.00050 0.00048 0.00044 0.00038 0.00027 diberikan

KESIMPULAN Hasil penghitungan dari elemen kubik untuk masalah potensial masih cukup akurat, bahkan untuk daerah dengan batas berbentuk garis lurus lebih baik dari elemen kuadratik. Untuk daerah yang berbentuk melingkar, hasil hampirannya relatif lebih buruk dibandingkan dengan hasil yang diperoleh dari elemen kuadratik.

REFERENSI 1. Ames, W.F, Numerical Methods for Partial Differential Equations, Third

Edition Academic Press, Inc., 1992. 2. Becker, A.A., The Boundary Element Method in Engineering. A Complete Course, McGraw-Hill Book Company, New York 1992. 3. Brebbia, C.A. and Dominguez, J. Boundary Elements. An Introductary Course, Second Edition, McGraw-Hill Book Company, New York 1992. 4. Hayakawa, K. and Iso, Y., High-order Uniform Convergence Estimation of Boundary Solution for Laplace’s Equation, Publ. RIMS, Kyoto Univ, 27 (1991), pp. 333-345.


99

5. Sloan, I.H, Boundary Element Methods, Applied mathematical report,

AMR95/9, The University of New South Wales, Autralia. 6. van Daalen,E.F.G., Numerical and Theoretical Studies of Water Waves and Floating Bodies, Disertation, Universiteit Twente.

Tulus: Jurusan Matematika, FMIPA, Universitas Sumatera Utara, Jl.



Biote-

100

Bulletin of Mathematics (Indonesian Mathematical Society Aceh - Nord Sumatera) c/o Department of Mathematics, University of Sumatera Utara Jalan Bioteknologi 1, Kampus USU, Medan 20155, Indonesia E-mail:[email protected] Website: http://indoms-nadsumut.org

Editor-in-Chief: Herman Mawengkang, Dept. of Math., University of Sumatera Utara, Medan 20155, Indonesia Managing Editor: Hizir Sofyan, Dept. of Math., Unsyiah University, Banda Aceh 23111, Indonesia Associate Editors: Algebra and Geometry: Irawaty, Dept. of Math., Institut Teknologi Bandung, Bandung 40132, Indonesia Pangeran Sianipar, Dept. of Math., University of Sumatera Utara, Medan 20155, Indonesia Analysis: Oki Neswan, Dept. of Math., Institut Teknologi Bandung, Bandung 40132, Indonesia Maslina Darus, School of Math. Sciences, Universiti Kebangsaan Malaysia, Bangi, 43600 Malaysia. Applied Mathematics: Tarmizi Usman, Dept. of Math., Unsyiah University, Banda Aceh 23111, Indonesia Imran, Dept. of Math., Riau University, Pekanbaru, Indonesia Discrete Mathematics & Combinatorics: Edy T. Baskoro, Dept. of Math., Institut Teknologi Bandung, Bandung 40132, Indonesia Saib Suwilo, Dept. of Math., University of Sumatera Utara, Medan 20155, Indonesia Statistics & Probability Theory: Soebanar, Dept. of Math., Gadjah Mada University, Yogyakarta 55281, Indonesia Sutarman, Dept. of Math., University of Sumatera Utara, Medan 20155, Indonesia

Instruction to Authors and Others Information Manuscript submission. Manuscript should be submitted, in triplicate, to the Editor-in-Chief by e-mail: [email protected]. It is understood the submitted paper is not being considered for publication in other journals. If a paper is accepted for publication in our journal, the author is assumed to have transferred the copyright to the IndoMS Aceh-Sumut. Papers published in this journal are reviewed by Reviewer of this Journal. Form of manuscript. Authors are encouraged to write their paper in English or in Bahasa Indonesia. Manuscripts should be written by using software LaTeX. The template can be download from our website. Make sure that your paper has an abstract, key words and phrases, and the current Mathematics Subject Classification. Avoid complicated formulas in the title as well as in the abstract. Abstract. Even if the paper is not written in English, the abstract must be presented in English. It should summarize the main result(s) and, possibly, the method(s) used, in at most 250 words. If the paper is a survey paper, it must be indicated so in the abstract. Figures. If any, figures should be supplied in a form suitable for photographic reproduction after the printer has inserted any necessary lettering in the correct type. They should be printed or drawn in black ink on good white paper, at approximately the size that they will be printed. References. References should be listed in alphabetical order according to the surnames of the authors at the end of the paper. Title of the journals should be abbreviated in accordance with the Abbreviations of the Names of Serials in the Mathematical Reviews Annual Index (AMS). The following reference styles should be used: 1. D.R. Shier and Ch. Witzgall (1980). ”Arc Tolerance in Shortest Path and Network Flow Problems”, Networks, Vol. 10, pp. 227 - 291. 2. B. Liu. Generalized exponents of primitive simple graphs. Acta Math. Appl. Sin, 9:1 (1993), 36-43. References should be cited in the text as, for example, [1] or [5, Theorem 6]. Reprints and page charge. The first author of an accepted paper will receive 25 reprints or no more than 250 pages in total. To cover the expences of the Journal, authors are encouraged to contribute through their research grants or their institutions, at IDR 25,000 or USD 4 per page of the article. Subscription rate. Annual subscription rate is IDR 100,000 or USD 30 for individuals, or IDR 150,000 or USD 45 for institutions. Members of the Indonesian Mathematical Society are entitled to a 20

Bulletin of Mathematics Volume 01 Number 01, January 2009 Invers Z-Matriks P. Sianipar

1

Digraph jarang dengan eksponen besar Mardiningsih dan Astri Syafrianti

15

Menentukan titik awal penyelesaian quadratic asisgnment problem Elly Rosmaini

25

Developing sensitivity analysis method in minimum cost flow problem Irvan and Herman Mawengkang 39 Teori linear pembangkit gelombang multi arah tipe snake untuk flap tunggal Marwan 49 Prediksi pasang surut laut di selat Malaka dengan menggunakan Hamsom model Taufiq Iskandar 59 Vertex exponent of two-cycles Saib Suwilo and Masna Meisaroh

73

Metode elemen batas untuk masalah potensial dengan hampiran elemen kubik Tulus 79

INVERS Z-MATRIKS. Pangeran Sianipar

Recommend Documents