trapezoid <- function(ftn, a, b, n = 10) {
h <- (b - a) / n
x <- seq(a, b, by = h)
fx <- sapply(x, ftn)
I <- (h/2) * (fx[1] + 2*sum(fx[2:n]) + fx[n+1])
return(I)
}Praktikum Pemrograman Statistika
Menyusun Luas dari Titik-Titik: Eksperimen Integrasi Numerik di R
Pendahuluan
Dalam berbagai persoalan kuantitatif, sering kali diperlukan perhitungan luas di bawah kurva dari suatu fungsi yang tidak dapat diintegralkan secara analitik. Hal ini terjadi ketika fungsi tersebut terlalu kompleks, tidak memiliki bentuk antiturunan yang sederhana, atau bahkan hanya tersedia dalam bentuk data empiris. Untuk situasi seperti itu, digunakan pendekatan integrasi numerik (numerical integration) untuk memperoleh nilai hampiran integral melalui pendekatan komputasional.
Integrasi numerik memegang peranan penting dalam statistika dan sains data, terutama dalam menghitung nilai harapan, probabilitas kumulatif, atau area di bawah fungsi kepadatan peluang yang tidak dapat dievaluasi secara langsung. Dengan memanfaatkan algoritma numerik, integral kontinu dapat direpresentasikan sebagai penjumlahan diskret dari nilai fungsi pada sejumlah titik tertentu.
Beberapa pendekatan umum yang digunakan dalam integrasi numerik meliputi:
- Metode Trapezoidal, yang mendekati fungsi dengan segmen garis lurus pada setiap subinterval;
- Metode Simpson, yang memperbaiki ketepatan pendekatan dengan polinomial kuadratik;
- Metode Gauss (Adaptive) Quadrature, yang memanfaatkan titik dan bobot tidak seragam berdasarkan teori polinomial ortogonal; serta
- Metode Monte Carlo Integration, yang mengestimasi integral melalui rata-rata nilai fungsi dari bilangan acak.
Masing-masing metode tersebut memiliki karakteristik tersendiri dalam hal ketelitian, efisiensi, dan kompleksitas komputasi. Pemahaman yang baik terhadap prinsip dan perilaku tiap metode menjadi dasar penting dalam pengembangan teknik komputasi yang andal untuk pemodelan dan analisis data berbasis integral.
Konsep Dasar Integral Tentu
Gambar di atas memperlihatkan representasi geometris dari integral tentu
\[ \int_a^b f(x)\,dx, \]
yang menggambarkan luas daerah di bawah kurva fungsi \(f(x)\) pada interval \([a,b]\).
Konsep ini menjadi dasar bagi berbagai metode integrasi numerik yang bertujuan menghampiri nilai luas tersebut ketika fungsi \(f(x)\) tidak dapat diintegralkan secara analitik.
Secara intuitif, nilai integral tentu merepresentasikan akumulasi total perubahan dari suatu fungsi pada selang tertentu. Dalam praktik komputasi, pendekatan integral dilakukan dengan membagi area tersebut menjadi bagian-bagian kecil yang bentuknya sederhana—seperti trapesium, parabola, atau himpunan titik acak—yang kemudian dijumlahkan untuk mendekati nilai sebenarnya dari integral.
Metode Trapezoidal
Metode Trapezoidal merupakan pendekatan paling sederhana dalam integrasi numerik.
Ide dasarnya adalah menghampiri fungsi \(f(x)\) dengan garis lurus di setiap subinterval \([x_i, x_{i+1}]\), kemudian menghitung luas di bawah garis tersebut sebagai luas trapesium.
Luas satu trapesium dengan panjang sisi sejajar \(f(x_i)\) dan \(f(x_{i+1})\), serta lebar alas \(h=x_{i+1}-x_i\), diberikan oleh:
\[ L_i = \frac{1}{2} \, h \, [f(x_i) + f(x_{i+1})] \]
Jika interval \([a,b]\) dibagi menjadi \(n\) bagian yang sama besar \((x_0 = a, x_n = b)\), maka luas total seluruh trapesium dapat dituliskan sebagai:
\[ \int_a^b f(x)\,dx \approx \frac{h}{2}\Big[f(x_0) + 2\sum_{i=1}^{n-1} f(x_i) + f(x_n)\Big]. \]
Semakin besar jumlah subinterval \(n\), semakin halus pendekatan garis lurus terhadap kurva, dan semakin kecil pula kesalahan (error) integrasi yang dihasilkan.
Implementasi metode ini dapat dibuat dengan fungsi sederhana sebagai berikut:
Fungsi di atas melakukan langkah-langkah berikut:
- Menentukan lebar subinterval \(h\);
- Membentuk vektor titik \(x_i\) antara \(a\) dan \(b\);
- Menghitung nilai fungsi \(f(x_i)\) pada setiap titik;
- Menjumlahkan luas seluruh trapesium sesuai rumus di atas.
Sebagai contoh, integral dari fungsi \(f(x)=e^{-x}\cos(x)\) pada interval \([0,\pi]\) dapat dihitung secara eksak menggunakan pendekatan analitik sebagai berikut:
\[ \int_0^{\pi} e^{-x}\cos(x)\,dx = \frac{e^{-x}}{2}\,[\sin(x) + \cos(x)] \Big|_0^{\pi} = \frac{1}{2}\,(1 + e^{-\pi}) \approx 0.521607. \]
Nilai tersebut merupakan hasil eksak integral fungsi \(f(x)\).
Selanjutnya, hasil tersebut dapat dibandingkan dengan pendekatan numerik menggunakan metode Trapezoidal.
f <- function(x) exp(-x) * cos(x)
trapezoid(f, 0, pi)[1] 0.5302151Hasil perhitungan memberikan nilai hampiran sekitar 0.5302151. Akurasi metode ini akan meningkat apabila jumlah pembagian n diperbesar.
trapezoid(f, 0, pi, n = 100)[1] 0.5216928trapezoid(f, 0, pi, n = 1000)[1] 0.5216078trapezoid(f, 0, pi, n = 10000)[1] 0.521607trapezoid(f, 0, pi, n = 100000)[1] 0.521607trapezoid(f, 0, pi, n = 1000000)[1] 0.521607Nilai hasil integrasi numerik akan semakin mendekati 0.521607 seiring meningkatnya jumlah subinterval \(n\), menunjukkan konvergensi metode Trapezoidal terhadap hasil eksak.
Metode Simpson
“Simpson’s Family Picture” — tidak ada hubungannya dengan integral, tapi cukup untuk mengingatkan bahwa yang akan dibahas berikutnya adalah Simpson’s Rule 😄
Metode Simpson merupakan penyempurnaan dari metode Trapezoidal. Pendekatan ini menghampiri fungsi \(f(x)\) bukan dengan garis lurus, melainkan dengan polinomial kuadratik (parabola) pada setiap dua subinterval. Dengan demikian, fungsi di antara tiga titik berurutan \((x_{i-1}, f(x_{i-1}))\), \((x_i, f(x_i))\), dan \((x_{i+1}, f(x_{i+1}))\) dianggap mengikuti bentuk parabola.
Berbeda dengan metode Trapezoidal yang menggunakan garis lurus, metode Simpson menggunakan parabola untuk menghampiri bentuk kurva asli.
Visualisasi berikut memperlihatkan perbedaan pendekatan kedua metode tersebut.
Sumber: Marta Caligaris et al. (2015). Procedia - Social and Behavioral Sciences, 176, 270–275.
Gambar di atas memperlihatkan bahwa pendekatan Simpson lebih baik dalam menyesuaikan bentuk kurva asli dibandingkan pendekatan garis lurus pada metode Trapezoidal. Area bayangan pada Simpson’s Rule lebih mendekati area sebenarnya di bawah kurva, yang menyebabkan hasil integralnya lebih akurat walaupun menggunakan jumlah titik evaluasi yang hampir sama.
Jika interval \([a,b]\) dibagi menjadi \(n\) bagian yang sama besar (dengan \(n\) harus genap), maka rumus Simpson’s Rule dapat dinyatakan sebagai:
\[ \int_a^b f(x)\,dx \approx \frac{h}{3}\Big[f(x_0) + 4\sum_{i=1,3,5,\ldots}^{n-1} f(x_i) + 2\sum_{i=2,4,6,\ldots}^{n-2} f(x_i) + f(x_n)\Big] \]
dengan \(h = \frac{b - a}{n}\) sebagai lebar subinterval.
Bobot \(4\) dan \(2\) pada masing-masing suku muncul karena fungsi diaproksimasi oleh parabola yang melalui tiga titik, di mana titik tengah (ganjil) memiliki pengaruh lebih besar terhadap bentuk kurva.
Semakin besar nilai \(n\), semakin halus parabola yang digunakan untuk mendekati kurva asli, sehingga menghasilkan nilai integral yang lebih akurat dibandingkan metode Trapezoidal.
Implementasi metode Simpson dapat ditulis dalam fungsi berikut:
simpson_n <- function(ftn, a, b, n = 10) {
# Pastikan n genap
n <- max(c(2 * (n %/% 2), 4))
h <- (b - a) / n
# Titik ganjil dan genap
x.odd <- seq(a + h, b - h, by = 2 * h)
x.even <- seq(a + 2 * h, b - 2 * h, by = 2 * h)
# Evaluasi fungsi
fx.0 <- ftn(a)
fx.n <- ftn(b)
fx.odd <- sapply(x.odd, ftn)
fx.even <- sapply(x.even, ftn)
# Rumus Simpson
I <- (h / 3) * (fx.0 + fx.n + 4 * sum(fx.odd) + 2 * sum(fx.even))
return(I)
}Langkah-langkah yang dilakukan fungsi di atas adalah:
- Menentukan jumlah pembagian \(n\) yang genap;
- Menghitung lebar subinterval \(h\);
- Menentukan titik-titik ganjil dan genap di dalam interval;
- Menghitung nilai fungsi pada setiap titik;
- Menjumlahkan seluruh komponen sesuai bobot (4 untuk ganjil, 2 untuk genap) dalam rumus Simpson.
Sebagai contoh, dihitung integral dari fungsi \(f(x)=e^{-x}\cos(x)\) pada interval \([0,\pi]\) seperti pada metode trapezoidal dengan hasil eksak \(\approx 0.521607\). Selanjutnya, dilakukan perbandingan hasil pendekatan numerik dengan metode Simpson untuk berbagai jumlah pembagian \(n\).
simpson_n(f, 0, pi)[1] 0.5214968simpson_n(f, 0, pi, n = 100)[1] 0.5216069simpson_n(f, 0, pi, n = 1000)[1] 0.521607simpson_n(f, 0, pi, n = 10000)[1] 0.521607Hasil perhitungan menunjukkan bahwa nilai hampiran dengan metode Simpson sangat cepat mendekati nilai eksak 0.521607, bahkan dengan jumlah subinterval yang relatif kecil.
Hal ini menunjukkan bahwa metode Simpson memiliki orde ketelitian yang lebih tinggi dibandingkan metode Trapezoidal. Secara teoritis, kesalahan pendekatan Simpson berorde \(O(h^4)\), sedangkan Trapezoidal berorde \(O(h^2)\). Dengan demikian, setiap kali jumlah subinterval digandakan, galat (error) metode Simpson menurun jauh lebih cepat dibandingkan metode Trapezoidal.
Metode Gauss (Adaptive) Quadrature
“Carl Friedrich Gauss” — tidak sedang mengajarkan integral di kelas, tapi tanpa idenya kita tidak akan punya metode Gauss Quadrature 😄
Metode Gauss Quadrature merupakan pengembangan dari metode integrasi numerik yang bertujuan untuk mencapai ketelitian tinggi tanpa harus menambah banyak titik evaluasi seperti pada metode Trapezoidal atau Simpson.
Ide utamanya adalah memilih titik evaluasi (\(x_i\)) dan bobot (\(w_i\)) yang tidak seragam, sehingga kesalahan integral menjadi sekecil mungkin untuk kelas fungsi polinomial hingga orde tertentu.
Berbeda dari metode sebelumnya yang menggunakan titik-titik berjarak sama, metode Gauss memilih titik \(x_i\) sebagai akar-akar dari polinomial ortogonal tertentu—biasanya Legendre polynomial untuk interval standar \([-1,1]\).
Dengan demikian, integral pada interval \([a,b]\) dapat didekati sebagai:
\[ \int_a^b f(x)\,dx \approx \sum_{i=1}^{n} w_i f(x_i), \]
dengan \(x_i\) adalah titik kuadratur (nodes) dan \(w_i\) adalah bobot (weights) yang ditentukan dari tabel atau hasil perhitungan polinomial ortogonal.
Untuk mengubah integral dari interval \([a,b]\) ke bentuk standar \([-1,1]\), digunakan transformasi:
\[ x = \frac{b-a}{2}t + \frac{b+a}{2}, \quad dx = \frac{b-a}{2}dt, \]
sehingga integral dapat ditulis kembali sebagai:
\[ \int_a^b f(x)\,dx = \frac{b-a}{2} \int_{-1}^1 f\!\left(\frac{b-a}{2}t + \frac{b+a}{2}\right) dt. \]
Sebelum mengimplementasikan metode Gauss Quadrature di R, perlu dipahami bahwa setiap orde (\(n\)) memiliki pasangan titik (\(x_i\)) dan bobot (\(w_i\)) yang berbeda. Pasangan ini diturunkan dari akar-akar polinomial Legendre \(P_n(x)\) dan digunakan untuk menghitung integral dengan ketelitian optimal.
Tabel berikut menunjukkan contoh nilai titik dan bobot untuk beberapa orde \(n\) pertama pada interval standar \([-1,1]\).
| n | \(x_i\) (Titik Gauss) | \(w_i\) (Bobot) |
|---|---|---|
| 2 | \(x_1 = -0.5773503\), \(x_2 = 0.5773503\) | \(w_1 = w_2 = 1\) |
| 3 | \(x = \{-0.7745967, 0, 0.7745967\}\) | \(w = \{0.555556, 0.888889, 0.555556\}\) |
| 4 | \(x = \{-0.8611363, -0.3399810, 0.3399810, 0.8611363\}\) | \(w = \{0.347855, 0.652145, 0.652145, 0.347855\}\) |
| 5 | \(x = \{-0.906180, -0.538469, 0, 0.538469, 0.906180\}\) | \(w = \{0.236927, 0.478629, 0.568889, 0.478629, 0.236927\}\) |
Titik-titik \(x_i\) simetris terhadap nol, sedangkan bobot \(w_i\) bernilai positif dan jumlahnya sama dengan 2.
Dalam R, pasangan nilai ini dapat diakses secara langsung menggunakan library gaussquad.
Fungsi legendre.quadrature.rules() menghasilkan daftar aturan kuadratur hingga orde tertentu, masing-masing berisi vektor \(x_i\) dan \(w_i\).
Contoh untuk mengambil titik dan bobot orde 4:
library(gaussquad)Loading required package: orthopolynom# Ambil aturan Legendre Quadrature hingga orde 4
rules <- legendre.quadrature.rules(4)
# Lihat isi aturan untuk n = 4
rules[[4]]Nilai-nilai ini dapat langsung digunakan untuk menghitung integral dengan rumus Gauss Quadrature.
Implementasi Gauss–Legendre Quadrature di R dapat dilakukan menggunakan library gaussquad sebagai berikut:
library(gaussquad)
# Definisikan fungsi
f <- function(x) exp(-x) * cos(x)
# Gunakan Legendre Quadrature orde 4
Lq <- legendre.quadrature.rules(4)[[4]]
xi <- Lq$x
wi <- Lq$w
# Transformasi dari [-1,1] ke [0,pi]
a <- 0; b <- pi
x.trans <- ((b - a) / 2) * xi + (b + a) / 2
I <- ((b - a) / 2) * sum(wi * f(x.trans))
I[1] 0.5216002Kode di atas menjalankan langkah-langkah berikut:
- Mendefinisikan fungsi yang akan diintegralkan;
- Mengambil titik (\(x_i\)) dan bobot (\(w_i\)) dari Legendre polynomial orde 4;
- Mentransformasi interval \([-1,1]\) ke \([a,b]\);
- Menghitung nilai hampiran integral sesuai rumus kuadratur Gauss.
Hasil perhitungan memberikan nilai integral mendekati 0.521607, yang merupakan hasil eksak dari integral
\[ \int_0^{\pi} e^{-x}\cos(x)\,dx = \frac{1}{2}(1 + e^{-\pi}) \approx 0.521607. \]
Meskipun hanya menggunakan empat titik evaluasi, metode Gauss Quadrature sudah menghasilkan nilai yang sangat akurat, bahkan lebih efisien dibandingkan metode Trapezoidal atau Simpson dengan ratusan subinterval.
Pendekatan ini bekerja karena titik-titik kuadratur dan bobotnya dipilih secara optimal agar kesalahan integrasi minimum untuk polinomial hingga orde tertentu. Hal ini membuat metode Gauss sangat efisien untuk fungsi halus, terutama ketika bentuk eksplisit fungsi diketahui dan dapat dievaluasi dengan presisi tinggi.
Konsep Adaptive Quadrature
Metode Adaptive Quadrature memperluas ide Gauss Quadrature dengan menyesuaikan jumlah titik integrasi di setiap subinterval berdasarkan kompleksitas fungsi \(f(x)\). Jika fungsi berubah secara tajam atau berosilasi pada bagian tertentu, maka interval tersebut akan dibagi menjadi bagian yang lebih kecil. Sebaliknya, untuk bagian fungsi yang halus, jumlah titik integrasi dapat dikurangi tanpa mengurangi akurasi.
Dengan pendekatan adaptif ini, algoritma secara otomatis meningkatkan kepadatan titik integrasi di area yang sulit, dan menghemat komputasi di area yang sederhana. Hasilnya adalah keseimbangan yang baik antara akurasi tinggi dan efisiensi komputasi. Pendekatan ini sering digunakan dalam perangkat lunak komputasi ilmiah, misalnya fungsi integrate() di R atau quad() di Python (SciPy).
Bagaimana Menentukan Orde Gauss Quadrature yang Tepat
Salah satu aspek penting dalam penerapan metode Gauss Quadrature adalah menentukan orde (\(n\)) yang sesuai.
Nilai \(n\) menentukan jumlah titik evaluasi (\(x_i\)) dan bobot (\(w_i\)) yang digunakan dalam perhitungan integral.
Semakin besar \(n\), semakin banyak titik evaluasi yang digunakan, dan umumnya semakin tinggi pula ketelitian hasilnya.
Namun, peningkatan orde juga berarti meningkatnya beban komputasi. Oleh karena itu, pemilihan orde yang tepat menjadi kunci efisiensi.
Prinsip Ketepatan Orde
Keunggulan metode Gauss–Legendre Quadrature adalah bahwa dengan \(n\) titik kuadratur, metode ini menghasilkan hasil yang eksak untuk semua polinomial hingga orde \(2n-1\).
Sebagai contoh:
| Orde (\(n\)) | Tepat untuk polinomial hingga orde |
|---|---|
| 2 | 3 |
| 3 | 5 |
| 4 | 7 |
| 5 | 9 |
| 6 | 11 |
Hal ini berarti, semakin besar nilai \(n\), semakin tinggi derajat fungsi polinomial yang dapat diintegralkan tanpa kesalahan.
Panduan Praktis Pemilihan Orde
| Tipe Fungsi | Karakteristik | Rekomendasi Orde (\(n\)) |
|---|---|---|
| Fungsi halus dan sederhana (misal \(e^{-x}\), \(\sin(x)\), \(\cos(x)\)) | Perubahan perlahan, tanpa singularitas | 3–5 |
| Fungsi berosilasi ringan (misal \(e^{-x}\cos(x)\)) | Ada osilasi tapi kontinu | 4–6 |
| Fungsi kompleks atau sangat berosilasi (misal \(\sin(10x)\), \(\ln(x)\) dekat 0) | Perlu lebih banyak titik | 8–12 |
| Fungsi tak halus atau berbasis data empiris | Tidak cocok untuk Gauss murni → gunakan Adaptive atau Monte Carlo | – |
Eksperimen di R
Untuk memastikan kestabilan hasil terhadap perubahan orde, dapat dilakukan perbandingan nilai integral untuk berbagai nilai \(n\) sebagai berikut:
library(gaussquad)
f <- function(x) exp(-x) * cos(x)
a <- 0; b <- pi
for (n in c(2, 3, 4, 5, 6, 8, 10)) {
Lq <- legendre.quadrature.rules(n)[[n]]
xi <- Lq$x
wi <- Lq$w
x.trans <- ((b - a) / 2) * xi + (b + a) / 2
I <- ((b - a) / 2) * sum(wi * f(x.trans))
print(paste("n =", n, "→ integral =", round(I, 6)))
}[1] "n = 2 → integral = 0.533096"
[1] "n = 3 → integral = 0.524073"
[1] "n = 4 → integral = 0.5216"
[1] "n = 5 → integral = 0.521606"
[1] "n = 6 → integral = 0.521607"
[1] "n = 8 → integral = 0.521607"
[1] "n = 10 → integral = 0.521607"Kode di atas menghitung hasil integral untuk beberapa nilai orde \(n\). Jika hasilnya tidak banyak berubah setelah nilai tertentu, berarti ordo tersebut sudah cukup akurat untuk fungsi yang digunakan.
gauss_legendre <- function(ftn, a, b, n = 4, verbose = TRUE) {
if (!requireNamespace("gaussquad", quietly = TRUE)) {
stop("Paket 'gaussquad' belum terpasang. Instal dengan install.packages('gaussquad')")
}
rules <- gaussquad::legendre.quadrature.rules(n)
Lq <- rules[[n]]
xi <- Lq$x
wi <- Lq$w
x.trans <- ((b - a) / 2) * xi + (b + a) / 2
I <- ((b - a) / 2) * sum(wi * ftn(x.trans))
if (verbose) {
message("Titik evaluasi (xi): ", paste(round(x.trans, 6), collapse = ", "))
message("Bobot (wi): ", paste(round(wi, 6), collapse = ", "))
}
return(I)
}Fungsi gauss_legendre() di atas merupakan fungsi user-defined yang menggeneralisasi penerapan metode Gauss–Legendre Quadrature untuk berbagai fungsi dan batas integral, dengan jumlah titik \(n\) yang dapat ditentukan pengguna.
# Definisikan fungsi yang akan diintegralkan
f <- function(x) exp(-x) * cos(x)
# Hitung integral pada interval [0, π]
gauss_legendre(f, 0, pi, n = 4)Titik evaluasi (xi): 2.923466, 2.104837, 1.036755, 0.218127Bobot (wi): 0.347855, 0.652145, 0.652145, 0.347855[1] 0.5216002Metode Monte Carlo Integration
“Monte Carlo” — bukan kasino di Monako, tetapi metode stokastik yang memanfaatkan bilangan acak untuk menghitung integral 😄
Setelah memahami pendekatan deterministik seperti Trapezoidal, Simpson, dan Gauss Quadrature, kini kita beralih ke pendekatan yang stokastik, yaitu Metode Monte Carlo Integration. Berbeda dengan metode sebelumnya yang menggunakan titik-titik evaluasi berjarak tetap atau terpilih optimal, Monte Carlo Integration menggunakan bilangan acak untuk memperkirakan nilai integral.
Metode ini sangat berguna ketika fungsi yang ingin diintegralkan:
- tidak memiliki bentuk analitik yang sederhana,
- berada dalam dimensi tinggi, atau
- hanya dapat dievaluasi pada titik-titik tertentu (misalnya hasil simulasi atau data empiris).
Untuk integral satu variabel pada interval \([a,b]\):
\[ I = \int_a^b f(x)\,dx \]
Metode Monte Carlo mendekatinya dengan:
\[ I \approx (b - a) \cdot \frac{1}{N} \sum_{i=1}^{N} f(x_i), \]
di mana:
- \(x_i\) diambil secara acak dari distribusi Uniform(\(a,b\)),
- \(N\) adalah jumlah sampel acak yang digunakan.
Secara intuitif, kita menghitung rata-rata nilai fungsi pada sejumlah titik acak, lalu mengalikannya dengan lebar interval.
Hubungan ini juga dapat dilihat dari sifat ekspektasi:
\[ \mathbb{E}[f(X)] = \frac{1}{b-a} \int_a^b f(x)\,dx, \]
sehingga estimasi integralnya adalah:
\[ I \approx (b - a)\, \mathbb{E}[f(X)]. \]
Implementasi fungsi Monte Carlo Integration dapat dibuat secara sederhana sebagai berikut:
montecarlo_integral <- function(ftn, a, b, N = 10000, verbose = TRUE) {
# Bangkitkan N bilangan acak uniform di [a,b]
x <- runif(N, a, b)
# Evaluasi fungsi pada titik-titik acak tersebut
fx <- ftn(x)
# Estimasi integral
I <- (b - a) * mean(fx)
return(I)
}Langkah-langkah dalam fungsi di atas:
- Membentuk \(N\) titik acak seragam antara \(a\) dan \(b\);
- Menghitung nilai fungsi \(f(x_i)\) pada titik-titik tersebut;
- Mengambil rata-ratanya;
- Mengalikan dengan panjang interval \((b - a)\) untuk mendapatkan nilai integral hampiran.
Sebagai contoh, kita gunakan fungsi yang sama seperti pada metode sebelumnya, yaitu \(f(x)=e^{-x}\cos(x)\) pada interval \([0,\pi]\), dengan hasil eksak \(\approx0.521607\).
# Definisikan fungsi
f <- function(x) exp(-x) * cos(x)
# Jalankan metode Monte Carlo
montecarlo_integral(f, 0, pi, N = 10000)[1] 0.5206168Hasilnya mendekati nilai eksak \(0.521607\). Namun, karena menggunakan bilangan acak, setiap kali fungsi dijalankan, hasilnya akan sedikit berbeda. Semakin besar jumlah sampel \(N\), semakin stabil dan akurat hasilnya.
Gambaran intuitif metode Monte Carlo dapat divisualisasikan dengan titik-titik acak di bawah kurva fungsi.
library(ggplot2)Warning: package 'ggplot2' was built under R version 4.4.3set.seed(123)
N <- 1000
x <- runif(N, 0, pi)
y <- runif(N, 0, 1)
# Titik di bawah kurva
under_curve <- y <= exp(-x) * cos(x) + abs(min(exp(-x) * cos(x)))
# Buat data frame
df <- data.frame(
x = x,
y = y,
under_curve = under_curve
)
# Data kurva fungsi
curve_df <- data.frame(
x = seq(0, pi, length.out = 200)
)
curve_df$y <- exp(-curve_df$x) * cos(curve_df$x) + abs(min(exp(-curve_df$x) * cos(curve_df$x)))
# Plot dengan ggplot2
ggplot(df, aes(x, y)) +
geom_point(aes(fill = under_curve), shape = 21, color = "black", size = 3) +
scale_fill_manual(values = c("gray80", "skyblue"),
labels = c("Di atas kurva", "Di bawah kurva"),
name = "") +
geom_line(data = curve_df, aes(x, y), color = "orange", size = 1) +
labs(
title = "Visualisasi Monte Carlo Integration",
x = "x", y = "y"
) +
theme_minimal(base_size = 13) +
theme(
legend.position = "top",
plot.title = element_text(face = "bold", size = 14, hjust = 0.5)
)Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
Visualisasi di atas memperlihatkan bahwa luas di bawah kurva diestimasi dari proporsi titik acak yang jatuh di area tersebut.
Metode Monte Carlo memiliki sifat konvergensi yang khas:
\[ \text{Error} \propto \frac{1}{\sqrt{N}}, \]
yang berarti untuk meningkatkan ketelitian dua kali lipat, jumlah sampel harus ditingkatkan empat kali lipat. Meskipun konvergensinya relatif lambat dibandingkan metode deterministik seperti Simpson atau Gauss, Monte Carlo sangat kuat untuk integrasi berdimensi tinggi, di mana metode deterministik menjadi tidak efisien.