JURNAL MATEMATIKA DAN KOMPUTER Vol. 5. No. 1, 14 - 24, April 2002, ISSN : 1410-8518 __________________________________________________________________ METODE ELEMEN HINGGA UNTUK MASALAH SYARAT BATAS DARI OPERATOR DIFERENSIAL POSITIF Sutrima Jurusan matematika FMIPA UNS Abstract The purpose of this research is to investigate numerical solutions of boundary value problems of positive differential operator by finite-element methods. The results showed that the finite-element methods can be implemented to obtained numerical solutions of boundary value problems of second-order positive differential operator. Moreover, rate of convergence of Ritz finite-element with Lagrange quadratic interpolation be α = 3 in the the L2 -norm and α = 2 in the energy norm. Keywords : boundary problems, positive differential operator finite-element methods.
1. PENDAHULUAN Banyak permasalahan-permasalahan praktis dalam kehidupan nyata yang dapat diformulasikan ke dalam model matematika, khususnya Persamaan Diferensial (PD). Misalnya masalah yang dapat dimodelkan menjadi PD
−
d du p +q u = f , dx dx
0< x< L
(1)
dengan p fungi yang mempunyai turunan kontinu, q fungsi kontinu, dan f fungsi yang diketahui. Misalkan dari masalah (1) dilengkapi masalah syarat batas Dirichlet atau campuran di 0 dan L. Pada kenyataannya banyak model dari masalah syarat batas (1) tidak dapat diselesaikan secara analitis. Oleh karena itu metode penyelesaian secara numerik akan berperan dalam kasus ini. Metode variasional merupakan salah satu metode yang efektif dalam analisis fungsional untuk menyelesaikan masalah yang berkaitan dengan operator positif pada ruang Hilbert. Metode variasional ini didasarkan kepada minimisasi dari fungsional kuadratik terkait dari masalah yang
14
Metode Elemen Hingga … ( Sutrima ) __________________________________________________________________ diberikan. Dari metode variasional inilah metode lemen hingga dikembangkan. Salah satu metode variasional yang terkenal adalah metode variasional Ritz. Berdasarkan algoritma metode variasional Ritz (Reddy, 1986), untuk kasus masalah kontinu Au = f pada Ω
(2)
dengan A operator linear yang definit positif pada DA dalam ruang Hilbert H .. Penyelesaian persamaan (2) yaitu u ∈ DA dapat diperoleh dengan meminimumkan fungsional kuadratik J (u ) =< Au , u > -2 < f , u > = B (u , u ) - 2 < f , u > H atau oleh Bader (2001) dapat dinyatakan dengan formulasi variasional. B (v, u ) =< v, f > H , ∀v ∈ H
(3)
dengan B : H × H → R fungsi bilinear yang simetris kontinu dan definit positif sehingga B (u , v) = l (v), untuk setiap v didalam H .
untuk suatu fungsional linear kontinu l pada H . Widianto (2001) telah menggunakan metode elemen hingga yang merupakan pengembangan dari metode variasional untuk menyelesaikan masalah persamaan (1) dengan fungsi hampiran untuk tiap-tiap elemen menggunakan interpolasi Lagrange linear. Secara umum prinsip dari metode elemen hingga untuk persamaan (1) adalah membagi domain Ω menjadi beberapa subdomain, dalam setiap subdomain digunakan formulasi variasional sesuai dengan metode variasional yang digunakan. Untuk kasus-kasus tertentu seperti yang telah diberikan oleh Reddy (1984) yaitu persamaan diferensial order empat, penggunaan basis dari interpolasi Lagrange linear lebih sulit untuk digunakan. Sehingga perlu cara lain yang lebih mudah untuk menyelesaikannya. Cara tersebut antara lain adalah menggunakan interpolasi Lagrange kuadratik. Dengan menggunakan interpolasi Lagrange kuadratik diharapkan hasil yang diperoleh lebih baik dari interpolasi Lagrange linear.
15
JURNAL MATEMATIKA DAN KOMPUTER Vol. 5. No. 1, 14 - 24, April 2002, ISSN : 1410-8518 __________________________________________________________________ Untuk mengukur besar error dalam elemen hingga dapat digunakan norma kuadrat rata-rata atau norma energi, lihat (Reddy, 1984, 1986) dan (Akin, 2000). Hasil pendekatan elemen hingga U H dapat dikatakan konvergen dalam norma energi ke u jika berlaku
u −U H
m
≤ cH α untuk α > 0
(4)
dengan c adalah konstanta yang besarnya tidak bergantung pada u dan U H (Reddy, 1984). Konstanta α disebut laju konvergensi. Eksperimen numerik menunjukkan bahwa laju konvergensi elemen hingga Galerkin yang menggunakan interpolasi Lagrange linear adalah α = 2 terhadap norma kuadrat rata-rata dan α = 1 terhadap norma energi, lihat Widianto (2001), Reddy dan Volpi (1992). Dengan mendefinisikan operator
A≡−
d d p + q , maka dapat dx dx
dibuktikan operator ini adalah operator linear simetri dan definit positif pada DA . Dari fakta-fakta ini layak untuk diteliti tentang aplikasi dari metode elemen hingga untuk menyelesaikan masalah syarat batas (1).
2. HASIL DAN PEMBAHASAN Misal Ω = Ω ∪ Γ domain dari (1) dengan Ω = (0, L) dan Γ himpunan titik-titik limit dari Ω . Misal domain Ω
dibagi menjadi E subdomain
Ωe , e = 1, 2,..., E . Jika digunakan interpolasi Lagrange kudratik sebagai fungsi
hampirannya,
maka
setiap
subdomain
Ωe , e = 1, 2,..., E mempunyai
sifat
Ωe = x1e , x3e dengan x11 = 0 , x3E = L, dan x3e −1 = x1e untuk indeks e = 2,..., E . Panjang elemen ke-e adalah h1e + h2e dengan h1e = x2e − x1e dan h2e = x3e − x2e . Karena dalam penelitian ini diambil subdomain yang seragam, maka h e cukup ditulis dengan h . Untuk lebih jelasnya lihat pada Gambar 1.
16
Metode Elemen Hingga … ( Sutrima ) __________________________________________________________________
Elemen ke-1 • x =0
•
1 1
x
1 2
Elemen ke-2
• x = x12
•
1 3
x
2 2
Elemen ke-E
• x = x13
E −1 3
2 3
……… x
• = x1E
• x
E 2
• x =L E 3
Gambar 1. Pembagian subdomain untuk tiap-tiap elemen. Dengan demikian untuk subdomain ke-e berlaku −
d du pe + qe u = f e , dx dx
x1e ≤ x ≤ x3e
(5)
dengan pe , qe dan f e secara berurutan adalah fungsi p , q dan f yang terdefinisi pada domain Ωe . Dengan algoritma metode variasional Ritz diperoleh
0 =< Au − f , v > x3e
= ∫[− x1e
d du pe + qe u − f e ]v dx dx dx x3e
x3e
du dv du = ∫ ( pe + qe u v − f e v)dx + v − pe . dx dx dx x1e x1e
du du e Dengan mengambil Q1e = − pe dan −Q2 = − pe , maka hasil terakhir dx x1e dx x3e menjadi x3e
0 = ∫ ( pe x1e
du dv + qe u v − f e v)dx − Q1e v( x1e ) − Q2e v( x3e ) dx dx
atau dapat ditulis dalam bentuk 0 = Be (u , v) − le (v)
(6)
dengan xe3
xe
1
1
3 du dv + qe u v)dx dan le (v) = ∫ f e v dx + Q1e v( x1e ) + Q2e v( x3e ) . Be (u , v) = ∫ ( pe dx dx xe xe
Mudah ditunjukkan bahwa Be (u , v) suatu fungsi bilinear yang simetris.
17
JURNAL MATEMATIKA DAN KOMPUTER Vol. 5. No. 1, 14 - 24, April 2002, ISSN : 1410-8518 __________________________________________________________________ n
Dikonstruksikan penyelesaian Ritz •e ( x) = ∑ • jeψ ej ( x) pada domain Ωe . j =1
Jika diambil v = ψ ie , sehingga dari persamaan (6) diperoleh n
0 = ∑ kijeν ej − f je , i = 1,2,…, n j =1
dengan x3 e dψ ie dψ j k = ∫ ( pe + qeψ ieψ ej )dx dan f i e = ∫ f eψ ie dx + Q1eψ ie ( x1e ) + Q2eψ ie ( x3e ) . dx dx x1e x1e x3e
e
e ij
Untuk selanjutnya bentuk tersebut dapat disajikan dalam persamaan matriks K eu e = f e
(7)
dengan K e = kije , f e = f i e dan u e = • je .
Jika digunakan ide interpolasi Lagrange kuadratik maka fungsi •e ( x) menjadi 3
•e ( x) = ∑ •i eψ ie
(8)
i =1
3
dengan basis ψ = ∏ e i
j =1 j ≠i
x − x ej xie − x ej
. Anggota basis ψ ie ini sering juga disebut dengan
0 untuk i ≠ j . Selanjutnya, jika diambil fungsi bentuk, dan bersifat ψ i ( xi ) = 1 untuk i = j transformasi x = x1e + ξ dengan he = x2e − x1e = x3e − x2e dan 0 < ξ < 2he , dan disubtitusikan kedalam basis akan diperoleh entri-entri dari K e dan f e ,
k = e ij
2 he
∫ 0
dan
( pˆ e
2 he 2 he e dψ ie dψ j e e e e e e ˆ ˆ + qe ψ i ψ j )dξ , f1 = ∫ f e ψ 1 dξ + Q1 , f 2 = ∫ fˆe ψ 2e dξ , dξ dξ 0 0
f = e 3
2 he
∫
fˆe ψ 3e dξ + Q2e
0
dengan pˆ e ≡ pe ( xe + ξ ) , qˆe ≡ qe ( xe + ξ ) dan fˆe ≡ f e ( xe + ξ ) .
18
(9)
Metode Elemen Hingga … ( Sutrima ) __________________________________________________________________ Perakitan dari semua matriks K eu e = f e , e = 1, 2,..., E diperoleh
k111 k121 k131 0 0 0 1 1 1 k23 0 0 0 k21 k22 1 k31 k321 k331 + k112 k122 k132 0 2 2 2 k21 k22 k23 0 0 0 2 2 2 3 0 0 k31 k32 k33 + k11 k123 0 0 k213 k223 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
⋯ ⋯
0 0
0 0
⋯ ⋯
0 ⋮
0 ⋮
⋯ k33E−1 + k11E k12E ⋯ k21E k22E ⋯
k31E
k32E
0 U1 f11 0 U2 f21 0 U3 f31 + f12 0 U4 f22 0 U5 f32 + f13 = 0 U6 f23 ⋮ ⋮ ⋮ k13E U2E−1 f3E−1 + f1E k23E U2E f2E k33E U2E+1 f3E
atau dapat ditulis Ku = f
(10)
yang berlaku pada domain Ω = [ 0, L ] . Matriks K dalam persamaan (10) sering disebut sebagai matriks kekakuan. Dengan menerapkan syarat batasnya terhadap (10) akan diperoleh penyelesaian hampiran elemen hingga. 3. IMPLEMENTASI KASUS Kasus 1. Perhatikan masalah syarat batas
−
d 2 du ( x + 1) + 12 u = 18 x, dx dx
dengan syarat batas Dirichlet masalah ini adalah u ( x) =
0 < x <1
(11)
u (0) = 0 dan u (1) = 1 . Penyelesaian eksak dari
3 1 x − x 3 , ( Reddy, 1986). 2 2
Dalam hal ini DA = {v ∈ C 2 [ 0,1] : v(0) = 0 dan v(1) = 1} . Pada artikel ini pembagian domain seragam, artinya panjang karakteristik H
dari semua
subdomain sama. Misal banyaknya elemen adalah E . Dapat ditentukan 1 e −1 e subdomain ke-e adalah Ωe = . Jelas , , e = 1, 2,.., E dan h = h e = 2E E E panjang karakteristik dari tiap subdomain adalah H = 2h e =
1 . E
19
JURNAL MATEMATIKA DAN KOMPUTER Vol. 5. No. 1, 14 - 24, April 2002, ISSN : 1410-8518 __________________________________________________________________ Untuk mencari penyelesaiannya digunakan persamaan (9) dan (10) dengan syarat batas terkait. Misal banyaknya elemen adalah 2, sehingga dari proses 1 1 1 pembagian subdomain diperoleh Ω1 = 0, , Ω 2 = ,1 dan h = . Dengan 4 2 2 menerapkan persamaan (9) dan (10) diperoleh
224 15 − 29 5 0 Dengan
menyelesaikan
−
29 5
202 15 −
107 15
2 U 3 2 107 3 − U3 = . 5 15 409 304 U 4 30 15 0
persamaan
(11)
ini
(11)
diperoleh
U 2 = 0.367363,
U 3 = 0.687234, U 4 = 0.914586 dan dari syarat batas diketahui U1 = 0 dan U5 = 1.
Dalam
hal
ini
telah
ditentukan
bahwa
•11 = U1 , •21 = U 2 , •31 = U 3 = •1 2 , •22 = U 4 , dan •3 2 = U 5 . Sehingga penyelesaian hampiran adalah 1 2 ,0 < x ≤ 1.5644 x − 0.3799 x 2 U H ( x) = . −0.1933 + 2.3288 x − 1.1355 x 2 , 1 < x < 1 2 Untuk banyaknya elemen E =2, E =4, E =8, E =16, E = 32 hasilnya dapat dilihat pada Tabel 1 dan besar error dapat dilihat pada Tabel 2.
Tabel 1. u dan U H dengan E =2 , E = 4, E = 8, E = 16, dan E = 32 untuk Kasus 1.
20
x
E=2
E=4
E=8
E=16
E=32
u ( x)
0
0
0
0
0
0
0
0.25 0.5 0.75
0.367363 0.687234 0.914586
0.367176 0.687483 0.914049
0.367187 0.687499 0.914062
0.367187 0.687500 0.914062
0.367187 0.687500 0.914062
0.367188 0.687500 0.914063
1
1
1
1
1
1
1
Metode Elemen Hingga … ( Sutrima ) __________________________________________________________________ Tabel 2. Error U H untuk Kasus 1. E
u −UH
u −U H
2 4 8 16
0.002174 0.000270 0.000033 4.212 × 10−6
0.028140 0.007001 0.001748 0.000437
32
5.265 × 10 −7
0.000109
0
1
Untuk mencari laju konvergensi error dari Kasus 1 dapat digunakan ide dari Reddy dan Volpi [1992], sehingga dapat ditentukan Y1 = Log u − U H
0
dan
Y2 = Log u − U H 1 . Dari Gambar 2 terlihat bahwa Y1 dan Y2 bentuknya mendekati garis lurus. Dengan data error dari Tabel 3 dan menggunakan metode pendekatan kuadrat terkecil dapat diperoleh bahwa
Y1 = Log u − U H
0
≈ 3.0027 Log H − 4.0502
dan
Y2 = Log u − U H
1
≈ 2.0022 Log H − 2.1846.
Y
(12)
Y2
-4
Y1
-6
-3.5
-2.5
-2
-1.5
-1
Log H
-10 -12 -14
Gambar 2. Plot Y1 dan Y2 terhadap Log H untuk Kasus 1. Berdasarkan persamaan (12) dapat diketahui bahwa gradien Y1 terhadap Log H mendekati 3 dan Y2 terhadap Log H mendekati 2. Dengan kata lain konvergensi
21
JURNAL MATEMATIKA DAN KOMPUTER Vol. 5. No. 1, 14 - 24, April 2002, ISSN : 1410-8518 __________________________________________________________________ error untuk Kasus 1 adalah α = 3 terhadap norma kuadrat rata-rata dan α = 2 terhadap norma energi. Kasus 2. Perhatikan masalah syarat batas −
d du x ( x + 1) + x u = −2e − 2 x, dx dx
0 < x <1
(13)
dengan syarat batas campuran u (0) = 0 dan u '(1) = e − e −1 . Penyelesaian eksak dari (13) adalah u ( x) = e x + e− x − 2 , [Reddy, 1986]. Proses pembagian domain dilakukan sama seperti pada Kasus 1.
1 e −1 e Subdomain ke-e adalah Ωe = dan E , , e = 1, 2,.., E dengan h = h e = 2E E E banyaknya elemen. Fungsi hampiran untuk setiap elemen dipilih sesuai persamaan (8). Selanjutnya dapat diterapkan persamaan (9) dan (10) untuk mencari penyelesaiannya. Untuk nilai E =2, E =4, E = 8, E =16, dan E = 32 hasilnya dapat dilihat pada Tabel 3. Sedangkan besar errornya disajikan pada Tabel 4. Tabel 3. u dan U H dengan E =2 , E = 4, E = 8, E = 16, dan E = 32 untuk Kasus 2. x
E=2
E=4
E=8
E=16
E=32
u ( x)
0
0
0
0
0
0
0
0.25 0.5 0.75
0.062637 0.255224 0.58897
0.06286 0.25525 0.589363
0.062826 0.255252 0.589366
0.062826 0.255252 0.589367
0.062826 0.255252 0.589367
0.062826 0.255252 0.589367
1
1.08606
1.08606
1.08606
1.08606
1.08606
1.08606
Tabel 4. Error U H untuk Kasus 2.
22
E
u −UH
2
0,0008875
0.0114916
4 8 16
0,0001137 0,0000143 1.7964 × 10 −6
0.0029463 0.0007412 0.0001856
32
2.2211 × 10−7
0.0000464
0
u −U H
1
Metode Elemen Hingga … ( Sutrima ) __________________________________________________________________
Y
Y2 -6 -3.5
Y1 -2.5
-2
-1.5
-1
Log H
-10 -12 -14 Gambar 3. Plot Y1 dan Y2 terhadap Log H untuk Kasus 2. Dengan data error dari Tabel 4 disimpulkan bahwa
Y1 = Log u − U H dan
Y2 = Log u − U H
1
0
≈ 2.992 Log H − 4.943
≈ 1.99 Log H − 3.008.
Dari persamaan (14) dapat diketahui bahwa gradien Y1
(14) terhadap Log H
mendekati 3 dan gradien Y2 terhadap Log H mendekati 2, lihat Gambar 3. Dengan kata lain laju konvergensi error untuk Kasus 2 ini adalah α = 3 terhadap norma kuadrat rata-rata dan α = 2 terhadap norma energi. 4. KESIMPULAN Berdasarkan hasil dari yang telah diuraikan bahwa : 1. Metode elemen hingga Ritz dengan interpolasi Lagrange kuadratik dapat diterapkan untuk menyelesaikan masalah syarat batas (1). 2. Berdasarkan eksperimen numerik dalam implementasi kasus, laju konvergensi elemen hingga Ritz dengan interpolasi Lagrange kuadratik adalah α = 3 terhadap norma kuadrat rata-rata dan α = 2 terhadap norma energi.
5. UCAPAN TERIMA KASIH Artikel ini adalah sebagain dari hasil penelitian hibah dari Proyek DUE UNS Tahun 2001. Untuk itu peneliti mengucapkan terima kasih kepada Proyek DUE UNS atas dana yang dipercayakan kepada peneliti.
23
JURNAL MATEMATIKA DAN KOMPUTER Vol. 5. No. 1, 14 - 24, April 2002, ISSN : 1410-8518 __________________________________________________________________ DAFTAR PUSTAKA 1. Akin J. E, Finite Element Analysis with Error Estimation, Rice, University, Department of Mechanical Enginering and Materials Science, Houston, 2000. 2. Bader, Mathias, Energy Minimization Methods for Feature Displacements in Map Generalization, Dissertation zur Erlanggung der naturwissenschaftlichen Doktorwurde vorgelegtde Mathemathisch - naturwissenschaftlichen Fakultat der Universitat Zurich, Zurich, 2001. 3. Dwi Widianto, Analisis Elemen Hingga Galerkin untuk Persamaan Aliran Laminer, Skripsi Mahasiswa FMIPA UNS, Universitas Sebelas Maret Surakarta, 2001. 4. Reddy B.D, Volpi, M. B, Mixed Finite Element Methods for the Circular Arch Problem, Computer Methods in Applied Mechanics and Zengineering Elsevier Science Publishers, North-Holland, 1992, 97 : 125 - 145. 5. Reddy J. N, An Introduction to The Finite Element Method, Mc Graw-Hill Inc, New York, 1984. 6. Reddy J. N, Applied Functional Analysis and Variational Method in Engineering, Mc Graw-Hill Inc, Singapore, 1986.
24
