Ilmu berikut menggunakan teknik penyelesaian numerik ketika menyelesaikan permasalahan

Chapter ini memberikan pengantar bagi pembaca untuk mengenal terlebih dahulu mengenai metode numerik. Pada chapter ini akan dibahas mengenai apa itu metode numerik, perbedaannya dengan metode analitik, dan analisis error.

Metode numerik merupakan teknik penyelesaian permsalahn yang diformulasikan secara matematis dengan menggunakan operasi hitungan (aritmatik) yaitu operasi tambah, kurang, kali, dan bagi. Metode ini digunakan karena banyak permasalahan matematis tidak dapat diselesaikan menggunakan metode analitik. Jikapun terdapat penyelesaiannya secara analitik, proses penyelesaiaannya sering kali cukup rumit dan memakan banyak waktu sehingga tidak efisien.

Terdapat keuntungan dan kerugian terkait penggunaan metode numerik. Keuntungan dari metode ini antara lain:

1. Solusi persoalan selalu dapat diperoleh.
2. Dengan bantuan komputer, perhitungan dapat dilakukan dengan cepat serta hasil yang diperoleh dapat dibuat sedekat mungkin dengan nilai sesungguhnya.
3. Tampilan hasil perhitungan dapat disimulasikan.

Adapun kelemahan metode ini antara lain:

1. Nilai yang diperoleh berupa pendekatan atau hampiran.
2. Tanpa bantuan komputer, proses perhitungan akan berlangsung lama dan berulang-ulang.

Perbedaan antara metode numerik dan metode analitik dapat dijelaskan sebagai berikut:

1. Solusi metode numerik selalu berbentuk angka, sedangkan solusi metode analitik dapat berbentuk fungsi matematik yang selanjutnya dapat dievaluasi untuk menghasilkan nilai dalam bentuk angka.

2. Solusi dari metode numerik berupa hampiran, sedangkan metode analitik berupa solusi sejati. Kondisi ini berakibat pada nilai error metode analitik adalah 0, sedangkan metode numerik \(\neq 0\).

3. Metode analitik cocok untuk permasalahan dengan model terbatas dan sederhana, sedangkan metode numerik cocok dengan semua jenis permasalahan.

Terdapat beberapa tahapan dalam menyelesaikan suatu permasalahan dengan metode numerik. Tahapan-tahapan tersebut antara lain:

Persoalan dunia nyata dimodelkan ke dalam persamaan matematika. Persamaan matematika yang terbentuk dapat berupa persamaan linier, non-linier, dan sebagainya sesuai dengan persoalan yang dihadapi.

Model matematika yang dihasilkan dari tahap 1 mungkin saja terlalu kompleks. Semakin kompleks suatu model, semakin rumit penyelesaiaannya, sehingga model perlu disederhanakan.

Seberapa sederhana model yang akan kita buat? tergantung pada permasalahan apa yang hendak pembaca selesaikan. Model yang terlalu sederhana akan tidak cocok digunakan untuk digunakan sebagai pendekatan sistem nyata atau lingkungan yang begitu kompleks. Penyederhanaan dapat berupa asumsi sejumlah variabel yang terlibat tidak signifikan, atau asumsi kondisi reaktor (steady atau non-steady).

Setelah model matematika sederhana diperoleh, tahap selanjutnya adalah memformulasikan model matematika secara numerik. Tahapan ini terdiri atas: + menentukan metode numerik yang akan dipakai bersama-sama dengan analisis galat (error) awal. + menyusun algoritma dari metode numerik yang dipilih.

Tahap selanjutnya adalah menerjemahkan algoritma ke dalam program komputer. Pada tahapan ini pembaca bisa memilih bahasa pemrograman yang pembaca kuasai.

Dalam buku ini kita hanya akan berfokus pada bahasa pemrograman R. Pembaca dapat menggunakan bahasa pemrograman lain selain dari buku ini. Pembaca hanya perlu memperhatikan bagaimana penulis membangun algoritma penyelesaian dan memtransfernya menjadi bentuk sintaks R. Dari sintaks tersebut pembaca dapat melihat bagaimana meletakakkan tiap tahapan algoritma menjadi sintaks pada bahasa pemrograman.

Sebelum digunakan dengan data sesungguhnya, program komputer perlu dilakukan uji coba dengan data simulasi dan dievaluasi hasilnya. jika hasil keluaran diyakini sudah sesuai, baru dioperasikan dengan data yang sesungguhnya.

Bila program sudah selesai dijalankan dengan data yang sesungguhnya, maka hasil yang diperoleh dilakukan interpretasi, meliputi analisis hasil keluaran dan membandingkannya dengan prinsip dasar dan hasil-hasil empriik untuk menaksir kualitas soluasi numerik termasuk keputusan untuk menjalankan kembali progrma dengan memperoleh hasil yang lebih baik.

Perhatikan Gambar 5.1 berikut:

Gambar 5.1: 4 ilustrasi terkait akurasi dan presisi

Pada Gambar 5.1 terdapat 4 buah kondisi ketika kita menembakkan beberapa perluru pada sebuah sasaran. Tujuan kita disini adalah untuk menembak bagian tengah sasaran tersebut.

Pada Gambar (a) dan (c) pada Gambar 5.1 merupakan gambar yang menunjukkan seseorang telah berhasil mengenai bagian tengah sasaran tersebut dapat kita katakan pula tembakan pada kedua gambar tersebut akurat. Akurat dalam hal ini dapat diartikan suatu kondisi dimana kedekatan lubang peluru dengan pusat sasaran. Secara umum akurasi diartikan sebagai tingkat kedekatan pengukuran kuantitas terhadap nilai sebenarnya.

Terdapat dua buah cara untuk mengukur akurasi. Metode pengukuran akurasi antara lain: error absolut dan error relatif.

Error absolut merupakan nilai absolut dari selisih antara nilai sebenarnya \(x\) dengan nilai observasi \(x'\). Error absolut dapat dituliskan menggunakan Persamaan (5.1).

\[\begin{equation} \epsilon_A=\left|x-x'\right| \tag{5.1} \end{equation}\]

Pengukuran lain yang sering digunakan untuk mengukur akurasi adalah error relatif. Berbeda dengan error absolut, error relatif membagi selisih antara nilai sebenarnya \(x\) dan nilai observasi \(x'\) dengan nilai sebenarnya. Hasil yang diperoleh merupakan nilai tanpa satuan. Persamaan error relatif disajikan pada Persamaan (5.2).

\[\begin{equation} \epsilon_R=\left|\frac{x-x'}{x}\right| \tag{5.2} \end{equation}\]

Dalam suatu pengukuran, hal lain yang perlu diperhatikan selain akurasi adalah presisi. Presisi adalah sejauh mana pengulangan pengukuran dalam kondisi yang tidak berubah mendapat hasil yang sama. Berdasarkan Gambar 5.1, Gambar (a) dan (b) menunjukkan kepresisian yang tinggi. Hal ini terlihat dari jarak antara lubang peluru yang saling berdekatan dan mengelompok.

Berdasarkan Gambar 5.1 dapat kita simpulkan bahwa dalam suatu sistem pengukuran akan terdapat 4 buah kondisi. Pengukuran akurat dan presisi (Gambar (a)), tidak akurat namun presisi (Gambar (b)), akurat namun tidak presisi (Gambar (c)), dan tidak akurat serta tidak presisi (Gambar (d)).

Dari kondisi-kondisi tersebut, akan meuncul yang dinamakan error. Dalam analisa numerik error atau kesalahan menjadi hal yang perlu diperhatikan.

Kesalahan numerik merupakan error atau kesalahan yang timbul akibat adanya proses pendekatan atau hampiran. Kesalahan numerik terjadi karena tiga hal, antara lain:

• Kesalahan bawaan (inherent error), merupakan kesalahan data yang timbul akibat adanya pengkuran, human error seperti kesalahan pencatatan, atau tidak memahami hukum-hukum fisik dari data yang diukur.

• Kesalahan pembulatan (round-off error), adalah kesalahan yang terjadi karena adanya pembulatan. Contoh: 3,142857143… menjadi 3,14.

• Kesalahan pemotongan (truncation error), adalah kesalahan yang ditimbulkan pada saat dilakukan pengurangan jumlah angka signifikan.

Kesalahan atau error dapat diukur menggunakan Persamaan (5.1) dan Persamaan (5.2) yang telah penulis jelaskan pada Chapter 5.2.

1. Howard, J.P. 2017. Computational Methods for Numerical Analysis with R. CRC Press.
2. Sidiq, M. Tanpa Tahun. Materi Kuliah Metode Numerik. Repository Universitas Dian Nuswantoro.
3. Subakti, I. 2006. Metode Numerik. Institut Teknologi Sepuluh Nopember.
4. Sutarno,H., Rachmatin,D. 2008. Hands Out Metode Numerik. Universitas Pendidikan Indonesia.

Pada chapter ini penulis akan menjelaskan mengenai cara untuk menyelesaikan sistem persamaan linier. Adapun yang akan dibahas pada chapter ini antara lain:

• operasi Vektor dan matriks
• Metode Eliminasi Gauss
• Metode Dekomposisi matriks
• Studi Kasus

Pada Chapter 2.7 dan Chapter 2.8 telah dijelaskan sekilas bagaimana cara melakukan operasi pada vektor dan matriks. Pada chapter ini, penulis akan menambahkan operasi-operasi lain yang dapat dilakukan pada vektor dan matriks. Dasar-dasar operasi ini selanjutnya akan digunakan sebagai dasar menyusun algoritma penyelesaian sistem persamaan linier.

Misalkan saja diberikan vektor \(u\) dan \(v\) yang ditunjukkan pada Persamaan (6.1).

\[\begin{equation} u = \begin{bmatrix} u_1 \\[0.3em] u_2 \\[0.3em] \vdots \\[0.3em] u_n \end{bmatrix} dan\ v\ = \begin{bmatrix} v_1 \\[0.3em] v_2 \\[0.3em] \vdots \\[0.3em] v_n \end{bmatrix} \tag{6.1} \end{equation}\]

Jika kita menambahkan atau mengurangkan nilai elemen vektor dengan suatu skalar (konstanta yang hanya memiliki besaran), maka operasi penjumlahan/pengurangan akan dilakukan pada setiap elemen vektor.

\[\begin{equation} u \pm x = \begin{bmatrix} u_1 \\[0.3em] u_2 \\[0.3em] \vdots \\[0.3em] u_n \end{bmatrix} \pm x = \begin{bmatrix} u_1 \pm x \\[0.3em] u_2 \pm x \\[0.3em] \vdots \\[0.3em] u_n \pm x \end{bmatrix} \tag{6.2} \end{equation}\]

Jika kita melakukan penjumlahan pada vektor \(u\) dan \(v\), maka operasi akan terjadi pada masing-masing elemen dengan indeks yang sama.

\[\begin{equation} u \pm v = \begin{bmatrix} u_1 \\[0.3em] u_2 \\[0.3em] \vdots \\[0.3em] u_n \end{bmatrix} \pm v\ = \begin{bmatrix} v_1 \\[0.3em] v_2 \\[0.3em] \vdots \\[0.3em] v_n \end{bmatrix} = \begin{bmatrix} u_1 \pm v_1 \\[0.3em] u_2 \pm v_2 \\[0.3em] \vdots \\[0.3em] u_n \pm v_n \end{bmatrix} \tag{6.3} \end{equation}\]

Untuk lebih memahami operasi tersebut, berikut penulis berikan contoh penerapannya pada R:

u <- seq(1,5) v <- seq(6,10) # penjumlahan u+v

Bagaimana jika kita melakukan operasi dua vektor, dimaana salah satu vektor memiliki penjang yang berbeda?. Untuk memnjawab hal tersebut, perhatikan sintaks berikut:

Berdasarkan contoh tersebut, R akan mengeluarkan peringatan yang menunjukkan operasi dilakukan pada vektor dengan panjang berbeda. R akan tetap melakukan perhitungan dengan menjumlahkan kembali vektor \(u\) yang belum dijumlahkan dengan vektor \(x\) sampai seluruh elemen vektor \(u\) dilakukan operasi penjumlahan.

Operasi lain yang dapat dilakukan pada vektor adalah menghitung inner product dan panjang vektor. Inner product dihitung menggunakan Persamaan (6.4).

\[\begin{equation} u.v=\sum_{i=1}^nu_1v_1+u_2v_2+\dots+u_nv_n \tag{6.4} \end{equation}\]

Panjang vektor atau vektor yang telah dinormalisasi dihitung menggunakan Persamaan (6.5)

\[\begin{equation} \left|u\right|=\sqrt{u_1^2+u_2^2+\dots+u_n^2} \tag{6.5} \end{equation}\]

Berikut adalah contoh bagaimana cara menghitung inner product dan panjang vektor menggunakan R:

# panjang vektor u sqrt(sum(u*u))

Misalkan kita memiliki 2 buah matriks \(A\) dan \(B\).

\[\begin{equation} A = \begin{bmatrix} a_{1.1} & a_{1.2} &\cdots& a_{1.n} \\[0.3em] a_{2.1} & a_{2.2} &\cdots& a_{2.n} \\[0.3em] \vdots & \vdots &\ddots& \vdots \\[0.3em] a_{m.1} & a_{m.2} &\cdots& a_{m.n} \end{bmatrix} dan\ B = \begin{bmatrix} b_{1.1} & b_{1.2} &\cdots& b_{1.n} \\[0.3em] b_{2.1} & b_{2.2} &\cdots& b_{2.n} \\[0.3em] \vdots & \vdots &\ddots& \vdots \\[0.3em] b_{m.1} & b_{m.2} &\cdots& b_{m.n} \end{bmatrix} \tag{6.6} \end{equation}\]

Jika salah satu matriks tersebut dijumlahkan atau dikurangkan dengan skalar.

\[\begin{equation} A \pm x = \begin{bmatrix} a_{1.1} & a_{1.2} &\cdots& a_{1.n} \\[0.3em] a_{2.1} & a_{2.2} &\cdots& a_{2.n} \\[0.3em] \vdots & \vdots &\ddots& \vdots \\[0.3em] a_{m.1} & a_{m.2} &\cdots& a_{m.n} \end{bmatrix} \pm x = \begin{bmatrix} a_{1.1}\pm x & a_{1.2}\pm x &\cdots& a_{1.n}\pm x \\[0.3em] a_{2.1}\pm x & a_{2.2}\pm x &\cdots& a_{2.n}\pm x \\[0.3em] \vdots & \vdots &\ddots& \vdots \\[0.3em] a_{m.1}\pm x & a_{m.2}\pm x &\cdots& a_{m.n}\pm x \end{bmatrix} \tag{6.7} \end{equation}\]

Jika kedua matriks \(A\) dan \(B\) saling dijumlahkan atau dikurangkan. Perlu diperhatikan bahwa penjumlahan dua buah matriks hanya dapat dilakukan pada matriks dengan ukuran yang seragam.

\[\begin{equation} \begin{split} A \pm B & = \begin{bmatrix} a_{1.1} & a_{1.2} &\cdots& a_{1.n} \\[0.3em] a_{2.1} & a_{2.2} &\cdots& a_{2.n} \\[0.3em] \vdots & \vdots &\ddots& \vdots \\[0.3em] a_{m.1} & a_{m.2} &\cdots& a_{m.n} \end{bmatrix} \pm \begin{bmatrix} b_{1.1} & b_{1.2} &\cdots& b_{1.n} \\[0.3em] b_{2.1} & b_{2.2} &\cdots& b_{2.n} \\[0.3em] \vdots & \vdots &\ddots& \vdots \\[0.3em] b_{m.1} & b_{m.2} &\cdots& b_{m.n} \end{bmatrix} \\ & = \begin{bmatrix} a_{1.1}\pm b_{1.1} & a_{1.2}\pm b_{1.2} &\cdots& a_{1.n}\pm b_{1.n} \\[0.3em] a_{2.1}\pm b_{2.1} & a_{2.2}\pm b_{2.2} &\cdots& a_{2.n}\pm b_{2.n} \\[0.3em] \vdots & \vdots &\ddots& \vdots \\[0.3em] a_{m.1}\pm b_{m.1} & a_{m.2}\pm b_{m.2} &\cdots& a_{m.n}\pm b_{m.n} \end{bmatrix} \end{split} \tag{6.8} \end{equation}\]

Untuk lebih memahaminya, berikut disajikan contoh operasi penjumlahan pada matriks:

Operasi pehitungan lain yang penting pada matriks adalah operasi perkalian matriks. Perlu diperhatikan bahwa untuk perkalian matriks, jumlah kolom matriks sebelah kiri harus sama dengan jumlah baris pada matriks sebelah kanan. Perkalian antara dua matriks disajikan pada Persamaan (6.9).

\[\begin{equation} A_{m.n}\times B_{n.r}=AB_{m.r} \tag{6.9} \end{equation}\]

Pada R perkalian matriks dilakukan menggunakan operator %*%. Berikut adalah contoh perkalian matriks pada R:

# Perkalian matriks A%*%B

Terdapat tiga buah operasi dasar pada baris matriksoperasi baris elementer. Ketiga operasi ini akan menjadi dasar operasi sub-chapter selanjutnya. Ketiga operasi dasar tersebut antara lain:

1. Row Scalling. Mengalikan baris matriks dengan konstanta bukan nol.
2. Row Swaping. Menukar urutan baris pada sebuah matriks (contoh: menukar baris 1 dengan baris 2 dan sebaliknya).
3. Row Replacement. Baris matriks diganti dengan hasil penjumlahan atau pengurangan baris matriks tersebut dengan baris matriks lainnya, dimana baris matriks lainnya yang akan dijumlahkan/dikurangkan dengan matriks tersebut telah dilakukan proses row scalling. Luaran yang diperoleh pada umumnya adalah nilai nol pada baris matriks awal atau akhir.

Ketiga proses tersebut akan terjadi secara berulang, khusunya jika kita hendak mengerjakan sistem persamaan linier menggunakan algoritma eliminasi Gauss. Untuk mempermudah proses tersebut, kita dapat membuat masing-masing fungsi untuk masing-masing operasi tersebut. Algoritma fungsi-fungsi tersebut selanjutnya menjadi dasar penyusunan algoritma fungsi-fungsi eliminasi Gauss dan dekomposisi matriks yang akan dijelaskan pada chapter selanjutnya.

Fungsi row scalling pada R dapat dituliskan pada sintaks berikut:

scale_row <- function(m, row, k){ m[row, ] <- m[row, ]*k return(m) }

Berikut adalah contoh penerapannya:

# membuat matriks A (A <- matrix(1:15, nrow=5))

# lakukan scaling pada row 2 dengan nilai 10 scale_row(m=A, row=2, 10)

Catatan: Untuk menyimpan hasil perhitungan, simpan proses perhitungan dalam sebuah objek (lihat Chapter 2.5).

Row swapping merupakan proses yang berulang, kita perlu menyimpan terlebih dahulu baris matriks pertama kedalam sebuah objek. Baris matriks pertama selanjutnya diganti dengan baris matriks kedua, sedangkan baris matriks kedua selanjutnya akan diganti dengan baris matriks pertama yang telah terlebih dahulu disimpan dalam sebuah objek. Fungsi row swapping pada R dapat dituliskan pada sintaks berikut:

Berikut merupakan contoh penerapan fungsi swap_row():

# pertukarkan baris 2 dengan baris 5 swap_row(m=A, row1=2, row2=5)

Pada proses row replacement, proses perhitungan dilakukan dengan melakukan penjumlahan suatu baris matriks dengan baris matriks lainnya dengan terlebih dahulu melakukan row scalling terhadap matriks lainnya. Berikut adalah fungsi replace_row() yang ditulis pada R:

Berikut adalah contoh penerapan fungsi replace_row():

replace_row(m=A, row1=1, row2=3, k=-3)

Pada sub-chapter ini kita akan menggunakan operasi baris elementer yang telah dijelaskan pada Chapter 2.5. Terdapat dua topik yang akan dibahas pada sub-chapter ini, yaitu: row echelon form termasuk reduced row echelon form dan matriks tridiagonal.

Eliminasi Gauss merupakan sebuah cara untuk mencari penyelesaian sistem persamaan linier. Ide dasar dari eliminasi Gauss adalah melakukan operasi matematika pada baris matriks (lihat Chapter 2.5) dan melanjutkannya sampai hanya tersisa satu variabel saja. Kita dapat melakukan lebih dari satu operasi baris elementer pada proses elmininasi ini (contoh: mengalikan sebuah baris dengan konstanta dan menjumlahkan hasilnya pada baris lain).

Sebuah matriks merupakan row echelon form jika matriks tersebut memenuhi beberapa kondisi:

1. Angka bukan nol pertama dari kiri (leading coefficient) selalu di sebelah kanan angka bukan nol pertama pada baris di atasnya.
2. Baris yang terdiri dari semua nol ada di bagian bawah matriks.

Misalkan terdapat persamaan linier seperti yang ditunjukkan pada Persamaan (6.10).

\[\begin{equation} \begin{matrix} a_{1.1}x_1+a_{1.2}x_2+a_{1.3}x_3+\cdots+a_{1.n}x_n=b_1 \\ a_{2.1}x_1+a_{2.2}x_2+a_{2.3}x_3+\cdots+a_{2.n}x_n=b_2 \\ a_{3.1}x_1+a_{3.2}x_2+a_{3.3}x_3+\cdots+a_{3.n}x_n=b_3 \\ \cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots \\ a_{m.1}x_1+a_{m.2}x_2+a_{m.3}x_3+\cdots+a_{m.n}x_n=b_n \end{matrix} \tag{6.10} \end{equation}\]

dimana \(a_{i.j}\) untuk \(i=1\) sampai dengan \(m\) dan \(j=1\) sampai dengan \(n\) merupakan koefisien persamaan linier. \(x_i\) untuk \(i=1\) sampai dengan \(n\) merupakan variabel bebas pada sistem persamaan linier.

Persamaan linier pada Persamaan (6.10) dapat dinyatakan ke dalam bentuk matriks pada Persamaan (6.11).

\[\begin{equation} \begin{bmatrix} a_{1.1} & a_{1.2} & a_{1.3} &\cdots& a_{1.n} \\[0.3em] a_{2.1} & a_{2.2} & a_{2.3} &\cdots& a_{2.n} \\[0.3em] a_{3.1} & a_{3.2} & a_{3.3} &\cdots& a_{3.n} \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] a_{m.1} & a_{m.2} & a_{m.3} &\cdots& a_{m.n} \end{bmatrix} \begin{bmatrix} x_1 \\[0.3em] x_2 \\[0.3em] x_3 \\[0.3em] \cdots \\[0.3em] x_n \end{bmatrix} = \begin{bmatrix} b_1 \\[0.3em] b_2 \\[0.3em] b_3 \\[0.3em] \cdots \\[0.3em] b_n \end{bmatrix} \tag{6.11} \end{equation}\]

\[\begin{equation} AX=B \tag{6.12} \end{equation}\]

dimana:

• matriks A merupakan matriks koefisien / Jacobian
• vaktor X merupakan vaktor variabel
• vektor B merupakan vektor konstanta

matriks pada Persamaan (6.11) dapat diubah menjadi augmented matrix, yaitu: perluasan matriks A dengan menambahkan vektor B pada kolom terakhirnya.

\[\begin{equation} \begin{bmatrix} a_{1.1} & a_{1.2} & a_{1.3} &\cdots& a_{1.n} & b_1 \\[0.3em] a_{2.1} & a_{2.2} & a_{2.3} &\cdots& a_{2.n} & b_2 \\[0.3em] a_{3.1} & a_{3.2} & a_{3.3} &\cdots& a_{3.n} & b_3 \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] a_{m.1} & a_{m.2} & a_{m.3} &\cdots& a_{m.n} & b_n \end{bmatrix} \tag{6.13} \end{equation}\]

\[\begin{equation} A=\left[A|B\right] \tag{6.14} \end{equation}\]

Teorema 6.1 (spltheorem) Suatu sistem persamaan linier mempunyai penyelesaian tunggal bila memenuhi syarat-syarat sebagai berikut:

• ukuran persamaan linier simultan bujursangkar (jumlah persamaan sama dengan jumlah variabel bebas).
• sistem persamaan linier non-homogen di mana minimal ada satu nilai vektor konstanta \(B\) tidak nol atau terdapat \(b_{n}\neq 0\).
• Determinan dari matriks koefisiensistem persamaan linier tidak sama dengan nol.

Untuk memperoleh penyelesaian sistem persamaan linier, Persamaan (6.13) perlu dilakukan operasi baris elementer. Hasil operasi baris dasar akan menghasilkan matriks row echelon form yang disajikan pada Persamaan (6.15).

\[\begin{equation} \begin{bmatrix} a_{1.1} & a_{1.2} & a_{1.3} &\cdots& a_{1.n} & b_1 \\[0.3em] a_{2.1} & a_{2.2} & a_{2.3} &\cdots& a_{2.n} & b_2 \\[0.3em] a_{3.1} & a_{3.2} & a_{3.3} &\cdots& a_{3.n} & b_3 \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] a_{m.1} & a_{m.2} & a_{m.3} &\cdots& a_{m.n} & b_n \end{bmatrix} \implies \begin{bmatrix} c_{1.1} & c_{1.2} & c_{1.3} &\cdots& c_{1.n} & d_1 \\[0.3em] 0 & c_{2.2} & c_{2.3} &\cdots& c_{2.n} & d_2 \\[0.3em] 0 & 0 & c_{3.3} &\cdots& c_{3.n} & d_3 \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] 0 & 0 & 0 &\cdots& c_{m.n} & d_n \end{bmatrix} \tag{6.15} \end{equation}\]

Sehingga penyelesaian sistem persamaan linier dapat diperoleh menggunakan Persamaan (6.16).

\[\begin{equation} \begin{matrix} x_n=\frac{d_n}{c_{m.n}} \\ x_{n-1}=\frac{1}{C_{m-1.n-1}}\left(d_{n-1}-c_{m-1.n}x_n\right) \\ \cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots \\ x_2=\frac{1}{c_{2.2}}\left(d_2-c_{2.3}x_3-c_{2.4}x_4-\dots-c_{2.n}x_n\right) \\ x_1=\frac{1}{c_{1.1}}\left(d_1-c_{1.2}x_2-c_{1.3}x_3-\dots-c_{1.n}x_n\right) \end{matrix} \tag{6.16} \end{equation}\]

Contoh 6.1 Selesaikan sistem persamaan berikut:

\[ \begin{matrix} x_1+x_2+x_3=6 \\ x_1+2x_2-x_3=2 \\ 2x_1+x_2+2x_3=10 \\ \end{matrix} \]

Jawab:

Augmented matrix sistem persamaan linier tersebut adalah sebagai berikut:

\[ \begin{bmatrix} 1 & 1 & 1 & 6 \\[0.3em] 1 & 2 & -1 & 2 \\[0.3em] 2 & 1 & 2 & 10 \end{bmatrix} \]

Operasi baris elementer selanjutnya dilakukan pada matriks tersebut. Pada langkah pertama, baris ke-2 dikurangkan dengan baris ke-1 (\(B_2-B_1\)) dan baris ke-3 dikurangkan oleh dua kali baris ke-1 (\(B_3-2B_1\)).

\[\begin{equation*} \begin{bmatrix} 1 & 1 & 1 & 6 \\[0.3em] 1 & 2 & -1 & 2 \\[0.3em] 2 & 1 & 2 & 10 \end{bmatrix} \begin{matrix} B_2-B_1 \\ \implies \\ B_3-2B_1 \end{matrix} \begin{bmatrix} 1 & 1 & 1 & 6 \\[0.3em] 0 & 1 & -2 & -4 \\[0.3em] 0 & -1 & 0 & -2 \end{bmatrix} \end{equation*}\]

Hasil dari langkah pertama tersebut, selanjutnya menjadi input dari langkah selanjutnya. Pada langkah selanjutnya operasi baris elementer kembali dilanjutkan. Baris ke-3 dikurangkan denganbaris ke-2 (\(B_3-B_2\)).

\[\begin{equation*} \begin{bmatrix} 1 & 1 & 1 & 6 \\[0.3em] 0 & 1 & -2 & -4 \\[0.3em] 0 & -1 & 0 & -2 \end{bmatrix} \begin{matrix} B_3-B_2 \\ \implies \end{matrix} \begin{bmatrix} 1 & 1 & 1 & 6 \\[0.3em] 0 & 1 & -2 & -4 \\[0.3em] 0 & 0 & -2 & -6 \end{bmatrix} \end{equation*}\]

Setelah diperoleh matriks row echelon form selanjutnya penyelesaian persamaan dapat dikerjakan menggunakan Persamaan (6.16).

\[\begin{equation*} \begin{matrix} x_3=\frac{-6}{-2}=3 \\ x_2=\frac{1}{1}\left(-4-\left(2\right)3\right)=2 \\ x_1=\frac{1}{1}\left(6-2-3\right)=1 \end{matrix} \end{equation*}\]

Algoritma Row Echelon Form

1. Masukkan matriks \(A\), dan vektor \(B\) beserta ukurannya \(n\)
2. Buat augmented matrix \(\left[A|B\right]\) namakan dengan \(A\)
3. Untuk baris ke-\(i\) dimana \(i=1\) s/d \(n\), perhatikan apakah nilai \(a_{i,j}\) sama dengan nol. a) Bila iya, lakukan row swapping antara baris ke-\(i\) dan baris ke-\(i+k\leq n\), dimana \(a_{i+k,j}\) tidak sama dengan nol. Bila tidak ada berarti perhitungan tidak bisa dilanjutkan dan proses dihentikan dengan tanpa penyelesaian, b) Bila tidak, lanjutkan.
4. Untuk baris ke-\(j\), dimana \(j=i+1\) s/d \(n\), lakukan operasi baris elementer:a) Hitung \(c=\frac{a_{j,i}}{a_{i,i}}\), b) untuk kolom \(k\), dimana \(k=1\) s/d \(n+1\), hitung \(a_{j,k}=a_{j,k}-c.a_{i,k}\).
5. Hitung akar, untuk \(i=n\) s/d 1 (bergerak dari baris pertama) menggunakan Persamaan (6.16).

Berdasarkan algoritma tersebut, kita dapat menyusun fungsi pada R untuk menyelesaikan sistem persamaan linier menggunakan matriks row echelon form. Fungsi yang akan dibentuk hanya sampai pada algoritma ke-4. Proses substitusi akan dilakukan secara manual. Berikut adalah sintaks yang digunakan:

Dengan menggunakan fungsi ref_matrix(), kita dapat membentuk matriks row echelon form pada Contoh 6.1.

am <- c(1,1,2, 1,2,1, 1,-1,2, 6,2,10) (m <- matrix(am, nrow=3))

matriks yang diperoleh selanjutnya dapat diselesaikan menggunakan Persamaan (6.16).

Contoh 6.2 Dengan menggunakan fungsi ref_matrix(), buatlah matriks row echelon form dari sistem persamaan linier berikut:

\[ \begin{matrix} 2x_1+x_2-x_3=1 \\ 3x_1+2x_2-2x_3=1 \\ x_1-5x_2+4x_3=3 \\ \end{matrix} \]

Jawab:

Augmented matrix dari sistem persamaan tersebut adalah sebagai berikut:

\[ \begin{bmatrix} 2 & 1 & -1 & 1 \\[0.3em] 3 & 2 & -2 & 1 \\[0.3em] 1 & -5 & 4 & 3 \end{bmatrix} \]

Penyelesaian matriks tersebut adalah sebagai berikut:

(m <- matrix(c(2,3,1, 1,2,-5, -1,-2,4, 1,1,3), nrow=3))

Proses lebih lanjut akan menghasilkan penyelesaian sebagai berikut:

\[ \begin{matrix} x_1=1 \\[0.3em] x_2=2 \\[0.3em] x_3=3 \end{matrix} \]

Berbeda dengan metode eliminasi Gauss yang telah dijelaskan pada Chapter 6.3.1, metode eliminasi Gauss-Jordan membentuk matriks menjadi bentuk reduced row echelon form. Metode ini merupakan pengembangan metode eliminasi Gauss, dimana matriks sebelah kiri augmented matrix diubah menjadi matriks diagonal (lihat Persamaan (6.17)).

\[\begin{equation} \begin{bmatrix} a_{1,1} & a_{1,2} & a_{1,3} &\cdots& a_{1,n} & b_1 \\[0.3em] a_{2,1} & a_{2,2} & a_{2,3} &\cdots& a_{2,n} & b_2 \\[0.3em] a_{3,1} & a_{3,2} & a_{3,3} &\cdots& a_{3,n} & b_3 \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] a_{m,1} & a_{m,2} & a_{m,3} &\cdots& a_{m,n} & b_n \end{bmatrix} \implies \begin{bmatrix} 1 & 0 & 0 &\cdots& 0 & d_1 \\[0.3em] 0 & 1 & 0 &\cdots& 0 & d_2 \\[0.3em] 0 & 0 & 1 &\cdots& 0 & d_3 \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] 0 & 0 & 0 &\cdots& 1 & d_n \end{bmatrix} \tag{6.17} \end{equation}\]

Sehingga penyelesaian persamaan linier tersebut adalah nilai \(d_1,d_2,d3,\dots,d_n\) dan atau:

\[\begin{equation} \begin{matrix} x_1=d_1 \\ x_2=d_2 \\ x_3=d_3 \\ \cdots\cdots\cdots\cdots\cdots\cdots\cdots\cdots \\ x_n=d_n \end{matrix} \tag{6.18} \end{equation}\]

Contoh 6.3 Selesaikan sistem persamaan berikut:

\[ \begin{matrix} x_1+x_2=3 \\ 2x_1+4x_2=8 \\ \end{matrix} \]

Jawab:

Augmented matrix dari persamaan linier tersebut adalah sebagai berikut:

\[ \begin{bmatrix} 1 & 1 & 3 \\[0.3em] 2 & 4 & 8 \end{bmatrix} \]

Operasi baris elementer selanjutnya dilakukan pada matriks tersebut.

\[\begin{equation*} \begin{bmatrix} 1 & 1 & 3 \\[0.3em] 2 & 4 & 8 \end{bmatrix} \begin{matrix} B_2-2B_1 \\ \implies \end{matrix} \begin{bmatrix} 1 & 1 & 3 \\[0.3em] 0 & 2 & 2 \end{bmatrix} \end{equation*}\]

\[\begin{equation*} \begin{bmatrix} 1 & 1 & 3 \\[0.3em] 0 & 2 & 2 \end{bmatrix} \begin{matrix} \frac{B_2}{2} \\ \implies \end{matrix} \begin{bmatrix} 1 & 1 & 3 \\[0.3em] 0 & 1 & 1 \end{bmatrix} \end{equation*}\]

\[\begin{equation*} \begin{bmatrix} 1 & 1 & 3 \\[0.3em] 0 & 1 & 1 \end{bmatrix} \begin{matrix} B_1-B_2 \\ \implies \end{matrix} \begin{bmatrix} 1 & 0 & 2 \\[0.3em] 0 & 1 & 1 \end{bmatrix} \end{equation*}\]

Penyelesaian persamaan linier tersebut adalah sebagai berikut:

\[ x_1=2\ dan\ x_2=1 \]

Algoritma Metode Eliminasi Gauss-Jordan

1. Masukkan matriks \(A\) dan vektor \(B\) beserta ukurannya \(n\)
2. Buat augmented matrix \(\left[A|B\right]\) namakan dengan \(A\)
3. Untuk baris ke-\(i\) dimana \(i=1\) s/d \(n\)

• Perhatikan apakah nilai \(a_{i.i}\) sama dengan nol:

  • Bila ya: pertukakan baris ke-\(i\) dan baris ke-\(i+k\le n\), dimana \(a_{i+k.i}\) tidak sama dengan nol, bila tidak ada berarti perhitungan tidak bisa dilanjutkan dan proses dihentikan dengan tanpa penyelesaian.
  • Bila tidak: lanjutkan

• Jadikan nilai diagonalnya menjadi satu dengan cara untuk setiap kolom \(k\) dimana \(k=1\) s/d \(n+1\), hitung \(a_{i.k}=\frac{a_{i.k}}{a_{i.i}}\)

1. Untuk baris ke-\(j\), dimana \(j=i+1\) s/d \(n\). Lakukan operasi baris elementer untuk kolom \(k\) dimana \(k=1\) s/d \(n\).

• Hitung \(c=a_{j.i}\)
• Hitung \(a_{j.k}=a_{j.k}-c.a_{i.k}\)

1. Penyelesaian untuk \(i=n\) s/d 1 disajikan pada Persamaan (6.18).

Dari algoritma tersebut, kita dapat membangun sebuah fungsi menggunakan R. Fungsi tersebut adalah sebagai berikut:

gauss_jordan <- function (a){ m <- nrow (a) n <- ncol (a) piv <- 1 # cek elemen diagonal utama apakah bernilai nol for(row_curr in 1:m){ if(piv <= n){ i <- row_curr while(a[i, piv] == 0 && i < m){ i <- i + 1 if(i > m){ i <- row_curr piv <- piv + 1 if(piv > n) return (a) } } # jika diagonal utama bernilai nol,lakukan row swapping if(i != row_curr) a <- swap_row(a, i, row_curr) # proses pembentukan matriks reduced row echelon form piv_val <- a[row_curr , piv] a <- scale_row (a, row_curr , 1 / piv_val) for(j in 1: m){ if(j != row_curr){ k <- a[j, piv]/a[row_curr, piv] a <- replace_row (a, row_curr, j, -k) } } piv <- piv + 1 } } return (a) }

Dengan menggunakan fungsi gauss_jordan(), sistem persamaan linier pada Contoh 6.3:

(m <- matrix(c(1,2,1,4,3,8), nrow=2))

Contoh 6.4 Dengan menggunakan fungsi gauss_jordan(), carilah penyelesaian sistem persamaan linier pada Contoh 6.1 dan Contoh 6.2:

Jawab:

Untuk Contoh 6.1:

am <- c(1,1,2, 1,2,1, 1,-1,2, 6,2,10) m <- matrix(am, nrow=3) gauss_jordan(m)

Untuk Contoh 6.2:

m <- matrix(c(2,3,1,1,2,-5, -1,-2,4,1,1,3), nrow=3) gauss_jordan(m)

Metode eliminasi Gauss merupakan metode yang sederhana untuk digunakan khususnya jika semua koefisien bukan nol berkumpul pada diagonal utama dan beberapa diagonal sekitarnya. Suatu sistem yang bersifat demikian disebut sebagai banded dan banyaknya diagonal yang memuat koefisien bukan nol disebut sebagai bandwidth. Contoh khusus yang sering dijumpai adalah matriks tridiagonal yang memiliki bandwidth tiga.

Proses eliminasi untuk matriks tridiagonal bersifat trivial karena dengan membentuk sebuah subdiagonal tambahan, proses substitusi mundur segera dapat dilakukan. Bentuk matriks tridiagonal disajikan pada Persamaan (6.19).

\[\begin{equation} \begin{bmatrix} a_{1,1} & a_{1,2} & 0 &\cdots& 0 \\[0.3em] a_{2,1} & a_{2,2} & a_{2,3} &\cdots& 0 \\[0.3em] 0 & a_{3,2} & a_{3,3} &\cdots& 0 \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] 0 & 0 & 0 &\cdots& a_{m,n} \end{bmatrix} \begin{bmatrix} x_1 \\[0.3em] x_2 \\[0.3em] x_3 \\[0.3em] \cdots \\[0.3em] x_n \end{bmatrix} = \begin{bmatrix} b_1 \\[0.3em] b_2 \\[0.3em] b_3 \\[0.3em] \cdots \\[0.3em] b_n \end{bmatrix} \tag{6.19} \end{equation}\]

Penyelesaian persamaan tersebut disajikan pada Persamaan (6.20).

\[\begin{equation} x_n=\frac{b_n}{a_{m.n}};\ x_i=\frac{b_i-a_{i,j+1}x_{i+1}}{a_{i,j}} \tag{6.20} \end{equation}\]

dimana \(i=n-1,n-2,\dots,1\).

Pada beberapa textbook, diagonal matriks sering dilambangkan dengan \(l\)(diagonal bawah), \(d\)(diagonal tengah), dan \(u\) (diagonal atas). Bentuk matriksnya disajikan pada Persamaan (6.21).

\[\begin{equation} \begin{bmatrix} d_{1} & u_{2} & 0 &\cdots& 0 \\[0.3em] l_{2} & d_{2} & u_{3} &\cdots& 0 \\[0.3em] 0 & l_{3} & d_{3} &\cdots& 0 \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] 0 & 0 & 0 &\cdots& d_{n} \end{bmatrix} \begin{bmatrix} x_1 \\[0.3em] x_2 \\[0.3em] x_3 \\[0.3em] \cdots \\[0.3em] x_n \end{bmatrix} = \begin{bmatrix} b_1 \\[0.3em] b_2 \\[0.3em] b_3 \\[0.3em] \cdots \\[0.3em] b_n \end{bmatrix} \tag{6.21} \end{equation}\]

Algoritma Penyelesaian Matrik Tridiagonal

1. Bentuk sistem persamaan linier menjadi matriks pada Persamaan (6.21).
2. Lakukan foward sweep. Setiap elemen diagonal \(l\) dieliminasi menggunakan reduksi baris.

• Untuk \(i=1\)

  • Hitung \(u_1=\frac{u_1}{d_1}\)
  • Hitung \(b_1=\frac{b_1}{d_1}\)

• Untuk \(i=2\) s/d \(n-1\)

  • Hitung \(u_i=\frac{u_i}{d_i-l_i\times u_{i-1}}\)
  • Hitung \(b_i=\frac{b_i-l_i\times u_{i-1}}{d_i-l_i\times u_{i-1}}\)

• Hitung \(b_n=\frac{b_n-l_n\times u_{n-1}}{d_n-l_n\times u_{n-1}}\)

1. Lakukan backward sweep. Setiap elemen diagonal \(u\) dilakukan eliminasi.

• Untuk \(i=n-1\) s/d \(1\)

  • Hitung \(x_n=b_i-u_i\times x_{i+1}\)

• Hitung \(x_n=b_n\)

Berdasarkan algoritma tersebut, kita dapat membangun sebuah fungsi pada R. Fungsi penyelesaian matriks tridiagonal disajikan sebagai berikut:

Contoh 6.5 Selesaikan sistem persamaan berikut menggunakan fungsi tridiagmatrix() dan fungsi gauss_jordan()!

\[ \begin{matrix} 3x_1+4x_2=20 \\ 4x_1+5x_2-2x_3=28 \\ 2x_2+5x_3-3x_4=18 \\ 3x_3+5x_4=18 \end{matrix} \]

Jawab:

Langkah pertama untuk menyelesaikannya, kita harus merubah persamaan tersebut kedalam bentuk matriks

\[\begin{equation*} \begin{bmatrix} 3 & 4 & 0 & 0 \\[0.3em] 4 & 5 & 2 & 0 \\[0.3em] 0 & 2 & 5 & 3 \\[0.3em] 0 & 0 & 3 & 5 \end{bmatrix} x = \begin{bmatrix} 20 \\[0.3em] 28 \\[0.3em] 18 \\[0.3em] 18 \end{bmatrix} \end{equation*}\]

Untuk menyelesaikan persamaan tersebut menggunakan fungsi tridiagmatrix(), kita perlu membentuk vektor diagonal \(l\), \(d\), \(u\), dan \(b\).

l <- u <- c(4, 2, 3); d <- c(3, 5, 5, 5) b <- c(20, 28, 18, 18)

Setelah terbentuk, vektor tersebut dapat langsung dimasukkan ke dalam fungsi tridiagmatrix().

tridiagmatrix(L=l, D=d, U=u, b=b)

Untuk menyelesaikannya menggunakan fungsi gauss_jordan(), kita perlu membentuk augmented matrix-nya terlebih dahulu.

m <- matrix(c(3,4,0,0,4,5,2,0, 0,2,5,3,0,0,3,5, 20,28,18,18), nrow=4) gauss_jordan(m)

R menyediakan fungsi bawaan solve() untuk menyelesaiakan sistem persamaan linier. Format fungsi solve() adalah sebagai berikut:

Catatan:

• a: matriks koefisien atau matriks segiempat
• b: vektor konstanta

Berikut adalah contoh penerapan fungsi solve() pada sistem persamaan linier yang disajikan pada Contoh 6.2:

# memecah matriks m menjadi matriks koefisien dan vektor konstanta a <- matrix(c(2,3,1,1,2,-5,-1,-2,4),nrow=3) b <- c(1,1,3) solve(a,b)

Jika kita hanya memasukkan matriks persegi, maka output yang akan dihasilkan adalah invers dari matriks yang kita masukkan.

Jika kita mengalikan invers dengan matriks semula, maka akan dihasilkan output berupa matriks identitas.

Penyelesaian matriks tridiagonal selain menggunakan fungsi solve(), juga dapat menggunakan fungsi Solve.tridiag() dari Paket limSolve. Untuk menginstall dan mengaktifkan Paket tersebut, jalankan sintaks berikut:

install.packages("limSolve")

Fungsi Solve.tridiag() memiliki format sebagai berikut:

Solve.tridiag ( diam1, dia, diap1, B=rep(0,times=length(dia)))

Catatan:

• diam1: vektor bukan nol di bawah diagonal matriks
• dia: vektor bukan nol pada diagonal matriks
• diap1: vektor bukan nol di atas diagonal matriks
• B: vektor konstanta

Untuk memahami penerapannya, kita akan menggunakan kembali matriks yang ada pada Contoh 6.5.

l <- u <- c(4, 2, 3); d <- c(3, 5, 5, 5) b <- c(20, 28, 18, 18) Solve.tridiag(diam1=l, dia=d, diap1=u, B=b)

Seringkali kita diminta untuk memperoleh nilai penyelesaian suatu persamaan linier \(Ax=B\), dimana nilai vektor \(B\) yang selalu berubah-ubah. Penggunaan metode eliminasi Gauss mengharuskan untuk menyelesaikan sistem persamaan linier \(Ax=B\) secara terpisah untuk setiap perubahan vektor \(B\). Untuk menghindari pekerjaan eliminasi yang selalu berulang-ulang, faktorisasi menjadi suatu hal yang dapat dilakukan untuk mempersingkat prosesnya. Faktorisasi atau dekomposisi matriks merupakan suatu algoritma untuk memecah matriks \(A\), hasil pemecahan ini selanjutnya digunakan untuk memperoleh penyelesaian sistem persamaan linier melalui perkalian antara vektor \(B\) dan hasil faktorisasi matriks \(A\).

Misalkan kita memiliki persamaan linier seperti yang ditunjukkan oleh Persamaan (6.12). Pada metode dekomposisi LU, matriks \(A\) difaktorkan menjadi matriks \(L\) dan matriks \(U\), dimana ukuran kedua matriks tersebut harus sama dengan ukuran matriks \(A\) atau dapat kita tuliskan bahwa hasil perkalian kedua matriks tersebut akan menghasilkan matriks \(A\).

\[\begin{equation} A=LU \tag{6.22} \end{equation}\]

Sehingga Persamaan (6.12) akan menjadi Persamaan (6.23).

\[\begin{equation} LUx=b \tag{6.23} \end{equation}\]

Langkah penyelesaian sistem persamaan linier, diawali dengan menghadirkan vektor \(t\) yang ditunjukkan pada Persamaan (6.24).

\[\begin{equation} Ux=t \tag{6.24} \end{equation}\]

Langkah pada Persamaan (6.24) tidak dimaksudkan untuk menghitung vektor \(t\), melainkan untuk menghitung vektor \(x\). Vektor \(t\) diperoleh dengan menggunakan Persamaan (6.25).

\[\begin{equation} Lx=t \tag{6.25} \end{equation}\]

Kita dapat menyelesaikan sistem persamaan yang ditunjukkan pada Persamaan (6.24) dan Persamaan (6.25) menggunakan berbagai algoritma penyelesaian yang telah dibahas sebelumnya. Namun, karena matriks \(L\) merupakan matriks segitiga bawah dengan nilai nol berada pada bagian atas diagonal utama, penyelesaian \(t\) mengambil langkah yang lebih sedikit. Kondisi ini sama dengan kondisi penyelesaian matriks tridiagonal, dimana kita memanfaatkan sejumlah jalan pintas penyelesaiaannya guna mempercepat komputasi. Matriks segitia bawah \(L\) akan berupa matriks persegi dengan ukuran \(m\), di mana \(m\) merupakan jumlah baris matriks \(A\). Persamaan (6.25) dalam bentuk matriks akan terlihat seperti Persamaan (6.26).

\[\begin{equation} Lt= \begin{bmatrix} 1 & 0 & 0 &\cdots& 0 \\[0.3em] l_{2,1} & 1 & 0 &\cdots& 0 \\[0.3em] l_{3,1} & l_{3,2} & 1 &\cdots& 0 \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] l_{m,1} & l_{m,2} & l_{m,3} &\cdots& 1 \end{bmatrix} \begin{bmatrix} t_1 \\[0.3em] t_2 \\[0.3em] t_3 \\[0.3em] \cdots \\[0.3em] t_n \end{bmatrix} = \begin{bmatrix} b_1 \\[0.3em] b_2 \\[0.3em] b_3 \\[0.3em] \cdots \\[0.3em] b_n \end{bmatrix} \tag{6.26} \end{equation}\]

Berdasarkan Persamaan (6.26), diketahui nilai \(t_1=b_1\). Nilai ini selanjutnya dapat digunakan untuk melakukan proses substitusi guna memperoleh seluruh nilai vektor \(t\). Proses ini disebut sebagai foward substitution. Proses substitusi dapat dituliskan menggunakan Persamaan (6.27).

\[\begin{equation} t_i=b_i-\sum_{j=1}^{i-1}l_{i,j}t_i \tag{6.27} \end{equation}\]

Seteleh nilai vektor \(t\) dihitung, kita dapat menghitung nilai \(x\) pada Persamaan (6.28).

\[\begin{equation} Ux= \begin{bmatrix} u_{1,1} & u_{1,2} & u_{1,3} &\cdots& u_{1,n} \\[0.3em] 0 & u_{2,2} & u_{2,3} &\cdots& u_{2,n} \\[0.3em] 0 & 0 & u_{3,3} &\cdots& u_{3,n} \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] 0 & 0 & 0 &\cdots& u_{m,n} \end{bmatrix} \begin{bmatrix} x_1 \\[0.3em] x_2 \\[0.3em] x_3 \\[0.3em] \cdots \\[0.3em] x_n \end{bmatrix} = \begin{bmatrix} t_1 \\[0.3em] t_2 \\[0.3em] t_3 \\[0.3em] \cdots \\[0.3em] t_n \end{bmatrix} \tag{6.28} \end{equation}\]

Jika diperhatikan, kita dapat mengetahui mengetahui nilai \(x_n=\frac{t_n}{u_{m,n}}\). Nilai tersebut selanjutnya dapat digunakan untuk melakukan proses susbtitusi pada nilai lainnya. Proses substitusi ini disebut sebagai backward substitution. Proses dekomposisi atau faktorisasi LU digambarkan pada Gambar 6.1.

Gambar 6.1: Tahapan dekomposisi LU.

Dekomposisi LU didasarkan pada operasi baris elementer. Pertama, kita perlu menemukan matriks segitiga atas yang sesuai dengan matriks \(A\). Solusi untuk melakukan dekomposisi bisa jadi tak terhingga, namun solusi yang paling sederhana adalah mengubah matriks \(A\) menjadi matriks row echelon form. Kedua, \(L\) harus menjadi matriks segitiga bawah yang mereduksi ke-\(l\) dengan mengikuti operasi baris yang sama yag menghasilkan \(U\). Kita dapat menggunakan algoritma Doolittle untuk menghasilkan \(L\), di mana nilai setiap entri dalam matriks segitiga bawah merupakan pengali yang digunakan untuk menghilangkan entri yang sesuai untuk setiap proses row replacement.

Pada praktiknya, proses eliminasi Gauss untuk memperoleh matriks \(U\) kadang menghasilkan nol di kolom pivotnya. Kondisi tersebut mengharuskan kita untuk melakukan proses row swapping atau pertukaran baris (biasanya dengan baris bawahnya) untuk pivot bukan nol. Jika proses tersebut berhasil dilakukan bisa jadi matriks \(A\) mungkin setara dengan matriks LU, tetapi tidak sama dalam hal urutan nilai pada tiap barisnya. Agar kita dapat memperoleh hasil yang sama (matriks A sama dengan matriks LU), diperlukan matriks ketiga, \(P\). Matriks ini merupakan matriks identitas dengan ukuran sama dengan matriks \(A\). Jika pertukaran baris dilakukan selama proses pembentukan matriks \(U\), maka pertukaran baris yang sama juga akan diimplemenntasikan pada matriks \(P\). oleh karena itu, dalam praktiknya matriks \(A=PLU\) dan perkalian dengan matriks \(P\) berfungsi untuk mengembalikan urutan baris.

Contoh 6.6 Selesaikan sistem persamaan linier berikut menggunakan faktorisasi LU

\[ \begin{matrix} x_1+x_2+3x_4=4 \\ 2x_1+x_2-x_3+x_4=1 \\ 3x_1-x_2-x_3+2x_4=-3 \\ -x_1-2x_2+3x_3-x_4=4 \end{matrix} \]

Jawab:

Nayatakan sistem persamaan tersebut ke dalam bentuk matriks \(Ax=b\).

\[\begin{equation*} Ux= \begin{bmatrix} 1 & 1 & 0 & 3 \\[0.3em] 2 & 1 & -1 & 1 \\[0.3em] 3 & -1 & -1 & 2 \\[0.3em] -1 & 2 & 3 & -1 \end{bmatrix} \begin{bmatrix} x_1 \\[0.3em] x_2 \\[0.3em] x_3 \\[0.3em] x_4 \end{bmatrix} = \begin{bmatrix} 4 \\[0.3em] 1 \\[0.3em] -3 \\[0.3em] 4 \end{bmatrix} \end{equation*}\]

Lakukan operasi baris elementer pada matriks \(A\) untuk memperoleh matriks \(U\). Urutan operasi baris elementer yang dilakukan adalah sebagai berikut:

• \(\left(B_2-2B_1\right)\to B_2 \to l_{2,1}=2\),
• \(\left(B_3-3B_1\right)\to B_3 \to l_{3,1}=3\),
• \(\left(B_4+B_1\right)\to B_4 \to l_{4,1}=-1\),
• \(\left(B_3-4B_2\right)\to B_3\to l_{3,2}=4\),
• \(\left(B_4+3B_2\right)\to B_4 \to l_{4,2}=-3\),
• \(l_{4,3}=0\)

Simpan pengali tiap tahapan pada masing-masing elemen matriks \(L\). Hasil operasi tersebut akan menghasilkan matriks triangular \(U\).

\[ U= \begin{bmatrix} 1 & 1 & 0 & 3 \\[0.3em] 0 & -1 & -1 & -5 \\[0.3em] 0 & 0 & 3 & 13 \\[0.3em] 0 & 0 & 0 & -13 \end{bmatrix} \]

Untuk matriks \(L\) sebagai berikut:

\[ L= \begin{bmatrix} 1 & 0 & 0 & 0 \\[0.3em] 2 & 1 & 0 & 0 \\[0.3em] 3 & 4 & 1 & 0 \\[0.3em] -1 & -3 & 0 & 1 \end{bmatrix} \]

Karena pada proses operasi baris elementer tidak terdapat operasi pertukaran baris, maka matriks \(P\) tidak mengalami perubahan:

\[ P= \begin{bmatrix} 1 & 0 & 0 & 0 \\[0.3em] 0 & 1 & 0 & 0 \\[0.3em] 0 & 0 & 1 & 0 \\[0.3em] 0 & 0 & 0 & 1 \end{bmatrix} \]

Lakukan operasi forward substitution menggunakan Persamaan (6.26).

\[\begin{equation*} \begin{bmatrix} 1 & 0 & 0 & 0 \\[0.3em] 2 & 1 & 0 & 0 \\[0.3em] 3 & 4 & 1 & 0 \\[0.3em] -1 & -3 & 0 & 1 \end{bmatrix} \begin{bmatrix} t_1 \\[0.3em] t_2 \\[0.3em] t_3 \\[0.3em] t_4 \end{bmatrix} = \begin{bmatrix} 4 \\[0.3em] 1 \\[0.3em] -3 \\[0.3em] 4 \end{bmatrix} \end{equation*}\]

Berdasarkan hasil perhitungan diperoleh nilai vektor \(t\).

\[ t_1=4, t_2=1, t_3=-3, t_4=4 \]

Operasi terakhir yang perlu dilakukan untuk memperoleh nilai \(x\) adalah dengan melakukan backward substitution menggunakan nilai vektor \(t\) yang telah dihitung.

\[\begin{equation*} \begin{bmatrix} 1 & 1 & 0 & 3 \\[0.3em] 0 & -1 & -1 & -5 \\[0.3em] 0 & 0 & 3 & 13 \\[0.3em] 0 & 0 & 0 & -13 \end{bmatrix} \begin{bmatrix} x_1 \\[0.3em] x_2 \\[0.3em] x_3 \\[0.3em] x_4 \end{bmatrix} = \begin{bmatrix} 4 \\[0.3em] 7 \\[0.3em] 13 \\[0.3em] -13 \end{bmatrix} \end{equation*}\]

Berdasarkan hasil perhitungan diperoleh nilai \(x\) sebagai berikut:

\[ x_1=-1, x_2=2, x_3=0, x_4=1 \]

Algoritma Dekomposisi LU

1. Masukkan matriks \(A\), dan vektor \(B\) beserta ukurannya \(n\)
2. Lakukan langkah poin ke-4 s/d poin 5 untuk meperoleh matriks \(U\).
3. Untuk baris ke-\(i\) di mana \(i=1\) s/d \(n\), perhatikan apakah nilai \(a_{i,j}\) sama dengan nol.

• Bila iya, lakukan row swapping antara baris ke-\(i\) dan baris ke-\(i+k\leq n\), dimana \(a_{i+k,j}\) tidak sama dengan nol. Bila tidak ada berarti perhitungan tidak bisa dilanjutkan dan proses dihentikan dengan tanpa penyelesaian.
• Bila tidak, lanjutkan.

1. Untuk baris ke-\(j\), dimana \(j=i+1\) s/d \(n\), lakukan operasi baris elementer:

• Hitung \(c=\frac{a_{j,i}}{a_{i,i}}\)
• untuk kolom \(k\), dimana \(k=1\) s/d \(n+1\), hitung \(a_{j,k}=a_{j,k}-c_i.a_{i,k}\)

1. Lakukan langkah poin ke-7 s/d poin 9 untuk memperoleh matriks \(L\)
2. Untuk diagonal matriks \(L\) isikan dengan nilai 1 dan elemen di atas diagonal dengan nilai nol.
3. Untuk elemen di bawah diagonal isikan dengan faktor pengali operasi baris elementer matriks \(U\).
4. Lakukan proses forward substitution menggunakan Persamaan (6.27) untuk memperoleh nilai vektor \(t\).
5. Lakukan backward substituion menggunakan Persamaan (6.16).

Berdasarkan algoritma tersebut, kita dapat menyusun algoritma faktorisasi LU menggunakan R. Berikut adalah sintaks yang digunakan:

lu_solve <- function(a, b=NULL){ m <- nrow(a) n <- ncol(a) piv <- 1 # membentuk matriks identitas P dan L P <- L <- diag(n) # cek elemen diagonal utama apakah bernilai nol for(row_curr in 1:m){ if(piv <= n){ i <- row_curr while(a[i, piv] == 0 && i < m){ i <- i + 1 if(i > m){ i <- row_curr piv <- piv + 1 if(piv > n) return(list(P = P, L = L, U = a)) } } # jika elemen diagonal utama bernilai nol,lakukan row swapping if(i != row_curr){ a <- swap_row(a, i, row_curr) P <- swap_row(P, i, row_curr) } # pembentukan matriks L dan U for(j in row_curr:m) if(j != row_curr){ k <- a[j, piv]/a[row_curr, piv] # matriks U a <- replace_row(a, row_curr, j, -k) # pengisian elemen matriks L L[j, piv] <- k } piv <- piv + 1 } } # penyelesaian persamaan linier if(is.null(b)){ return(list(P = P, L = L, U = a)) }else{ # forward substitution t <- forwardsolve(L, b) # backward substitution x <- backsolve(a, t) return(list(P = P, L = L, U = a, result=x)) } }

Kita dapat menyelesaikan sistem persamaan linier pada Contoh 6.6 menggunakan fungsi yang telah kita buat.

# membuat matriks a dan vektor b a <- matrix(c(1,2,3,-1,1,1,-1,2, 0,-1,-1,3,3,1,2,-1), nrow=4) b <- c(4,1,-3,4) # penyelesaian decomp<-lu_solve(a,b)

Untuk membentuk kembali matriks \(A\), kita dapat mengalikan matriks \(L\), \(U\), dan \(P\).

decomp$L%*%decomp$U%*%decomp$P

Contoh 6.7 Lakukan dekomposisi LU pada matriks berikut dan lakukan pengecekan apakah perkalian hasil dekomposisi matriks akan menghasilkan matriks semula!

\[ \begin{bmatrix} 0 & 1 & -1 \\[0.3em] 1 & 5 & 9 \\[0.3em] 7 & -1 & -5 \end{bmatrix} \]

Jawab:

Lakukan proses dekomposisi menggunakan fungsi lu_solve().

# membentuk matriks a (A <- matrix(c(0, 1, 7, 1, 5, -1, -2, 9, -5), 3))

# dekomposisi lu decomp<-lu_solve(A)

Lakukan pengecekan apakah matriks hasil dekomposisi akan menghasilkan matriks \(A\).

decomp$P %*% decomp$L %*% decomp$U

Fungsi lu() pada Paket Matrix dapat digunakan untuk melakukan dekomposisi LU. Untuk meggunakan fungsi tersebut, kita harus menginstall dan mengaktifkan Paket Matrix.

install.packages("Matrix")

Untuk dapat menggunakannya kita hanya perlu menginputkan matriks kedalam fungsi tersebut. Berikut adalah contoh penerapannya:

Untuk menampilkan hasil dekomposisi, jalankan fungsi expand().

decomp <- Matrix::expand(lum) decomp

Dekomposisi Cholesky memberikan faktorisasi matriks alternatif sehingga \(A = LL^T\), di mana \(L^T\) merupakan transpose konjugat dari matriks \(L\). Dalam kasus ini, penulis hanya bekerja dengan matriks rill dengan nilai rill dan bagian imajiner nol. Jadi untuk tujuan sub-chapter ini, matriks \(L^T\) hanyalah transpose dari matriks \(L\).

Seperti dekomposisi LU, dekomposisi Cholesky dapat digunakan untuk menyelesaikan sistem persamaan linier. Kelebihannya, Menemukan dekomposisi Cholesky jauh lebih cepat daripada dekomposisi LU. Namun, dekomposisi ini hanya terbatas pada matriks tertentu saja. Dekomposisi Cholesky hanya dapat digunakan pada matriks definit positif dan simetris. Matriks simeteris merupakan matriks yang nilai di atas dan di bawah diagonalnya simetris atau sama; secara matematis, untuk semua \(i\) dan \(j\) pada matriks \(A\), \(a_{i;j}=a_{j;i}\). Definit positif berarti bahwa setiap entri pivot (nilai elemen diagonal utama) selelu bernilai positif. Selain itu, untuk matriks definit positif, hubungan \(xAx>0\) untuk semua vektor, \(x\).

Karena \(L^∗\) transpose dari matriks \(L\), maka \(l^{T}_{i,j} = l_{j,i}\) untuk semua nilai \(i\) dan \(j\). Tanpa kendala (constraint) ini, dekomposisi Cholesky akan mirip dekomposisi LU. Tetapi dengan kendala ini, nilai elemen matriks \(L\) dan \(L^T\) harus dipilih dengan cermat sehingga hubungan \(A = LL^T\) berlaku. Bentuk dekomposisi Cholesky disajikan pada Persamaan (6.29).

\[\begin{equation} \begin{bmatrix} a_{1,1} & a_{1,2} & a_{1,3} &\cdots& a_{1,m} \\[0.3em] a_{2,1} & a_{2,2} & a_{2,3} &\cdots& a_{2,m} \\[0.3em] a_{3,1} & a_{3,2} & a_{3,3} &\cdots& a_{3,m} \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] a_{m,1} & a_{m,2} & a_{m,3} &\cdots& a_{m,m} \end{bmatrix} = \begin{bmatrix} l_{1,1} & 0 & 0 &\cdots& 0 \\[0.3em] l_{2,1} & l_{2,2} & 0 &\cdots& 0 \\[0.3em] l_{3,1} & l_{3,2} & l_{3,3} &\cdots& 0 \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] l_{m,1} & l_{m,2} & l_{m,3} &\cdots& l_{m,m} \end{bmatrix} \begin{bmatrix} l_{1,1} & l_{1,2} & l_{1,3} &\cdots& l_{1,m} \\[0.3em] 0 & l_{2,2} & l_{2,3} &\cdots& l_{2,m} \\[0.3em] 0 & 0 & l_{3,3} &\cdots& l_{3,m} \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] 0 & 0 & 0 &\cdots& l_{m,m} \end{bmatrix} \tag{6.29} \end{equation}\]

Untuk setiap elemen matriks \(A\) memiliki hubungan yang dituliskan pada Persamaan (6.30).

\[\begin{equation} a_{i,j}=\sum_{k=1}^mL_{i,k}L_{k,j} \tag{6.30} \end{equation}\]

Berdasarkan Persamaan (6.29), sejumlah nilai elemen \(L_{i,k}\) dan \(L_{k,j}\) adalah nol. Nilai tiap elemen diagonal utama yang tidak bernilai nol dihitung menggunakan Persamaan (6.31).

\[\begin{equation} l_{i,i}=\sqrt{\left(a_{i,i}-\sum_{k=1}^{i-1}l_{i,k}^2\right)} \tag{6.31} \end{equation}\]

Elemen diagonal dihitung menggunakan Persamaan (6.32)

\[\begin{equation} l_{i,j}=\frac{1}{l_{i,i}}\left(a_{i,j}-\sum_{k=1}^{i-1}l_{i,k}l_{j,k}\right) \tag{6.32} \end{equation}\]

Algoritma Dekomposisi Cholesky

1. Masukkan matriks \(A\), dan vektor \(B\) beserta ukurannya \(n\).
2. Untuk elemen matriks \(L\), hitung menggunakan Persamaan (6.32).
3. Untuk nilai diagonal utama matriks \(L\), hitung menggunakan Persamaan (6.31).
4. Untuk memperoleh matriks \(L^T\), lakukan transpose pada matriks \(L\).
5. Untuk memperoleh nilai \(x\),

• Hitung vektor \(t\) menggunakan Persamaan (6.25).
• Hitung vektor \(x\) menggunakan Persamaan (6.24), dimana matriks \(U=L^T\).

Berdasarkan algoritma tersebut, kita dapat menyusun fungsi pada R untuk melakukan dekomposisi Cholesky. Fungsi tersebut disajikan pada sintaks berikut:

cholesky_solve <- function(a, b=NULL){ m <- nrow(a) # membentuk matriks L dengan elemen nol L = diag(0,m) # Perhitungan elemen matriks L for(i in 1:m){ for(k in 1:i){ p_sum <- 0 for(j in 1:k) p_sum <- p_sum + L[j,i]*L[j,k] # Pehitungan elemen diagonal utama if(i==k) L[k,i]<-sqrt(a[i,i]-p_sum) else L[k,i]<-(a[k,i]-p_sum)/L[k,k] } } # Perhitungan elemn matriks L* tL <- t(L) # penyelesaian persamaan linier if(is.null(b)){ return(list(L = L, tL = tL, a = a)) }else{ # forward substitution t <- forwardsolve(L, b) # backward substitution x <- backsolve(tL, t) return(list(L = L, tL = tL, a = a, result=x)) } }

Contoh 6.8 Dengan menggunakan fungsi cholesky_solve(), lakukan dekomposisi pada matriks berikut! Lakukan pengecekan pada hasil dekomposisi apakah hasil kali matriks dekomposisi akan menghasilkan matriks semula!

\[ \begin{bmatrix} 9 & -3 & 6 \\[0.3em] -3 & 17 & -10 \\[0.3em] 6 & -10 & 12 \end{bmatrix} \]

Jawab:

Dekomposisi Cholesky menggunakan fungsi cholesky_solve(), disajikan pada sintaks berikut:

# mengecek hasil dekomposisi decomp$tL %*% decomp$L

Fungsi lain yang dapat digunakan untuk melakukan dekomposisi Cholesky adalah menggunakan fungsi chol() pada Paket Matrix. Pada fungsi tersebut, kita hanya perlu menginputkan objek matrik kedalamnya. Berikut adalah contoh penerapan fungsi tersebut menggunakan matriks pada Contoh 6.8.

Penting!!!

Fungsi chol() hanya menampilkan matriks \(L^T\). Untuk menampilkan matriks \(L\), kita perlu melakukan transpose

Terdapat beberapa algoritma lain yang telah dikembangkan untuk melakukan dekomposisi matriks. Pada buku ini hanya akan dijelaskan secara singkat terkait fungsi yang digunakan dalam melakukan dekomposisi matriks. Algoritma yang akan dijelaskan pada sub-chapter ini antara lain: QR, singular value decomposition (SVD), dan dekomposisi eigen. Untuk algoritma lainnya, pembaca dapat membaca buku terkait atau mengecek dokumentasinya pada Paket base.

Dekomposisi QR merupakan dekomposisi yang penting dalam menyelesaikan sistem persamaan linier. Dekomposisi ini juga berperan penting untuk menghitung koefisien regresi dan pengaplikasian algoritma Newton-Raphson.

Untuk memperoleh informasi terkait dekomposisi ini, pembaca dapat mengetikkan sintaks berikut pada R:

Berikut merupakan contoh penerapan fungsi qr() untuk menyelesaikan sistem persamaan linier:

# memperoleh penyelesaian SPL qr.solve(A,b)

Singular value decomposition (SVD) merupakan algoritma faktorisasi matriks yang mendekomposisi matriks segiempat menjadi matriks \(UDV_H\), dimana \(D\) merupakan mmatriks diagonal non negatif, \(U\) dan \(V\) merupakan matriks unitary, dan \(V_H\) merupakan matriks tanspose konjugat dari matriks \(V\). Algoritma ini banyak digunakan dalam analisis principal component.

Pada R, SVD dapat dilakukan menggunakan fungsi svd() dari Paket base. Berikut adalah sintaks untuk memperoleh informasi terkait fungsi tersebut:

Berikut adalah contoh penerapan fungsi svd():

# dekomposisi matriks A svd(A)

Proses umum yang digunakan untuk menemukan nilai eigen dan vektor eigen suatu matriks segiempat dapat dilihat sebagai proses dari dekomposisi eigen. Proses ini akan mendekomposisi matriks menjadi \(VDV^{-1}\), dimana \(D\) merupakan matriks diagonal yang terbentuk dari nilai eigen, dan \(V\) merupakan vektor eigen. Proses dekomposisi ini akan berguna bagi pembaca yang ingin mempelajari principal component analysis.

Fungsi eigen() pada Paket base dapat digunakan untuk melakukan dekomposisi eigen. Untuk mempelajari lebih jauh terkait fungsi ini, pambaca dapat menjalankan sintaks berikut:

Berikut adalah contoh sintaks untuk melakukan dekomposisi eigen:

Pada Chapter 6.5 kita akan membahas penyelesaian persamaan linier dengan menggunakan metode iterasi. Terdapat dua metode iterasi yang akan dibahas yaitu iterasi Jacobi dan Gauss-Seidel.

Metode iterasi dimulai dengan estimasi nilai akhir. Setelah menerapkan beberapa perlakuan pada nilai estimasi, hasil perlakuan selanjutnya menjadi nilai estimasi untuk iterasi berikutnya. Proses tersebut akan berlangsung secara terus-menerus hingga ambang batas dipenuhi. Nilai ambang batas dapat berupa jumlah iterasi maksimum atau selisih antara nilai estimasi baru dan estimasi semula lebih kecil dari suatu nilai toleransi yang ditetapkan.

Jumlah kuadrat merupakan metode yang sering digunakan untuk mengecek apakah selisih nilai estimasi baru terhadap estimasi lama lebih kecil dari nilai toleransi yang ditetapkan. Persamaan (6.33) menampilkan hubungan antara jumlah kuadrat dan nilai toleransi pada proses iterasi.

\[\begin{equation} \sqrt{\sum_{i=1}^n\left(x_i^{n+1}-x_i^{n}\right)^2}<t_0 \tag{6.33} \end{equation}\]

dimana \(x^{n}\) merupakan iterasi ke-\(n\) dari algoritma dan \(t_0\) merupakan nilai toleransi maksimum yang diterima.

Untuk menyelesaikan matriks menggunakan metode iterasi, kita dapat mulai dengan premis terdapat matriks \(A\) dan vektor \(x\) dan b, sehingga \(Ax = b\). Dengan menggunakan metode Jacobi, pertama-tama kita dapat amati bahwa terdapat matriks \(R\) dan \(D\) yang memiliki hubungan \(A = R + D\). Berdasarkan kedua hubungan tersebut, dapat diturunkan operasi matriks melalui persamaan berikut:

\[\begin{equation} Ax=b \tag{6.34} \end{equation}\]

\[\begin{equation} Rx+Dx=b \tag{6.35} \end{equation}\]

\[\begin{equation} Dx=b-Rx \tag{6.36} \end{equation}\]

\[\begin{equation} x=D^{-1}\left(b-Rx\right) \tag{6.37} \end{equation}\]

Persamaan (6.37) merupakan persamaan yang dapat kita gunakan untuk memperoleh nilai \(x\). Jika kita menulis kembali persamaan tersebut, maka kita akan memperoleh persamaan yang digunakan sebagai acuan iterasi Jacobi.

\[\begin{equation} x^{n+1}=D^{-1}\left(b-Rx^{n}\right) \tag{6.38} \end{equation}\]

dimana \(D\) merupakan matriks diagonal dengan nilai elemen diagonal berupa diagonal utama matriks \(A\). Invers dari matriks \(D\) secara sederhana sebagai matriks diagonal sama dengan satu dibagi dengan elemen diagonal utama matriks \(A\). Matriks \(R\) identik dengan matriks \(A\). Namun, diagonal utamanya bernilai nol. Suatu iterasi dikatakan konvergen jika jumlah kuadrat dari vektor \(x^{\left(n+1\right)}\) dan vektor \(x^{\left(n\right)}\) semakin mengecil.

Suatu persamaan linier yang hendak diselesaikan dengan menggunakan metode iterasi Jacobi harus memenuhi syarat nilai elemen diagonal utama matriks harus lebih dominan. Maksudnya adalah nilai absolut diagonal utama matriks harus lebih besar dari jumlah nilai absolut elemen matriks lainnya pada satu kolom.

Contoh 6.9 Selesaikan sistem persamaan berikut menggunakan iterasi Jacobi!

\[\begin{equation*} \begin{bmatrix} 5 & 2 & 3 \\[0.3em] 2 & 7 & 4 \\[0.3em] 1 & 3 & 8 \end{bmatrix} x = \begin{bmatrix} 40 \\[0.3em] 39 \\[0.3em] 55 \end{bmatrix} \end{equation*}\]

Jawab:

Berdasarkan matriks \(A\) (matriks koefisien), kita dapat memastikan bahwa matriks tersebut memiliki nilai dominan pada elemen diagonal utama. Sebagai contoh:

\[ \left|5\right|>\left|2\right|+\left|1\right|\ \ \ \left(kolom\ 1\right) \]

\[ \left|7\right|>\left|2\right|+\left|3\right|\ \ \ \left(kolom\ 2\right) \] Untuk mempermudah proses iterasi, kita akan menggunakan bantuan R untuk melakukan komputasi. Langkah pertama yang perlu dilakukan adalah menyiapkan matriks \(A\), vektor \(b\), dan vektor \(x\) (nilai taksiran awal).

(A <- matrix(c(5,2,1,2,7,3,3,4,8), 3))

Langkah selanjutnya adalah memperoleh invers matriks \(D\).

(Dinv <- diag(1/diag(A)))

Persiapan terakhir sebelum iterasi dilakukan adalah menyiapkan matriks \(R\).

Iterasi selanjutnya dilakukan menggunakan Persamaan (6.38).

iterasi 1

(x1 <- Dinv %*% (b-R%*%x))

iterasi 2

(x2 <- Dinv %*% (b-R%*%x1))

iterasi 3

(x3 <- Dinv %*% (b-R%*%x2))

Selama proses iterasi,jumlah akar jumlah kuadrat dihitung. Sebagai contoh berikut disajikan akar jumlah kuadrat pada iterasi ke-3:

Selama proses iterasi nilai tersebut terus mengecil. Iterasi dihantikan jika nilai akar jumlah kuadrat tersebut lebih kecil dari nilai toleransi. Pada contoh ini digunakan nilai toleransi \(10^{-7}\).

Proses iterasi berlangsung sampai dengan iterasi ke-62 dengan nilai \(x\) akhir sebagai berikut:

\[ x = \begin{bmatrix} 4 \\[0.3em] 1 \\[0.3em] 6 \end{bmatrix} \]

Algoritma Iterasi Jacobi

1. Masukkan matriks \(A\), dan vektor \(B\) beserta ukurannya \(n\).
2. Hitung invers matriks \(D\), dimana nilai invernya merupakan matriks diagonal dari satu per diagonal utama matriks \(A\).
3. Hitung matriks \(R\), dimana \(R\) merupakan selisih matriks \(A\) dikurangi dengan matriks diagonal dengan entri dari diagonal utama matriks \(A\).
4. Tetapkan vektor \(x\) estimasi.
5. Tetapkan nilai toleransi maksimum yang dapat diterima.
6. Lakukan iterasi menggunakan Persamaan (6.38).
7. Hitung akar jumlah kuadrat dari vektor \(x^{n+1}\) dan vektor \(x^n\).
8. Jadikan nilai \(x^{n+1}\) sebagai nilai taksiran \(x\) untuk iterasi berikutnya.
9. Hentikan proses iterasi jika telah memenuhi syarat yang ditampilkan pada Persamaan (6.33).

Berdasarkan algoritma tersebut, kita dapat menyusun fungsi sebuah fungsi untuk melakukan iterasi Jacobi. Berikut sintaks yang digunakan:

Berikut adalah penerpan fungsi jacobi() tersebut:

Contoh 6.10 Selesaikan sistem persamaan berikut menggunakan fungsi jacobi()

\[ \begin{bmatrix} 27 & 6 & -1 \\[0.3em] 6 & 15 & 2 \\[0.3em] 1 & 1 & 54 \end{bmatrix} x = \begin{bmatrix} 85 \\[0.3em] 72 \\[0.3em] 110 \end{bmatrix} \]

Jawab:

Matriks \(A\) (matriks koefisien) berdasarkan sistem persamaan linier tersebut telah memenuhi syarat dari algoritma Jacobi (nilai diagonal utama dominan dibanding nilai lainnya pada satu kolom). Penyelesaian sistem persamaan tersebut, sebagai berikut:

A <- matrix(c(27,6,1,6,15,1,-1,2,54), 3) b <- c(85,72,110) jacobi(A,b)

Nilai vektor \(x\) sesungguhnya dapat diperoleh menggunakan fungsi solve().

Berdasarkan hasil perhitungan, vektor \(x\) hasil iterasi memiliki nilai identik dengan nilai penyelesaian yang sebenarnya.

Perlu diperhatikan dalam penggunaan fungsi jacobi() syarat utama matriks haruslah terpenuhi, seperti: nilai diagonal matriks \(A\) lebih besar dari nilai elemen lainnya pada satu kolom. Selain itu, nilai diagonal matriks \(D\) tidak boleh sama dengan nol agar inver matriks \(D\) dapat diperoleh. Jika syarat-syarat tersebut terpenuhi, maka metode Jacobi dapat diterapkan. Jika tidak terpenuhi, maka penyelesaian yang konvergen mungkin masih dapat diperoleh meskipun penulis tidak dapat menjamin hal tersebut dapat terjadi.

Metode iterasi Gauss-Seidel melakukan dekomposisi pada matriks \(A\) menjadi matriks segitiga atas \(U\) dan matriks segitiga bawah \(L\). Dekomposisi ini tidak sama dengan dekomposisi LU pada Chapter 6.4.1. Matriks \(U\) pada metode Gauss-Seidel merupakan elemen (entri) matriks \(A\) pada bagian atas diagonal utama, sedangkan matriks \(L\) merupakan elemen diagonal utama dan bagian bawah diagonal utama matriks \(A\). Elemen selain yang penulis sebutkan pada kedua matriks tersebut akan bernilai nol. Persamaan iterasi Gauss-Seidel ditampilkan pada Persamaan (6.39).

\[\begin{equation} x^{n+1}=L^{-1}\left(b-Ux^{n}\right) \tag{6.39} \end{equation}\]

Syarat agar suatu sistem persamaan linier dapat diselesaikan menggunakan metode Gauss-Seidel adalah matriks harus memiliki nilai diagonal utama yang dominan. Maksudnya, nilai absolut diagonal utama lebih besar dari jumlah nilai absolut elemen lainnya dalam satu kolom. Jika syarat ini tidak terpenuhi maka metode ini tidak akan memperoleh penyelesaian yang konvergen.

Contoh 6.11 Selesaikan sistem persamaan pada Contoh 6.10 menggunakan iterasi Gauss-Seidel!

Jawab:

Kita akan kembali menggunakan bantuan R untuk melakukan kalkulasi pada proses iterasi Gauss-Seidel. Kita telah melakukan pengecekan pada sistem persamaan linier pada contoh tersebut dan menghasilkan kesimpulan bahwa persamaan linier tersebut dapat diselesaikan dengan metode Gauss-Seidel. Langkah selanjutnya adalah membentuk matriks \(L\) dan \(U\).

# membentuk matriks U dan L dari matriks A (L <- U <- A)

Selanjutya lakukan invers terhadap matriks \(L\) menggunakn fungsi solve().

Tetapkan nilai estimasi awal dan nilai toleransi yang dikehendaki. Nilai toleransi pada proses ini ditetapkan sebesar \(10^-7\).

# tebakan awal nilai x (x <- rep(0, length(b)))

Lakukan iterasi menggunakan Persamaan (6.39).

Iterasi 1

(x1 <- Linv %*% (b - U %*% x))

# akar jumlah kuadrat sqrt(sum(x1-x)^2)

Iterasi 2

(x2 <- Linv %*% (b - U %*% x1))

# akar jumlah kuadrat sqrt(sum(x2-x1)^2)

Iterasi terus dilakukan sampai dengan nilai akar jumlah kuadrat lebih kecil dari nilai toleransi. Setelah iterasi ke-7 diperoleh nilai vektor \(x\) sebesar:

\[ x = \begin{bmatrix} 2,425476 \\[0.3em] 3,573016 \\[0.3em] 1,925954 \end{bmatrix} \]

Algoritma Iterasi Gauss-Seidel

1. Masukkan matriks \(A\), dan vektor \(B\) beserta ukurannya \(n\).
2. Lakukan dekomposisi LU, dimana matriks \(L\) merupakan matriks segitiga bawah dengan nilai entri diagonal utama matriks \(A\) dan bagian bawah diagonalnya dan matriks \(U\) merupakan matriks segitiga atas dengan entri berasal dari elemen atas diagonal utama matriks \(A\). Isi elemen lain yang tidak disebut pada kedua matriks tersebut dengan nol.
3. Tetapkan vektor \(x\) estimasi.
4. Tetapkan nilai toleransi maksimum yang dapat diterima.
5. Lakukan iterasi menggunakan Persamaan (6.39).
6. Hitung akar jumlah kuadrat dari vektor \(x^{n+1}\) dan vektor \(x^n\).
7. Jadikan nilai \(x^{n+1}\) sebagai nilai taksiran \(x\) untuk iterasi berikutnya.
8. Hentikan proses iterasi jika telah memenuhi syarat yang ditampilkan pada Persamaan (6.33).

Berdasarkan algoritma tersebut, kita dapat menyusun fungsi sebuah fungsi untuk melakukan iterasi Gauss-Seidel. Berikut sintaks yang digunakan:

Contoh 6.12 Selesaikan sistem persamaan pada Contoh 6.10 menggunakan fungsi gauss_seidel()!

Jawab:

Penyelesaiansistem persamaan linier tersebut menggunakan fungsi gauss_seidel() disajikan pada sintaks berikut:

Aljabar linier banyak diaplikasikan baik dalam bidang engineering, fisika, sampai dengan statistika. Pada sub-chapter ini penulis akan menjelaskan penerapan aljabar linier pada metode kuadrat terkecil dan aliran massa dalam reaktor. Untuk penerapan lainnya pembaca dapat membaca buku lainnya terkait aljabar linier.

Metode kuadrat terkecil merupakan salah satu aplikasi penerapan aljabar linier yang paling populer. Intuisi dibalik metode ini adalah bagaimana kita meminimalkan jarak antara sejumlah titik dengan garis regresi. Misalkan kita menggambarkan scatterplot antara dua buah variabel. Pola yang terbentuk dari plot tersebut adalah terjadi korelasi positif antara variabel pada sumbu \(x\) dan sumbu \(y\). Kita ingin menggambarkan garis regresi terbaik yang dapat menangkap seluruh pola tersebut. Garis regresi terbaik terjadi ketika jumlah kuadrat jarak antara titik observasi dan garis regresi yang terbentuk seminimal mungkin.

Untuk lebih memahaminya kita akan melakukan latihan menggunakan dataset trees yang berisi data hasil pengukuran kayu dari pohon yang ditebang. Pada dataset ini terdapat 31 observasi dan 3 buah kolom. Keterangan dari ketiga buah kolom tersebut adalah sebagai berikut:

• Girth: diameter pohon dalam satuan inch.
• Height: tinggi pohon dalam satuan feet.
• Volume: volume kayu dalam satuan cubic feet.

Untuk mengecek 6 observasi pertama dan struktur data, jalankan sintaks berikut:

Scatterplot matriks sangat bagus untuk mengecek korelasi antar variabel dalam dataset tersebut. Berikut adalah sintaks untuk membuatnya:

(#fig:trees, LUfig)Scatterplot matriks dataset trees

Kita ingin membuat sebuah model linier untuk memprediksi Volume kayu berdasarkan variabel Girth dan Heiht atau volume sebagia fungsi dari variabel Girth dan Heiht. Kita dapat menuliskan relasi antara variabel volume sebagai fungsi dari variabel Girth dan Heiht menggunakan Persamaan (6.40).

\[\begin{equation} Volume=\beta_{girth} Girth+\beta_{height}Height+\beta_0 \tag{6.40} \end{equation}\]

dimana \(\beta_0\) merupakan intersep persamaan regresi linier dan nilai \(\beta\) lainnya merupakan koefisien dari variabel Girth dan Heiht. Variabel Volume disebut sebagai variabel respon, sedangkan variabel Girth dan Heiht disebut sebagai variabel prediktor.

Metode kuadrat terkecil berusaha memperoleh seluruh koefisien variabel dan intersep dari persamaan regresi linier. Berdasarkan yang telah penulis jelaskan garis regresi terbaik adalah garis yang memiliki nilai kuadrat terkecil jarak antara titik observasi dan garis regresi. Dasar dari metode kuadrat terkecil merupakan persamaan yang relatif sederhana yang ditunjukkan pada Persamaan (6.41).

\[\begin{equation} A^{T}A=A^{T}b \tag{6.41} \end{equation}\]

dimana \(b\) merupakan vektor dari variabel respon (Volume) dan matrik \(A\) merupakan matriks variabel prediktor (variabel Girth dan Heiht).

Untuk menginputkan intercept kedalam persamaan linier kita perlu menmabhakan satu kolom di awal matriks \(A\) yang berisi nilai 1. Berikut adalah sintaks yang digunakan untuk membentuk matriks \(A\):

Langkah selanjutnya adalah membentuk matriks \(b\). Berikut adalah sintaks yang digunakan:

resp<- trees$Volume head(resp)

Untuk memperoleh koefisien \(\beta\), kita dapat mencarinya dengan cara menyelesaikan Persamaan (6.41). Berikut adalah sintaks yang digunakan:

Berdasarkan hasil yang diperoleh, persamaan linier yang terbentuk disajikan pada Persamaan (6.42).

\[\begin{equation} Volume=4.7081605 Girth + 0.3392512 Height -57.9876589 \tag{6.42} \end{equation}\]

Pembaca juga dapat menggunakan fungsi lain untuk memperoleh nilai koefisien tersebut, seperti: lu_solve()dansolve(). untuk fungsi jacobi() dan gauss_seidel(), kita harus pastikan syarat-syarat terkait metode tersebut. Berikut adalah contoh penyelesaian menggunakan sintaks lainnya:

# metode LU lu_solve(A,b)

# fungsi solve() solve(A,b)

R juga menyediakan fungsi untuk membentuk model regresi linier. Fungsi yang digunakan adalah lm(). Berikut sintaks yang digunakan untuk membentuk model linier menggunakan fungsi lm():

lm(Volume~Girth+Height, data=trees)

Pada sub-chapter ini penulis akan memberikan penerapan aljabar linier untuk menghitung konsentrasi suatu zat atau parameter lingkungan dalam reaktor yang saling terhubung. Pada contoh kasus kali ini diasumsikan terdapat lima buah reaktor yang saling terhubung satu sama lain sesuai Gambar 6.2. Debit air (\(\frac{m^{3}}{detik}\)) dan konsentrasi zat pencemar (\(\frac{mg}{m^3}\)) disajikan pula diagram alir tersebut. Diasumsikan kelima buah reaktor tersebut dalam kondisi steady dan volume reaktor diasumsikan sama. Kesetimbangan massa persatuan waktu dalam kondisi steady disajikan pada Persamaan (6.43).

Gambar 6.2: Aliran massa dalam reaktor.

\[\begin{equation} m_in=m_out \tag{6.43} \end{equation}\]

\[\begin{equation} Q_{in}C_{in}= Q_{out}C_{out} \tag{6.44} \end{equation}\]

Berdasarkan Gambar 6.2, dapat dibentuk lima buah sistem persamaan linier. Persamaan linier yang terbentuk disajikan sebagai berikut:

\[ \begin{matrix} 6c_{1}-c_{3}=50 \\ -3c_{1}+3c_{2}=0 \\ -c_{2}+9c_{3}=160 \\ -c_{2}-8c_{3}+11c_{4}-2c_{5}=0 \\ -3c_{1}-c_{2}+4c_{5}=0 \end{matrix} \]

Untuk menyelesaiakan sistem persamaan linier tersebut dan memperoleh nilai \(c\) dari masing-masing reaktor, kiat perlu mengubahnya dulu kedalam bentuk matriks \(Ax=b\). Berikut adalah matriks yang terbentuk:

\[\begin{equation*} \begin{bmatrix} 6 & 0 & -1 & 0 & 0 \\[0.3em] -3 & 3 & 0 & 0 & 0 \\[0.3em] 0 & -1 & 9 & 0 & 0 \\[0.3em] 0 & -1 & -8 & 11 & -2 \\[0.3em] -3 & -1 & 0 & 0 & 4 \end{bmatrix} c = \begin{bmatrix} 50 \\[0.3em] 0 \\[0.3em] 160 \\[0.3em] 0 \\[0.3em] 0 \end{bmatrix} \end{equation*}\]

Kita akan menyelesaikannya dengan menggunakan metode elminasi Gauss-Jordan, dekomposisi LU, iterasi Jacobi, dan iterasi Gauss-Seidel. Untuk dapat menyelesaikannya menggunakan metode-metode tersebut pada R, kita perlu membentuk matriksnya terlebih dahulu:

(A <- matrix(c(6,-3,0,0,-3, 0,3,-1,-1,-1, -1,0,9,-8,0, 0,0,0,11,0, 0,0,0,-2,4),nrow=5))

Metode Eliminasi Gauss-Jordan

Metode Dekomposisi LU

Metode Iterasi Jacobi

Metode Iterasi Gauss-Seidel

gauss_seidel(A,b, maxiter=200)

Berdasarkan seluruh metode tersebut, diperoleh konsentrasi zat pencemar pada masing-masing reaktor adalah sebagai berikut:

\[ \begin{matrix} c_{1}=11,50943 \frac{mg}{m^3} \\ c_{2}=11,50943 \frac{mg}{m^3} \\ c_{3}=19,05660 \frac{mg}{m^3} \\ c_{4}=16,99828 \frac{mg}{m^3} \\ c_{5}=11.50943 \frac{mg}{m^3} \end{matrix} \]

1. Bloomfield, V.A. 2014. Using R for Numerical Analysis in Science and Engineering. CRC Press
2. Howard, J.P. 2017. Computational Methods for Numerical Analysis with R. CRC Press.
3. Kreyszig, E. 2011. Advanced Engineering Mathematics, 10th Edition. John Wiley & Sons.
4. Primartha, R. 2018. Belajar Machine Learning Teori dan Praktik. Penerbit Informatika : Bandung.
5. Sanjaya, M. 2015. Metode Numerik Berbasis Phython. Penerbit Gava Media: Yogyakarta.
6. Suparno, S. 2008. Komputasi untuk Sains dan Teknik Edisi II. Departemen Fisika-FMIPA Universitas Indonesia.

1. Selesaikan sistem persamaan linier berikut menggunakan eliminasi Gauss!

\[ \begin{matrix} -4x+4y=-1 \\ -2x+2y-3z=-3 \\ 3x+1y-3z=-3 \end{matrix} \]

1. Carilah penyelesaian dari sistem persamaan linier soal no.1 menggunakan algoritma dekomposisi LU!
2. Tunjukan 5 iterasi pertama sistem persamaan linier berikut menggunakan algoritma Jacobi dan Gauss-Seidel!

\[ \begin{matrix} 3x+2y-1z=-3 \\ -3x-3y-3z=9 \\ 1y-1z=-1 \end{matrix} \]

1. Gunakan fungsi jacobi() dan gauss_seidel() untuk menyelesaikan sistem persamaan linier pada soal no.4 dan tentukan metode mana yang paling cepat memperoleh penyelesaian? (petunjuk: gunakan fungsi system.time() dan jumlah iterasi yang diperlukan untuk memperoleh hasil yang konvergen)
2. Apakah yang terjadi jika kita menginputkan matriks segiempat \(A\) kedalam fungsi solve() dan apa yang akan terjadi jika selanjutnya argumen pada fungsi tersebut juga menyertakan vektor \(b\)?
3. Dengan menggunakan dataset mtcars buatlah persamaan linier variabel mpg sebagai fungsi dari variabel wt, hp, dan qsec menggunakan algoritma dekomposisi LU?

Persamaan non-linier dapat diartikan sebagai persamaan yang tidak mengandung syarat seperti persamaan linier, sehingga persamaan non-linier dapat merupakan:

1. Persamaan yang memiliki pangkat selain satu (misal: \(x^2\))
2. Persamaan yang mempunyai produk dua variabel (misal: \(xy\))

Dalam penyelesaian persamaan non-linier diperlukan akar-akar persamaan non-linier, dimana akar sebuah persamaan non-linier \(f\left(x\right)=0\) merupakan nilai \(x\) yang menyebabkan nilai \(f\left(x\right)\) sama dengan nol. Dalam hal ini dapat disimpulkan bahwa akar-akar penyelesaian persamaan non-linier merupakan titik potong antara kurva \(f\left(x\right)\) dengan sumbu \(x\). Ilustrasi penjelasan tersebut ditampilkan pada Gambar 7.1.

Gambar 7.1: Penyelesaian persamaan non-linier.

Contoh sederhana dari penentuan akar persamaan non-linier adalah penentuan akar persamaan kuadratik. Secara analitik penentuan akar persamaan kuadratik dapat dilakukan menggunakan Persamaan (7.1).

\[\begin{equation} x_{1,2}=\frac{-b\pm\sqrt{b^2-4a}}{2a} \tag{7.1} \end{equation}\]

Untuk masalah yang lebih rumit, penyelesaian analitik sudah tidak mungkin dilakukan. Metode numerik dapat digunakan untuk menyelesaikan masalah yang lebih kompleks. Untuk mengetahui apakah suatu persamaan non-linier memiliki akar-akar penyelesaian atau tidak, diperlukan analisa menggunakan Teorema berikut:

Teorema 7.1 (root) Suatu range x=[a,b] mempunyai akar bila f(a) dan f(b) berlawanan tanda atau memenuhi f(a).f(b)<0

Untuk memahami teorema tersebut perhatikan ilustrasi pada Gambar 7.2.

Gambar 7.2: Ilustrasi teorema Bolzano.

Pada Chapter 7 ini, akan dilakukan sejumlah pembahasan antara lain:

• penentuan akar persamaan dengan metode tertutup
• penentuan akar persamaan dengan metode terbuka
• fungsi-fungsi R untuk mementukan akar persamaan non-linier
• studi kasus

Metode tertutup disebut juga metode bracketing. Disebut sebagai metode tertutup karena dalam pencarian akar-akar persamaan non-linier dilakukan dalam suatu selang \(\left[a,b \right]\).

Penyelesaian persamaan non-linier menggunakan metode tabel dilakukan dengan membagi persamaan menjadi beberapa area, dimana untuk \(x=\left[a,b \right]\) dibagi sebanyak \(N\) bagian dan pada masing-masing bagian dihitung nilai \(f\left(x \right)\) sehingga diperoleh nilai \(f\left(x \right)\) pada setian \(N\) bagian.

Bila nilai \(f\left(x_k \right)=0\) atau mendekati nol, dimana \(a \le k \le b\), maka dikatakan bahwa \(x_k\) adalah penyelesaian persamaan \(f\left(x \right)\). Bila tidak ditemukan, dicari nilai \(f\left(x_k \right)\) dan \(f\left(x_{k+1} \right)\) yang berlawanan tanda. Bila tidak ditemukan, maka persamaan tersebut dapat dikatakan tidak mempunyai akar untuk rentang \(\left[a,b \right]\).

Bila akar persamaan tidak ditemukan, maka ada dua kemungkinan untuk menentukan akar persamaan, yaitu:

1. Akar persamaan ditentukan oleh nilai mana yang lebih dekat. Bila \(f\left(x_k\right)\le f\left(x_{k+1}\right)\), maka akarnya \(x_k\). Bila \(f\left(x_{k+1}\right)\le f\left(x_{k}\right)\), maka akarnya \(x_{k+1}\).
2. Perlu dicari lagi menggunakan rentang \(x=\left[x_{k}, x_{k+1} \right]\).

Secara grafis penyelesaian persamaan non-linier menggunakan metode table disajikan pada Gambar 7.3.

Gambar 7.3: Ilustrasi metode tabel.

Algoritma Metode Tabel

1. Definisikan fungsi \(f\left(x \right)\)
2. Tentukan rentang untuk \(x\) yang berupa batas bawah \(a\) dan batas atas \(b\).
3. Tentukan jumlah pembagi \(N\)
4. Hitung step pembagi

\[\begin{equation} h=\frac{b+a}{N} \tag{7.2} \end{equation}\]

1. Untuk \(i=0\) s/d \(N\), hitung:

\[\begin{equation} x_i=a+i.h \tag{7.3} \end{equation}\]

\[\begin{equation} y_i=f\left(x_i \right) \tag{7.4} \end{equation}\]

1. Untuk \(i=0\) s/d \(N\), dimana

• Bila \(f\left(x \right)=0\), maka akarnya \(x_k\)

• Bila \(f\left(a \right) f\left(b \right) <0\), maka:

  • \(f\left(x_k\right)\le f\left(x_{k+1}\right)\), maka akarnya \(x_k\)
  • Bila tida, \(x_{k+1}\) adalah penyelesaian atau dapat dikatakan penyelesaian berada diantara \(x_k\) dan \(x_{k+1}\).

Kita dapat membuat suatu fungsi pada R untuk melakukan proses iterasi pada metode Tabel. Fungsi root_table() akan melakukan iterasi berdasarkan step algoritma 1 sampai 5. Berikut adalah sintaks yang digunakan:

Contoh 7.1 Carilah akar persamaan \(f\left(x \right)=x+e^{x}\) pada rentang \(x=\left[-1,0 \right]\)?

Jawab:

Sebagai permulaan, jumlah pembagi yang digunakan adalah \(N=10\). Dengan menggunakan fungsi root_table() diperoleh hasil yang disajikan pada Tabel 7.1.

tabel <- root_table(f=function(x){x+exp(x)}, a=-1, b=0, N=10)

Berdasarkan Tabel 7.1 diperoleh penyelesaian di antara \(-0,6\) dan \(-0,5\) dengan nilai \(f\left(x \right)\) masing-masing sebesar \(-0,0512\) dan \(-0,1065\), sehingga dapat diambil penyelesaian \(x=-0,6\). Kita dapat terus melakukan iterasi sampai memperoleh nilai \(f\left(x \right)\) < nilai toleransi dengan terus merubah rentang yang diberikan. Iterasi berikutnya dengan nilai pembagi sama dan rentang nilai \(x=\left[-0,6;-0,5\right]\) diperoleh nilai \(x=-0,57\) dan \(f\left(x \right)=0,00447\).

Untuk melihat gambaran lokasi akar, kita dapat pulang mengeplotkan data menggunakan fungsi plot. Berikut adalah fungsi yang digunakan:

Gambar 7.4: Plot fungsi x+exp(x) pada rentang -1 sampai 0.

Untuk mengetahui lokasi akar dengan lebih jelas, kita dapat memperkecil lagi rentang nilai yang dimasukkan dalam fungsi curve().

Metode tabel pada dasarnya memiliki kelemahan yaitu cukup sulit untuk memdapatkan error penyelesaian yang cukup kecil, sehingga metode ini jarang sekali digunakan untuk menyelesaikan persamaan non-linier. Namun, metode ini cukup baik digunakan dalam menentukan area penyelesaian sehingga dapat dijadikan acuan metode lain yang lebih baik.

Prinsip metode bagi dua adalah mengurung akar fungsi pada interval \(x=\left[a,b \right]\) atau pada nilai \(x\) batas bawah \(a\) dan batas atas \(b\). Selanjutnya interval tersebut terus menerus dibagi 2 hingga sekecil mungkin, sehingga nilai hampiran yang dicari dapat ditentukan dengan tingkat toleransi tertentu. Untuk lebih memahami metode biseksi, perhatikan visualisasi pada Gambar 7.5.

Gambar 7.5: Ilustrasi metode biseksi.

Metode biseksi merupakan metode yang paling mudah dan paling sederhana dibanding metode lainnya. Adapun sifat metode ini antara lain:

1. Konvergensi lambat
2. Caranya mudah
3. Tidak dapat digunakan untuk mencari akar imaginer
4. Hanya dapat mencari satu akar pada satu siklus.

Algoritma Metode Biseksi

1. Definisikan fungsi \(f\left(x \right)\)
2. Tentukan rentang untuk \(x\) yang berupa batas bawah \(a\) dan batas atas \(b\).
3. Tentukan nilai toleransi \(e\) dan iterasi maksimum \(N\)
4. Hitung \(f\left(a \right)\) dan \(f\left(b \right)\)
5. Hitung:

\[\begin{equation} x=\frac{a+b}{2} \tag{7.5} \end{equation}\]

1. Hitung \(f\left(x \right)\)
2. Bila \(f\left(x \right).f\left(a \right)<0\), maka \(b=x\) dan \(f\left(b \right)=f\left(x \right)\). Bila tidak, \(a=x\) dan \(f\left(a \right)=f\left(x \right)\)
3. Bila \(\left|b-a \right|<e\) atau iterasi maksimum maka proses dihentikan dan didapatkan akar=\(x\), dan bila tidak ulangi langkah 6.
4. Jika sudah diperoleh nilai dibawah nilai toleransi, nilai akar selanjutnya dihitung berdasarkan Persamaan (7.5) dengan nilai \(a\) dan \(b\) merupakan nilai baru yang diperoleh dari proses iterasi.

Berdasarkan algoritma tersebut, kita dapat menyusun suatu fungsi pada R yang dapat digunakan untuk melakukan iterasi tersebut. Fungsi root_bisection() merupakan fungsi yang telah penulis susun untuk melakukan iterasi menggunakan metode biseksi. Berikut adalah sintaks dari fungsi tersebut:

Contoh 7.2 Carilah akar persamaan \(f\left(x \right)=xe^{-x}+1\) pada rentang \(x=\left[-1,0 \right]\) dengan nilai toleransi sebesar \(10^-7\)?

Jawab:

Langkah pertama dalam penghitungan adalah menghitung nilai \(x\) menggunakan Persamaan (7.5).

\[ x=\frac{-1+0}{2}=-0,5 \]

Hitung nilai \(f\left(x \right)\) dan \(f\left(a \right)\).

\[ f\left(x \right)=-0,5.e^{0,5}+1=0,175639 \]

\[ f\left(a \right)=-1.e^{1}+1=-1,71828 \]

Berdasarkan hasil perhitungan diperoleh:

\[ f\left(x \right).f\left(a \right)<0 \]

Sehingga \(b=x\) dan \(f\left(b \right)=f\left(x \right)\). Iterasi dilakukan kembali dengan menggunakan nilai \(b\) tersebut.

Untuk mempersingkat waktu iterasi kita akan menggunakan fungsi root_bisection() pada R. Berikut adalah sintaks yang digunakan:

root_bisection(function(x){x*exp(-x)+1}, a=-1, b=0)

Berdasarkan hasil iterasi diperoleh akar persamaan \(x=-2.980232e-08\) dan iterasi yang diperlukan untuk memperolehnya sebanyak \(24\) iterasi.

Metode regula falsi merupakan metode yang menyerupai metode biseksi, dimana iterasi dilakukan dengan terus melakukan pembaharuan rentang untuk memperoleh akar persamaan. Hal yang membedakan metode ini dengan metode biseksi adalah pencarian akar didasarkan pada slope (kemiringan) dan selisih tinggi dari kedua titik rentang. Titik pendekatan pada metode regula-falsi disajikan pada Persamaan (7.6).

\[\begin{equation} x=\frac{f\left(b\right).a-f\left(a\right).b}{f\left(b\right)-f\left(a\right)} \tag{7.6} \end{equation}\]

Ilustrasi dari metode regula falsi disajikan pada Gambar 7.6.

Gambar 7.6: Ilustrasi metode regula falsi.

Algoritma Metode Regula Falsi

1. Definisikan fungsi \(f\left(x \right)\)
2. Tentukan rentang untuk \(x\) yang berupa batas bawah \(a\) dan batas atas \(b\).
3. Tentukan nilai toleransi \(e\) dan iterasi maksimum \(N\)
4. Hitung \(f\left(a \right)\) dan \(f\left(b \right)\)
5. Untuk iterasi \(i=1\) s/d \(N\)

• Hitung nilai \(x\) berdasarkan Persamaan (7.6)
• Hitung \(f\left(x \right)\)
• Hitung \(error=\left|f\left(x \right) \right|\)
• Jika \(f\left(x \right).f\left(a \right)<0\), maka \(b=x\) dan \(f\left(b \right)=f\left(x \right)\). Jika tidak, \(a=x\) dan \(f\left(a \right)=f\left(x \right)\).

1. Akar persamaan adalah \(x\)

Fungsi root_rf() didasarkan pada langkah-langkah di atas. Sintaks fungsi tersebut adalah sebagai berikut:

Contoh 7.3 Selesaikan persamaan non-linier pada Contoh 7.2 menggunakan metode regula falsi pada rentang \(x=\left[-1,0 \right]\) dengan nilai toleransi sebesar \(10^-7\)?

Jawab:

Langkah pertama penyelesaian dilakukan dengan mencari nilai \(f\left(a \right)\) dan \(f\left(b \right)\).

\[ f\left(a \right)=-1.e^{1}+1=-1,71828 \] \[ f\left(b \right)=0.e^{0}+1=1 \] Hitung nilai \(x\) dan \(f\left(x \right)\).

\[ x=\frac{\left(1.-1\right)-\left(-1,71828.0\right)}{1+1,71828}=-0.36788 \]

\[ f\left(x \right)=-0.36788.e^{0.36788}+1=0.468536 \]

Berdasarkan hasil perhitungan diperoleh:

\[ f\left(x \right).f\left(a \right)<0 \]

Sehingga \(b=x\) dan \(f\left(b \right)=f\left(x \right)\). Iterasi dilakukan kembali dengan menggunakan nilai \(b\) tersebut.

Untuk mempercepat proses iterasi, kita dapat pula menggunakan fungsi root_rf() pada R. Berikut adalah sintaks yang digunakan:

root_rf(function(x){x*exp(-x)+1}, a=-1, b=0)

Berdasarkan hasil perhitungan diperoleh nilai \(x=-0,5671433\) dan jumlah iterasi yang diperlukan adalah \(15\). Jumlah ini lebih sedikit dari jumlah iterasi yang diperlukan pada metode iterasi biseksi yang juga menunjukkan metode ini lebih cepat memperoleh persamaan dibandingkan metode biseksi.

Metode terbuka merupakan metode yang menggunakan satu atau dua tebakan awal yang tidak memerlukan rentang sejumlah nilai. Metode terbuka terdiri dari beberapa jenis yaitu metode iterasi titik tetap, metode Newton-Raphson, dan metode Secant.

Metode iterasi titik tetap merupakan metode penyelesaian persamaan non-linier dengan cara menyelesaikan setiap variabel \(x\) yang ada dalam suatu persamaan dengan sebagian yang lain sehingga diperoleh \(x=g\left(x \right)\) untuk masing-masing variabel \(x\). Sebagai contoh, untuk menyelesaikan persamaan \(x+e^{x}=0\), maka persamaan tersebut perlu diubah menjadi \(x=e^x\) atau \(g\left(x \right)=e^x\). Secara grafis metode ini diilustrasikan seperti Gambar 7.7.

Gambar 7.7: Ilustrasi metode iterasi titik tetap.

Algoritma Metode Iterasi Titik Tetap

1. Definisikan \(f\left(x \right)\) dan \(g\left(x \right)\)
2. Tentukan nilai toleransi \(e\) dan iterasi masimum (N)
3. Tentukan tebakan awal \(x_0\)
4. Untuk iterasi \(i=1\) s/d \(N\) atau \(f\left(x_iterasi \right)\ge e \to x_i=g\left(x_{i-1} \right)\), Hitung \(f\left(x_i \right)\)
5. Akar persamaan adalah \(x\) terakhir yang diperoleh

FUngsi root_fpi() dapat digunakan untuk melakukan iterasi dengan argumen fungsi berupa persamaan non-linier, nilai tebakan awal, nilai toleransi, dan jumlah iterasi maksimum. Berikut adalah sintaks fungsi tersebut:

Contoh 7.4 Selesaikan persamaan non-linier pada Contoh 7.2 menggunakan metode iterasi titik tetap?

Jawab:

Untuk menyelesaikan persamaan non-linier tersebut kita perlu mentransformasi persamaan non-linier tersebut terlebih dahulu.

\[ xe^{-x}+1=0\ \to x=-\frac{1}{e^{-x}} \]

Untuk tebakan awal digunakan nilai \(x=-1\)

\[ x_1 = -\frac{1}{e^{1}}=-2,718282 \]

Nilai \(x\) tersebut selanjutnya dijadikan nilai input pada iterasi selanjutnya:

\[ x_2 = -\frac{1}{e^{2,718282}}=-0,06598802 \]

iterasi terus dilakukan sampai diperoleh \(\left| x_{i+1}-x_i \right|\le e\).

Untuk mempercepat proses iterasi kita dapat menggunakan bantuan fungsi root_fpi(). Berikut adalah sintaks yang digunakan:

root_fpi(function(x){-1/exp(-x)}, x0=-1)

Berdasarkan hasil iterasi diperoleh nilai \(x=-0,5671433\) dengan jumlah iterasi yang diperlukan sebanyak \(29\) kali. Jumlah iterasi akan bergantung dengan nilai tebakan awal yang kita berikan. Semakin dekat nilai tersebut dengan akar, semakin cepat nilai akar diperoleh.

Metode Newton-Raphson merupakan metode penyelesaian persamaan non-linier dengan menggunakan pendekatan satu titik awal dan mendekatinya dengan memperhatikan slope atau gradien. titik pendekatan dinyatakan pada Persamaan (7.7).

\[\begin{equation} x_{n+1}=x_n-\frac{f\left(x_n\right)}{f'\left(x_n\right)} \tag{7.7} \end{equation}\]

Ilustrasi metode Newton-Raphson disajikan pada Gambar 7.8.

Gambar 7.8: Ilustrasi metode Newton-Raphson.

Algoritma Metode Newton-Raphson

1. Definisikan \(f\left(x \right)\) dan \(f'\left(x \right)\)
2. Tentukan nilai toleransi \(e\) dan iterasi masimum (N)
3. Tentukan tebakan awal \(x_0\)
4. Hitung \(f\left(x_0 \right)\) dan \(f'\left(x_0 \right)\)
5. Untuk iterasi \(i=1\) s/d \(N\) atau \(\left|f\left(x \right) \right|\ge e\), hitung \(x\) menggunakan Persamaan (7.7)
6. Akar persamaan merupakan nilai \(x_i\) terakhir yang diperoleh.

Fungsi root_newton() merupakan fungsi yang dibuat menggunakan algoritma di atas. Fungsi tersebut dituliskan pada sintaks berikut:

Contoh 7.5 Selesaikan persamaan non-linier \(x-e^{-x}=0\) menggunakan metode Newton-Raphson?

Jawab:

Untuk dapat menggunakan metode Newton-Raphson, terlebih dahulu kita perlu memperoleh turunan pertama dari persamaan tersebut.

\[ f\left(x\right)=x-e^{-x}\to f'\left(x\right)=1+e^{-x} \]

Tebakan awal yang digunakan adalah \(x=0\).

\[ f\left(x_0\right)=0-e^{-0}=-1 \] \[ f'\left(x_0\right)=1+e^{-0}=2 \]

Hitung nilai \(x\) baru:

\[ x_1=x_0-\frac{f\left(x_0\right)}{f'\left(x_0\right)}=0-\frac{-1}{2}=0,5 \]

Untuk mempercepat proses iterasi, kita dapat menggunakan fungsi root_newton(). Berikut adalah sintaks yang digunakan:

root_newton(function(x){x-exp(-x)}, function(x){1+exp(-x)}, x0=0)

Berdasarkan hasil iterasi diperoleh akar penyelesaian persamaan non-linier adalah \(x=0,5671433\) dengan jumlah iterasi yang diperlukan adalah \(5\) iterasi.

Dalam penerapannya metode Newton-Raphson dapat mengalami kendala. Kendala yang dihadapi adalah sebagai berikut:

1. titik pendekatan tidak dapat digunakan jika merupakan titik ekstrim atau titik puncak. Hal ini disebabkan pada titik ini nilai \(f'\left(x \right)=0\). Untuk memahaminya perhatikan ilustasi yang disajikan pada Gambar 7.9. Untuk menatasi kendala ini biasanya titik pendekatan akan digeser.

Gambar 7.9: Ilustrasi titik pendekatan di titik puncak.

1. Sulit memperoleh penyelesaian ketika titik pendekatan berada diantara 2 titik stasioner. Untuk memahami kendala ini perhatikan Gambar 7.10. Untuk menghindarinya, penentuan titik pendekatan dapat menggunakan bantuan metode tabel.

Gambar 7.10: Ilustrasi titik pendekatan diantara 2 titik stasioner.

1. Turunan persamaan sering kali sulit untuk diperoleh (tidak dapat dikerjakan dengan metode analitik).

Metode Secant merupakan perbaikan dari metode regula-falsi dan Newton Raphson, dimana kemiringan dua titik dinyatakan secara diskrit dengan mengambil bentuk garis lurus yang melalui satu titik. Persamaan yang dihasilkan disajikan pada Persamaan (7.8).

\[\begin{equation} y-y_0=m\left(x-x_0\right) \tag{7.8} \end{equation}\]

Nilai \(m\) merupakan transformasi persamaan tersebut.

\[\begin{equation} m_n=\frac{f\left(x_n\right)-f\left(x_{n-1}\right)}{x_n-x_{n-1}} \tag{7.9} \end{equation}\]

Bila \(y=f\left(x\right)\) dan \(y_n\) dan \(x_n\) diketahui, maka titik ke \(n+1\) adalah:

\[\begin{equation} y_{n+1}-y_n=m_n\left(x_{n+1}-x_n\right) \tag{7.10} \end{equation}\]

Bila titik \(x_{n+1}\) dianggap akar persamaan maka nilai \(y_{n+1}=0\), sehingga diperoleh:

\[\begin{equation} -y_n=m_n\left(x_{n+1}-x_n\right) \tag{7.11} \end{equation}\]

\[\begin{equation} \frac{m_nx_n-y_n}{m_n}=x_{n+1} \tag{7.12} \end{equation}\]

atau

\[\begin{equation} x_{n+1}=x_n-y_n\frac{1}{m_n} \tag{7.13} \end{equation}\]

\[\begin{equation} x_{n+1}=x_n-f\left(x_n\right)\frac{x_n-x_{n+1}}{f\left(x_n\right)-f\left(x_{n+1}\right)} \tag{7.14} \end{equation}\]

Berdasarkan Persamaan (7.14) diketahui bahwa untuk memperoleh akar persamaan diperlukan 2 buah titik pendekatan. Dalam buku ini akan digunakan titik pendekatan kedua merupakan titik pendekatan pertama ditambah sepuluh kali nilai toleransi.

\[\begin{equation} x_1=x_0+10*tol \tag{7.15} \end{equation}\]

Algoritma Metode Secant

1. Definisikan \(f\left(x \right)\) dan \(f'\left(x \right)\)
2. Tentukan nilai toleransi \(e\) dan iterasi masimum (N)
3. Tentukan tebakan awal \(x_0\) dan \(x_1\)
4. Hitung \(f\left(x_0 \right)\) dan \(f\left(x_1 \right)\)
5. Untuk iterasi \(i=1\) s/d \(N\) atau \(\left|f\left(x \right) \right|\ge e\), hitung \(x\) menggunakan Persamaan (7.14)
6. Akar persamaan adalah nilai x yang terakhir.

Fungsi root_secant() merupakan fungsi yang penulis buat untuk melakukan iterasi menggunakan metode Secant. Berikut merupakan sintaks dari fungsi tersebut:

Contoh 7.6 Selesaikan persamaan non-linier pada Contoh 7.5 menggunakan metode Secant?

Jawab:

Untuk menyelesaikan persamaan tersebut digunakan nilai pendekatan awal \(x_0=0\) dan \(x_1=0+10*10^{-7}=10^{-6}\).

\[ f\left(x_0 \right)=0-e^{-0}=-1 \]

\[ f\left(x_1 \right)=10^{-6}-e^{-10^{-6}}=-0,999998 \]

Hitung nilai \(x_2\) dan \(f\left(x_2 \right)\).

\[ x_2=0+0,999998\frac{10^{-6}-0}{-0,999998+1}=0,499999 \]

Untuk mempercepat proses iterasi kita dapat menggunakan fungsi root_secant() pada R. Berikut sintaks yang digunakan:

root_secant(function(x){x-exp(-x)}, x=0)

Berdasarkan hasil iterasi diperoleh nilai akar penyelesaian adalah \(x=0,5671433\) dengan iterasi dilakukan sebanyak \(6\) kali.

Secara umum metode Secant menawarkan sejumlah keuntungan dibanding metode lainnya. Pertama, seperti metode Newton-Raphson dan tidak seperti metode tertutup lainnya, metode ini tidak memerlukan rentang pencarian akar penyelesaian. Kedua, tidak seperti metode Newton-Raphson, metode ini tidak memerlukan pencarian turunan pertama persamaan non-linier secara analitik, dimana tidak dapat dilakukan otomasi pada setiap kasus.

Adapun kerugian dari metode ini adalah berpotensi menghasilkan hasil yang tidak konvergen sama seperti metode terbuka lainnya. Selain itu, kecepatan konvergensinya lebih lambat dibanding metode Newton-Raphson.

Paket base pada R menyediakan fungsi uniroot() untuk mencari akar persamaan suatu fungsi pada rentang spesifik. Fungsi ini menggunakan metode Brent yaitu kombinasi antara root bracketing, biseksi, dan interpolasi invers kuadrat. Format fungsi tersebut secara sederhana adalah sebagai berikut:

uniroot(f, interval, tol=.Machine$double.eps^0.25, maxiter=1000)

Catatan:

• f: persamaan non-linier
• interval: vektor interval batas bawah dan atas
• tol: nilai toleransi
• maxiter: iterasi maksimum

Berikut adalah contoh penerapan fungsi uniroot():

uniroot(function(x){x*exp(-x)+1}, interval=c(-1,0), tol=1e-7)

Berdasarkan hasil iterasi diperoleh akar persamaan tersebut adalah \(-0,5671433\) dengan jumlah iterasi sebanyak \(7\) iterasi dan tingkat presisi sebesar \(5e-08\).

Fungsi lain yang dapat digunakan untuk mencari akar persamaan adalah uniroot.all() dari paket rootSolve. Fungsi ini mengatasi kelemahan dari uniroot(), dimana uniroot() tidak bekerja jika fungsi hanya menyentuh dan tidak melewati sumbu nol \(y=0\). Untuk memahaminya perhatikan contoh berikut:

uniroot(function(x){sin(x)+1}, c(-pi,0))

Bandingkan dengan sintaks berikut:

uniroot(function(x){sin(x)+1}, c(-pi,-pi/2))

Untuk menggunakan fungsi uniroot.all(), jalankan sintaks berikut:

Jalankan kembali fungsi dan rentang di mana uniroot() tidak dapat bekerja:

uniroot.all(function(x){sin(x)+1}, c(-pi,0))

Fungsi polyroot() pada paket base dapat digunakan untuk memperoleh akar dari suatu polinomial. Algortima yang digunakan dalam fungsi tersebut adalah algoritma Jenkins dan Traub.

Untuk dapat menggunakannya kita hanya perlu memasukkan vektor koefisien dari polinomial. Pengisian elemen dalam vektor dimulai dari variabel dengan pangkat tertinggi menuju variabel dengan pangkat terendah. Berikut adalah contoh bagaimana fungsi polyroot() digunakan untuk mencari akar polinomial \(f\left(x\right)=x^2+1\):

Contoh lainnya adalah mencari akar polinomial \(f\left(x\right)=4x^2+5x+6\):

Pembaca dapat mencoba membuktikan hasil yang diperoleh tersebut menggunakan metode analitik.

Penerapan penyelesaian sistem persamaan non-linier banyak dijumpai dalam berbagai kasus di bidang lingkungan. Pada bagian ini penulis tidak akan menjelaskan seluruhnya. Penulis hanya akan menjelaskan penerapannya pada sebuah persamaan yaitu Hukum Bernoulli.

Misalkan terdapat sebuah saluran dengan penampang sesuai dengan Gambar 7.11.

Gambar 7.11: Aliran fluida pada sebuah pipa.

Berdasarkan hukum Bernoulli, maka diperoleh persamaan berikut:

\[\begin{equation} \frac{Q^2}{2gb^2h_0^2}+h_0=\frac{Q^2}{2gb^2h^2}+h+H \tag{7.16} \end{equation}\]

Persamaan tersebut dapat dilakukan transformasi menjadi persamaan berikut:

\[\begin{equation} f\left(h\right)=h^3+\left(H-\frac{Q^2}{2gb^2h_0^2}-h_0\right)h^2+\frac{Q^2}{2gb^2}=0 \tag{7.17} \end{equation}\]

Data-data terkait saluran tersebut adalah sebagai berikut:

• \(Q=1,2\ \frac{m^3}{\det\ }\) = volume aliran fluida tiap satuan waktu
• \(g=9,81\ \frac{m}{s^2}\) =percepatan gravitasi
• \(b=1,8\ m\) =lebar pipa
• \(h_0=0,6\ m\) =ketinggian air maksimum
• \(H=0,075\ m\) =tinggi pelebaran pipa
• \(h\) = ketinggian air

Kita dapat menggunakan pendekatan numerik untuk menentukan \(h\). Pada studi kasus ini tidak dijelaskan lokasi dimana akar penyelesaian berada, sehingga metode terbuka seperti Secant cukup sesuai untuk menyelesaikannya:

Berikut adalah persamaan yang baru setelah seluruh data dimasukkan kedalam tiap variabelnya:

\[ f\left(h\right)=h^3+\left(0,075-\frac{1,2^2}{2\times9,81\times1,8^2\times0,6^2}-0,6\right)h^2+\frac{1,2^2}{2\times9,81\times1,8^2}=0 \]

Untuk penyelesaiannya penulis akan memberikan tebakan awal nilai \(h=h_0=0,6\). Berikut adalah sintaks penyelesaian menggunakan metode secant:

Berdasarkan hasil perhitungan diperoleh nilai \(h=-0,2870309\) atau ketinggian air sekitar \(0,3\ m\) dengan jumlah iterasi sebanyak \(26\) kali.

Pembaca dapat mencoba menggunakan metode lain seperti metode tertutup. Untuk dapat melakukannya, pembaca perlu memperoleh rentang lokasi akar persamaan tersebut berada menggunakan metode tabel.

1. Atmika, I.K.A. 2016. Diktat Mata Kuliah: Metode Numerik. Jurusan Teknik Mesin Universitas Udayana.
2. Bloomfield, V.A. 2014. Using R for Numerical Analysis in Science and Engineering. CRC Press
3. Howard, J.P. 2017. Computational Methods for Numerical Analysis with R. CRC Press.
4. Jones, O. Maillardet, R. Robinson, A. 2014. Introduction to Scientific Programming and Simulation Using R. CRC Press
5. Kreyszig, E. 2011. Advanced Engineering Mathematics, 10th Edition. John Wiley & Sons.
6. Sanjaya, M. 2015. Metode Numerik Berbasis Phython. Penerbit Gava Media: Yogyakarta.
7. Sudiadi dan Teguh R. 2015. Metode Numerik. STMIK

1. Temukan akar persamaan dari persamaan non-linier \(f\left(x\right)=x^3-2x+2\) menggunakan metode terbuka dengan \(x_0=0\) dan \(x_0=\frac{1}{2}\)!
2. Apakah kelebihan dari metode tertutup (contoh: metode biseksi) dibanding metode terbuka (contoh: Newton-Raphson)? (catatan: pembaca dapat pula mencari dari referensi lainnya)
3. Temukan akar persamaan dari persamaan \(f\left(x\right)=\frac{sin\left(x\right)}{x}\) dengan rentang pencarian \(x=0,5\) dan \(x=1\)!
4. Pada kondisi apakah metode Secant lebih dipilih dibanding metode Newton-Raphson?
5. Modifikasilah fungsi root_bisection() dan root_rf() sehingga kita tidak perlu memasukkan argumen a dan b dan hanya perlu memasukkan satu vektor interval kedalam fungsi tersebut! (contoh: interval=c(a,b))

Pada dunia nyata, data sering kali tidak tersaji secara lengkap. Seringkali terdapat nilai data yang hilang (missing value). Terdapat banyak penyebab dari kondisi tersebut, baik akibat kesalahan manusianya maupun keterbatasan kemampuan alat ukur.

Kondisi lain yang muncul dari data yang kita miliki adalah adanya outlier atau nilai yang berbeda jauh dengan mayoritas data yang kita miliki. Nilai tersebut akan menentukan hasil analisis atau uji statistik yang kita lakukan, terlebih lagi jika uji statistik yang kita lakukan menggunakan metode parametrik.

Terdapat banyak cara untuk menangani kondisi-kondisi tersebut. Sejumlah peneliti memilih untuk menghapus data tersebut. Hal ini dapat dilakukan jika jumlah data yang kita miliki cukup besar. Bagaimana jika data yang kita miliki sedikit dan pengukuran ulang cukup mahal atau cukup sulit dilakukan?. Salah satu cara yang dapat dilakukan adalah dengan melakukan interpolasi terhadap data.

Interpolasi dan ekstrapolasi adalah proses “menebak” nilai data dengan memperhatikan data lain yang kita miliki. Interpolasi merupakan teknik untuk mencari nilai suatu variabel yang hilang pada rentang data yang diketahui, sedangkan ektrapolasi merupakan teknik menemukan nilai suatu variabel diluar rentang data yang telah diketahui. Data lain yang kita miliki seringkali memiliki sejumlah pola. Pola yang terbentuk dapat berupa polinomial atau mengelompok. Tiap pola akan memiliki metode pendekatan yang berbeda-beda. Terdapat kemungkinan tak terbatas dari pola data tersebut. Penilaian profesional atau ahli diperlukan untuk menentukan metode mana yang sesuai berdasarkan riwayat penelitian atau pekerjaan yang pernah dilakukan sebelumnya.

Pada Chapter 8 penulis akan menjelaskan teknik-teknik interpolasi yang dapat kita lakukan. Adapun yang akan dibahas pada Chapter ini adalah sebagai berikut:

• Teknik interpolasi polinomial
• Teknik interpolasi piecewise
• Studi kasus penerapan teknik interpolasi

Interpolasi polinomial merupakan teknik interpolasi dengan mengasumsikan pola data yang kita miliki mengikuti pola polinomial baik berderajat satu (linier) maupun berderajat tinggi. Interpolasi dengan metode ini dilakukan dengan terlebih dahulu membentuk persamaan polinomial. Persamaan polinomial yang terbentuk selanjutnya digunakan untuk melakukan interpolasi dari nilai yang diketahui atau ekstrapolasi (prediksi) dari nilai diluar rentang data yang diketahui.

Pada Chapter 8.1 pembahasan akan dibagi menjadi 3 bagian. Bagian pertama kita akan mengulang kembali teknik evaluasi polinomial, sedangkan dua bagian selanjutnya akan membahas teknik interpolasi linier dan polinomial orde tinggi dengan menjadikan pembahasan bagian pertama sebagai dasar pada dua bagian berikutnya.

Pada Chapter ini pembaca akan mempelajari teknik untuk melakukan substitusi nilai \(x\) pada persamaan polinomial untuk memperoleh nilai \(y\). Terdapat berbagai pendekatan dalam melakukan proses tersebut, mulai dari metode naive maupun metode Horner. Kedua metode akan menghasilkan hasil yang sama namun dengan proses komputasi yang berbeda. Metode naive cenderung lambat dalam proses komputasi karena jumlah proses yang dilakukan dalam sekali proses lebih banyak dari pada metode Horner.

Untuk memahami metode-metode evaluasi polinomial yang telah disebutkan tersebut, secara umum persamaan polinomial disajikan pada Persamaan (8.1).

\[\begin{equation} f\left(x\right)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0 \tag{8.1} \end{equation}\]

dimana \(a\) merupakan koefisien polinomial, \(x\) merupakan variabel, dan \(n\) merupakan indeks dan pangkat polinomial.

Pada metode naive kita melakukan evaluasi polinomial sama dengan cara kita melakukan evaluasi polinomial saat kita SMA. Nilai \(x\) akan disubstitusikan pada masing-masing elemen persamaan polinomial. Masing-masing elemen polinomial selanjutnya dijumkahkan untuk menghitung \(y\).

Pada R kita dapat menuliskan sebuah fungsi untuk melakukan evaluasi polinomial menggunakan metode naive tersebut. Pada fungsi tersebut, koefisien polinomial akan disimpan kedalam sebuah vektor dengan urutan pengisian mulai dari koefisien dengan pangkat \(x\) terendah ke tertinggi.

Contoh 8.1 Hitung nilai \(y\) pada persamaan \(f\left(x\right)=x^4+3x^3-15x^2-19x+30\), jika diketahui nilai \(x\) adalah -1, 0, dan 1!

Jawab:

Untuk dapat menghitung nilai \(y\) menggunakan fungsi naive_poly() pada persamaan tersebut dengan nilai \(x\) yang diketahui, kita perlu merubah koefisien persamaan tersebut dan nilai \(x\) yang diketahui menjadi vektor:

x <- c(-1,0,1) coeff <- c(30,-19,-15,3,1)

Masukkan vektor-vektor yang telah terbentuk tersebut kedalam fungsi naive_poly().

Berdasarkan hasil perhitungan diperoleh nilai \(y\) masing-masing sebesar 32, 30, dan 0. Pembaca dapat mengeceknya sendiri hasil perhitungan tersebut menggunakan cara manual.

Kita dapat meningkatkan efisiensi proses perhitungan pada fungsi naive_poly() tersebut. Sebagai contoh, setiap kali kita melakukan loop untuk menghitung nilai \(y\) pada polinomial, kita dapat memperoleh eksponensial dari \(x\). Namun untuk setiap koefisien \(a_i\), nilai eksponensial yang terkait \(x_i\) merupakan seri produk.

\[\begin{equation} x^i=\prod_{j=0}^ix \tag{8.2} \end{equation}\]

Untuk \(x^i\), terdapat sebanyak \(i\) koefisien \(x\) dikalikan bersamaan. Namun, untuk setiap \(x^{i-1}\) terdapat lebih sedikit 1 perkalian dibanding koefisien \(x\) dengan pangkat yang lebih besar dan seterusnya.

Berdasarkan ilustrasi tersebut, kita dapat membentuk fungsi better_poly() sebagai perbaikan dari fungsi naive_poly(). Berikut adalah sintaks yang digunakan:

Contoh 8.2 Hitung nilai \(y\) pada persamaan yang disajikan pada Contoh 8.1 menggunakan nilai \(x\) yang telah diketahui pada soal tersebut menggunakan fungsi better_poly()?

Jawab:

Sejauh ini kita telah membentuk 2 fungsi yaitu naive_poly() dan better_poly(). Kedua fungsi tersebut memiliki perbedaan proses menghitung yang mempengaruhi efisiensinya masing-masing. Sebagai contoh jika diberikan polinomial berderajat 10, fungsi naive_poly() akan mengejakan perkalian dalam proses loop sebanyak 55 kali, sedangkan fungsi better_poly() akan melakukannya sebanyak 20 kali (\(2n\) perkalian).

Metode lain yang lebih efisien dalam melakukan evaluasi polinomial adalah metode Horner. Metode ini oleh William Horner pada abad ke-18. Dalam metode Horner, bentuk polinomial pada Persamaan (8.1) akan ditransformasi menjadi Persamaan (8.3).

\[\begin{equation} \begin{split} f\left(x\right)&=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0\\ &= a_0+a_1x+\cdots+a_{n-1}x^{n-1}+a_nx^n\\ &= a_0+x\left(a_1+\cdots+a_{n-1}x^{n-2}+a_nx^{n-1}\right)\\ &= a_0+x\left(a_1+\cdots+x\left(a_{n-1}+x\left(a_n\right)\right)\cdots\right) \end{split} \tag{8.3} \end{equation}\]

Berdasarkan Persamaan (8.3), jika kita melakukan perhitungan pada persamaan polinomial berderajat 10, kita dapat mereduksi perhitungan menjadi 10 perkalian dan 10 penjumlahan. Jumlah tersebut sangat kecil dibandingkan kedua metode sebelumnya dan dapat dikatakan lebih efisien dibandingkan metode lainnya.

FUngsi horner_poly() merupakan fungsi yang dibentuk bedasarkan Persamaan (8.3). Sintaks fungsi tersebut adalah sebagai berikut:

Contoh 8.3 Kerjakan kembali Contoh 8.2 fungsi horner_poly()?

Jawab:

Misalkan kita memiliki 3 buah data dengan dua buah variabel kita misalkan variabel \(x\) dan variabel \(y\). Pada salah satu data terdapat data yang hilang pada. Agar ketiga data tersebut tetap dapat digunakan dalam iterasi diperlukan interpolasi untuk “menebak” nilai dari data yang hilang.

Berdasarkan pengukuran yang sebelumnya pernah dilakukan diketahui bahwa pola data yang terbentuk variabel \(x\) dan \(y\) divisualisasikan menggunakan scatterplot adalah pola linier. Berdasarkan hal tersebut interpolasi dilakukan dengan menggunakan metode linier.

Interpolasi linier dilakukan dengan terlebih dahulu membentuk fungsi linier. Dengan kata lain kita perlu mencari nilai slope \(m\) dan intercept \(b\). Nilai \(m\) dihitung sebagai rasio selisih jarak dua titik pada sumbu \(y\) dan sumbu \(x\) yang dapat dituliskan melalui Persamaan (8.4).

\[\begin{equation} m=\frac{y_2-y_1}{x_2-x_1} \tag{8.4} \end{equation}\]

Nilai intercept (titik potong pada sumbu \(y\)) dihitung menggunakan Persamaan (8.5).

\[\begin{equation} b=y_2-mx_2 \tag{8.5} \end{equation}\]

Algoritma Interpolasi Linier

1. Tentukan dua buah titik \(\left(x,y\right)\) sebagai dasar pembentukan persaman linier.
2. Hitung \(m\) menggunakan Persamaan (8.4)
3. Hitung \(b\) menggunakan Persamaan (8.5)
4. Definiskan fungsi linier berdasarkan nilai \(m\) dan \(b\)
5. Hitung \(y\) dengan cara substitusi nilai \(x\) pada persamaan linier untuk melakukan interpolasi atau ekstrapolasi nilai \(y\) yang ingin dicari.

Algoritma poin 1 sampai 4 tersebut, kita dapat membentuk fungsi pembentuk persamaan linier dari 2 titik yang diketahui. Fungsi tersebut disajikan pada sintaks berikut:

Contoh 8.4 Diketahui koordinat 2 buah titik yaitu (0,-1) dan (2,3). Jika diketahui titik ketiga memiliki koordinat sumbu \(x\) sebesar 1. Lakukan interpolasi untuk menentukan koordinat sumbu \(y\) titik ketiga tersebut!

Jawab:

Berdasarkan data-data yang terdapat pada soal terserbut, kita dapat menghitung nilai \(m\) dan \(b\). Nilai \(m\) dapat dihitung sebagai berikut:

\[ m=\frac{-1-3}{0-2}=2 \]

Dengan menggunakan nilai \(m\) tersebut kita dapat menghitung nilai \(b\).

\[ b=3-2*3=-1 \]

Berdasarkan hasil perhitungan diperoleh persamaan linier yang terbentuk adalah sebagai berikut:

\[ y=2x-1 \]

Berdasarkan persamaan linier tersebut nilai \(y\) dapat dihitung.

\[ y=2*1-1=1 \]

Kita dapat pula membentuk persamaan linier menggunakan fungsi linear_inter(). Berikut adalah sintaks yang digunakan:

x <- c(0,2) y <- c(-1,3) (coeff <- linear_inter(x,y))

Setelah diperoleh koefisien persamaan linier berdasarkan dua titik tersebut, kita akan menggunakan fungsi horner_poly() untuk memperoleh nilai \(y\). Berikut adalah sintaks yang digunakan:

Hasil interpolasi yang diperoleh dapat dikatakan sesuai dengan lokasi kedua titik data yang ditunjukkan pada Gambar 8.1. Berdasarkan hal tersebut, kita dapat yakin bahwa hasil interpolasi yang telah kita lakukan telah sesuai.

Gambar 8.1: Interpolasi linier dua titik (Sumber:Howard, 2017).

Metode interpolasi linier dapat dibilang merupakan metode interpolasi yang sangat sederhana. Disamping kemudahannya, metode ini memiliki potensi error numerik jika jarak antara kedua titik cukup berdekatan terlebih lagi jika selisih penyebut (\(x_2-x_1\)) sangat kecil sehingga akan menghasilkan nilai \(y\) yang sangat besar.

Disamping adanya potensi error numerik tersebut, metode ini menjadi dasar bagi metode interpolasi lain yang lebih kompleks. Metode selanjutnya merupakan pengembangan dari metode interpolasi ini.

Dengan menggunakan dua titik, kita dapat membentuk garis lurus (linier) yang tepat pada dua titik tersebut. Masalah timbul jika selisih nilai \(x\) kedua titik tersebut sangat kecil atau kedua titik tersebut memiliki nilai \(x\) yang sama. Hal ini akan menyebabkan slope yang dihasilkan menjadi tidak terhingga atau garis yang terbentuk adalah garis vertikal tegak lurus.

Bagaimana jika terdapat tiga buah titik? apakah kita masih bisa menggunakan interpolasi linie?. Ya, asalkan ketiga titik tersebut membentuk pola linier atau terletak pada satu garis yang sama. Pada kenyatannya kondisi tersebut jarang terjadi, sehingga pendekatan menggunakan polinomial orde lebih tinggi diperlukan. Persamaan kuadratik (polinomial orde dua) dapat digunakan untuk membentuk persamaan polinomial pada ketiga titik tersebut, sehingga iterasi dapat dilakukan. Untuk 4 buah titik data, polinomial orde tiga dapat digunakan untuk melakukan interpolasi. Secara umum berdasarkan penjelasan tersebut, untuk \(n\) titik data interpolasi dapat dilakukan menggunakan persamaan polinomial orde \(n-1\).

Diberikan set data berpasangan yang telah diurutkan \(\left(x_i,y_i\right)\), fungsi interpolasi harus memenuhi persyaratan berikut:

\[\begin{equation} p\left(x_i\right)=y_i \tag{8.6} \end{equation}\]

Untuk setiap \(i\). Sebagai tambahan, fungsi interpolasi berupa fungsi polinomial dengan bentuk umum sebagai berikut:

\[\begin{equation} y_i=\beta_nx_i^n+\beta_{n-1}x_i^{n-1}+\cdots+\beta_1x_i+\beta_0 \tag{8.7} \end{equation}\]

Persamaan (8.7) dapat dituliskan kedalam bentuk matriks yang ditampilkan pada Persamaan (8.8).

\[\begin{equation} \begin{bmatrix} x_1^n & x_1^{n-1} & \cdots & x_1 & 1 \\[0.3em] x_2^n & x_2^{n-1} & \cdots & x_2 & 1 \\[0.3em] \vdots & \vdots & \vdots &\ddots& \vdots \\[0.3em] x_n^n & x_n^{n-1} & \cdots & x_n & 1 \end{bmatrix} \begin{bmatrix} \beta_n \\[0.3em] \beta_{n-1} \\[0.3em] \vdots \\[0.3em] \beta_0 \end{bmatrix} = \begin{bmatrix} y_1 \\[0.3em] y_2 \\[0.3em] \vdots \\[0.3em] y_n \end{bmatrix} \tag{8.8} \end{equation}\]

Persamaan matrik tersebut dapat dituliskan sebagai \(X\beta = y\). Untuk menyelesaikan persamaan tersebut (memperoleh nilai \(\beta\)), pembaca dapat membaca kembali Chapter 6. Matriks \(X\) disebut sebagai matriks Vandemonde dan matriks tersebut mengandung sejumlah nilai \(x\) dengan pangkat sampai dengan \(n\).

Algoritma Interpolasi Polinomial Orde Tinggi

1. Tentukan set titik berpasangan \(\left(x,y\right)\) yang telah diurutkan.
2. Bentuk matriks Vandermonde sesuai dengan Persamaan (8.8).
3. Definiskan persamaan matriks \(X\beta=y\)
4. Selesaikan persamaan matriks pada poin 3 untuk memperoleh nilai \(\beta\)
5. Definisikan persamaan polinomial berdasarkan koefisien \(\beta\) yang diperoleh
6. Lakukan substitusi \(x\) persamaan polinomial pada poin 5 untuk memperoleh nilai \(y\)

Berdasarkan algoritma poin 1 sampai 5, kita dapat membentuk suatu fungsi untuk membentuk persamaan polinomial berdasarkan Persamaan (8.8). Berikut adalah sintaks fungsi tersebut:

Contoh 8.5 Diketahui koordinat 3 buah titik yaitu (-1,-2), (1,2) dan (0,1). Jika diketahui titik keempat memiliki koordinat sumbu \(x\) sebesar -2. Lakukan interpolasi untuk menentukan koordinat sumbu \(y\) titik keempat tersebut!

Jawab:

Untuk menyelesaiakn contoh soal tersebut, kita perlu terlebih dahulu membentuk matriks sesuai dengan Persamaan (8.8). Berdasarkan soal tersebut, terdapat tiga buah titik data yang diketahui, sehingga polinomial yang hendak dibentuk selanjutnya adalah polinomial berderajat 2.

\[\begin{equation*} \begin{bmatrix} -1^2 & -1^{1} & 1 \\[0.3em] 1^2 & 1^{1} & 1 \\[0.3em] 0^2 & 0^{1} & 1 \end{bmatrix} \begin{bmatrix} \beta_2 \\[0.3em] \beta_1 \\[0.3em] \beta_0 \end{bmatrix} = \begin{bmatrix} -2 \\[0.3em] 2 \\[0.3em] 1 \end{bmatrix} \end{equation*}\]

Setelah matriks tersebut terbentuk, pembaca dapat menyelesaikannya menggunakan berbagai metode yang telah penulis jelaskan pada Chapter 6 untuk memperoleh nilai \(\beta\).

Untuk menyelesaikan contoh soal tersebut pada R, kiat perlu membentuk matriks \(x\) dan \(y\) terlebih dahulu.

x <- c(-1, 1, 0) y <- c(-2, 2, -1)

Koefisien persamaan polinomial dihitung menggunakan fungsi poly_inter().

(coeff <- poly_inter(x, y))

Berdasarkan hasil perhitungan, diperoleh nilai \(\beta\). Nilai tersebut selanjutnya digunakan untuk membentuk persamaan polinomial. Berikut merupakan persamaan polinomial yang terbentuk:

\[ f\left(x\right)=1x^2+2x-1 \]

Fungsi horner_poly() selanjutnya digunakan untuk mengevaluasi polinomial tersebut. Berikut adalah hasil substitusi \(x\) pada persamaan tersebut:

Hasil yang diperoleh terlihat cukup sesuai jika kita perhatikan visualisasi ketiga titik tersebut pada Gambar 8.2.

Gambar 8.2: Interpolasi kuadratik tiga titik (Sumber:Howard, 2017).

Contoh 8.6 Dengan menggunakan data pada Contoh 8.4, lakukan proses perhitungan untuk membentuk persamaan polinomial menggunakan fungsi poly_inter()

Jawab:

Berdasarkan data pada Contoh 8.4, polinomial yang terbentuk merupakan polinomial derajat 1 (linier). Berikut adalah nilai koefisien yang dihasilkan dari perhitungan menggunakan fungsi poly_inter().

x <- c(2, 0) y <- c(3, -1) poly_inter(x, y)

Meskipun proses perhitungan menggunakan fungsi poly_inter() lebih rumit dibandingkan dengan fungsi linear_poly(), hasil perhitungan keduanya menggunakan data pada Contoh 8.4 menghasilkan hasil yang sama.

Interpolasi dengan polinomial sering memberikan hasil yang tidak dapat diterima. Interpolasi polinomial yang dihasilkan dari sejumlah besar data titik biasanya berderajat tinggi. Polinomial berderajat tinggi pada umumnya bersifat osilatif (grafiknya naik turun secara cepat). Akibatnya, perubahan data pada interval kecil dapat menyebabkan fluktuasi besar pada keseluruhan interval. Karena alasan ini, biasanya interpolasi hanya menggunakan polinomial berderajat rendah.

Interpolasi piecewise menawarkan alternatif lain. Pada interpolasi piecewise, pada titik yang berbeda sepanjang kurva, nilai fungsi lebih mungkin lebih baik didekati menggunakan dua atau lebih interpolasi. Pada metode ini kita akan membuat fungsi interpolasi ditiap antara dua titik observasi.

Pada sub-Chapter ini akan dijelaskan 2 buah metode interpolasi piecewise, yaitu: interpolasi linier piecewise dan interpolasi kubik spline. Interpolasi pertama dilakukan menggunakan persamaan linier, sehingga kurva yang terbentuk bukan merupakan kurva kontinu. Interpolasi selanjutnya dilakukan menggunakan persamaan polinomial berderajat tinggi sehingga kurva yang dihasilkan lebih halus (tidak ada sudut siku pada setiap titik).

Interpolasi linier piecewise merupakan interpolasi yang menggunakan pendekatan interpolasi linier. Fungsi linier akan dibentuk pada setiap dua titik observasi. Untuk lebih memahaminya perhatikan kembali Gambar 8.2. Pada gambar tersebut sebelumnya kita telah membentuk persamaan kudratik untuk menghubungkan titik-titik tersebut. Dibanding menggunakan persamaan polinomial seperti kuadratik tersebut, interpolasi piecewise akan menghubungkan tiap dua titik observasi tersebut dengan garis lurus.

Algoritma Interpolasi Linier Piecewise

1. Tentukan set titik berpasangan \(\left(x,y\right)\) yang telah diurutkan berdasarkan nilai sumbu \(x\).
2. Hitung \(m\) pada setiap dua titik berdekatan menggunakan Persamaan (8.4)
3. Hitung \(b\) pada setiap dua titik berdekatan menggunakan Persamaan (8.5)
4. Definiskan fungsi linier berdasarkan nilai \(m\) dan \(b\)
5. Hitung \(y\) dengan cara substitusi nilai \(x\) pada persamaan linier untuk melakukan interpolasi nilai \(y\) yang ingin dicari.
6. Untuk melakukan ekstrapolasi dengan titik observasi diluar rentang titik diketahui, gunakan persamaan linier yang berada pada bagian ujung terdekat dengan nilai \(x\) yang hendak dicari nilai \(y\)-nya.

Berdasarkan algoritma tersebut, kita dapat menyusun fungsi pada R untuk membentuk persamaan linier piecewise. Berikut adalah sintaks yang digunakan:

Contoh 8.7 Tentukan persamaan-persamaan linier yang dihasilkan dari titik observasi yang ditampilkan pada Contoh 8.5?

Jawab:

Untuk menentukan persamaan-persamaan linier yang menghubungkan setiap titik, kita akan menggunakan fungsi pwise_linterp() yang telah kita buat. Berikut adalah sintaks yang digunakan:

x <- c(-1, 1, 0, -2) y <- c(-2, 2, -1, -1) pwise_linterp(x, y)

Jika persamaan-persamaan yang terbentuk tersebut divisualisasikan akan terlihat seperti pada Gambar 8.3.

Gambar 8.3: Interpolasi linier piecewise empat titik (Sumber:Howard, 2017).

R juga menyediakan fungsi untuk melakukan interpolasi linier piecewise. Fungsi approxfun() dapat digunakan untuk melakukan interpolasi tersebut. Fungsi approxfun() hanya memerlukan dua input yaitu vektor \(x\) dan vektor \(x\). Untuk lebih memahaminya, kita akan menggunakan kembali data yang disajikan pada Contoh 8.5 untuk membuat fungsi linier piecewise. Berikut adalah sintaks yang digunakan:

f <- approxfun(x, y) # tentukan nilai y jika x= 0 f(0)

# tentukan nilai y jika x = 0.5 f(0.5)

Jika menggunakan interpolasi polinomial berderajat satu (sebuah garis) lebih dari beberapa interval merupakan peningkatan dari satu baris interpolasi, dan jika menggunakan polinomial berderajat tinggi juga merupakan peningkatan dari satu garis interpolasi tunggal, maka dapat disimpulkan bahwa penggunaan polinomial berderajat tinggi pada selang beberapa interval juga akan menjadi peningkatan dalam proses interpolasi. Dalam beberapa kasus, hal tersebut benar, tetapi kita masih menghadapi sudut tajam di mana masing-masing kurva interpolasi tergabung (interpolasi linier piecewise). Sudut tajam ini mencegah diferensiasi dan pada prakteknya tidak dapat digunakan untuk memodelkan beberapa fungsi di dunia nyata, seperti roller coaster span.

Interpolasi spline kubik memecahkan masalah ini. Interpolasi ini akan memberikan kurva tergabung yang halus. Hal tersebut juga membuat spline terintegrasi. Karena setiap bagian individu diwakili oleh kurva kubik (polinomial derajat 3), maka masing-masing bagian individu juga dapat dianalisis sebagai kurva kubik. Dengan asumsi ada \(n\) titik data untuk interpolasi, kita akan mendefinisikan \(S_i\) sebagai fungsi polinomial kubik yang mewakili kurva pada domain \(\left[x_i; x_{i + 1}\right]\). Kemudian untuk \(n\) titik observasi, ada \(n - 1\) interpolasi polinomial kubik.

Bentuk umum seri polinomial dituliskan pada Persamaan (8.9).

\[\begin{equation} S_i=d_i\left(x-x_i\right)^3+c_i\left(x-x_i\right)^2+b_i\left(x-x_i\right)+a_i \tag{8.9} \end{equation}\]

Ini mengarah ke 4 nilai yang tidak diketahui tidak diketahui, yaitu: \(a_i\), \(b_i\), \(c_i\), dan \(d_i\) pada setiap Persamaan (8.9). Oleh karena itu, ada \(4 \left(n - 1\right) = 4n - 4\) nilai yang tidak diketahui. Karena kita ingin spline membentuk garis kontinu dan dapat didiferensiasi, ada satu set persamaan yang menentukan spline kubik:

\[\begin{equation} S_i\left(x_i\right)=y_i,\ \ \ \ \ \ \ \ \ \ i=1,\dots\dots,n-1 \tag{8.10} \end{equation}\]

\[\begin{equation} S_i\left(x_{i+1}\right)=y_{i+1},\ \ \ \ \ \ \ \ \ \ i=1,\dots\dots,n-1 \tag{8.11} \end{equation}\]

\[\begin{equation} S_i'\left(x_{i+1}\right)=S'_{i+1}\left(x_i\right),\ \ \ \ \ \ \ \ \ \ i=1,\dots\dots,n-2 \tag{8.12} \end{equation}\]

\[\begin{equation} S_i''\left(x_{i+1}\right)=S''_{i+1}\left(x_i\right),\ \ \ \ \ \ \ \ \ \ i=1,\dots\dots,n-2 \tag{8.13} \end{equation}\]

Persamaan (8.10) dan (8.11) sudah cukup jelas. Persyaratan ini memastikan bahwa jika kita mengevaluasi spline di salah satu node internal, hasil yang akan kita peroleh merupakan jawaban yang telah ditentukan, yaitu spline yang dievaluasi pada \(x_i\) untuk beberapa \(i\) adalah \(y_i\), dan setiap komponen spline bergabung dengan rapi. Persamaan (8.12) memastikan bahwa kita memiliki turunan pertama yang berkelanjutan di setiap simpul internal. Ini mencegah terbentuknya sudut tajam pada tiap node. Persamaan (8.13) memastikan turunan kedua juga kontinu, dimana kondisi ini menguntungkan karena itu berarti turunan pertama itu sendiri dapat didiferensiasi juga.

Kondisi ini menyebabkan ada \(4n − 6\) kondisi yang harus kita penuhi. Jika \(4n − 4\) kita yang tidak diketahui dipecahkan sebagai sebuah matriks, dan akhirnya matriks tersebut akan terpecahkan, matriks akan menjadi kurang ditentukan. Kita dapat menyelesaikan kondisi tersebut dengan memasukkan dua ketentuan tambahan. Dengan splines kubik, secara normal adalah menentukan akhir di kedua ujung untuk mencapai dua kondisi tambahan. Untuk contoh ini, dua kondisi yang akan kita tambahkan adalah \(S_i''\left(x_{i}\right)= 0\) dan \(S_i''\left(x_{n}\right)= 0\). Kedua kondisi ini memastikan bahwa pada titik akhir, turunan pertamanya linier dan oleh karena itu fungsi spline berlanjut ke arah yang sudah berjalan. Interpolasi spline kubik ini disebut juga sebagai “natural spline”.

Awalnya, kita dapat melihat bahwa setiap polinomial kubik \(S_i\) digeser ke kanan oleh unit \(x_i\) atau bergeser ke kiri jika \(x_i\) negatif. Pada \(x_i\) nilai fungsi adalah \(a_i\) yang berarti \(a_i = y_i\) untuk setiap nilai \(i\). Menyelesaikan sisa koefisien yang ada akan lebih kompleks, tetapi sekarang \(3n\) tidak diketahui dengan \(3n\) kondisi. Derivasi penuh tersedia dari berbagai sumber, tetapi secara garis besar dari Persamaan (8.12) kita mendapatkan \(S_i'\left(x_{i+1}\right)=S'_{i+1}\), kemudian \(S'_{i + 1} - S'_i \left(x\right) = 0\), dan kita dapat mensubstitusi Persamaan (8.9) ke dalam kedua komponen, pemecahan untuk \(d_i\) dalam hal ini \(x_i\), \(y_i\), dan \(c_i\). Proses yang sama dapat direplikasi dengan Persamaan (8.13) dan \(b_i\). Hasilnya adalah matriks tridiagonal, seperti yang ada pada Chapter 6.3.3, yang dapat dipecahkan untuk menemukan koefisien. Terdapat sebuah matriks, A, sedemikian rupa sehingga,

\[\begin{equation} AC=V \tag{8.14} \end{equation}\]

dimana \(A\) merupakan matriks tridiagonal. Pada matriks tridiagonal ini \(\left(u_1,u_2,\dots,u_{n-1}\right)\) merupakan vektor \(U\), dan vektor \(M\), \(L\), dan \(V\) ditentukan dengan cara serupa. Matriks ini sedemikian rupa,

\[\begin{equation} u_i=l_i=x_{i+1}-x_i \tag{8.15} \end{equation}\]

kecuali pada \(u_0=l_n\). Lebih jauh, diagonal utama \(D\) ditentukan menggunakan Persamaan (8.16).

\[\begin{equation} m_i=2\left(x_{i+1}-x_i+x_i-x_{i-1}\right) \tag{8.16} \end{equation}\]

kecuali pada \(d_0=d_n=1\). Akhirnya vektor \(V\) ditentukan dengan Persamaan (8.17).

\[\begin{equation} v_i=3\left(\frac{y_{i+1}-y_i}{x_{i+1}-x_i}-\frac{y_i-y_{i-1}}{x_i-x_{i-1}}\right) \tag{8.17} \end{equation}\]

kecuali pada \(v_0=v_n=0\). Sehingga,

\[\begin{equation} \begin{bmatrix} m_{1} & u_{1} & 0 & 0 & 0 \\[0.3em] l_{1} & m_{2} & u_{2} & 0 & 0 \\[0.3em] 0 & l_{2} & \ddots & \ddots & 0 \\[0.3em] 0 & 0 & \ddots & \ddots & u_{n-1} \\[0.3em] 0 & 0 & 0 &l_{n-1} & m_{n} \end{bmatrix} C = \begin{bmatrix} v_1 \\[0.3em] v_2 \\[0.3em] \vdots \\[0.3em] \vdots \\[0.3em] v_n \end{bmatrix} \tag{8.18} \end{equation}\]

Penyelesaian Persamaan (8.18) akan menghasilkan vektor koefisien \(c\), sehingga koefisien \(b\) dapat dihitung menggunakan Persamaan (8.19).

\[\begin{equation} b_i=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}-\frac{x_{i+1}-x_i}{3}\left(2c_i+c_{i+1}\right) \tag{8.19} \end{equation}\]

Dengan menggunakan vektor \(C\) yang sudah diketahui, koefisien \(d\) dapat dihitung menggunakan Persamaan (8.20).

\[\begin{equation} d_i=\frac{c_{i+1}-c_i}{3\left(x_{i+1}-x_i\right)} \tag{8.20} \end{equation}\]

Algoritma Interpolasi Spline Kubik

1. Tentukan set titik berpasangan \(\left(x,y\right)\).
2. Tentukan koefisien \(a_i\) menggunakan Persamaan (8.10), dimana \(a_i = y_i\).
3. Hitung elemen diagonal bawah (\(l_i\)) dan elemen diagonal atas (\(u_i\)) matriks tridiagonal menggunakan Persamaan (8.15), dimana \(u_0=l_n=0\).
4. Hitung elemen diagonal utama (\(m_i\)) menggunakan Persamaan (8.16), dimana \(d_0=d_n=1\).
5. Hitung elemen vektor \(V\) menggunakan Persamaan (8.17), dimana \(v_0=v_n=0\).
6. Susunlah elemen \(l_i\), \(u_i\), dan \(m_i\) menjadi matriks tridiagonal \(A\).
7. Deifinisikan persamaan linier seperti pada Persamaan (8.14)
8. Selesaikan sistem persamaan linier pada Persamaan (8.14) sehingga diperoleh vektor \(C\) yang merupakan kumpulan koefisien \(c\).
9. Hitung koefisien \(b_i\) menggunakan Persamaan (8.19)
10. Hitung koefisien \(d_i\) menggunakan Persamaan (8.20).
11. Bentuk seri persamaan polinomial menggunakan elemen koefisien \(a\), \(b\), \(c\), dan \(d\) yang telah dihitung.
12. Untuk melakukan ekstrapolasi dengan titik observasi diluar rentang titik diketahui, gunakan persamaan polinomial yang berada pada bagian ujung terdekat dengan nilai \(x\) yang hendak dicari nilai \(y\)-nya.

Berdasarkan algoritma tersebut, kita dapat menyusun suatu fungsi pada R untuk mencari seri persamaan polinomial derajat tiga. Fungsi tersebut adalah sebagai berikut:

cubic_spline <- function(x, y){ n <- length(x) d_vec <- b_vec <- a_vec <- rep(0, n-1) vec <- rep(0, n) delta_x <- delta_y <- rep(0, n-1) ## Menghitung nilai selisih dan vektor A for(i in 1:(n-1)){ a_vec[i] <- y[i] delta_x[i] <- x[i+1] - x[i] delta_y[i] <- y[i+1] - y[i] } ## Menyusun matriks tridiagona Au <- c(0, delta_x[2:(n-1)]) Am <- c(1, 2*(delta_x[1:(n-2)]+delta_x[2:(n-1)]), 1) Al <- c(delta_x[1:(n-2)], 0) vec[0] <- vec[n] <- 0 for(i in 2:(n-1)) vec[i] <- 3 * (delta_y[i]/delta_x[i] - delta_y[i-1]/delta_x[i-1]) ## penyelesaian tridiagonal matriks nm <- length(Am) Al <- c(NA , Al) ### forward sweep Au[1] <- Au[1] / Am[1] vec[1] <- vec[1] / Am[1] for(i in 2:(n - 1)){ Au[i] <- Au[i] / (Am[i] - Al[i] * Au[i - 1]) vec[i] <- (vec[i] - Al[i] * vec[i - 1]) / (Am[i] - Al[i] * Au[i - 1]) } vec[n] <- (vec[n] - Al[n] * vec[n - 1])/ (Am[n] - Al[n] * Au[n - 1]) ### backward sweep c_vec <- rep.int (0, n) c_vec[n] <- vec[n] for(i in (n - 1) :1) c_vec[i] <- vec[i] - Au[i] * c_vec[i + 1] ## Hitung vektor B dan D dari vektor C for(i in 1:(n-1)){ b_vec[i] <- (delta_y[i]/delta_x[i])- (delta_x[i]/3)*(2*c_vec[i]+c_vec[i+1]) d_vec[i] <- (c_vec[i+1]-c_vec[i]) / (3*delta_x[i]) } return(list(a = a_vec, b = b_vec, c = c_vec[-n], d = d_vec)) }

Contoh 8.8 Tentukan persamaan-persamaan spline kubik yang dihasilkan dari titik observasi yang ditampilkan pada Contoh 8.5?

Jawab:

Untuk menentukan persamaan-persamaan linier yang menghubungkan setiap titik, kita akan menggunakan fungsi cubic_spline() yang telah kita buat. Berikut adalah sintaks yang digunakan:

x <- c(-2, -1, 0, 1) y <- c(-1, -2, -1, 2) cubic_spline(x, y)

Sebagai contoh untuk domain \(\left[0,1\right]\), persamaan spline kubiknya adalah \(-0,4x^3+1,2x^2+2,2x-1\). Jika persamaan-persamaan yang terbentuk tersebut divisualisasikan akan terlihat seperti pada Gambar 8.4.

Gambar 8.4: Interpolasi spline kubik empat titik (Sumber:Howard, 2017).

Pada R juga terdapat fungsi splinefun() untuk melakukan interpolasi spline. Format fungsi tersebut adalah sebagai berikut:

Catatan:

• x,y : vektor titik yang akan dilakukan interpolasi
• method : spesifikasi jenis spline yang digunakan. Nilai yang mungkin antara lain: "fmm", "natural", "periodic", "monoH.FC" dan "hyman".
• ties : metode untuk menangani nilai imbang pada titik yang akan diinterpolasi.

Untuk melakukan interpolasi spline kubik, metode yang digunakan adalah "natural". Berikut adalah contoh interpolasi menggunakan kembali data pada Contoh 8.5:

FUngsi splinefun() menghasilkan nilai yang sama persis dengan proses interpolasi yang telah kita lakukan.

Pada studi kasus kali ini, kita akan membahas teknik mengisi nilai yang hilang pada data runtun waktu. Terdapat banyak teknik untuk yang dapat digunakan untuk mengisi data yang hilang. Salah satu teknik yang dapat digunakan adalah dengan melakukan interpolasi.

Pengisian data hilang pada data runtun waktu (time series) dapat dilakukan dengan berbagai cara sesuai dengan situasi yang dihadapi. Pengisian dapat menggunakan nilai rata-rata jika data memiliki pola white noise, observasi terakhir atau observasi dimasa mendatang, dan interpolasi linier.

Data yang digunakan pada contoh kasus kali ini adalah data airquality. Dataset tersebut merupakan data kualitas udara bulan Mei sampai September 1973 yang ada di New York. Berikut adalah ringkasan data airquality tersebut:

Pada contoh kasus kali ini kita akan mencoba melakukan pengisian data hilang pada data Solar.R pada dataset airquality. Langkah pertama yang perlu dilakukan adalah membuat objek data runtun waktu pada data tersebut.

Solar <- ts(airquality[,"Solar.R"])

Visualisasi data tersebut disajikan pada Gambar 8.5.

Gambar 8.5: Visualisasi variabel radiasi matahari pada dataset airquality.

Berdasarkan visualisasi tersebut terdapat garis yang terputus yang menunjukkan data yang hilang. Agar garis tersebut dapat tersambung, kita perlu melakukan pengisian nilai yang hilang pada data tersebut. Berikut adalah sintaks yang digunakan untuk melakukan pengisian nilai hilang tersebut:

Berikut adalah visualisasi menggunakan metode nilai rata-rata:

Gambar 8.6: Visualisasi data radiasi matahari menggunakan metode nilai rata-rata.

Secara berturut-turut berikut adalah visualisasi dari metode locf, nocb, dan interpolasi linier:

Gambar 8.7: Visualisasi data radiasi matahari menggunakan metode locf.

Gambar 8.8: Visualisasi data radiasi matahari menggunakan metode nocb.

Gambar 8.9: Visualisasi data radiasi matahari menggunakan metode interpolasi linier.

Pemilihan metode interpolasi mana yang sesuai akan berbeda pada setiap situasi dan jenis data yang akan dilakukan interpolasi. Interpolasi data runtun waktu pada bidang lingkungan umumnya menggunakan metode interpolasi nilai rata-rata dan linier dan tidak menutup kemungkinan interpolasi dengan metode lain yang telah dijelaskan pada buku ini dapat pula digunakan.

Pada situasi dimana data membentuk pola white noises (pola acak disekitar nilai rata-rata dan memiliki varians yang konstan) seperti yang ditunjukkan variabel Solar.R, interpolaasi dengan nilai rata-rata cukup sesuai untuk digunakan untuk mengisi nilai hilang (missing value) pada data runtun waktu tersebut.

1. Howard, J.P. 2017. Computational Methods for Numerical Analysis with R. CRC Press.
2. Kreyszig, E. 2011. Advanced Engineering Mathematics, 10th Edition. John Wiley & Sons.
3. Suparno, S. 2008. Komputasi untuk Sains dan Teknik Edisi II. Departemen Fisika-FMIPA Universitas Indonesia.

1. Diberikan data titik (3,5), (0,-2), dan (4,1). Tentukan persamaan polinomial untuk melakukan interpolasi pada ketiga titik tersebut!
2. Jika diberikan 13 titik observasi, apakah saudara cenderung akan menggunakan interpolasi polinomial atau spline? jelaskan alasan saudara?
3. Bentuklah kembali fungsi poly_inter() dengan menambahkan metode evaluasi polinomial ke dalam fungsi tersebut!
4. Pada Latihan No.1, bentuklah persamaan spline kubik menggunakan titik observasi tersebut!

Pada Chapter 9, penulis akan menjabarkan mengenai metode numerik untuk melakukan diferensiasi dan integrasi pada suatu fungsi. Adapun yang akan dibahas pada Chapter ini antara lain:

• Metode Beda Hingga
• Metode Integrasi Newton-Cotes
• Metode Integrasi Kudratur Gauss
• Metode Integrasi Adaptif
• Metode Integrasi Romberg
• Metode Integrasi Monte Carlo
• Studi Kasus

Diferensiasi merupakan proses mencari slope suatu garis pada titik yang diberikan. Secara umum proses diferensiasi dinyatakan melalui Persamaan (9.1).

\[\begin{equation} f'\left(x\right) \approx \frac{f\left(x+h\right)-f\left(x\right)}{h} \tag{9.1} \end{equation}\]

Kita dapat menyatakan secara formal proses diferensiasi sebagai limit Persamaan (9.1) dimana \(h\) mendekati nol. Jadi kita ingin membuat nilai \(h\) sekecil mungkin untuk memperoleh pendekatan terbaik terhadap nilai turunan suatu fungsi. Kita membatasi nilai \(h\) pada sejumlah nilai yang masuk akal untuk mencegah pembagian dengan nilai yang tidak biasa. Kita juga harus memastikan \(f\left(x\right)\) dan \(f\left(x+h\right)\) terpisah cukup jauh untuk mencegah floating point round off error mempengaruhi proses substraksi.

Terdapat 3 buah metode untuk memperoleh turunan pertama suatu fungsi dengan menggunakan metode numerik, yaitu: metode selisih maju, metode selisih mundur, dan metode selisih tengah. Error pada ketiga metode numerik tersebut ditaksir menggunakan deret Taylor. Persamaan (9.2) dan Persamaan (9.3) menunjukkan persamaan untuk memperoleh turunan pertama dan taksiran error menggunakan metode selisih maju dan metode selisih mundur.

\[\begin{equation} f'\left(x\right) = \frac{f\left(x+h\right)-f\left(x\right)}{h} - \frac{h}{2} f''\left(c\right) \tag{9.2} \end{equation}\]

\[\begin{equation} f'\left(x\right) = \frac{f\left(x\right)-f\left(x-h\right)}{h} - \frac{h}{2} f''\left(c\right) \tag{9.3} \end{equation}\]

Metode nilai tengah menggunakan ukuran langkah \(h\) dua kali dibandingkan dengan 2 metode lainnya. Error yang dihasilkan juga berbeda dengan kedua metode sebelumnya, dimana error dihasilkan dari pemotongan turunan ketiga pada deret Taylor. Secara umum metode selisih tengah memiliki akurasi yang lebih baik dibandingkan kedua metode sebelumnya karena metode ini mempertimbangkan dua sisi untuk memeriksa nilai \(x\). Persamaan (9.4) merupakan persamaan untuk memperoleh nilai turunan pertama suatu fungsi dan estimasi error menggunakan deret Taylor.

\[\begin{equation} f'\left(x\right) = \frac{f\left(x+h\right)-f\left(x-h\right)}{2h} - \frac{h^2}{6} f'''\left(c\right) \tag{9.4} \end{equation}\]

Bagaimana menentukan \(h\)? beberapa literatur menggunakan pendekatan machine error \(\epsilon\) berdasarkan program yang digunakan untuk melakukan proses perhitungan. Metode selisih maju dan selisih mundur menggunakan pendekatan yang ditunjukkan pada Persamaan (9.5).

\[\begin{equation} h^{\ast}=x\sqrt{\epsilon} \tag{9.5} \end{equation}\]

Untuk metode selisih tengah pendekatan nilai \(h\) menggunakan Persamaan (9.6).

\[\begin{equation} h^{\ast}=x\sqrt[3]{\epsilon} \tag{9.6} \end{equation}\]

Kita dapat menggunakan Persamaan (9.2) sampai Persamaan (9.4) untuk membentuk sebuah program yang digunakan untuk menghitung turunan pertama suatu fungsi. Sintaks yang digunakan adalah sebagai berikut:

Contoh 9.1 Hitunglah turunan pertama persamaan berikut menggunakan metode selisih titik tengah pada x =1 dan nilai h=0,05!

\[ f\left(x\right)=e^{-x}\sin\left(2x\right)+1 \]

Jawab:

Untuk menghitung turunan pertama menggunakan metode selisih tengah, kita dapat menggunakan Persamaan (9.4). Berikut adalah proses perhitungannya:

\[ f'\left(1\right) = \frac{f\left(1+0.05\right)-f\left(1-0,05\right)}{2\times0.05}=--0.6390352 \]

Dengan menggunakan fungsi findiff(), hasil yang diperoleh adalah sebagai berikut:

findiff(function(x) exp(-x)*sin(2*x)+1, x=1, h=0.05, method="central")

Kita dapat memperkecil nilai \(h\) untuk memperoleh akurasi yang lebih baik berdasarkan pendekatan Persamaan (9.6).

Penyelesaian persamaan matematik dalam bidang Teknik Lingkungan pada umumnya tidak hanya melibatkan turunan pertama, pada penyelesaian persamaan difusi umumnya menggunakan turunan kedua. Persamaan (9.7) merupakan pendekatan numerik untuk memperoleh nilai turunan kedua suatu persamaan dengan pendekatan deret Taylor.

\[\begin{equation} f''\left(x\right)=\frac{f\left(x+h\right)-2f\left(x\right)+f\left(x-h\right)}{h^2}-\frac{h^2}{12}f^{\left(4\right)}\left(c\right) \tag{9.7} \end{equation}\]

Fungsi findiff2() merupakan fungsi yang digunakan untuk menghitung turunan kedua suatu persamaan yang didasarkan pada Persamaan (9.7).

findiff2 <- function(f, x, h){ return((f(x+h)-2*f(x)+f(x-h))/(h^2)) }

Kita dapat menghitung kembali turunan kedua fungsi pada Contoh 9.1 menggunakan fungsi findiff2(). Berikut adalah sintaks yang digunakan:

findiff2(function(x){ exp(-x)*sin(2*x)+1 }, x=1, h=0.05)

Terdapat sejumlah fungsi R yang dapat digunakan untuk menghitung turunan suatu persamaan matematik. Fungsi-fungsi tersebut tersedia dalam sejumlah Paket, baik base package maupun yang berasal dari Paket lainnya.

Paket standar yang sering digunakan untuk melakukan taksiran numerik turunan suatu fungsi adalah Paket numDeriv. Pada Paket tersebut terdapat fungsi grad() yang digunakan untuk menaksir turunan pertama suatu persamaan. Format fungsi tersebut adalah sebagai berikut:

grad(func, x, method="Richardson", method.args=list(), ...)

Catatan:

• func: Fungsi persamaan matematik yang akan dicari turunannya.
• x: Lokasi atau titik yang akan dicari gradiennya
• method: Metode estimasi yang digunakan. Metode yang dapat digunakan antara lain:
• “simple”: metode selisih maju dengan \(h = 10^{-4}\)
• “Richardson”: metode interpolasi Richardson
• “complex”: complex-step derivative approach dan dapat digunakan untuk persamaan complex-differentiable.
• method.args: argumen tambahan yang digunakan bersama dengan argumen method. Jika metode yang digunakan adalah “simple”, nilai \(h\) dapat dispesifikasikan pada method.args jika diinginkan nilai \(h\) lainnya.

Berikut adalah contoh penerapan fungsi grad() untuk memperoleh turunan pertama persamaan matematik pada Contoh 9.1:

Terdapat sejumlah fungsi pada Paket pracma yang dapat digunakan untuk melakukan diferensiasi suatu persamaan matematik. FUngsi-fungsi yang dapat digunakan untuk melakukan diferensiasi sederhana antara lain: fderiv(), numderiv(), numdiff(), dan grad().

Fungsi fderiv() dapat digunakan untuk melakukan diferensiasi orde pertama sampai dengan orde tinggi. Perlu dicatat bahwa diferensiasi cenderung kurang akurat jika orde diferensiasi semakin tinggi. Format yang digunakan untuk melakukan diferensiasi menggunakan fungsi fderiv() adalah sebagai berikut:

Catatan:

• f: Fungsi persamaan matematik yang akan dicari turunannya.
• x: Lokasi atau titik yang akan dicari gradiennya
• n: Orde diferensiasi yang digunakan. Orde diferensiasi yang dapat digunakan adalah 1 sampai 8.
• method: Metode estimasi yang digunakan. Metode yang dapat digunakan antara lain:
• “central”: metode titik pusat
• “forward”: metode selisih maju
• “bacward”: metode selisih mundur.
• …: argumen tambahan fungsi f.

Berikut adalah contoh penerapan fungsi fderiv() untuk memperoleh turunan pertama dan kedua persamaan matematik pada Contoh 9.1:

Fungsi numderiv() menggunakan ekstrapolasi Richardson untuk melakukan taksiran turunan suatu persamaan matematik. Berbeda dengan fungsi lainnya, fungsi numderiv() tidak hanya menampilkan hasil diferensiasi, fungsi ini juga menampilkan error absolut, error relatif, dan jumlah iterasi yang berlangsung. Berikut adalah format fungsi yang digunakan:

numderiv(f, x0, maxiter = 16, h = 1/2, ..., tol = .Machine$double.eps)

Catatan:

• f: Fungsi persamaan matematik yang akan dicari turunannya.
• x0: Lokasi atau titik yang akan dicari gradiennya
• maxiter: Iterasi maksimum yang digunakan.
• h: step size yang digunakan
• tol: toleransi error yang digunakan.

Berikut adalah contoh penerapan fungsi numderiv() untuk memperoleh turunan pertama persamaan matematik pada Contoh 9.1:

Fungsi numderiv() memiliki keterbatasan dalam penggunaannya. Argumen x0 yang digunakan haruslah angka numerik tunggal. Fungsi numdiff() mengatasi keterbatasan tersebut. Fungsi ini dapat menerima input berupa vektor, sehingga dapat digunakan untuk mencari nilai turunan pada sejumlah titik. Selain itu, output fungsi ini lebih sederhana, dimana hanya menampilkan hasil diferensiasinya saja. Berikut adalah format fungsi numdiff():

numdiff(f, x, maxiter = 16, h = 1/2, ..., tol = .Machine$double.eps)

Catatan:

• f: Fungsi persamaan matematik yang akan dicari turunannya.
• x: Vektor titik yang akan dicari gradiennya
• maxiter: Iterasi maksimum yang digunakan.
• h: step size yang digunakan
• tol: toleransi error yang digunakan.

Berikut adalah contoh penerapan fungsi numdiff() untuk memperoleh turunan pertama persamaan matematik pada Contoh 9.1:

Fungsi grad() pada Paket pracma berbeda dengan yang digunakan pada Paket numDeriv. Perbedaan utama fungsi pada kedua Paket tersebut adalah metode estimasi yang digunakan untuk menghitung turunan pertama suatu persamaan matematik. Pada Paket pracma, metode yang digunakan adalah metode titik pusat, sedangkan pada Paket numDeriv metode yang digunakan adalah metode selisih maju, ekstrapolasi Richardson, dan complex. Format fungsi grad() pada Paket pracma adalah sebagai berikut:

grad(f, x0, heps = .Machine$double.eps^(1/3), ...)

Catatan:

• f: Fungsi persamaan matematik yang akan dicari turunannya.
• x0: Titik yang akan dicari gradiennya
• heps: step size yang digunakan
• …: Argumen lain yang digunakan pada fungsi f.

Berikut adalah contoh penerapan fungsi grad() untuk memperoleh turunan pertama persamaan matematik pada Contoh 9.1:

grad(function(x){exp(-x)*sin(2*x)+1}, x0 = 1)

Metode integrasi Newton-Cotes secara umum merupakan metode integrasi yang dilakukan dengan membagi area di bawah kurva suatu fungsi menjadi beberapa panel dengan terlebih dahulu menetapkan batas atas dan batas bawah interval. Integral atau luas area di bawah kurva ditentukan berdasarkan jumlah luas panel yang digunakan untuk mendekati luas area di bawah kurva.

Terdapat beberapa metode yang akan penulis jelaskan pada sub-Chapter ini. Metode-metode tersebut antara lain:

• Metode integral Riemann
• Metode trapezoida
• Metode Simpson 1/3
• Metode Simpson 3/8

Metode integral Riemann dilakukan dengan membagi interval di bawah kurva suatu fungsi matematik sebanyak \(m\) subinterval sama besar. Pada setiap subinterval dibentuk persegi panjang setinggi kurva pada setiap titik tengah persegi panjang tersebut. Area setiap subinterval diperoleh dengan mengalikan panjang dan lebar masing-masing persegi panjang. Jumlah masing-masing area tersebut digunakan untuk menaksir interval integral suatu fungsi dengan interval tertentu. Fungsi proses integrasi menggunakan metode titik tengah dapat dituliskan pada Persamaan (9.8).

\[\begin{equation} \int_a^bf\left(x\right)dx\ \approx\sum_{i=1}^mf\left(i\ \frac{\left|b-a\right|}{m}-\frac{\left|b-a\right|}{2m}\right)\ \frac{\left|b-a\right|}{m} \tag{9.8} \end{equation}\]

dimana \(b\) dan \(a\) masing-masing merupakan batas atas dan bawah interval kurva yang hendak dihitung integralnya.

Error dari metode ini dapat diestimasi menggunakan Persamaan (9.9).

\[\begin{equation} \int_a^bh\left(x\right)dx=-\ \frac{\left(b-a\right)^3}{24m^2}f^{\left(2\right)}\left(\xi\right) \tag{9.9} \end{equation}\]

dimana \(\xi\) merupakan nilai antara \(a\) dan \(b\).

Contoh 9.2 Hitunglah intergral fungsi di bawah ini menggunakan metode integral Reimann dengan interval 0 sampai 1 dan jumlah panel 2 dan 4!

\[ \int_0^1 x^2 dx \]

Jawab:

Fungsi pada Contoh 9.2 dapat diselesaikan menggunakan metode analitik. Penyelesaian analitik fungsi tersebut adalah sebagai berikut:

\[ \int_0^1 x^2 dx = \left[\frac{x^3}{3}\right]_{_0}^{^1}=\frac{1^3}{3}-\frac{0^3}{3}=0,333... \]

Penyelesaian numerik menggunakan metode titik tengah dengan jumlah panel 2 dapat dilakukan dengan menentukan lokasi titik tengah kedua panel. Berdasarkan interval fungsi dapat kita tentukan titik tengah kedua panel berada pada \(x=0,25\) dan \(x=0,75\). Perhitungan dilakukan seperti berikut:

\[ \int_0^1 x^2 dx \approx \left(f\left(0,25\right)+f\left(0,75\right)\right)\ \frac{1-0}{2}=\frac{0,25^2+0,75^2}{2}=0,3125 \] Untuk meningkatkan akurasi dari nilai yang dihasilkan, jumlah panel dapat ditingkatkan. Untuk jumlah panel 4, titik tengah berada pada \(x=\left\{0,125;0,375;0,625;0,875\right\}\).

\[ \int_0^1 x^2 dx \approx \left(f\left(0,125\right)+f\left(0,375\right)+f\left(0,625\right)+f\left(0,875\right)\right)\ \frac{1-0}{4}=0,328125 \]

Visualisasi proses integrasi dengan metode Riemann dapat dilihat pada Gambar 9.1.

Gambar 9.1: Visualisasi integral Riemann dengan 2 panel dan 4 panel (sumber:Howard, 2017).

Berdasarkan Persamaan (9.8), kita dapat mengembangkan fungsi R yang dapat digunakan untuk melakukan perhitungan integral Riemann. Sintaks fungsi tersebut adalah sebagai berikut:

Kita akan menghitung kembali fungsi pada Contoh 9.2 dengan menggunakan jumlah panel 2, 4 dan 100. Berikut adalah sintaks yang digunakan:

# m=2 riemann(function(x) x^2, a=0, b=1, m=2)

# m=4 riemann(function(x) x^2, a=0, b=1, m=4)

# m=100 riemann(function(x) x^2, a=0, b=1)

Berdasarkan teori yang telah dipaparkan sebelumnya, kita ketahui bahwa untuk memperoleh nilai pendekatan integral yang sebenarnya kita dapat meningkatkan jumlah panel yang digunakan. Untuk mengetahui jumlah panel minimum yang diperlukan untuk memperoleh hasil integrasi yang stabil, kita akan melakukan simulasi menggunakan data yang disajikan pada Contoh 9.2 dengan memvariasikan jumlah panel yang akan digunakan. Pada simulasi yang akan dilakukan kita akan coba memvariasikan jumlah panel dari 2 hingga 100. Berikut adalah sintaks yang digunakan:

Gambar 9.2: Visualisasi simulasi pemilihan jumlah panel minimum metode integrasi Riemann.

Berdasarkan hasil simulasi dapat disimpulkan jumlah panel minimum yang diperlukan untuk memperoleh hasil integrasi yang stabil kira-kira sebesar \(m=40\).

Pendekatan trapezoida dilakukan dengan melakukan pendekatan area dibawah kurva fungsi \(y=f\left(x\right)\) dengan subinterval \(\left[x_i,x_{i+1}\right]\) menggunakan trapesium. Untuk memahami pendekatan yang digunakan pembaca dapat memperhatikan Gambar 9.3.

Gambar 9.3: Visualisasi integragrasi numerik menggunakan metode tapezoida (sumber: Jones et.al., 2014).

Fungsi proses integrasi menggunakan metode trapezoida dapat dituliskan pada Persamaan (9.10).

\[\begin{equation} \int_a^bf\left(x\right)dx\ = \lim\limits_{m \to \infty} \sum_{i=1}^m\frac{\left(c_{i+1}-c_i\right)\times\left(f\left(c_{i+1}\right)+f\left(c_i\right)\right)}{m} \tag{9.10} \end{equation}\]

dimana

\[ c_i=a+\frac{\left(b-a\right)}{n}i \] \(n\) merupakan nilai subinterval dan \(m\) merupakan jumlah panel trapesium yang digunakan.

Error dari metode ini dapat diestimasi menggunakan Persamaan (9.11).

\[\begin{equation} \int_a^bh\left(x\right)dx=-\ \frac{\left(b-a\right)^3}{12m^2}f^{\left(2\right)}\left(\xi\right) \tag{9.11} \end{equation}\]

dimana \(\xi\) merupakan nilai antara \(a\) dan \(b\).

Contoh 9.3 Hitung kembali nilai intergrasi persamaan pada Contoh 9.2 menggunakan metode trapezoida dengan jumlah panel m=2!

Jawab:

Penyelesaian numerik menggunakan trapezoida dengan jumlah panel 2 dapat dilakukan dengan menentukan lokasi titik evaluasi. Berdasarkan Gambar 9.4, terdapat 3 batas subinterval yaitu pada \(a\),\(\frac{\left(b-a\right)}{2}\), dan \(b\). Perhitungan intergral menggunakan ketiga titik evaluasi tersebut adalah sebagai berikut:

\[ \int_0^1 x^2 dx \approx \frac{\left(0,5\right)\left(0,25+1,25\right)}{2}=0,375 \]

Gambar 9.4: Visualisasi integrasi metode trapezoida dengan 2 panel (sumber:Howard, 2017).

Berdasarkan Persamaan (9.10), kita dapat mengembangkan fungsi R yang dapat digunakan untuk melakukan perhitungan integral metode trapezoida. Sintaks fungsi tersebut adalah sebagai berikut:

Kita dapat menghitung kembali intergral persamaan pada Contoh 9.2 menggunakan fungsi trap() yang telah dibuat. Berikut adalah sintaks yang digunakan:

trap(function(x)x^2, a=0, b=1, m=2)

Untuk mengetahui jumlah panel minimum yang diperlukan untuk memperoleh hasil integrasi yang stabil pada persamaan tersebut, kita akan kembali melakukan simulasi menggunakan variasi jumlah panel yang digunakan. Dalam simulasi variasi jumlah panel yang digunakan adalah 2 sampai 100. Berikut adalah sintaks yang digunakan:

Gambar 9.5: Visualisasi simulasi pemilihan jumlah panel minimum metode integrasi trapezoida.

Berdasarkan hasil simulasi diperoleh nilai panel minimum sebesar \(m=20\). Hasil yang diperoleh tersebut menujukkan bahwa metode trapezoida lebih efisien dalam proses komputasi dibandingkan metode Riemann.

Metode Simpson membagi subinterval \(\left[a,b\right]\) menjadi \(n\) subinterval, dimana \(n\) merupakan bilangan genap. Untuk setiap pasang subinterval, luas area di bawah fungsi \(f\left(x\right)\) ditaksir menggunakan polinomial berderajat 2.

Misalkan \(u<v<w\) merupakan titik sembarang pada suatu fungsi yang akan dicari integralnya yang terpisah sejauh \(h\). Untuk \(x\in\left[u,w\right]\) kita ingin menaksir \(f\left(x\right)\) menggunakan parabola yang melalui titik \(\left(u, f\left(u\right)\right)\), \(\left(v, f\left(v\right)\right)\), dan \(\left(w, f\left(w\right)\right)\). Terdapat tepat 1 parabola \(p\left(x\right)\) yang dapat dibentuk dari ketiga titik koordinat tersebut yang ditunjukkan melalui Persamaan (9.12).

\[\begin{equation} p\left(x\right)=f\left(u\right)\frac{\left(x-v\right)\left(x-w\right)}{\left(u-v\right)\left(u-w\right)}+f\left(v\right)\frac{\left(x-u\right)\left(x-w\right)}{\left(v-u\right)\left(v-w\right)}+f\left(w\right)\frac{\left(x-u\right)\left(x-v\right)}{\left(w-u\right)\left(w-v\right)} \tag{9.12} \end{equation}\]

Sebagai taksiran luas di bawah kurva \(y=f\left(x\right)\) digunakan \(\int_{w}^u p\left(x\right)dx\). Hasil integrasi kurva Persamaan (9.12) disajikan pada Persamaan (9.13).

\[\begin{equation} \int_w^up\left(x\right)dx=\frac{h}{3}\left(f\left(u\right)+4f\left(v\right)+f\left(w\right)\right) \tag{9.13} \end{equation}\]

Sekarang asumsikan \(n\) merupakan bilangan genap, maka kita perlu menambahkan taksiran untuk subinterval \(\left[x_{2i},x_{2i+2}\right]\) untuk memperoleh taksiran \(S\) pada integral \(\int_{a}^b f\left(x\right)dx\) yang disajikan pada Persamaan (9.14).

\[\begin{equation} S\approx\frac{h}{3}\left(f_0+4\sum_{i=1,3,5.\dots}^{n-1}f_i+2\sum_{i=2,4,6,\dots}^{n-2}f_i+f_n\right) \tag{9.14} \end{equation}\]

Persamaan (9.14) disebut sebagai kaidah Simpson 1/3 karena terdapat koefisien 1/3 pada bagian depan persamaan tersebut. Persamaan tersebut juga mudah diingat mengingat pola koefisien persamaan tersebut adalah \(1,4,2,4,2,\dots ,2,4,1\). Namun penggunaan kaidah 1/3 Simpson mengharuskan jumlah subinterval \(n\) genap. Kondisi tersebut jelas berbeda dengan metode trapezoida yang tidak mensyaratkan jumlah selang.

Error dari metode Simpson 1/3 dapat dihitung menggunakan Persamaan (9.14).

\[\begin{equation} \int_a^bh\left(x\right)dx=-\ \frac{\left(b-a\right)^5}{180m^4}f^{\left(4\right)}\left(\xi\right) \tag{9.15} \end{equation}\]

dimana \(\xi\) merupakan nilai antara \(a\) dan \(b\).

Algoritma Metode Simpson

1. Tentukan fungsi \(f\left(x\right)\) dan selang integrasinya \(\left[a,b\right]\).
2. Tentukan jumlah subinterval \(n\).
3. Hitung nilai selang subinterval \(h\), \(h=\frac{b-a}{n}\).
4. Tentukan awal integrasi \(x_0=a\) dan akhir \(x_n=b\) dan hitung nilai \(f\left(a\right)\) dan \(f\left(b\right)\).
5. Untuk \(x=1,2,\dots,n-1\),

• jika ganjil, hitung: \(4\times f\left(x\right)\)
• jika genap, hitung: \(2\times f\left(x\right)\)

1. Jumlahkan nilai-nilai taksiran tersebut menggunakan Persamaan (9.14).

Berdasarkan algoritma tersebut, kita dapat membentuk fungsi simpson() yang dapat digunakan untuk melakukan integrasi menggunakan metode Simpson 1/3. Berikut adalah sintaks yang digunakan:

Contoh 9.4 Hitung integral persamaan di bawah ini dengan menggunakan jumlah panel m=8!

\[ \int_{0}^1 \frac{1}{1+x}dx \]

Jawab:

lebar selang \(h\) dapat dihitung seperti berikut:

\[ h = \frac{1-0}{8}=0,125 \]

Integral persamaan tersebut selanjutnya dapat dihitung menggunakan Persamaan (9.14):

\[ S\approx\frac{0,125}{3}\left(f\left(0\right)+4f\left(0,125\right)+2f\left(0,25\right)\dots+f\left(1\right)\right)\approx 0,69412 \]

Fungsi simpson() juga menghasilkan nilai yang serupa dengan perhitungan manual yang telah dilakukan. Berikut adalah sintaks yang digunakan:

simpson(function(x)1/(1+x), a=0, b=1, m=8)

Jika pada metode Simpson 1/3 digunakan pendekatan polinomial berderajat 2 untuk mencari luas dibawah kurva, pada metode Simpson 3/8 digunakan pendekatan polinomial berderajat 3 untuk memperoleh hasil yang lebih baik. Bentuk umum integrasi yang digunakan disajikan pada Persamaan (9.16).

\[\begin{equation} \int_a^bf\left(x\right)dx\approx\frac{3h}{8}\left(f_0+3\sum_{i=1;i\ne3,6,9,\dots}^{n-1}f_i+2\sum_{i=3,6,9,\dots}^{n-3}f_i+f_n\right) \tag{9.16} \end{equation}\]

Error dari metode ini dapat diestimasi menggunakan Persamaan (9.17).

\[\begin{equation} \int_a^bh\left(x\right)dx=-\ \frac{\left(b-a\right)^5}{80m^4}f^{\left(5\right)}\left(\xi\right) \tag{9.17} \end{equation}\]

Algoritma Metode Simpson 3/8

1. Tentukan fungsi \(f\left(x\right)\) dan selang integrasinya \(\left[a,b\right]\).
2. Tentukan jumlah subinterval \(n\).
3. Hitung nilai selang subinterval \(h\), \(h=\frac{b-a}{n}\).
4. Tentukan awal integrasi \(x_0=a\) dan akhir \(x_n=b\) dan hitung nilai \(f\left(a\right)\) dan \(f\left(b\right)\).
5. Untuk \(x=1,2,\dots,n-1\),

• jika bukan kelipatan 3, hitung: \(3\times f\left(x\right)\)
• jika kelipatan 3, hitung: \(2\times f\left(x\right)\)

1. Jumlahkan nilai-nilai taksiran tersebut menggunakan Persamaan (9.16).

Berdasarkan algoritma tersebut, fungsi R dapat disusun untuk melakukan komputasi metode simpson 3/8. Berikut adalah sintaks fungsi tersebut:

Contoh 9.5 Hitung kembali integral persamaan yang disajikan pada Contoh 9.4 menggunakan fungsi simpson38() yang telah dibuat sebelumnya!

Jawab:

Berikut adalah sintaks yang digunakan untuk melakukan proses integrasi Simpson 3/8 menggunakan fungsi simpson38():

simpson38(function(x)1/(1+x), a=0, b=1, m=9)

Terdapat sejumlah Paket yang menyediakan fungsi yang dapat digunakan untuk melakukan proses integrasi suatu fungsi pada R. Paket pracma menyediakan fungsi-fungsi seperti trapz() dan cotes().

Fungsi trapz() merupakan fungsi yang digunakan untuk melakukan integrasi dengan pendekatan trapesium. Penggunaan fungsi tersebut untuk melakukan perhitungan integral tidak dapat dilakukan secara langsung, kita perlu membuat terlebih dahulu seri koordinat x dan koordinat y. Program selanjutnya akan melakukan interpolasi linier terhadap selang yang saling berdekatan. Luas masing-masing panel selang selanjutnya dihitung dan dijumlahkan untuk memperoleh nilai integral. Format fungsi tersebut adalah sebagai berikut:

Catatan:

• x: vektor sumbu x.
• y: vektor sumbu y.

Untuk lebih memahami penerapannya berikut disajikan contoh untuk mencari integral persamaan pada Contoh 9.2:

Fungsi lain yang dapat digunakan untuk menghitung integral suatu fungsi menggunakan metode Newton-Cotes adalah cotes(). Pada fungsi tersebut kita perlu menyatakan jumlah subinterval yang digunakan dan jumlah nodes interpolasi yang digunakan. Jumlah nodes akan menentukan fungsi polinomial pendekatan yang digunakan untuk menghitung luas di bawah suatu fungsi. Jika jumlah nodes diatur menjadi dua, maka fungsi pendekatannya berupa garis linier (metode trapezoida). Jika nodes diatur menjadi 3 maka pendekatannya berupa fungsi kuadrat (metode Simpson 1/3). Secara sederhana derajat polinomial \(n\) memerlukan \(n+1\) nodes. Format fungsi cotes() disajikan sebagai berikut:

cotes(f, a, b, n, nodes, ...)

Catatan:

• f: fungsi yang akan dicari integralnya
• a: batas atas.
• b: batas bawah.
• n: jumlah subinterval atau panel.
• nodes: jumlah nodes yang digunakan untuk interpolasi fungsi pada tiap subinterval.

Untuk lebih memahami penerapannya berikut adalah contoh perhitungan intergral menggunakan persamaan pada Contoh 9.4.

# metode Simpson 1/3 cotes(f, a, b, n=8, nodes=3)

# metode Simpson 3/8 cotes(f, a, b, n=9, nodes=4)

Metode Newton-Cotes sangat powerful, tetapi metode tersebut memiliki dua fitur yang kurang diinginkan. Pertama, kita harus menggunakan evaluasi fungsi sejumlah \(n +1\) untuk hasil presisi dan pada polinomial berderajat \(n\). Hal tersebut mungkin tampak seperti rasio yang baik, tetapi dalam praktiknya, jumlah titik evaluasi akan sering ditingkatkan untuk memperoleh akurasi yang lebih tinggi, namun hasil yang diperoleh juga tidak presisi. Sebagai contoh pembaca dapat melakukan simulasi dengan memvariasikan jumlah penel yang digunakan untuk memperoleh nilai integral sebuah fungsi. Jika kita plotkan hasil yang kita peroleh, grafik yang muncul berupa garis yang berosilasi khusunya pada penggunaan polinomial berderajat tinggi yang menunjukkan hasil yang diperoleh menjadi kurang presisi.

Kelemahan kedua, metode Newton-Cotes memerlukan fungsi terintegrasi untuk dievaluasi pada node yang berjarak sama. Ini benar terlepas dari fungsi yang digunakan. Setiap panel membutuhkan node dengan jarak yang sama di dalamnya. Ini bisa menjadi masalah dengan fungsi periodik, di mana diskontinuitas periodik dapat secara kebetulan mendarat di titik evaluasi. Jika panel Newton-Cotes dapat menyebabkan masalah, kita dapat menggunakan integrasi Gaussian untuk menyelesaikan integral.

Secara umum integrasi Gauss berusaha memperoleh pendekatan luas dibawah kurva fungsi dengan memecah fungsi tersebut menjadi faktor bobot \(c\) dan \(f\left(x\right)\) yang merupakan polinomial pendekatannya. Integral diperoleh melalui hasil kali dari bobot dan fungsi polinomial. Jumlah bobot dan fungsi yang digunakan bergantung pada orde \(n\) polinomial yang akan digunakan untuk mengestimasi integral suatu fungsi. Bentuk umum dari kaidah Gauss tersebut ditampilkan pada Persamaan (9.18).

\[\begin{equation} \int_{-1}^1f\left(x\right)dt\approx c_1f\left(x_1\right)+c_2f\left(x_2\right)+\dots+c_nf\left(x_n\right) \tag{9.18} \end{equation}\]

dimana \(c\) merupakan faktor bobot dan \(x\) merupakan titik evaluasi.

Nilai faktor bobot dan nilai masing-masing titik evaluasi disajikan pada Gambar 9.6.

Gambar 9.6: Tabulasi faktor bobot, titik evaluasi dan galat pemotongan.

Fungsi yang akan dicari nilai integrannya pada umumnya tidak hanya memiliki daerah batas \(\left[-1,1\right]\), sehingga pendekatan kuadratur Gauss tidak dapat digunakan secara langsung pada fungsi yang tidak memiliki batas tersebut. Agar kuadratur Gauss tetap dapat digunakan, fungsi tersebut perlu dilakukan transformasi. Proses transformasi dituliskan pada Persamaan (9.19).

\[\begin{equation} \int_a^bf\left(x\right)dx=\frac{b-a}{2}\int_{-1}^1f\left(\frac{a+b+\left(b-a\right)t}{2}\right)dt \tag{9.19} \end{equation}\]

Algoritma Kuadratur Gauss-Legendre

1. Tentukan fungsi \(f\left(x\right)\) dan selang integrasinya \(\left[a,b\right]\).
2. Lakukan transformasi fungsi tersebut hingga diperoleh fungsi dengan selang \(\left[-1,1\right]\) menggunakan Persamaan (9.19).
3. Tentukan orde polinomial \(n\) yang akan digunakan.
4. Lakukan proses integrasi dengan mengalikan faktor bobot \(c\) dengan \(f\left(x_i\right)\) seperti yang ditujukkan pada Persamaan (9.18).

Kita dapat membangun suatu fungsi pada R untuk melakukan integrasi Gauss berdasarkan algoritma tersebut. Sintaks fungsi tersebut adalah sebagai berikut:

gauss_legendre <- function(f, n){ if(n==2){ c <- c(1,1) x <- c(-0.577350269, 0.577350269) integral <- 0 for(i in 1:n){ integral <- integral + c[i]*f(x[i]) } }else if(n==3){ c <- c(0.555555556,0.888888889,0.555555556) x <- c(-0.774596669,0,0.774596669) integral <- 0 for(i in 1:n){ integral <- integral + c[i]*f(x[i]) } }else if(n==4){ c <- c(0.347854845,0.652145155,0.652145155, 0.347854845) x <- c(-0.861136312,-0.339981044,0.339981044, 0.861136312) integral <- 0 for(i in 1:n){ integral <- integral + c[i]*f(x[i]) } }else if(n==5){ c <- c(0.236926885,0.478628670,0.568888889, 0.478628670,0.236926885) x <- c(-0.906179846,-0.538469310,0, 0.538469310, 0.906179846) integral <- 0 for(i in 1:n){ integral <- integral + c[i]*f(x[i]) } }else if(n==6){ c <- c(0.171324492,0.360761573,0.467913935, 0.467913935,0.360761573,0.171324492) x <- c(-0.932469514,-0.661209386,-0.238619186, 0.238619186, 0.661209386,0.932469514) integral <- 0 for(i in 1:n){ integral <- integral + c[i]*f(x[i]) } }else{ stop("n harus ditentukan") } return(integral) }

Contoh 9.6 Hitunglah integral fungsi berikut menggunakan metode Gauss-Legendre 2 titik!

\[ \int_1^2 \left(x^2+1\right)dx \]

Jawab:

Agar fungsi tersebut dapat dicari nilai integralnya menggunakan metode Gauss-Legendre, fungsi tersebut perlu ditransformasi terlebih dahulu:

\[ \int_1^2 \left(x^2+1\right)dx=0,5\int_{-1}^1 \left[\left(1,5+0,5t\right)^2+1\right]dt \]

Jadi dalam hal ini

\[ f\left(t\right)=\left(1,5+0,5t\right)^2+1 \]

maka

\[ \int_1^2 \left(x^2+1\right)dx=0,5\left[1\times f\left(\frac{1}{\sqrt{3}}\right)+1\times f\left(\frac{1}{\sqrt{3}}\right)\right]\approx 3,333333333 \]

Kita juga dapat menggunakan fungsi gauss_legendre() untuk melakukan integrasi Gauss-Legendre pada dua titik. Berikut adalah sintaks yang digunakan:

0.5*gauss_legendre(function(x)((1.5+0.5*x)^2)+1, n=2)

Fungsi legendre.quadrature() dari Paket gaussquand dapat dijadikan alternatif untuk menghitung integral suatu fungsi menggunakan metode Gauss-Lagendre. Fortmat fungsi tersebut adalah sebagai berikut:

legendre.quadrature(functn, rule, lower=-1, upper=1, weighted = TRUE, ...)

Catatan:

• functn: fungsi yang akan dicari integralnya
• rule: data frame yang terdiri atas orde n aturan kuadratur legendre
• lower: batas atas.
• upper: batas bawah.
• weighted: nilai boolean yang menyatakan apakah bobot fungsi disertakan dalam integran
• …: argumen tambahan functn

Untuk menentukan rule pada fungsi legendre.quadrature(), diperlukan fungsi lain yang tersedia pada Paket gaussquad. Fungsi tersebut adalah legendre.quadrature.rules(). Fungsi tersebut akan menampilkan list data frame berdasarkan orde polinomial yang digunakan sebagai taksiran yang terdiri atas faktor bobot dan titik evaluasi. Format fungsi tersebut adalah sebagai berikut:

legendre.quadrature.rules(n,normalized=FALSE)

Catatan:

• n: orde tertinggi polinomial yang akan ditampilkan
• normalized: nilai boolean. Jika bernilai TRUE, aturan digunakan untuk polinomial ortogonal.

Untuk memahami penerapan kedua fungsi tersebut berikut disajikan contoh penerapannya:

Integrasi adaptif menyediakan pendekatan yang berbeda untuk memperoleh nilai intergral suatu fungsi. Salah satu prinsip utama dari analisis numerik adalah bahwa kita harus berkomitmen pada semacam analisis manusia terhadap suatu masalah sebelum mencoba menyelesaikannya secara algoritmik. Metode analisis numerik umumnya tidak dapat menyelesaikan semua masalah dengan sangat baik. Jadi pengetahuan terhadap masalah yang hendak diselesaikan dapat memungkinkan kita memilih metode numerik yang lebih baik sesuai dengan masalah. Misalnya, dalam konteks integrasi numerik, diskontinuitas pada titik akhir tidak akan cocok untuk solusi Newton-Cotes yang bersifat tertutup.

Tentu saja, akan lebih baik jika kita bisa memprogram komputer untuk mempelajari sesuatu tentang masalah, daripada kita melakukan itu sendiri. Metode adaptif memberikan pendekatan untuk melakukan hal ini. Metode integrasi adaptif memeriksa integral yang mereka operasikan dan mengubah parameter mereka sendiri untuk meningkatkan kualitas integrasi. Algoritma adaptif yang paling sederhana memberikan pendekatan brute force untuk peningkatan kualitas dengan memeriksa error integrasi. Di sisi lain, jika kita tahu error-nya, secara teoritis kita bisa memperbaiki estimasi. Di situlah batas error pada algoritma Newton-Cotes dapat membantu.

Jika kita dapat menemukan sesuatu tentang error tersebut, kita dapat menggunakan informasi itu untuk memperbaiki estimasi pada proses integrasi. Bayangkan kita sedang mengintegrasikan suatu fungsi, \(f \left(x\right)\) pada batas \(\left[a, b\right]\). Jika kita menggunakan metode titik tengah (integral Riemann). Kita telah menetahui error maksimum yang mungkin terjadi pada metode tersebut melalui Persamaan (9.9). Dua pengamatan segera menjadi jelas. Terlepas dari apa fungsi \(f \left(x\right)\) itu, atau turunan keduanya, dua perubahan dapat dilakukan pada integrasi untuk meningkatkan kualitasnya. Pertama, error adalah proporsi terhadap kubik dari panjang domain integrasi. Mengurangi panjang meningkatkan kualitas dan memotong panjang menjadi dua memberikan peningkatan kualitas delapan kali lipat. Kedua, error berbanding terbalik dengan jumlah panel m yang digunakan . Meningkatkan panel mengurangi error dan menggandakan jumlah panel yang disediakan dapat meeningkatkan kualitas empat kali lipat.

Kita dapat merancang algoritma di sekitar pengamatan ini. Pertama, kita dapat memperkirakan nilai integral \(Q_1\), menggunakan aturan titik tengah 1-point. Kedua, kita bisa memperkirakannya lagi menggunakan aturan titik tengah 2-point \(Q_2\). Karena kita telah menggandakan jumlah titik dalam aturan, kita sekarang tahu bahwa perbedaan maksimum antara \(Q_1\) dan \(Q\), nilai sebenarnya dari integral, tidak lebih dari 4 kali lebih besar dari perbedaan maksimum antara \(Q_2\) dan \(Q\). Jika \(Q_2 −Q\) kurang dari toleransi tertentu, maka perbedaan antara \(Q_1\) dan \(Q_2\) harus kurang dari tiga kali toleransi yang sama.

Kita peeriksa perbedaan antara dua perkiraan. Jika perbedaannya lebih besar dari toleransi, kita mungkin masih berada dalam toleransi, tetapi kita belum pasti. Jadi proses membagi wilayah integrasi menjadi dua, dan menerapkan integrator adaptif untuk kedua bagian, secara terpisah, menjumlahkan hasilnya akan terus dilakukan.

Algoritma Metode Riemann Adaptif

1. Tentukan fungsi \(f\left(x\right)\) dan selang integrasinya \(\left[a,b\right]\).
2. Tentukan jumlah subinterval \(m\).
3. Jika \(m=1\), hitung luas area di bawah kurva dengan pendekatan metode Riemann dengan \(m=2\).
4. Jika \(n>1\),

• Hitung \(Q_1\) dengan pendekatan metode Riemann dan \(m=1\)
• Hitung \(Q_2\) dengan pendekatan metode Riemann dan \(m=2\)

1. Jika \(Q_1-Q_2>3\times\text{nilai toleransi}\),

• Perkecil \(m\) sebanyak 1, \(m-1\)
• Perkecil \(\text{nilai toleransi}\) menjadi setengahnya
• Bagi area integrasi menjadi dua bagian dengan menetapkan \(c\) sebagai batas, sehingga terdapat dua batas yaitu: \(\left[a, c\right]\) dan \(\left[c, b\right]\).
• Lakukan perhitungan kembali integral pada masing-masing batas tersebut menggunakan metode Riemann dan cek apakah \(Q_1-Q_2>3\times\text{nilai toleransi}\).

1. Jika \(Q_1-Q_2<3\times\text{nilai toleransi}\), luas integral = \(Q_2\).

Kita dapat membangun sebuah fungsi integral adaptif menggunakan algoritma tersebut. Sintaks fungsi tersebut adalah sebagai berikut:

riemann_adaptint <- function(f, a, b, m=10, tol=1e-8){ if(m<1){ stop("m harus >= 1") }else if(m==1){ m <- 2 n_width <- (b-a)/m x <- seq(a, b-n_width, length.out = m) + n_width/2 y <- f(x) area <- sum(y)*abs(b-a)/m }else{ m1 <- 1 n_width1 <- (b-a)/m1 x1 <- seq(a, b-n_width1, length.out = m1) + n_width1/2 y1 <- f(x1) q1 <- sum(y1)*abs(b-a)/m1 m2 <- 2 n_width2 <- (b-a)/m2 x2 <- seq(a, b-n_width2, length.out = m2) + n_width2/2 y2 <- f(x2) q2 <- sum(y2)*abs(b-a)/m2 if(abs(q1-q2)>3*tol){ m <- m-1 tol <- tol/2 c <- (a+b)/2 lt <- riemann_adaptint(f, a, c, m=m, tol=tol) rt <- riemann_adaptint(f, c, b, m=m, tol=tol) area <- lt+rt }else{ area <- q2 } } return(area) }

Contoh 9.7 Hitung integral fungsi di bawah ini dengan menggunakan integral Riemann adaptif dengan jumlah panel yang digunakan m=100!

\[ \int_{1}^{10} \sin\left(x\right)^2+\log\left(x\right)dx \]

Jawab:

Kita akan menghitung integral dari fungsi tersebut menggunakan fungsi R yang telah dibuat. Berikut adalah sintaks yang digunakan:

riemann_adaptint(function(x)sin(x)^2 + log(x), a=1, b=10, m=100)

Terdapat dua buah fungsi yang hendak penulis kenalkan pada pembaca yang berfungsi untuk melakukan integrasi adaptif pada R. Fungsi-fungsi tersebut antara lain: integrate() dari Paket base dan integral() dari Paket pracma.

Fungsi integrate() merupakan fungsi yang akan melakukan integrasi numerik menggunakan metode kudratur adaptif untuk sebuah variabel dengan selang terbatas (finite) maupun tidak terbatas (infinite). Format fungsi tersebut secara umum adalah sebagai berikut:

Catatan:

• f: fungsi yang akan dicari integralnya
• lower: batas bawah.
• upper: batas atas.
• …: argumen tambahan functn
• subdivision: jumlah subinterval atau panel yang akan digunakan.
• rel.tol: nilai akurasi relatif yang hendak dicapai
• abs.tol: nilai akurasi absolut yang hendak dicapai

Contoh penerapan fungsi integrate() adalah sebagai berikut:

integrate(function(x)sin(x)^2 + log(x), lower=1, upper=10, rel.tol = 1e-8)

Fungsi lainnya yang dapat digunakan untuk melakukan komputasi integral adaptif adalah fungsi integral() dari Paket pracma. Terdapat dua buah metode integrasi adaptif yang dapat digunakan pada fungsi tersebut yaitu: Gauss-Konrod dan Simpson. Metode Clenshaw-Curtis yang tersedia masih belum dapat melakukan integrasi adaptif melalui fungsi tersebut. Format umum fungsi integral() adalah sebagai berikut:

Catatan:

• fun: fungsi yang akan dicari integralnya
• xmin: batas bawah.
• xmax: batas atas.
• method: metode imtegrasi yang digunakan.
• …: argumen tambahan fun.
• no_intervals: jumlah subinterval atau panel yang akan digunakan.
• reltol: nilai akurasi relatif yang hendak dicapai
• abstol: nilai akurasi absolut yang hendak dicapai

Berikut merupakan contoh penerapan fungsi integral():

Seperti halnya algoritma integrasi adaptif, integrasi Romberg adalah perluasan yang relatif mudah dari keluarga algoritma Newton-Cotes. Keduanya bekerja dengan menggunakan iterasi yang disempurnakan dari beberapa metode Newton-Cotes yang mendasarinya untuk memberikan perkiraan nilai integral yang lebih akurat. Tidak seperti proses komputasi fungsi riemann_adapint(), integrasi Romberg bukanlah pendekatan adaptif terhadap integrasi. Hal tersebut berarti metode Romberg tidak mengubah perilakunya sendiri berdasarkan perilaku fungsi yang akan diintegrasikan. Sebaliknya, kita mengeksploitasi perilaku fungsi trapesium pada batas untuk menghasilkan estimasi integral.

Untuk memahami integrasi Romberg, kita harus mulai dengan implementasi rekursif dari aturan trapesium. Jika kita mulai dengan suatu fungsi, \(T\left(f, m\right)\) di mana \(T\) adalah fungsi trapesium, \(f\) adalah fungsi yang akan diintegrasikan, dan \(m\) adalah jumlah panel untuk diintegrasikan, maka,

\[\begin{equation} S\left(f, m\right)=\frac{4T\left(f, m\right)-T\left(f, m/2\right)}{3} \tag{9.20} \end{equation}\]

di mana \(S\) adalah aturan Simpson. Kemudian, jika kita mendefinisikan \(T \left(f, 0\right) = \left(b − a\right) \left(f \left(b\right) + f \left(a\right)\right) = 2\), maka fungsi rekursif kita selesai, karena berdasarkan hubungan ini, fraksi yang diberikan dalam Persamaan (9.20) juga merupakan perkiraan untuk integral.

Secara umum,

\[\begin{equation} I_{j,k}=\frac{4^k I_{j,k-1}-I_{j-1,k-1}}{4^k-1} \tag{9.21} \end{equation}\]

di mana \(I_{0, 0}\) adalah aturan trapesium satu panel dan \(I_{j, 0}\) adalah aturan trapesium dengan panel \(2^j\). Dengan menggunakan fungsi-fungsi dasar ini, \(I_{j, k}\) dapat ditemukan secara iteratif sebagai matriks segitiga-bawah di mana masing-masing nilai di kolom yang bukan paling kiri adalah fungsi dari nilai di sebelah kiri dan entri di atasnya.

Definisi rekursif ini muncul dari ekstrapolasi Richardson. Ketika diterapkan pada algoritma trapesium, yang konvergen menuju nilai sebenarnya dari integral sebagai \(m\) (jumlah panel) meningkat, hubungan dalam Persamaan (9.21) muncul. Penting untuk dipahami bahwa pada batas ketika \(k\) mendekati tak terhingga, nilai \(I_{j, k}\) adalah nilai sejati integral. Untuk nilai yang lebih kecil dari \(k\), integral Romberg masih hanya perkiraan, meskipun hasil yang diperoleh sangat bagus.

Algoritma Metode Integrasi Romberg

1. Tentukan fungsi \(f\left(x\right)\) dan selang integrasinya \(\left[a,b\right]\).
2. Tentukan jumlah subinterval \(m\).
3. Bentuk matrik \(R\) dengan ukuran \(m\times m\) yang akan menampung hasil perhitungan.
4. Untuk \(R_{1,1}\) hitung integral fungsi menggunakan metode trapezoida dengan \(m=1\).
5. Untuk \(j=2,\dots,m\) dan \(k=1\), hitung integral dengan jumlah panel \(m=2^{j-1}\)
6. Untuk \(j=2,\dots,m\) dan \(k=2,\dots,m\) hitung nilai perbaikan nilai integrasi menggunakan Persamaan (9.21).
7. Solusi integrasi diperoleh pada \(R_{m,m}\).

Berdasarkan algoritma tersebut, kita akan menyusun suatu fungsi pada R untuk melakukan proses komputasi integrasi dengan metode Romberg. Berikut adalah sintaks fungsi yang dibuat:

romberg <- function(f, a, b, m, tab=FALSE){ R <- matrix(NA, nrow=m, ncol=m) R[1,1] <- trap(f, a, b, m=1) for(j in 2:m){ R[j,1] <- trap(f, a, b, m=2^(j-1)) for(k in 2:j){ k4 <- 4^(k-1) R[j,k] <- (k4*R[j,k-1]-R[j-1,k-1])/(k4-1) } } if(tab==TRUE){ return(R) }else{ return(R[m,m]) } }

Contoh 9.8 Hitung integral fungsi yang ditampilkan pada Contoh 9.7 dengan m=10!

Jawab:

Kita dapat menggunakan fungsi romberg() untuk melakukan proses integrasi menggunakan metode Romberg. Berikut adalah sintaks yang digunakan:

Berdasarkan hasil perhitungan nilai integral fungsi tersebut adalah 18.5249.

Fungsi romberg() pada Paket pracma dapat digunakan untuk melakukan integrasi metode Romberg. Format umum fungsi tersebut adalah sebagai berikut:

romberg(f, a, b, maxit = 25, tol = 1e-12, ...)

Catatan:

• f: fungsi yang akan dicari integralnya
• a: batas bawah.
• b: batas atas.
• …: argumen tambahan functn
• maxit: jumlah iterasi maksimum.
• tol: nilai akurasi yang hendak dicapai

Berikut adalah contoh penerapan fungsi romberg():

pracma::romberg(function(x)sin(x)^2 + log(x), a=1, b=10)

Nama Monte Carlo berasal dari daerah di Monako, yang terkenal karena aktikvitas kasino dan perjudiannya. Jelas, permainan kasino yang baik tergantung pada keacakan, seperti juga metode Monte Carlo. Nama ini menggambarkan pentingnya keacakan dalam proses karena algoritma Monte Carlo menggunakan generator angka acak untuk membedakan input ke suatu fungsi.

Angka acak harus berasal dari domain fungsi yang diharapkan. Selanjutnya, fungsi itu sendiri bersifat deterministik karena untuk diberikan dua input dari domain fungsi \(x_1\) dan \(x_2\). Jika \(x_1 = x_2\), maka \(f \left(x_1\right) = f \left(x_2\right)\). Generator angka acak digunakan untuk menghasilkan sejumlah besar input dan fungsinya dijalankan pada setiap input. Akhirnya, hasil yang diperoleh dikumpulkan sesuai dengan model logika yang sesuai dengan analisis yang dilakukan.

Metode Monte Carlo dapat digunakan untuk integrasi numerik dalam jumlah dimensi apa pun yang diberikan. Pendekatan mendasar metode ini adalah menempatkan beberapa titik \(m\) secara acak di atas domain untuk diintegrasikan. Jika titik terletak “di bawah” garis fungsi, maka titik tersebut dianggap dalam area integrasi. Jika titiknya “di atas” garis fungsi, maka titik tersebut bukan berada diluar garis integrasi. Area di bawah perkiraan kurva adalah persentase titik di bawah garis.

Beberapa algoritma Monte Carlo yang paling awal digunakan untuk menemukan area di bawah kurva atau untuk memperkirakan nilai \(\pi\) sebuah hobi favorit matematikawan sejak dahulu. Satu pendekatan menciptakan seperempat lingkaran, menggunakan fungsi \(f \left(x\right) = \sqrt{1-x^2}\). Melalui domain \(\left[0, 1\right]\), ini adalah fungsi dan merupakan hasil dari penyelesaian persamaan standar untuk lingkaran, \(x^2 + y^2 = r\) untuk \(y\) di mana \(r = 1\).

Gambar 9.7 menunjukkan plot fungsi ini. Selain itu, 20 titik acak dipilih. Jika titik di bawah kurva dilambangkan dengan lingkaran hitam terisi dan titik-titik kurva dilambangkan dengan titik bulat kosong. Dalam contoh ini, terdapat 15 titik berada di bawah kurva, mengarah ke estimasi area luas area 0,75. Karena kurva mewakili seperempat lingkaran, estimasi untuk \(\pi\) adalah 3. Meningkatkan jumlah tes titik acak meningkatkan ketepatan estimasi dan akurasi.

Gambar 9.7: Visualisasi metode Monte-Carlo untuk fungsi setengah lingkaran dengan jumlah bilangan acak 20 (sumber: Jones et.al., 2014).

Bentuk umum metode Monte-Carlo disajikan pada Persamaan (9.22).

\[\begin{equation} \int_a^bf\left(x\right)dx\approx\left(b-a\right)\frac{1}{N}\sum_{i=1}^Nf\left(x_i\right) \tag{9.22} \end{equation}\]

dimana \(N\) merupakan jumlah titik yang akan dievaluasi.

Algoritma Metode Monte Carlo

1. Tentukan fungsi \(f\left(x\right)\) dan selang integrasinya \(\left[a,b\right]\).
2. Tentukan jumlah titik acak yang akan digunakan \(N\).
3. Lakukan produksi titik acak \(x\) dengan selang \(\left[a,b\right]\) sejumlah \(N\)
4. Hitung \(f\left(x_i\right)\)
5. Hitung estimasi area menggunakan Persamaan (9.22)

Berdasarkan algoritma tersebut, fungsi R dapat dibangun untuk melakukan integrasi numerik menggunakan metode Monte Carlo. Berikut adalah sintaks yang digunakan:

Contoh 9.9 Hitung integral fungsi yang ditampilkan pada Contoh 9.7 menggunakan metode Monte-Carlo dengan m=1e6!

Jawab:

Integrasi Monte Carlo menggunakan menggunakan fungsi monte_int() disajikan pada sintaks berikut:

monte_int(function(x)sin(x)^2 + log(x), a=1, b=10)

Hasil yang diperoleh sedikit berbeda dengan yang dihasilkan oleh metode lainnya. Hal ini disebabkan oleh penggunaan bilangan acak pada proses integrasi. Selain itu, metode ini juga menghasilkan kualitas hasil yang rendah dengan tingkat komputasi yang tinggi dibandingkan dengan metode Newton-Cotes.

Keunggulan metode Monte Carlo dibandingkan metode sebelumnya adalah kemampuan untuk menangani proses integrasi berganda. Berikut adalah bentuk umum proses integrasi bivariat menggunakan metode Monte Carlo:

\[\begin{equation} \int \int f\left(x,y\right)dx\approx V\frac{1}{N}\sum_{i=1}^Nf\left(x_i,y_i\right) \tag{9.23} \end{equation}\]

dimana \(V\) merupakan area perpotongan \(x-y\) dimana fungsi \(f\left(x,y\right)\) diintegrasikan.

Algoritma Metode Monte Carlo Bivariat

1. Tentukan fungsi \(f\left(x\right)\) dan domain integrasinya pada masing-masing sumbu \(x\) dan \(y\)
2. Tentukan jumlah titik acak yang akan digunakan \(N\).
3. Lakukan produksi titik acak \(x\) dan \(y\) masing-masing domain sumbunya sejumlah \(N\)
4. Hitung \(V=\left( x_{max}-x_{min}\right) \times \left(y_{max}-y_{max}\right)\)
5. Hitung estimasi volume menggunakan Persamaan (9.23)

Contoh 9.10 Hitung volume melalui integrasi persamaan berikut menggunakan metode Monte Carlo dengan domain x = [0,1] dan y=[0,1]!

\[ \int_{0}^1 \int_{0}^1 x^2y \ dy \ dx \]

Jawab:

Berikut adalah sintaks yang digunakan untuk melakukan integrasi persamaan tersebut menggunakan metode Monte Carlo bivariat:

monte_int2(function(x,y)x^2*y, xdom=c(0,1), ydom=c(0,1))

Karena metode Monte Carlo tidak deterministik, error integrasi Monte Carlo tidak dibatasi dalam pengertian yang telah kita lihat sejauh ini. Namun, kita dapat memperkirakan varians dari estimasi yang dihasilkan, yang berkurang dengan meningkatnya jumlah poin:

\[\begin{equation} \text{Var}\frac{1}{N}\sum f\left(\right)=\frac{\sigma^2}{N} \tag{9.24} \end{equation}\]

dimana \(\sigma^2=\text{Var}\ f\left(\right)\). Definisi ini juga digunakan untuk proses integrasi dengan dimensi yang lebih tinggi.

Perlu diketahui pula bahwa metode Monte Carlo hanya dapat digunakan jika nilai terendah dari variabel bebas yang digunakan tidak negatif. Hal ini dilakukan untuk mencegah adanya pengurangan dengan nilai negatif sehingga hasil integrasi jauh lebih besar dari yang seharusnya.

Pada studi kasus kali ini, penulis akan memberikan contoh penerapan integrasi numerik dalam menganalisa jarak jatuh seorang penerjung payung yang melompat dari sebuah pesawat. Kecepatan penerjun payung dapat dituliskan ke dalam sebuah fungsi dari waktu,

\[\begin{equation} v\left(t\right)=\frac{gm}{c}\left(1-e^{-\left(c/m\right)t}\right) \tag{9.25} \end{equation}\]

dimana \(v\) adalah kecepatan penerjun payung dalam \(m/dt\), \(g\) adalah percepatan gravitasi sebesar \(9,8 \ m/dt^2\), \(m\) adalah massa penerjun payung sebesar \(68,1 \ kg\), dan \(c\) adalah koefisien tahanan udara sebesar \(12,5\ kg/dt\).

Misalkan kita ingin mengetahui seberapa jauh penerjun telah jatuh setelah waktu tertentu \(t\). Karena kecepatan merupakan turunan pertama dari fungsi jarak, maka jarak penerjun dari titik terjun \(\left(t=0\right)\) adalah:

\[\begin{equation} d=\int_{0}^{t}v\left(t\right)dt=\int_{0}^{t}\frac{gm}{c}\left(1-e^{-\left(c/m\right)t}\right)dt \tag{9.26} \end{equation}\]

Jika kita ingin mengetahui jarak yang telah ditempuh saat \(t=10\), kita dapat melakukan integrasi pada persamaan tersebut dengan domain \(t=\left[0,10\right]\). Persamaan (9.26) dapat dinyatakan menjadi Persamaan (9.27) dengan memasukkan semua komponen yang telah diketahui sebelumnya.

\[\begin{equation} d=\int_{0}^{10}\frac{9,8\times 68,1}{12,5}\left(1-e^{-\left(12,5/68,1\right)^t}\right)dt \tag{9.27} \end{equation}\]

Kita dapat menyelesaikan Persamaan (9.27) dengan menggunakan metode-metode integrasi numerik yang telah dijabarkan sebelumnya. Berikut adalah sintaks untuk masing-masing metode tersebut:

f <- function(x)((9.8*68.1)/12.5)*(1-exp(-(12.5/68.1)*x)) a <- 0; b <- 10

Newton-Cotes

# Metode Riemann riemann(f, a, b)

# Metode Trapezoida trap(f, a, b)

# Metode Simpson 1/3 simpson(f, a, b)

# Metode Simpson 3/8 simpson38(f, a, b)

Metode Adaptif

riemann_adaptint(f, a, b, m=100)

Metode Romberg

Metode Monte Carlo

1. Bloomfield, V.A. 2014. Using R for Numerical Analysis in Science and Engineering. CRC Press.
2. Chapra, S.C. Canale, R.P. 2015. Numerical Methods For Engineers, Seventh Edition. Mc Graw Hill.
3. Howard, J.P. 2017. Computational Methods for Numerical Analysis with R. CRC Press.
4. Kreyszig, E. 2011. Advanced Engineering Mathematics, 10th Edition. John Wiley & Sons.
5. Sanjaya, M. 2015. Metode Numerik Berbasis Phython. Penerbit Gava Media: Yogyakarta.
6. Suparno, S. 2008. Komputasi untuk Sains dan Teknik Edisi II. Departemen Fisika-FMIPA Universitas Indonesia.

1. Hitung integral fungsi \(f\left(x\right)=\sin^2\left(x\right)\) pada domain \(x \in \left[0,\pi\right]\) !
2. Tuliskan fungsi R yang dapat melakukan integrasi Riemann dengan aturan titik kiri !
3. Buatlah sebuah fungsi R yang dapat melakukan integrasi adaptif menggunakan metode Simpson 1/3 !
4. Kerjakan kembali soal 3 dengan menggunakan metode Simpson 3/8! (Note: pembaca dapat melakukan pencarian algoritma di internet dan mentrasformasikannya menjadi sintaks R)
5. Fungsi monte_int() hanya mampu melakukan integrasi pada domain positif. Buatlah algoritma baru sehingga metode ini dapat melakukan integrasi pada domain negatif !

Persamaan diferensial merupakan persoalan matematis yang sering dijumpai dalam bidang teknik lingkungan. Sering kali suatu persamaan diferensial tidak dapat diselesaikan secara analitik sehingga diperlukan metode numerik untuk menyelesaikannya. Pada Chapter 10, kita akan membahas masalah-masalah dalam persamaan diferensial dan metode penyelesaiannya. Adapun yang akan dibahas pada Chapter 10 kali ini antara lain:

• Initial value problems
• Sistem persamaan diferensial
• Persamaan diferensial parsial

Initial value problems merupakan permasalahan yang sering ditemukan pada proses dekomposisi zat kimia atau polutan dalam reaktor. Penyelesaiaan persamaan diferensial biasanya dipersulit dengan tidak tersedianya informasi yang cukup untuk menyelesaikannya.Sebuah persamaan diferensial \(f'\left(x,\dots\right)\) merupakan hasil diferensiasi beberapa fungsi \(f\left(x,\dots\right)\). Proses penyelesaian persamaan diferensial, dan menemukan nilai \(f \left(x,\dots\right)\) untuk beberapa nilai x,… tidak dimungkinkan karena integral dari \(f'\left(x,\dots\right)\) hanya digambarkan bentuk umum. Pergeseran vertikal, atas atau bawah, tidak diketahui. Pergeseran vertikal ini menghasilkan konstanta integrasi.

Selama proses diferensiasi, nilai apapun dari proses pergeseran vertikal (integrasi) akan hilang sebagai akibat dari eliminasi konstanta yang memiliki turunan 0. Kita biasa melakukannya ketika mengintegrasikan fungsi dengan menambahkan konstanta \(+ C\) pada proses integrasi ke integral yang tidak terbatas. Hal ini terkadang bukan menjadi permasalahan sebab jika menemukan nilai integrasi pada suatu batas tertentu syarat \(+ C\) dibatalkan dan konstanta integrasi tidak diperlukan.

Untuk persamaan diferensial biasa, tidak ada pembatalan yang nyaman, yang mengarah ke initial value problems. Initial value problems memberikan nilai \(f\left(x_0,\dots\right)\), di mana \(x_0\) biasanya bernilai 0, meski tidak diharuskan. Nilai awal ini memberikan informasi yang cukup untuk menyelesaikan persamaan dan menemukan nilai aktual dari \(f\left(x,\dots\right)\) untuk sejumlah nilai \(x\). Terdapat beberapa metode yang akan dibahas pada Chapter 10.1, antara lain:

• Metode Euler
• Metode Heun
• Metode Titik Tengah
• Metode Runge-Kutta Orde 4
• Metode multistep linier

Metode Euler merupakan metode paling sederhana yang diturunkan dari deret Taylor. Penyelesaian initial value problems menggunakan metode Euler dilakukan melalui Persamaan (10.1).

\[\begin{equation} y_{i+1}=y_i+f\left(x_i,y_i\right)h \tag{10.1} \end{equation}\]

dimana \(i\) merupakan tahapan iterasi.

Algoritma Metode Euler

1. Tentukan titik awal integrasi \(x_0\) dan \(y_0\).
2. Tentukan jumlah iterasi \(n\) dan step size \(h\) yang digunakan.
3. Lakukan integrasi menggunakan Persamaan (10.1).

Algoritma tersebut, selanjutnya dapat disusun ke dalam sebuah fungsi R. Fungsi tersebut adalah sebagai berikut:

Contoh 10.1 Selesaikan persamaan diferensial di bawah ini, jika diketahui f(0)=1 menggunakan h=0,05 dan n=100!

\[ f'\left(x,y\right)=\frac{y}{2x+1} \]

Jawab:

Penyelesaian secara analitik persamaan tersebut untuk nilai \(f\left(0\right)=1\) sebagai berikut:

\[ f\left(x\right)=\sqrt{2x+1} \]

Secara numerik persamaan tersebut dapat diselesaikan sebagai berikut:

iterasi 1

\[ y\left(0,5\right)=1+0,05\times\frac{1}{\left(2\cdot 0\right)+1}=1,05 \]

iterasi 2

\[ y\left(0,1\right)=1,05+0,05\times\frac{1}{\left(2\cdot0,05\right)+1}=1.097727 \]

Kita dapat juga menggunakan fungsi euler() untuk menyelesaikan persamaan tersebut secara numerik. Hasil yang diperoleh selanjutnya diplotkan dengan hasil yang diperoleh menggunakan metode analitik. Berikut adalah sintaks yang digunakan:

Gambar 10.1: Visualisasi integrasi numerik dengan metode Euler dan metode analitik

Berdasarkan hasil visualisasi dapat dilihat bahwa metode Euler dapat dengan baik memberikan pendekatan nilai integrasi persamaan. Pembaca dapat mencoba untuk melakukan simulasi kembali dengan nilai \(h\) yang lebih kecil.

Metode Heun merupakan salah satu peningkatan dari metode Euler. Metode ini melibatkan 2 buah persamaan. Persamaan pertama disebut sebagai persamaan prediktor yang digunakan untuk memprediksi nilai integrasi awal (Persamaan (10.2)). Persamaan kedua disebut sebagai persamaan korektor yang mengoreksi hasil integrasi awal (Persamaan (10.3)). Metode Heun pada Chapter ini merupakan metode prediktor-korektor satu tahapan. Akurasi integrasi dapat ditingkatkan dengan melakukan koreksi ulang terhadap nilai koreksi semula menggunakan persamaan kedua.

\[\begin{equation} y_{i+1}^{0}=y_i+f\left(x_i,y_i\right)h \tag{10.2} \end{equation}\]

\[\begin{equation} y_{i+1}=y_i+\frac{f\left(x_i,y_i\right)+f\left(x_{i+1},y_{i+1}^0\right)}{2}h \tag{10.3} \end{equation}\]

Algoritma Metode Heun

1. Tentukan titik awal integrasi \(x_0\) dan \(y_0\).
2. Tentukan jumlah iterasi \(n\) dan step size \(h\) yang digunakan.
3. Lakukan prediksi nilai awal dengan Persamaan (10.2).
4. Lakukan koreksi nilai awal menggunakan Persamaan (10.3).
5. Lakukan koreksi terhadap nilai koreksi yang dihasilkan sebelumnya menggunakan Persamaan (10.3).

Kita dapat membangun sebuah fungsi yang dapat melakukan proses integrasi menggunakan metode Heun. Berikut adalah sintaks yang digunakan:

heun <- function(f, x0, y0, h, n, iter=1){ x <- x0 y <- y0 for(i in 1:n){ ypred0 <- f(x0,y0) ypred1 <- y0 + h*ypred0 ypred2 <- f(x0+h,ypred1) ykor <- y0 + h*(ypred0+ypred2)/2 if(iter!=1){ for(i in 1:iter){ ykor <- y0 + h*(ypred0+f(x0+h,ykor))/2 } } y0 <- ykor x0 <- x0 + h x <- c(x, x0) y <- c(y, y0) } return(data.frame(x=x,y=y)) }

Contoh 10.2 Selesaikan kembali persamaan yang ditampilkan pada Contoh 10.1 menggunakan metode Heun!

Jawab:

Contoh perhitungan secara manual menggunakan metode Heun untuk sekali iterasi adalah sebagai berikut:

\[ f'\left(0;1\right)=\frac{1}{\left(2\cdot 0\right)+1}=1 \]

\[ y_{1}^0=1+0,05\cdot1=1,05 \]

\[ y'_1=f'\left(0,05;1,05\right)=\frac{1,05}{\left(2\cdot 0,05\right)+1}=0,9545455 \]

\[ y_1=1+0,05\cdot\frac{1+0,9545455}{2}=1,047727 \]

Penyelesaian persamaan tersebut menggunakan fungsi heun() dengan iterasi pada nilai koreksi sebanyak 1 kali disajikan pada sintaks berikut:

num <- heun(f1, x0=0, y0=1, h=0.05, n=100)

Gambar 10.2: Visualisasi integrasi numerik dengan metode Heun dan metode analitik

Metode titik tengah menggunakan setengah step size pada metode Euler untuk melakukan estimasi terhadap integral suatu persamaan diferensial. Metode ini melakukan perhitungan melalui dua tahapan yaitu: menghitung nilai estimasi integral pada setengah step size(Persamaan (10.4)) dan menghitung nilai integral menggunkan hasil perhitungan setengah step size sebelumnya (Persamaan (10.5)).

\[\begin{equation} y_{i+\frac{1}{2}}=y_i+f\left(x_i,y_i\right)\frac{h}{2} \tag{10.4} \end{equation}\]

\[\begin{equation} y_{i+1}=y_i+f\left(x_{i+\frac{1}{2}},y_{i,\frac{1}{2}}\right)h \tag{10.5} \end{equation}\]

Algoritma Metode Tengah

1. Tentukan titik awal integrasi \(x_0\) dan \(y_0\).
2. Tentukan jumlah iterasi \(n\) dan step size \(h\) yang digunakan.
3. Lakukan integrasi pada setengah tahapan iterasi menggunakan Persamaan (10.4).
4. Lakukan iterasi pada setengah tahapan selanjutnya menggunakan Persamaan (10.5).

Berdasarkan algoritma tersebut, kita dapat membangun sebuah fungsi pada R yang adapat digunakan untuk menyelesaikan persamaan diferensial menggunakan metode titik tengah. Berikut adalah sintaks yang digunakan:

Contoh 10.3 Selesaikan kembali persamaan yang ditampilkan pada Contoh 10.1 menggunakan metode titik tengah!

Jawab:

Contoh perhitungan secara manual menggunakan metode titik tengah untuk sekali iterasi adalah sebagai berikut:

\[ y_{\frac{1}{2}}=1+\frac{1}{\left(2\cdot 0\right)+1}\cdot\frac{0,05}{2}=1,025 \]

\[ y_{1}=1+\frac{1,025}{\left(2\cdot 0,025\right)+1}\cdot0,05=1,0488 \]

Penyelesaian persamaan tersebut menggunakan fungsi midpt() dengan iterasi pada nilai koreksi sebanyak 1 kali disajikan pada sintaks berikut:

num <- midpt(f1, x0=0, y0=1, h=0.05, n=100)

Gambar 10.3: Visualisasi integrasi numerik dengan metode titik tengah dan metode analitik

Runge-Kutta orde 4 merupakan metode yang paling populer dalam penyelesaian persamaan diferensial. Metode ini dapat memperoleh akurasi deret Taylor tanpa memerlukan diferensiasi orde yang lebih tinggi. Metode Runge-Kutta orde 4 dituliskan ke dalam Persamaan (10.6).

\[\begin{equation} y_{i+1}=y_i+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)h \tag{10.6} \end{equation}\]

dimana

\[\begin{equation} k_1=f\left(x_i,y_i\right) \tag{10.7} \end{equation}\]

\[\begin{equation} k_2=f\left(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_1h\right) \tag{10.8} \end{equation}\]

\[\begin{equation} k_3=f\left(x_i+\frac{1}{2}h,y_i+\frac{1}{2}k_2h\right) \tag{10.9} \end{equation}\]

\[\begin{equation} k_4=f\left(x_i+h,y_i+k_3h\right) \tag{10.10} \end{equation}\]

Algoritma Metode Tengah

1. Tentukan titik awal integrasi \(x_0\) dan \(y_0\).
2. Tentukan jumlah iterasi \(n\) dan step size \(h\) yang digunakan.
3. Lakukan integrasi menggunakan Persamaan (10.6).

Berdasarkan algoritma tersebut, kita dapat membangun sebuah fungsi pada R yang adapat digunakan untuk menyelesaikan persamaan diferensial menggunakan metode Runge-Kutta orde 4. Berikut adalah sintaks yang digunakan:

rk4 <- function(f, x0, y0, h, n){ x <- x0 y <- y0 for(i in 1:n){ k1 <- f(x0,y0) k2 <- f(x0+0.5*h,y0+0.5*k1*h) k3 <- f(x0+0.5*h,y0+0.5*k2*h) k4 <- f(x0+h,y0+k3*h) y0 <- y0 + (1/6)*(k1+2*k2+2*k3+k4)*h x0 <- x0 + h x <- c(x, x0) y <- c(y, y0) } return(data.frame(x=x,y=y)) }

Contoh 10.4 Selesaikan kembali persamaan yang ditampilkan pada Contoh 10.1 menggunakan metode Runge-Kutta orde 4!

Jawab:

Contoh perhitungan secara manual menggunakan metode titik tengah untuk sekali iterasi adalah sebagai berikut:

\[ k_1=\frac{1}{\left(2\cdot 0\right)+1}=1 \]

\[ k_2=\frac{1+\frac{1}{2}\cdot1\cdot0,05}{\left(2\cdot 0\cdot\frac{1}{2}\cdot0,05\right)+1}=1,025 \]

\[ k_2=\frac{1+\frac{1}{2}\cdot1,025\cdot0,05}{\left(2\cdot 0\cdot\frac{1}{2}\cdot0,05\right)+1}=1,025625 \]

\[ k_2=\frac{1+1,025625\cdot0,05}{\left(2\cdot 0\cdot0,05\right)+1}=1,051281 \]

\[ y_1=1+\frac{1}{6}\left(1+2\cdot1,025+2\cdot1,025625+1,051281\right)0,05=1.051271 \]

1+(1/6)*(1+2*1.025+2*1.025625+1.051281)*0.05

Iterasi dapat pula dilakukan dengan menggunakan fungsi rk4(). Berikut adalah sintaks yang digunakan:

num <- rk4(f1, x0=0, y0=1, h=0.05, n=100)

Gambar 10.4: Visualisasi integrasi numerik dengan metode Runge-Kutta orde 4 dan metode analitik

Jika metode Runge-Kutta mengalami kesulitan karena terlalu banyak evaluasi fungsi yang digunakan, masuk akal untuk bertanya apakah kita dapat menggunakan kembali beberapa evaluasi fungsi sebelumnya, yang sudah kita buat. Sebagai contoh, kita ingin tahu apakah kita dapat menggunakan kembali estimasi \(f\left(0,1\right)\) dan \(f\left(0,2\right)\) untuk memperkirakan nilai \(f\left(0,3\right)\). Jika kita dapat menggunakan kembali perkiraan sebelumnya, kita dapat memperoleh akurasi tambahan tanpa menimbulkan penalti kinerja yang terkait dengan evaluasi fungsi tambahan. Metode multistep linier dikembangkan untuk mengatasi masalah ini.

Di satu sisi, metode multistep linier dasar untuk persamaan diferensial hanya mencakup satu titik \(x_i\), dalam perhitungan \(x_i + 1\). Ini persis bagaimana fungsi metode Euler dan metode Euler merupakan metode multistep linier dasar. Metode selanjutnya menggunakan \(x_i − 1\) dan \(x_i\) untuk menghitung \(x_i + 1\). Metode Adams-Bashforth menggunakan tambahan berbobot, termasuk bobot negatif, dari langkah dan poin untuk sampai pada langkah berikutnya. Seperti metode numerik lainnya, bobot muncul dari interpolasi polinomial titik yang tersedia.

Metode Adam-Bashforth orde 2 didasarkan pada Persamaan (10.11).

\[\begin{equation} y_{i+2}=y_{i+1}+\frac{h}{2}\left(3f\left(x_{i+1},y_{i+1}\right)-f\left(x_i,y_i\right)\right) \tag{10.11} \end{equation}\]

Pendekatan ini melakukan interpolasi antara titik sebelumnya untuk memperkirakan titik ketiga dalam grup. Titik ketiga ini menjadi titik tengah dari iterasi berikutnya saat seluruh proses berlanjut. Karena nilai sebelumnya disimpan dan digunakan kembali, evaluasi fungsi tambahan tidak diperlukan. Jika metode Runge-Kutta dapat dibandingkan dengan tip-toeing melalui bidang vektor, maka metode Adams-Bashforth dapat sama dibandingkan dengan menjalankan melalui bidang vektor. Namun, ini tidak berarti metode Adams-Bashforth lebih unggul.

Algoritma Metode Tengah

1. Tentukan titik awal integrasi \(x_0\) dan \(y_0\).
2. Tentukan jumlah iterasi \(n\) dan step size \(h\) yang digunakan.
3. Lakukan pendekatan pada iterasi ke-1 menggunakan metode Euler.
4. Lakukan integrasi ke-2 sampai n menggunakan Persamaan (10.11).

Berdasarkan algoritma tersebut, kita dapat membangun sebuah fungsi pada R yang adapat digunakan untuk menyelesaikan persamaan diferensial menggunakan metode multistep linier. Berikut adalah sintaks yang digunakan:

adambashforth <- function(f, x0, y0, h, n){ # pendekatan Euler untuk x1 dan y1 y1 <- y0 + h*f(x0,y0) x1 <- x0 + h x <- c(x0,x1) y <- c(y0,y1) n <- n-1 for(i in 1:n){ yn <- y1 + 1.5*h*f(x1,y1) - 0.5*h*f(x0,y0) xn <- x1 + h y0 <- y1 x0 <- x1 y1 <- yn x1 <- xn y <- c(y,y1) x <- c(x,x1) } return(data.frame(x=x,y=y)) }

Contoh 10.5 Selesaikan kembali persamaan yang ditampilkan pada Contoh 10.1 menggunakan metode Runge-Kutta orde 4!

Jawab:

Untuk melakukan iterasi, kita dapat menggunakan fungsi adambashforth(). Berikut adalah sintaks yang digunakan:

num <- rk4(f1, x0=0, y0=1, h=0.05, n=100)

Gambar 10.5: Visualisasi integrasi numerik dengan metode Adam-Bashfoth orde 2 dan metode analitik

Sistem persamaan diferensial akan sering pembaca temui dalam pemodelan sistem dinamik. Pada proses pemodelan sistem tersebut akan ditemukan aksi-interasi antar komponennya yang dinyatakan ke dalam suatu sistem persamaan diferensial.

Pada fungsi rk4sys() disusun sebuah fungsi pada R untuk melakukan iterasi terhadap sistem persamaan diferensial menggunakan metode Runge-Kutta orde 4. Berikut adalah sintaks yang digunakan:

rk4sys <- function(f, x0, y0, h, n){ x <- x0 y <- y0 values <- data.frame(x=x,t(y0)) for(i in 1:n){ k1 <- f(x0,y0) k2 <- f(x0+0.5*h,y0+0.5*k1*h) k3 <- f(x0+0.5*h,y0+0.5*k2*h) k4 <- f(x0+h,y0+k3*h) y0 <- y0 + (1/6)*(k1+2*k2+2*k3+k4)*h x0 <- x0 + h values <- rbind(values, data.frame(x=x0,t(y0))) } return(values) }

Contoh 10.6 Sebuah model populasi yang dibagi menjadi tiga kelompok umur: anak (0-12 tahun), melahirkan anak (13-40 tahun), dan usia (41 tahun atau lebih). Kelompok 0-12 meningkat berdasarkan kelahiran di kelompok 13-40, dan menurun dengan kematian dan lewat bagian ke kelompok 13-40. Kelompok 13–40 meningkat dengan perolehan dari kelompok 0-12, dan berkurang dengan kematian dan lewat bagian ke dalam kelompok> = 41. Kelompok> = 41 meningkat dengan keuntungan dari kelompok 13-40, dan menurun dengan kematian. Parameter dipilih untuk mewakili tingkat kelahiran dan kematian yang cukup tinggi seperti yang ditemukan di banyak masyarakat berkembang. Model interaksi tersebut dinyatakan ke dalam persamaan interkasi di bawah ini! Simulasikan dinamika populasi pada model tersebut pada 10 tahun ke depan, jika diketahui populasi semula masing-masing kelompok usia secara berurutan adalah 200,400, dan 400, serta nilai masing-masing koefisien kelahiran, kematian kelompok 1, kematian kelompok2, dan kematian kelompok 3 secara berurutan adalah 0,5;0,1;0,1;0,25!

\[ pop1' = b\cdot pop2'+\frac{11}{12}\cdot pop1'\left(1-d1\right) \]

\[ pop2' = \frac{1}{12}\cdot pop1'\cdot\left(1-d1\right)+\frac{26}{27}\cdot pop2'\cdot \left(1-d2\right) \]

\[ pop3, = \frac{1}{27}\cdot pop2,\cdot\left(1-d2\right)+ pop3'\cdot \left(1-d3\right) \]

Jawab:

Sistem persamaan diferensial perlu ditransformasi ke dalam bentuk sebuah fungsi pada R.

Nilai awal dituliskan seperti berikut:

# Nilai awal y = c(pop1=200,pop2=400,pop3=400)

Simulasi model ditampilkan pada sintaks berikut:

pop <- rk4sys(Population, x0=0, y0=y, h=0.1, n=100) head(pop)

Gambar 10.6: Visualisasi hasil simulasi model dinamika populasi

Fungsi ode() pada paket deSolve merupakan salah satu fungsi yang dapat digunakan untuk menyelesaikan persamaan dan sistem persamaan diferensial. Fungsi ini mudah digunakan serta menyediakan berbagai macam metode iterasi numerik untuk menyelesaikan persamaan diferensial. Format umum fungsi ini antara lain:

ode(y, times, func, parms,method , ...)

Catatan:

• y : nilai awal (kondisi awal) suatu persamaan diferensial.
• times : deret waktu yang terdiri dari kapan simulasi dimulai, kapan simulasi berkahir dan berapa step size yang digunakan.
• func : fungsi yang berisi persamaan diferensial. Return value pada fungsi haruslah berupa list.
• parms : list parameter yang diinputan kedalam func.

• method : sebuah string yang berupa metode intgrasi yang digunakan. Metode yang tersedia dan kegunaan metode tersebut antara lain:

  • “bdf”: menangani persamaan diferensial menggunakan formula diferensiasi mundur (bacward differensiation) dan cocok untuk menangani kondisi stiff. Metode ini setara dengan method="lsode".
  • “bdf_d”: menggunakan formula diferensiasi mundur yang memanfaatkan iterasi Jacobi-Newton (mengabaikan elemen diagonal Jacobian). Cocok digunakan persamaan atau sistem persamaan diferensial kondisi stiff. Metode ini setara dengan method="lsode",mf=23.
  • “adams”: metode Adams yang menggunakan iterasi fungsional (tanpa menggunakan Jacobian). Cocok digunakan untuk persamaan atau sistem persamaan non stiff. Metode ini setara dengan method="lsode", mf=10.
  • “impAdams”: metode Adams implisit yang menggunakan iterasi Newton-Raphson. Metode ini setara dengan method="lsode",mf=12.
  • “impAdams_d”: metode Adams implisit yang menggunakan iterasi Jacobi-Newton. Metode ini setara dengan method="lsode, mf=13".
  • “euler”: metode iterasi Euler
  • “rk4”: metode iterasi Runge-Kutta orde 4.
  • metode lain: “lsoda”, “lsode”, “lsodes”, “lsodar”, “vode”, “daspk”, “ode23”, “ode45”, “radau”.

• ...: argumen tambahan untuk integrator atau method.

Contoh penerapan fungsi tersebut menggunakan Contoh 10.6, sebagai berikut:

library(deSolve) # sistem persamaan diferensial Population = function(t,y,param) { y1 = y[1] # 0-12 group population y2 = y[2] # 13-40 group population y3 = y[3] # 41 and older population y1.new = b*y2 + 11/12*y1*(1 - d1) y2.new = 1/12*y1*(1 - d1) + 26/27*y2*(1 - d2) y3.new = 1/27*y2*(1 - d2) + y3*(1 - d3) return(list(c(y1.new, y2.new, y3.new))) } # parameter b = 0.5 # Birth rate in 13-40 group d1 = 0.1 # Death rate of 0-12 group d2 = 0.1 # Death rate of 13-40 group d3 = 0.25 # Death rate of 41 and older # nilai awal y = c(pop1=200, pop2=400, pop3=400) # waktu Tahun = seq(0,10,0.1) # simulasi out = ode(func = Population, y = y, times = Tahun, parms = c(b,d1,d2,d3), method = "rk4")

Gambar 10.7: Visualisasi hasil simulasi model dinamika populasi menggunakan paket desolve

Persamaan diferensial parsial (PDE) banyak dijumpai pada pemodelan transport polutan dalam bidang teknik lingkungan. Persamaan diferensial parsial merupakan persamaan diferensial yang melibatkan lebih dari satu variabel independen, biasanya variabel waktu dan satu atau lebih variabel posisi atau beberapa variabel spasial. PDE diklasifikasikan menjadi 3 jenis: parabolik (time-dependent dan difusif), hiperbolik (time-dependent dan gelombang), dan eliptik (time-independent). Dalam penyelesaian PDE pada umumnya kita menggunakan metode FTCS (forward in time, centered in space). Untuk memahami definisi tersebut, pembaca dapat membaca kembali Chapter 9.1.

Persamaan difusi merupakan contoh PDE parabolik. Persamaan difusi satu dimensi spasial ditampilkan pada Persamaan (10.12).

\[\begin{equation} \frac{\partial C}{\partial t}=D\frac{\partial^2C}{\partial x^2} \tag{10.12} \end{equation}\]

Persamaan tersebut mirip dengan persamaan konduksi panas. Pada persamaan konduksi panas, variabel konsentrasi \(C\) diganti dengan variabel temperatur \(T\), dan koefisien difusi D diganti dengan koefsien difusi termal K.

Untuk menyelesaikan Persamaan (10.12), persamaan tersebut diubah ke dalam bentuk metode beda hingga, dimana turunan waktu menggunakan pendekatan Euler (metode beda hingga maju) dan turunan variabel spasial diubah ke dalam bentuk pendekatan titik pusat. Proses diskretisasi Persamaan (10.12), ditampilkan pada Persamaan (10.13).

\[\begin{equation} \frac{C\left(i+1,j\right)-C\left(i,j\right)}{\Delta t}=D\frac{C\left(i,j+1\right)+C\left(i,j-1\right)-2C\left(i,j\right)}{\Delta x^2} \tag{10.13} \end{equation}\]

dimana \(i\) merupakan step untuk variabel waktu \(t\) dan \(j\) merupakan step untuk variabel spasial \(x\).

Persamaan (10.13) dapat disusun kembali sehingga menjadi Persamaan (10.14) yang menyatakan persamaan konsentrasi \(C\) pada saat \(i+1\) pada posisi \(j\).

\[\begin{equation} C\left(i+1,j\right)=C\left(i,j\right)+A\left[C\left(i,j+1\right)+C\left(i,j-1\right)-2C\left(i,j\right)\right] \tag{10.14} \end{equation}\]

dimana

\[\begin{equation} A=D\frac{\Delta t}{\Delta x^2} \tag{10.15} \end{equation}\]

Untuk stabilitas komputasi, pemilihan peningkatan waktu terhadap jarak yang dinyatakan pada nilai \(A\) harus \(\le\frac{1}{2}\).

Pada contoh berikut, kita akan melakukan simulasi menggunakan Persamaan (10.14). Parameter yang digunakan dan nilai awal yang digunakan dinyatakan pada sintaks berikut:

Langkah selanjutnya adalah inisiasi konsentrasi awal, dimana seluruh konsentrasi awal pada tiap grid adalah nol kecuali pada grid pusat.

Simulasi selanjutnya dilakukan dengan melakukan iterasi pada grid points konsentrasi yang telah dibuat.

Visualisasi dari hasil simulasi tersebut ditampilkan pada Gambar 10.8.

Gambar 10.8: Visualisasi 3D simulasi difusi partikulat

Persamaan gelombang 1 dimensi ditampilkan pada Persamaan (10.16).

\[\begin{equation} \frac{\partial^2W}{\partial t^2}=c^2\frac{\partial^2W}{\partial x^2} \tag{10.16} \end{equation}\]

dimana \(W\) merupakan pemindahan dan \(c\) merupakan kecepatan gelombang. Persamaan (10.16) merupakan bentuk PDE hiperbolik. Versi lebih sederhana dari Persamaan (10.16) adalah persamaan adveksi yang ditampilkan pada Persamaan (10.17).

\[\begin{equation} \frac{\partial y}{\partial t}=-c\frac{\partial y}{\partial x} \tag{10.17} \end{equation}\]

Persamaan (10.17) menggambarkan evolusi pada bidang skalar \(y\left(x,y\right)\) dibawa oleh gelombang dengan kecepatan konstan \(c\) dan bergerak dari kiri ke kanan jika \(c>0\). Persamaan adveksi merupakan contoh palin sederhana dari persamaan konservasi flux.

Diskretisasi Persamaan (10.17) ditampilkan pada Persamaan (10.18).

\[\begin{equation} y\left(i+1,j\right)=y\left(i,j\right)-\frac{c\Delta t}{2\Delta x}\left[y\left(i,j+1\right)-y\left(i,j-1\right)\right] \tag{10.18} \end{equation}\]

Pada contoh berikut, kita akan melakukan simulasi menggunakan Persamaan (10.17). Parameter yang digunakan dan nilai awal yang digunakan dinyatakan pada sintaks berikut:

Kondisi awal dan kondisi batas periodik dinyatakan sebagai berikut:

Langkah selanjutnya melakukan interasi untuk melihat perubahan nilai pada grid points.

Visualisasi dari hasil simulasi tersebut ditampilkan pada Gambar 10.9.

Gambar 10.9: Visualisasi simulasi adveksi

Persamaan Laplace dalam dua dimensi disajikan pada Persamaan (10.19).

\[\begin{equation} \frac{\partial V}{\partial x^2}+\frac{\partial V}{\partial y^2}=0 \tag{10.19} \end{equation}\]

Persamaan (10.19) merupakan contoh tipe ketiga PDE, persamaan elips. Persoalan ini sering muncul dalam bidang elektrostatik, gravitasi, dan bidang lain di mana potensi \(V\) harus dihitung sebagai fungsi posisi. Jika ada muatan atau massa dalam ruang, dan jika kita menggeneralisasi ke tiga dimensi, persamaannya menjadi persamaan Poisson

\[\begin{equation} \frac{\partial V}{\partial x^2}+\frac{\partial V}{\partial y^2}+\frac{\partial V}{\partial z^2}=0 \tag{10.20} \end{equation}\]

Bergantung pada geometri masalahnya, Persamaan (10.20) juga dapat ditulis dalam koordinat bola, silindris, atau lainnya.

Untuk menyelesaikan persamaan eliptik jenis ini, persamaan harus diberi syarat batas. Biasanya dengan menentukan bahwa titik, garis, atau permukaan tertentu dipertahankan pada nilai konstan potensial. Kemudian potensi di titik lain disesuaikan sampai mencapai perkiraan yang diinginkan. (Dalam kasus yang jarang terjadi, persamaan dengan kondisi batas dapat diselesaikan dengan tepat secara analitis; tetapi biasanya solusi perkiraan harus cukup.)

Ada banyak pendekatan untuk solusi numerik dari persamaan Laplace. Mungkin yang paling sederhana adalah Jacobi, di mana titik-titik interior secara berturut-turut didekati dengan rata-rata titik-titik di sekitarnya, sedangkan titik-batas dibatasi pada nilai-nilai tetap dan yang ditentukan. Kita asumsikan sebagai contoh bidang persegi, dibatasi oleh \(\left(0,1\right)\) dalam arah \(x\) dan \(y\), di mana ujung pada \(y = 1\) dipertahankan pada \(V = 1\) dan tiga tepi lainnya dipertahankan pada \(V = 0\). Kita buat tebakan awal yang paling sederhana untuk potensi di titik interior, tetapi ini akan disamakan saat solusinya menyatu.

Pada contoh berikut, kita akan melakukan simulasi persamaan Laplace menggunakan metode Jacobi. Parameter yang digunakan dalam simulasi adalah sebagai berikut:

Tebakan awal untuk profil voltase adalah sebagai berikut:

Kondisi batas di tetapkan seperti berikut:

V[1,] <- 0 V[n,] <- 0 V[,1] <- 0 V[,n] <- V0*rep(1,n)

Selanjutnya visualisasikan tebakan awal.

Gambar 10.10: Visualisasi tebakan awal solusi persaman Laplace

Proses iterasi menggunakan metode Jacobi ditampilkan pada sintaks berikut:

Hasil simulasi selanjutnya di visualisasikan kembali.

Gambar 10.11: Visualisasi hasil simulasi solusi persaman Lplace

## Penyelesaian Persamaan Diferensial Parsial Menggunakan PaketReacTran`

Package ReacTran memfasilitasi pemodelan transport reaktif dalam 1, 2, dan 3 dimensi. Paket ini "berisi rutinitas yang memungkinkan pengembangan model transportasi reaktif dalam sistem perairan (sungai, danau, lautan), media berpori (agregat flok, sedimen, …) dan bahkan organisme yang diidealkan.

Pada paket ReacTran terdapat sejumlah fasilitas fungsi, antara lain:

• Fungsi untuk menyiapkan kisi beda hingga (1D atau 2D)
• Fungsi untuk melampirkan parameter dan properti ke kisi ini (1D atau 2D)
• Fungsi untuk menghitung jangka waktu transport adveksi-difusi di atas grid (1D, 2D, 3D)
• Berbagai fungsi lainnya.

Saat paket ReacTran dimuat, paket ini juga memuat dua paket pendukung, yaitu: rootSolve dan deSolve. Paket rootSolve berguna untuk memecahkan persamaan diferensial untuk kondisi tunak baik persamaan diferensial uniform atau multikomponen (1D, 2D, dan 3D). Sedangkan deSolve berguna untuk menyediakan fungsi yang digunakan untuk memperoleh penyelesaian numerik persamaan diferensial biasa (ODE), persamaan diferensial parsial (PDE), persamaan aljabar diferensial (DAE) dan persamaan delay differential.

Fungsi setup.grid.1D() digunakan untuk membentuk kisi satu dimensi. Secara sederhana fungsi ini membagi ruang satu dimensi \(L\) antara \(x.up\) dan \(x.down\) menjadi sejumlah \(N\) kisi sebesar dx.1. Format umum fungsi tersebut, adalah sebagai berikut:

setup.grid.1D(x.up=0, x.down=NULL,L=NULL, N=NULL, dx.1=NULL, p.dx.1= rep(1,length(L)), max.dx.1=L, dx.N=NULL, p.dx.N=rep(1,length(L)), max.dx.N=L)

Catatan:

• x.up : posisi hilir.
• x.down : posisi hulu.
• L : x.down-x.up.
• N : jumlah kisi = L/dx.1.

Pada situasi yang lebih kompleks, ukuran sel atau kisi dapat bervariasi, atau dapat pula lebih dari satu zona. Kondisi ini dijelaskan lebih jauh pada laman bantuan fungsi.

Nilai yang dihasilkan dari fungsi setup.grid.1D() termasuk x.mid (vektor sepanjang \(N\) yang menyatakan posisi titik tengah) dan x.int (vektor sepanjang \(N+1\) yang menyatakan posisi antar muka antara kisi sel dimana flux diukur).

Fungsi plot untuk grid.1D() memvisualissikan baik posisi sel dan ketebalan kotak, menampilkan x.mid dan x.int. Contoh di halaman bantuan menunjukkan perilaku ini.

setup.grid.1D() berfungsi sebagai titik awal untuk setup.grid.2D, yang membuat kisi di atas domain persegi panjang yang ditentukan oleh dua kisi 1D ortogonal.

Banyak model transportasi akan melibatkan kisi-kisi dengan properti konstan. Tetapi jika beberapa properti yang mempengaruhi difusi atau adveksi bervariasi dengan posisi di grid, variasi dapat digabungkan dengan fungsi setup.prop.1D() (atau setup.prop.2D() dalam dua dimensi).

Diberikan fungsi matematis atau matriks data, fungsi setup.prop.1D() menghitung nilai properti yang diminati di tengah sel kisi dan pada antarmuka antar sel. Format fungsi tersebut adalah sebagai berikut:

Catatan:

• func : fungsi yang mengatur ketergantungan spasial properti.
• value : nilai konstan yang diberikan ke properti jika tidak ada ketergantungan spasial.
• xy : matriks data di mana kolom pertama memberikan posisi, dan kolom kedua memberikan nilai-nilai yang diinterpolasi di atas grid.
• iterpolate : metode interpolasi yang digunakan (spline atau linier).
• grid: objek yang dibentuk melalui fungsi setup.grid.1D().
• ...: argumen tambahan func.

Fungsi ini menghitung istilah transportasi — laju perubahan konsentrasi akibat difusi dan adveksi — dalam model cairan 1D (fraksi volume = 1) atau padatan berpori (fraksi volume dapat bervariasi dan <1).tran.1D() juga digunakan untuk masalah dalam geometri bola atau silinder, meskipun dalam kasus ini antarmuka sel jaringan akan memiliki area variabel. Format fungsi ini adalah sebagai berikut:

tran.1D(C, C.up = C[1], C.down = C[length(C)], flux.up = NULL, flux.down = NULL, a.bl.up = NULL, a.bl.down = NULL, D = 0, v = 0, AFDW = 1, VF = 1, A = 1, dx, full.check = FALSE, full.output = FALSE)

Catatan:

• C : vektor konsentrasi pada titik tengah kisi sel.
• C.up, C.down : konsentrasi pada hulu dan hilir batas.
• flux.up,flux.down : flux dari dan keluar sistem pada hulu dan hilir batas.
• Jika terdapat tranport konvektif sepanjang hulu dan hilir batas, a.bl.up dan a.bl.down merupakan koefisiennya.
• D: Koefisien difusi.
• v: koefisien adveksi.
• VF: fraksi volume.
• A: fraksi area.
• dx: ketebalan kisi, baik nilai konstan atau vektor.
• full.check, full.output: nilai logik untuk memeriksa konsistensi dan mengatur output dari perhitungan. Keduanya FALSE secara default.

Ketika full.output = FALSE, nilai-nilai yang dikembalikan oleh trans.1D() adalah dC, yaitu: laju perubahan C di pusat setiap sel kisi karena transportasi, dan flux.up dan flux.down, yatu: fluks masuk dan keluar dari model di batas hulu dan hilir.

ReacTran juga memiliki fungsi untuk memperkirakan istilah difusi dan adveksi dalam model dua dan tiga dimensi, dan dalam koordinat silinder dan kutub. Jumlah input tumbuh dengan dimensi, tetapi input pada dasarnya sama seperti pada kasus 1D. Lihat halaman bantuan untuk tran.2D(), tran.3D(), tran.cylindrical(), dan tran.polar().

Namun penyempurnaan lain adalah fungsi tran.volume.1D(), yang memperkirakan istilah transportasi volumetrik dalam model 1D. Berbeda dengan tran.1D(), yang menggunakan fluks (massa per satuan luas per satuan waktu), tran.volume.1D() menggunakan aliran (massa per satuan waktu). Ini berguna untuk memodelkan saluran yang area lintas bagiannya berubah, ketika perubahan area tidak perlu dimodelkan secara eksplisit. Ini juga memungkinkan input lateral dari saluran samping.

Setelah kisi telah diatur dan properti ditetapkan dan model transport telah diformulasikan dengan tran.1D() (atau analog 2D atau 3D-nya), maka ReacTran memanggil ode.1D() dari paket deSolve jika solusi time-dependent diperlukan, atau stable.1D() dari paket rootSolve jika solusi kondisi-stabil diinginkan. Sistem ODE yang dihasilkan dari metode pendekatan garis biasanya sparse dan non-stiff. Integrator dalam deSolve, seperti “lsoda” (metode default 1D) sangat cocok untuk menangani sistem persamaan tersebut. Jika sistem ODE non-stiff, maka “adams” umumnya merupakan metode yang baik.

Pada contoh kali ini, kita akan memodifikasi persamaan difusi satu dimensi yang telah dilakukan sebelumnya dengan menambahkan bagian adveksi serta diselesaikan menggunakan paket Reactran.

Siapkan kisi menggunakan fungsi setup.grid.1D() dan persiapkan nilai parameter persamaan.

Bentuk fungsi yang menyatakan persamaan adveksi-difusi.

Inisiasi konsentrasi awal pada kisi sel.

Lakukan perhitungandengn menggunakan time step 0,01.

Visualisasikan hasil perhitungan.

Gambar 10.12: Visualisasi hasil simulasi persamaan adveksi-difusi menggunakan paket ReacTran

Persamaan (10.16) dapat diselesaikan dengan cara yang sama dengan persamaan difusi dengan mengatur nilai \(c^2=D\), memisalkan \(W=u\) dan \(\frac{\partial u}{\partial t}=v\), dan menyelesaikannya dengan cara yang familiar untuk variabel berpasangan \(\left(u,v\right)\). Di sini kita misalkan persamaan gelombang 1D untuk senar yang dipetik, yang awalnya ditahan pada \(0\) amplitudo untuk \(x <−25\) dan \(x> 25\), dan direntangkan secara linear hingga maksimum pada \(x = 0\). ode.1D() digunakan untuk menyelesaikan set himpunan ODE simultan dengan \(c = 1\).

Langkah pertama yang harus dilakukan adalah melakukan penetapan parameter.

Menetapkan kondisi awal seperti ketinggian senar dan kecepatan.

Menetapkan deret waktu yang digunakan dalam simulasi.

times <- seq(from = 0, to = 50, by = 1)

Membangun fungsi yang akan diselesaikan secara numerik.

Selesaikan persamaan menggunakan fungsi ode.1D() dengan metode adams.

Visualisasikan hasil perhitungan.

Gambar 10.13: Visualisasi hasil simulasi persamaan gelombang menggunakan paket ReacTran

Kita akan kembali melakukan simulasi terhadap Persamaan (10.19) menggunakan metode lainnya. Pada contoh kali ini, gradien pada sumbu \(y\) bernilai -1. Gradien hanyalah fluks, \(D \left(\partial C = \partial x\right)\), dengan \(D\) set sama dengan 1. Fungsi penyelesaian yang digunakan adalah steady.2D(), karena tidak ada ketergantungan waktu dalam persamaan. Sebagai kondisi awal yang arbitrer, kami menggunakan nomor acak \(N_x \times N_y\) yang terdistribusi secara seragam. Kita juga harus menentukan \(nspec\), jumlah spesies dalam model (hanya satu, potensi, dalam hal ini), \(dimens\), vektor dengan 2 nilai dengan jumlah sel dalam arah \(x\) dan \(y\), dan lrw, panjang array. Lihat halaman bantuan untuk steady.2D() untuk informasi lebih detail.

Langkah pertama dalam melakukan simulasi adalah menetapkan parameter model.

Bentuk fungsi yang akan diselesaikan secara numerik.

laplace <- function(t, U, parms) { w = matrix(nrow = Nx, ncol = Ny, data = U) dw = tran.2D(C = w, C.x.up = 0, C.y.down = 0, flux.y.up = 0, flux.y.down = -1, D.x = 1, D.y = 1, dx = xgrid, dy = ygrid)$dC list(dw) }

Mulia dengan bilangan acak uniform sebagai kondisi awal, kemudian selesaikan untuk memperoleh nilai pada kondisi tunak dan buat plot kontur.

Gambar 10.14: Visualisasi hasil simulasi persamaan laplace pada kondisi tunak menggunakan paket ReacTran

Pada contoh kali ini, kita akan menyelesaikan persamaan Poisson untuk dipol yang ditampilkan pada Persamaan (10.21).

\[\begin{equation} \frac{\partial^2w}{\partial x^2}+\frac{\partial^2w}{\partial y^2}=-\frac{\rho}{\varepsilon_0} \tag{10.21} \end{equation}\]

untuk dipol yang terletak di tengah selembar persegi jika tidak pada 0 potensial. Untuk mempermudah, kita dapat mengatur semua faktor skala sama dengan satu. Dalam definisi fungsi poisson, nilai-nilai dalam matriks \(N_x \times N_y\) w adalah input melalui vektor data \(U\). Seperti dalam persamaan Laplace di atas, kita menetapkan nilai awal \(w\) pada sel-sel grid sama dengan angka acak uniform.

Langkah pertama dalam melakukan simulasi adalah menetapkan parameter model.

Cari nilai \(x\) dan \(y\) pada titik kisi yang mendekati \(\left(0,4;0,5\right)\) untuk muatan positif dan \(\left(0,6;0,5\right)\) untuk muatan negatif.

Bentuk fungsi Poisson yang akan diselesaikan secara numerik.

poisson <- function(t, U, parms) { w = matrix(nrow = Nx, ncol = Ny, data = U) dw = tran.2D(C = w, C.x.up = 0, C.y.down = 0, flux.y.up = 0, flux.y.down = 0, D.x = 1, D.y = 1, dx = xgrid, dy = ygrid)$dC dw[ipos,jpos] = dw[ipos,jpos] + 1 dw[ineg,jneg] = dw[ineg,jneg] - 1 list(dw) }

Selesaikan untuk kondisi tunak dan buat visualisasinya.

Gambar 10.15: Visualisasi hasil simulasi persamaan Poisson pada kondisi tunak menggunakan paket ReacTran

Model Lotka-Voltera merupakan model yang menggambarkan interaksi antara predator dan mangsa. Pada studi kasus kali ini model yang akan digunakan merupakan model dengan 3 bentuk interaksi yaitu: tanaman \(u\), herbivora \(v\), dan karnivora \(w\). Sistem persamaan diferensial sistem tersebut ditampilkan pada kumpulan persamaan berikut:

\[\begin{equation} \frac{du}{dt}=au-\alpha_1f_1\left(u,v\right) \tag{10.22} \end{equation}\]

\[\begin{equation} \frac{dv}{dt}=-bv+\alpha_1f_1\left(u,v\right)-\alpha_2f_2\left(v,w\right) \tag{10.23} \end{equation}\]

\[\begin{equation} \frac{dw}{dt}=-c\left(w-w\ast\right)+\alpha_2f_2\left(v,w\right) \tag{10.24} \end{equation}\]

dimana \(w\ast\) merupakan tingkat predasi minimum untuk menstabilkan populasi ketika populasi mangsa rendah. Interaksi antar komponen digambarkan ke dalam bentuk persamaan logistik.

\[\begin{equation} f_i\left(x,y\right)=\frac{xy}{1+k_ix} \tag{10.25} \end{equation}\]

Langkah pertama untuk menyelesaikan model adalah melakukan penentuan parameter model dan pembentukan fungsi model.

library(deSolve) f <- function(x,y,k){x*y/(1+k*x)} model <- function(t, xx, parms) { u = xx[1] # plant resource v = xx[2] # herbivore w = xx[3] # carnivore with(as.list(parms),{ du = a*u - alpha1*f(u, v, k1) dv = -b*v + alpha1*f(u, v, k1) - alpha2*f(v, w, k2) dw = -c*(w - wstar) + alpha2*f(v, w, k2) list(c(du, dv, dw)) })} times <- seq(0, 200, 0.1) parms <- c(a=1, b=1, c=10, alpha1=0.2, alpha2=1, k1<-0.05, k2=0, wstar=0.006) xstart <- c(u=10, v=5, w=0.1)

Proses iterasi selanjutnya akan menggunakan metode lsoda, dimana metode iterasi ini dapat berganti secara otomatis dalam menangani sistem persamaan diferensial stiff dan non-stiff. Persamaan diferensial dikatakan stiff apabila variabel dependen berubah berdasarkan 2 atau lebih variabel independen yang sangat berbeda skalanya.

out = ode(xstart, times, model, parms, method ="lsoda")

Ketiga variabel selanjutnya divisualisasikan.

Gambar 10.16: Visualisasi simulasi model Lotka-Voltera

1. Bloomfield, V.A. 2014. Using R for Numerical Analysis in Science and Engineering. CRC Press.
2. Chapra, S.C. Canale, R.P. 2015. Numerical Methods For Engineers, Seventh Edition. Mc Graw Hill.
3. Griffiths, G.W. 2016. Numerical analysis using R : solutions to ODEs and PDEs. Cambridge University Press.
4. Howard, J.P. 2017. Computational Methods for Numerical Analysis with R. CRC Press.
5. Kreyszig, E. 2011. Advanced Engineering Mathematics, 10th Edition. John Wiley & Sons.
6. Soetaert,K., Cash J., Mazzia F. 2012. Solving Differential Equations in R. Springer.
7. Suparno, S. 2008. Komputasi untuk Sains dan Teknik Edisi II. Departemen Fisika-FMIPA Universitas Indonesia.

1. Tunjukkan 10 hasil iterasi metode Euler untuk persamaan diferensial \(f'\left(x,y\right)=y\) dimana \(x_0=0\), \(y_0=2\) dan step size \(h=0,1\)!
2. Kerjakan lagi soal no.1 dengan menggunakan metode lainnya dan bandingkan seluruh metode tersebut menggunakan penyelesaian analitiknya!
3. Bentuk kembali fungsi rk4sys() menggunakan algoritma metode Adams-Bashforth dan namai fungsi tersebut adamsys()!

Pada Chapter 11, kita akan membahas mengenai cara melakukan analisis data pada R. Pada Chapter ini penulis akan memperkenalkan fungsi-fungsi yang ada pada R yang dapat membantu kita menganalisis data dan melakukan sejumlah uji statistik.

Pada Chapter ini kita tidak akan berfokus pada persamaan-persamaan matematika yang menjadi dasar suatu uji statistik. Chapter ini menitik beratkan pada bagaimana pembaca dapat melakukan sejumlah uji statistik pada R dan gambaran metode yang digunakan.

Kita dapat melakukan import data dalam berbagai format pada R. Namun, pada sub-chapter ini hanya akan dibahas bagaimana cara mengimport data dari file dengan format .csvdan .txt. Secara umum fungsi-fungsi yang digunakan untuk membaca data pada file dengan format tersebut adalah sebagai berikut:

read.table(file, header = FALSE, sep = "", dec = ".", stringsAsFactors = default.stringsAsFactors()) read.csv(file, header = TRUE, sep = ",", dec = ".") read.csv2(file, header = TRUE, sep = ";", dec = ",") read.delim(file, header = TRUE, sep = "\t", dec = ".") read.delim2(file, header = TRUE, sep = "\t", dec = ",")

Catatan:

• file : lokasi dan nama file yang akan dibaca diakhiri dengan format file. Secara default fungsi akan membaca file yang ada pada working directory. Untuk mengetahui lokasi working directory, jalankan fungsi getwd(). Salin file yang akan dibaca pada lokasi working directory.
• header : nilai logik yang menunjukkan apakah baris pertama pada file yang dibaca akan dibaca sebagai nama kolom.
• sep : simbol yang menujukkan pemisah antar data. Pemisah antar data dapat berupa “",”;“,”.", dll.
• dec : simbol yang menujukkan desimal. Pemisah desimal dapat berupa “.” atau “,”.
• stringsAsFactors : nilai logik yang menunjukkan apakah jenis data string akan dikonversi menjadi factor.

Kelima fungsi tersebut digunakan untuk membaca data tabular atau data yang disusun kedalam format tabel. Fungsi read.table() merupakan bentuk umum dari keempat fungsi lainnya. Fungsi tersebut dapat digunakan untuk membaca data dalam kedua format yang telah disebutkan sebelumnya. Fungsi lainnya lebih spesifi, dimana fungsi read.csv() dan read.csv2() digunakan untuk membaca data dengan ekstensi .csv, sedangkan read.delim() dan read.delim2() untuk membaca data dengan ekstensi .txt. Berikut adalah contoh bagaimana cara membaca data dengan nama data.csv yang ada pada working directory dengan pemisah antar data berupa ; dan tanda koma berupa ,:

data <- read.table(file="data.csv", sep=";", dec=",")

Terdapat sejumlah fungsi yang akan pembaca sering gunakan untuk mengecek dataset yang akan pembaca analisa. Fungsi-fungsi tersebut antara lain:

• head(): mengecek \(n\) (default 6) observasi teratas.
• tail(): mengecek \(n\) (default 6) observasi terbawah.
• str(): mengecek struktur data atau jenis data pada masing-masing kolom. Jenis data yang ada pada R dapat berupa num (numerik), int (integer), Factor(factor), date (tanggal), dan chr (karakter atau string).
• summary(): ringkasan data.

Berikut adalah contoh penerapan fungsi-fungsi tersebut pada dataset iris:

# cek 10 observasi teratas head(iris, 10)

# cek 10 observasi terbawah tail(iris, 10)

# cek struktur data str(iris)

# ringkasan data summary(iris)

Fungsi-fungsi lainnya yang dapat digunakan untuk melakukan analisis statistika deskriptif adalah sebagai berikut:

• mean() : menghitung nilai rata-rata variabel numerik.
• sd() : menghitung simpangan baku variabel numerik.
• var() : menghitung varians variabel numerik.
• median() : menghitung median suatu variabel numerik.
• range() : memperoleh nilai minimum dan maksimum suatu variabel numerik.
• IQR() : memperoleh nilai jarak antar kuartil.
• quantile() : memperoleh kuantil variabel numerik.

Berikut adalah contoh penerapan fungsi-fungsi tersebut:

attach(airquality) # rata-rata konsentrasi ozon mean(Ozone, na.rm = TRUE)

# median konsentrasi ozon median(Ozone, na.rm = TRUE)

# simpangan baku konsentrasi ozon sd(Ozone, na.rm = TRUE)

# varians konsentrasi ozon var(Ozone, na.rm = TRUE)

# range konsentrasi ozon range(Ozone, na.rm = TRUE)

# IQR konsentrasi ozon IQR(Ozone, na.rm = TRUE)

Pada analisis statistik inferensial khususnya pada pengujian hipotesis, asumsi normalitas merupakan sesuatu yang harus terpenuhi jika prosedur uji yang digunakan merupakan prosedur uji parametrik. Terdapat dua buah cara untuk melakukan uji tersebut, antara lain:

1. Metode grafis (qq-plot, ECDF, plot densitas, histogram, dan boxplot).
2. Metode matematis (Shapiro-Wilk, Cramer-von Mises, Shapiro-Francia, Anderson-Darling, Liliefors, Pearson Chi-square, dll).

Pada Chapter ini, kita akan berfokus pada uji matematis karena cara pengujian dengan menggunakan metode grafis telah penulis jabarkan pada Chapter visualisasi data.

Metode uji normalitas yang sering digunakan pada R adalah metode Shapiro-Wilk. Metode ini merupakan metode uji yang memiliki power yang besar khusunya untuk ukuran sampel yang relatif kecil. Versi awal metode ini terbatas dengan jumlah sampel 3 sampai 50 sampel. Versi selanjutnya mengalami modifikasi sehingga dapat menangani sampel sampai dengan 5000 sampel bahkan lebih.

Untuk melakukan uji SHapiro-Wilk pada R, pembaca dapat menggunakan fungsi shapiro.test(). Format fungsi tersebut adalah sebagai beriku:

Catatan:

Untuk lebih memahami impelementasi fungsi tersebut pada data, berikut adalah contoh penerapan fungsi tersebut untuk menguji normalitas distribusi konsentrasi ozon pada dataset airquality:

shapiro.test(airquality$Ozone)

Berdasarkan hasil uji diperoleh nilai p-value < 0,05, sehingga \(H_0\) ditolak dan dapat disimpulkan bahwa distribusi konsentrasi ozon tidak mengikuti distribusi normal. Untuk lebih memahami prosedur pengujian normalitas distribusi suatu data pembaca dapat membaca lebih lanjut pada tautan Environmental Data Modelin.

Uji rata-rata satu sampel merupakan uji statistik untuk menguji apakah rata-rata suatu sampel berasal dari suatu populasi yang telah diketahui nilai rata-ratanya. Sedangkan uji rata-rata untuk dua populasi dilakukan untuk menguji apakah kedua selisis rata-rata populasi tersebut bernilai nol yang menujukkan bahwa kedua populasi tersebut memiliki nilai rata-rata yang sama. Uji rata-rata dua populasi dapat dilakukan untuk sampel independen (contoh: uji rata-rata performa dua buah IPAL) dan berpasangan (contoh: uji rata-rata input dan output IPAL).

Untuk melakukan uji rata-rata pada R dapat digunakan fungsi t.test() untuk uji parametrik dan wilcox.test() untuk melakukan uji non-parametrik sign rank test. Format fungsi-fungsi tersebut adalah sebagai berikut:

t.test(x, y = NULL, alternative = c("two.sided", "less", "greater"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, ...) wilcox.test(x, y = NULL, alternative = c("two.sided", "less", "greater"), mu = 0, paired = FALSE, conf.level = 0.95, ...)

Catatan:

• x,y : vektor numerik. Jika argumen x dan y diisikan maka uji hipotesis dilakukan untuk dua buah populasi.
• alternative: digunakan untuk menentukan jenis uji hipotesis apakah satu sisi(“less” dan “greater”), atau dua sisi (“two.sided”).
• mu : nilai rata-rata populasi atau nilai rata-rata selisih antar populasi jika dilakukan uji hipotesis terhadap dua populasi. Secara default nilainya 0.
• paired : nilai logikal yang menentukan apakah uji dua populasi digunakan untuk sampel berpasangan (TRUE) atau tidak (FALSE).
• var.equal : nilai logikal yang menunjukkan apakah varians kedua populasi diasumsikan sama atau berbeda.
• conf.level : tingkat kepercayaan. Secara default tingkat kepercayaan yang digunakan adalah 95%.

Berikut adalah contoh penerapan fungsi tersebut untuk uji hipotesis satu dan dua populasi:

# Uji hipotesis dua populasi dni3 <- dimnames(iris3) ii <- data.frame(matrix(aperm(iris3, c(1,3,2)), ncol = 4, dimnames = list(NULL, sub(" L.",".Length", sub(" W.",".Width", dni3[[2]])))), Species = gl(3, 50, labels = sub("S", "s", sub("V", "v", dni3[[3]])))) # parametrik t.test(x=iris$Sepal.Length[iris$Species=="setosa"], y=ii$Sepal.Length[iris$Species=="versicolor"])

Fungsi t.test() akan menghasilkan output berupa nilai t uji, derajat kebebasan (df), nilai p-value, rentang estimasi nilai rata-rata berdasarkan tingkat kepercayaan yang digunakan, serta estimasi nilai rata-rata sampel. Fungsi wilcox.test() akan menghasilkan dua buah output yaitu nilai W dan p-value berdasarkan nilai W yang dihasilkan.

Pada sebuah analisa, kita sering kali tertarik untuk menganalisa hubungan atau korelasi antara satu variabel terhadap variabel lainnya. Pengamatan adanya korelasi antar variabel dapat dilakukan melalui visualisasi menggunakan scatterplot dan perhitungan matematis menggunakan metode Pearson untuk metode parametrik dan metode rangking Spearman dan Kendall untuk metode non-parametrik. Pada Chapter ini kita akan berfokus untuk melakukan uji korelasi menggunakan R menggunakan metode matematis.

Pada R uji korelasi dapat dilakukan dengan menggunakan fungsi cor.test(). Format fungsi tersebut adalah sebagai berikut:

Catatan:

• x,y : vektor numerik.
• alternative: digunakan untuk menentukan jenis uji hipotesis apakah satu sisi(“less” dan “greater”), atau dua sisi (“two.sided”).
• method : metode perhitungan korelasi yang digunakan.
• conf.level : tingkat kepercayaan. Secara default tingkat kepercayaan yang digunakan adalah 95%.

Berikut adalah penerapan fungsi cor.test() berdasarkan metode-metode yang telah disediakan pada fungsi tersebut:

Berdasarkan output yang dihasilkan, metode Pearson menghasilkan output berupa nilai t uji, derajat kebebasan, nilai p-value, rentang estimasi nilai korelasi berdasarkan tingkat kepercayaan, dan estimasi nilai korelasi. Metode Kendall dan Spearman disisi lai menghasilkan output berupa nilai z uji dan S untuk masing-masing metode serta nilai p-value berdasarkan nilai statistika uji dan estimasi koefisien korelasi.

Pada sub-Chapter sebelumnya penulis telah menjelaskan uji rata-rata untuk satu sampel dan dua sampel. Pada kenyataannya dalam sebuah percobaan laboratorium, kita tidak hanya membandingkan dua buah grup sampel saja, namun beberapa grup dan sejumlah faktor. Untuk menganalisa apakah variasi perlakuan pada kelompok sampel akan memberikan hasil yang berbeda-beda pada rata-rata tiap grup atau tidak diperlukan analisis varians untuk menganilisa variasi perlakuan atau faktor pada masing-masing grup. Analisis varians dapat dilakukan baik untuk satu faktor maupun dua faktor atau lebih. Untuk melakukannya pada R, kita dapat menggunakan fungsi aov() untuk analisis varians dengan metode parametrik dan kruskal.test() untuk analisis varians dengan menggunakan metode nonparametrik. Berikut adalah format kedua fungsi tersebut:

aov(formula, data = NULL) kruskal.test(formula, data)

Catatan:

• formula : formula model yang digunakan.
• data: dataset yang akan digunakan.

Berikut adalah contoh penerapan kedua fungsi tersebut untuk melihat apakah terdapat beda pada rata-rata konsentrasi bulanan ozon menggunakan dataset airquality:

summary(aov(Ozone~Month, airquality))

kruskal.test(Ozone~Month, airquality)

Berdasarkan hasil yang diperoleh diketahui bahwa rata-rata konsentrasi bulanan ozon tidak sama tiap bulannya atau minimal terdapat satu bulan dimana konsentrasi ozonnya berbeda secara signifikan dengan konsentrasi ozon pada bulan-bulan lainnya. Untuk lebih memahami terkait analisis varians pada R dan cara membaca output kedua fungsi tersebut, pembaca dapat membaca tulisan pada halaman situs sthda.

Analisis komponen utama menggunakan transformasi ortogonal (umumnya nilai singular atau dekomposisi nilai eigen) untuk mengubah seperangkat variabel pengamatan yang mungkin berkorelasi menjadi seperangkat variabel tidak berkorelasi (ortogonal) yang disebut komponen utama. Transformasi didefinisikan sedemikian rupa sehingga komponen utama pertama memiliki varians setinggi mungkin (menyumbang variabilitas pada data sebanyak mungkin), dan masing-masing komponen berikutnya pada gilirannya memiliki varians tertinggi yang mungkin di bawah kendala, dimana komponen tersebut menjadi ortogonal ke komponen sebelumnya.

Dalam R, analisis komponen utama umumnya dilakukan dengan fungsi prcomp (). Format fungsi tersebut adalah sebagai berikut:

prcomp(x, retx = TRUE, center = TRUE, scale. = FALSE, tol = NULL)

Catatan:

• x : data frame atau matriks kompleks numerik.
• retx : nilai logik yang mengindikasikan apakah variabel hasil rotasi perlu ditampilkan.
• center : nilai logik yang mengidikasikan apakah variabel perlu dilakukan pergeseran sehingga nilai rata-ratanya berpusat pada nilai nol.
• scale : nilai logik yang mengidikasikan apakah variabel perlu dilakukan penskalaan sebelum dilakukan analisis.
• tol : nilai toleransi yang menunjukkan batas nilai bagi komponen yang akan dipertahankan. Komponen yang dihilangkan jika simpangan bakunya kurang dari atau sama dengan tol x simpangan baku pc1.

Untuk memahami penerapan fungsi tersebut, kita akan melakukan simulasi menggunakan dataset iris. Output yang dihasilkan di bawah ini menunjukkan bagaimana empat variabel numerik ditransformasikan menjadi empat komponen utama. Penskalaan data mungkin tidak diperlukan dalam kasus ini, karena keempat pengukuran memiliki unit yang sama dan besarnya sama. Namun, penskalaan umumnya merupakan praktik yang baik.

Fungsi summary memberikan proporsi varians total yang dikaitkan dengan masing-masing komponen utama, dan proporsi kumulatif ketika masing-masing komponen ditambahkan. Kita melihat bahwa dua komponen pertama berperan lebih dari 95% dari total varians.

Histogram (hasil plot dalam analisis prcomp) secara grafis merekapitulasi proporsi varian yang disumbangkan oleh masing-masing komponen utama, sementara biplot menunjukkan bagaimana variabel awal diproyeksikan pada dua komponen utama pertama (Gambar 11.1). Ini juga menunjukkan koordinat dari masing-masing sampel dalam ruang (PC1, PC2). Satu spesies iris (yang berubahmenjadi setosa dari analisis kluster di bawah ini) secara jelas dipisahkan dari dua spesies lain dalam ruang koordinat ini.

Gambar 11.1: Analisis komponen utama data iris.

Untuk informasi lebih lanjut terkait metode analisis komponen utama, pembaca dapat membacanya pada laman Little Book of R for Multivariate Analysis.

Analisis cluster mencoba untuk mengurutkan satu set objek ke dalam kelompok (cluster) sedemikian rupa sehingga objek dalam cluster yang sama lebih mirip satu sama lain dibandingkan objek pada cluster lainnya. Ini digunakan untuk analisis eksplorasi melalui proses penambangan data (data mining) di banyak bidang, seperti bioinformatika, biologi evolusi, analisis gambar, lingkungan, dan pembelajaran mesin.

Menurut Wikipedia: “Analisis Cluster itu sendiri bukanlah salah satu algoritma spesifik, tetapi tugas umum yang harus dipecahkan. Ini dapat dicapai dengan berbagai algoritma yang berbeda secara signifikan dalam pengertian mereka tentang apa yang merupakan sebuah cluster dan bagaimana cara menemukannya secara efisien. Gagasan populer mengenai cluster termasuk kelompok dengan jarak rendah di antara anggota cluster, area padat ruang data, interval atau distribusi statistik tertentu. Algoritma pengelompokan dan pengaturan parameter yang sesuai (termasuk nilai-nilai seperti fungsi jarak yang akan digunakan, ambang kepadatan atau jumlah cluster yang diharapkan) tergantung pada dataset individual dan tujuan penggunaan hasil. Analisis cluster seperti itu bukan tugas otomatis, tetapi proses berulang penemuan pengetahuan yang melibatkan trial and error. Seringkali diperlukan untuk memodifikasi preprocessing dan parameter sampai hasilnya mencapai properti yang diinginkan.

Hierarchical clustering membangun hierarki cluster, di mana metrik hierarki adalah suatu ukuran ketidaksamaan antar cluster. Menurut halaman bantuan untuk hclust(), metode pengelompokan hierarkis aglomeratif, “Fungsi ini melakukan analisis hierarki cluster menggunakan seperangkat ketidaksamaan untuk n objek yang dikelompokkan. Awalnya, masing-masing objek ditugaskan ke cluster sendiri dan kemudian algoritma melanjutkan secara iteratif, pada setiap tahap bergabung dengan dua cluster yang paling mirip, terus sampai hanya ada satu cluster. Pada setiap tahap, jarak antar kluster dihitung ulang dengan formula pembaruan ketidaksamaan Lance Williams sesuai dengan metode pengelompokan tertentu yang digunakan.”Ada tujuh metode aglomerasi yang tersedia, dengan lengkap — yang mencari kluster kompak, bola — sebagai default. Format umum fungsi hclust() adalah sebagai berikut:

hclust(d, method = "complete", members = NULL)

Catatan:

• d : struktur ketidaksamaan yang dihasilkan dengan fungsi dist.
• method : metode alglomerasi yang digunakan. Metode yang dapat digunakan antara lain: “ward.D”, “ward.D2”, “single”, “complete”, “average” (= UPGMA), “mcquitty” (= WPGMA), “median” (= WPGMC) atau “centroid” (= UPGMC)
• members: NULL atau vektor dengan ukuran sama dengan d. Untuk info lebih lanjut jalankan sintaks ?hclust.

Untuk memahami penerapan fungsi hclust(), kita akan kembali menggunakan data iris_use. Berikut adalah sintaks yang digunakan:

Visualisasi cluster dibentuk melalui dendogram yang ditampilkan pada Gambar 11.2.

Gambar 11.2: Analisis pengelompokan hierarkis alglomeratif data iris.

Menurut halaman bantuan diana() (DIvisive ANAlysis Clustering) dalam library cluster, “Algoritma diana membangun hierarki pengelompokan, dimulai dengan satu kluster besar yang berisi semua n pengamatan. Cluster dibagi sampai masing-masing cluster hanya berisi satu pengamatan. Pada setiap tahap, cluster dengan diameter terbesar dipilih. Format fungsi diana() adalah sebagai berikut:

diana(x, metric = "euclidean", stand = FALSE)

Catatan:

• x : struktur ketidaksamaan yang dihasilkan dengan fungsi dist atau data frame.
• metric : karakter string yang menyatakan metode pengukuran jarak yang digunakan. Metode pengukuran jarak dapat berupa “euclidean” dan “manhattan”.
• stand: vektor logik yang menyatakan apakah data akan dilakukan standardisasi terlebih dahulu sebelum dilakukan analisis.

Gambar 11.3: Analisis pengelompokan divisif data iris.

k-means melakukan pengelompokan n pengamatan ke dalam k cluster di mana setiap pengamatan akan tergabung dengan pusat cluster terdekat. Pengguna harus menentukan jumlah pusat (cluster) yang diinginkan sebagai output. Untuk melakukan pengelompokan dengan algoritma k-means pada R dapat menggunakan fungsi kmeans(). Format fungsi tersebut secara umum adalah sebagai berikut:

Catatan:

• x : data frame.
• centers : jumlah cluster yang ingin di buat.
• iter.max: jumlah iterasi maksimum yang diijinkan.
• algorithm : algoritma pengelompokan yang digunakan. Untuk informasi lebih lanjut jalankan sintaks ?kmeans.

iris_kmeans <- kmeans(iris_use, centers = 3) iris_kmeans

pam mem-partisi data menjadi k cluster di sekitar medoid. Medoid dari set data yang terbatas merupakan titik data dengan nilai ketidaksamaan rata-rata untuk semua titik data adalah minimum. Hal tersebut menujukkan bahwa medoid merupakan pusat dari set cluster. Menurut halaman bantuan pam(), pendekatan k-medoid lebih kuat daripada pendekatan k-means “karena meminimalkan jumlah ketidaksamaan daripada jumlah jarak euclidean kuadrat”. Format umum fungsi pam() adalah sebagai berikut:

pam(x, k, diss = inherits(x, "dist"), metric = c("euclidean", "manhattan"))

Catatan:

• x : data frame.
• k : jumlah cluster yang ingin di buat.
• method: metode perhitungan jarak yang digunakan. Untuk informasi lebih lanjut jalankan sintaks ?pam.

library(cluster) iris_pam <- pam(iris_use, k=3) iris_pam

Gambar 11.4: Analisis pengelompokan metode pam data iris.

Pada Chapter 12, kita akan membahas cara membentuk model statistik menggunakan R. Terdapat 2 buah jenis model yang akan dibahas pada Chapter ini, yaitu: regresi dan klasifikasi. Untuk informasi terkait cara untuk melakukan inferensi berdasarkan hasil yang diperoleh dan cara untuk melakukan prediksi menggunakan model yang terbentuk tidak akan dijelaskan dalam buku ini. Pembaca dapat membaca lebih lanjut pada referensi berikut:

Regresi linier merupakan model sederhana yang paling sering dibahas dalam buku-buku statistika. Modelnya cukup sederhana dimana kita berusaha membentuk model dengan pendekatan garis linier dengan prinsip meminimalkan jumlah kuadrat residual pada data. Model yang tebentuk akan menghasilkan dua buah nilai yaitu nilai konstanta (titik potong sumbu y) dan nilai slope kurva. Model yang terbentuk secara umum haruslah memenuhi asumsi dasar model linier berikut:

1. Asumsi liniearitas: kurva relasi yang terbentuk antara variabel independen terhadap variabel dependen harus linier. Asumsi ini dapat dipelajari melalui plot residual terhadap nilai fitted value. Jika asumsi liniearitas terpenuhi, maka titik-titik residual yang di plotkan akan membentuk pola acak. Jika pada plot yang dihasilkan terbentuk pola tidak linear maka transformasi data pada variabel prediktor atau independen diperlukan.
2. Error atau residu berdristribusi normal: normalitas error di cek menggunakan qq-plot atau uji normalitas yang telah dibahas pada Chapter 11.4.
3. Outlier dan high influence point: kedua pengamatan tersebut dideteksi melalui qq-plot, plot residual terhadap nilai fitted value, dan plot residuals vs leverage. Jika outlier terjadi akibat adanya error selama pengukuran maka outlier dapat dihilangkan.
4. Error bersifat independen: independensi residual dapat dideteksi melaui plot korelasi serial dengan mengeplotkan \(r_i\) vs \(r_{i-1}\).
5. Varians bersifat konstan: Varians bersifat konstan dicek melalui plot square root standardize residual vs fitted value. Pada kasus dimana varians tidak bersifat konstan, kita dapat memberikan bobot pada model yang akan kita bentuk (weighted least square), dimana bobot yang diberikan proporsional dengan invers varians.
6. multikolinearitas: tidak ada variabel dependen yang saling berfkorelasi. Multikolinearitas dapat dideteksi melalui plot matriks korelasi. Pada model adanya kolinearitas ditunjukkan dari nilai variance inflation factor (VIF) yang tinggi. Secara umum nilai VIF terkecil sebesar 1 dan jika kolinearitas terjadi nilainya dapat lebih besar dari 5 atau 10. Untuk mengatasi kolinearitas pada model dapat dilakukan dengan dua cara, yaitu: mengeluarkan variabel dengan nilai VIF yang tinggi pada model atau menggabungkan dua variabel prediktor yang saling berkorelasi menjadi satu variabel baru.

Pembentukan model linier pada R dilakukan dengan menggunakan fungsi lm(). Format umum fungsi tersebut adalah sebagai berikut:

lm(formula, data, subset, weights)

Catatan:

• formula : formula model yang hendak dibentuk.
• data: data yang digunakan untuk membentuk model.
• subset : subset data yang akan digunakan dalam pembentukan model.
• weight : nilai pembobotan dalam pembentukan model.

Pada Chapter 12.1.1 akan diberikan contoh pembentukan model linier sederhana menggunakan dataset Boston dari library MASS dengan jumlah observasi sebesar 506 observasi. Pada contoh kali ini kita akan mencoba membentuk model dengan variabel dependen berupa medv (median harga rumah) dan variabel independen berupa lstat (persen rumah tangga dengan status ekonomi menengah ke bawah). Berikut adalh sintaks untuk membentuk model tersebut:

lm.fit <- lm(medv~lstat, data=Boston) anova(lm.fit)

Plot residual disajikan pada Gambar 12.1.

Gambar 12.1: Analisis residual model linier medv vs lstat pada dataset Boston.

Berdasarkan hasil plot dapat dilihat bahwa seluruh asumsi model linier tidak terpenuhi. Selain melalui plot residual, uji asumsi model linier dapat juga dilakukan secara matematis. Berikut adalah sintaks yang digunakan:

# varians bersifat konstan # (varians tidak konstan) library(lmtest)

Pada Chapter 12.1.2, kita akan membuat tiga buah model regresi linier. Model pertama akan menambahkan variabel age (usia bangunan) pada model sebelumnya, model kedua akan menggunakan seluruh ariabel yang ada, dan model ketiga akan melakukan pembaharuan dengan mengeluarkan variabel dengan VIF paling tinggi dari model kedua. Berikut adalah sintaks untuk membentuk ketiag model tersebut:

Berdasarkan hasil perhitungan diketahui nilai VIF dari model < 10, sehingga asumsi multikolinearitas terpenuhi. Untuk asumsi lainnya dapat dicek pada plot residual yang ditampilkan pada Gambar 12.2.

Gambar 12.2: Analisis residual model linier 1 pada dataset Boston.

# Model 2 lm.fit2 <- lm(medv~., data=Boston) anova(lm.fit2)

Berdasarkan hasil perhitungan diperoleh nilai VIF untuk seluruh varaibel prediktor dalam model < 10, sehingga asumsi multikolinearitas terpenuhi. Untuk asumsi lainnya dapat dicek pada plot residual yang ditampilkan pada Gambar 12.3.

Gambar 12.3: Analisis residual model linier 2 pada dataset Boston.

Pada model ketiga, kita akan mencoba untuk melakukan pembaharuan pada model kedua dengan melakukan drop variabel dengan vif yang paling tinggi. Pada hasil perhitungan sebelumnya, variabel tax (pajak) memiliki nilai VIF yang paling tinggi, sehingga pada model ketiga variabel tersebut tidak disertakan. Terdapat dua cara untuk melakukannya berikut adalah sintaks yang digunakan:

Plot residual ditampilkan pada Gambar 12.4.

Gambar 12.4: Analisis residual model linier 3 pada dataset Boston.

Interaksi antar variabel pada model linier dapat dengan mudah dimasukkan kedalam fungsi lm(). Terdapat dua buah cara untuk melakukannya. Cara pertama dengan menggunakan tanda : pada formula (contoh: \(y1 ~ x1+x2+x1:x2\)). Tanda : menyatakan formula persamaan linier memasukkan interaksi antar variabel prediktor di dalamnya. Cara kedua adalah dengan menggunakan tanda *. Cara ini lebih sederhana, dimana fungsi lm() akan secara otomatis menerjemahkannya sebagai serangkaian variabel tunggal dan interaksinya. Berikut adalah contoh penerapannya menggunakan kedua cara tersebut:

# Cara 2 lm.inter <- lm(medv~lstat*age, data=Boston) anova(lm.inter)

Plot residual ditampilkan pada Gambar 12.5.

Gambar 12.5: Analisis residual model dengan melibatkan interaksi antar variabel pada dataset Boston.

Fungsi lm() juga dapat melibatkan transformasi non-linier prediktor pada argumen formula-nya. Transformasi non-linier dilakukan dengan menambahkan fungsi identitas I(). Sebagai contoh model berikut melibatkan transformasi kuadrat pada variabel lstat:

lm.trans <- lm(medv~lstat+I(lstat^2), data=Boston) anova(lm.trans)

Cara yang lebih sederhana untuk melibatkan tranformasi polinomial kedalam model linier adalah dengan menggunakan fungsi poly(). Berikut adalah contoh penerapannya:

lm.trans <- lm(medv~poly(lstat,2), data=Boston) anova(lm.trans)

Plot residual ditampilkan pada Gambar 12.6.

Gambar 12.6: Analisis residual model dengan melibatkan transformasi non-linier variabel prediktor pada dataset Boston.

FUngsi trasnformasi lainnya juga dapat digunakan pada pembentukan model linier. Berikut adalah contoh penerapan transformasi logaritmik dan eksponensial pada model linier:

summary(lm(medv~log(lstat), data=Boston))

summary(lm(medv~exp(lstat), data=Boston))

Regresi linier dapat pula dilakukan dengan jenis variabel prediktor berupa variabel kategorikal. Untuk dapat melakukannya jenis data variabel kategorikal terlebih dahulu dirubah kedalam factor.

Pada contoh kali ini, kita akan menggunakan dataset utds dari library edtaR untuk memodelkan konsentrasi uranium dalam air tanah menggunakan variabel TDS pada berbagai kondisi kesadahan. Berikut adalah sintaks yang digunakan

lm.fit <- lm(Uranium~., data=utds) anova(lm.fit)

Plot residual yang ditampilkan pada Gambar 12.7.

Gambar 12.7: Analisis residual model dengan melibatkan variabel kategorikal pada dataset utds.

Pada pembahasan sebelumnya kita telah menyaksikan bahwa sebagian model yang telah terbentuk tidak memenuhi asumsi varians yang konstan. Untuk mengatasi hal tersebut, kita dapat membentuk regresi dengan memberikan bobot sebesar invers variansnya. Untuk melakukannya diperlukan beberapa tahapan, antara lain:

1. Membentuk model linier dari variabel dataset: fit <- lm(y~(variabel prediktor)).
2. Menghitung error absolut (abse <- abs(resid(fit))) dan nilai *fitted value dari model (yhat <- fitted(fit).
3. Membentuk kembali model menggunakan data residual absolut (efit <- lm(abse~poly(yhat,2))) dan menghitung residual fitted value (shat <- fitted(efit)).
4. Gunakan nilai bobot w <- 1/shat^2 untuk membentuk model regresi dengan pembobotan (fitw <-lm(y~(variabel prediktor), weights=w)).

Kita akan membentuk kembali model menggunakan dataset Boston dengan menggunakan seluruh variabel, namun pada model kali ini kita akan memberikan bobot pada model yang terbentuk. Berikut adalah sintaks yang digunakan:

Plot residual ditampilkan pada Gambar 12.8.

Gambar 12.8: Analisis residual model regresi linier dengan pembobotan pada dataset Boston.

Pada Chapter 12.2, kita telah membahas cara untuk membangun model dengan output berupa variabel dengan nilai numerik. Pada Chapter 12.2, kita akan belajar cara membentuk model regresi dengan 2 respons (0 dan 1). Pada regresi ini pembentukan model didasarkan oleh kurva logistik, dimana melalui kurva tersebut nilai yang dihasilkan akan memiliki rentang dari 0 sampai 1. Karena model yang dibuat bertujuan untuk memprediksi dua buah kemungkinan (0 atau 1), maka diperlukan suatu nilai ambang (y<0,5 = 0 dan y >= 0,5 = 1).

Fungsi glm() dapat digunakan untuk membentuk model regresi logistik. Format umum fungsi tersebut adalah sebagai berikut:

glm(formula, family = gaussian, data, weights, subset, )

Catatan:

• formula : formula model yang hendak dibentuk.
• family : distribusi yang digunakan. Untuk regresi logistik digunakan argumen family=binomial
• data: data yang digunakan untuk membentuk model.
• subset : subset data yang akan digunakan dalam pembentukan model.
• weight : nilai pembobotan dalam pembentukan model.

Pada contoh berikut, kita akan membuat model untuk memprediksi contamination rating menggunakan dataset contamination dari library edtaR. Berikut adalah sintaks yang digunakan:

1. Akritas, M. 2016. PROBABILITY & STATISTICS WITH R FOR ENGINEERS AND SCIENTISTS. Pearson.
2. Bloomfield, V.A. 2014. Using R for Numerical Analysis in Science and Engineering. CRC Press.
3. James, G., Witten, D., Hastie, T., Tibshirani, R. 2013. An Introduction to Statistical Learning. Springer.
4. Kerns, G.J., 2018. Introduction to Probability and Statistics Using R. Course notes for University of Auckland Paper STATS 330. http://ipsur.r-forge.r-project.org/book/download/IPSUR.pdf.
5. Lee, A., Ihaka, R., Triggs, C. 2012. ADVANCED STATISTICAL MODELLING.
6. Primartha, R. 2018. Belajar Machine Learning Teori dan Praktik. Penerbit Informatika : Bandung.
7. Rosadi,D. 2016. Analisis Statistika dengan R. Gadjah Mada University Press: Yogyakarta.
8. STHDA. <(http://www.sthda.com/english/>