Algoritma Iteratif, Root Finding, dan Integral Numerik

Annebel Diestya Clarissa

4/27/2021

Algoritma Iteratif

Algoritma adalah urutan langkah untuk menyelesaikan masalah dengan baik, yaitu untuk memperoleh keluaran (output) sesuai dengan masukan (input). Sedangkan proses iteratif merupakan suatu proses dalam algoritma yang melakukan perhitungan dengan hasil yang diingin menggunakan cara siklus berulang melalui sebuah operasi perulangan.

Iteratif menggambarkan situasi dimana urutan instruksi dapat dieksekusi berulang kali. Dalam sekali eksekusi biasa disebut juga 1 iterasi. Dan jika iterasi dilakukan berulang kali, proses tersebut biasa dikenal sebagai loop. Sehinga proses iterasi terjadi karena adanya loop.

Ciri-ciri algoritma iteratif :
* pemanggilan procedure beberapa kali atau hingga suatu kondisi tertentu terpenuhi
* Tidak ada variabel lokal baru
* Program tidak sederhana
* Menggunakan perulangan for dan while

Terdapat beberapa algoritme iteratif yang bisa diterapkan untuk melakukan pengurutan:

  • Bubble Sort

  • Selection Sort

  • Insertion Sort

  • Merge Sort

  • Quick Sort

Namun yang akan kita bahas selanjutnya adalah Bubble Sort yaitu metode pengurutan algoritma dengan cara melakukan penukaran data secara terus menerus sampai bisa dipastikan dalam suatu iterasi tertentu tidak ada lagi perubahan/penukaran.

Cara kerja dari iterasi bubble sort dilakukan dalam 2 tingkat yaitu :

  1. Iterasi untuk memastikan semua elemen diperiksa

  2. Iterasi untuk mengakses elemen-elemen vektor yang digunakan untuk penukaran antar elemen

Pada iterasi tingkat kedua nilai-nilai yang kecil akan dipindahkan ke depan sehingga selanjutnya elemen-elemen di bagian depan tidak perlu diperiksa lagi.
Misal pada saat iterasi tingkat pertama sama dengan 3 maka elemen yang akan dibandingkan hanya dari element 4 sampai 7 saja (ingat ada j+1). Semakin meningkat iterasi tingkat 1 maka banyaknya iterasi tingkat 2 semakin mengecil (pada kasus ini diakomodir dengan nx-i).

Pengurutan data : Loop dengan For

Algoritma

Proses pertukaran

Program R

x=c(15,18,10,8,20) #Data yang akan diurutkan

n=length(x)

for(i in 2:n){
if (x[i]>x[1]){
a=x[1]
x[1]=x[i]
x[i]=a
}
}
a=x[1]
a
## [1] 20

Pengurutan data : Loop dengan While

Algoritma

Program R

x = c(15,18,10,8,20) #Data yang akan diurutkan

n=length(x) 

while(is.unsorted(x)){ 
  for (i in 1:(n - 1)) { 
    j=i+1 
    if (x[i] > x[j]) 
      { a = x[i] 
      x[i] = x[j] 
      x[j] = a 
      } 
    } 
  } 

x
## [1]  8 10 15 18 20

Pendugaan Jumlah deret tak hingga

Contoh 1

\[\mathrm{z=1+\frac{5}{6}+\frac{7}{11}+\frac{9}{18}+....}\]

Identifikasi Bentuk deret

Bentuk Umum dari deret

\[\mathrm{Z=\sum_{i=1}^{\infty} \frac{2i+1}{i^2+2}}\]

Program R

e = 10
z0 = 0
i= 1

while(e>10**(5)){
    z1 = z0 + (2*i+1)/(i^2 + 2)
    e = abs(z1-z0)
    z0 = z1
    i = i +1
}

Contoh 2

\[\mathrm{z=\frac{2}{5}+\frac{3}{40}+\frac{6}{135}+\frac{11}{320}+\frac{18}{625}+\frac{27}{1080}+....}\]

Bentuk Umum dari deret

\[\mathrm{Z=\sum_{i=1}^{\infty} \frac{i^2+2n+3}{5n^3}}\]

Program R

# Menuliskan rumus fungsi dalam bentuk function
fz2 <- function(n){
  ((n^2) - 2*n + 3) / (5*n^3)
}

# Mendefinisikan nilai-nilai awal (initial value)
error <- 10 #Nilai awal kriteria berhenti
z0 <- 0 #Nilai awal deret
n <- 1 #Nilai awal iterasi

# Membuat iterasi untuk mencari jumlah deret
while (error>1e-6) {
  z <- z0 + fz2(n)
  error <- abs(z-z0)
  z0 <- z
  n <- n+1
}

# Hasil Akhir
z 
## [1] 2.61992
# Banyaknya Iterasi
n
## [1] 2e+05

Contoh 3

Berdasarkan jawaban Contoh 1 dan 2, buatlah fungsi yang bernama deret dengan ketentuan:

  1. Jika inputnya hanya fungsi deret bilangan, outputnya adalah jumlah deret dan banyaknya iterasi

  2. Jika inputnya fungsi deret bilangan dan bilangan n, outputnya adalah n dan suku ke-n dari deret tersebut

Jawaban a

Membuat fungsi deret

deret <- function(fx,n){
  if(missing(n)) {
    error <- 10 #Nilai awal kriteria berhenti
    z0 <- 0 #Nilai awal deret
    i <- 1 #Nilai awal iterasi
    
    # Membuat iterasi untuk mencari jumlah deret
    while (error>1e-6) {
      z <- z0 + fx(i)
      error <- abs(z-z0)
      z0 <- z
      i <- i+1
    }
    return(list(HasilAkhir=z,BanyakIterasi=i))
  } else {
    return(list(SukuKe=n, Nilai=fx(n)))
    }
}
# Contoh output fungsi deret (1 argumen input)
deret(fz2)
## $HasilAkhir
## [1] 2.61992
## 
## $BanyakIterasi
## [1] 2e+05

Jawaban b

# Contoh output fungsi deret (2 argumen input)
deret(fz2,n=5)
## $SukuKe
## [1] 5
## 
## $Nilai
## [1] 0.0288

Bilangan Fibonacci

Hitung penjumlahan dari barisan bilangan Fibonacci yang nilainya kurang dari 1000000 dan merupakan bilangan kelipatan 5!

# Membuat barisan Fibonacci
fibo <- c(0,1)
n <- length(fibo)
x <- fibo[n]+fibo[n-1]
while (x<1000000){
  n <- n+1
  fibo[n] <- x
  x <- fibo[n]+fibo[n-1]
}

# Barisan Fibonacci
fibo
##  [1]      0      1      1      2      3      5      8     13     21     34
## [11]     55     89    144    233    377    610    987   1597   2584   4181
## [21]   6765  10946  17711  28657  46368  75025 121393 196418 317811 514229
## [31] 832040
# Mencari suku barisan Fibonacci kelipatan 5
fibo[fibo %% 5 == 0]
## [1]      0      5     55    610   6765  75025 832040
# Menjumlahkan bilangan Fibonacci kelipatan 5
sum(fibo[fibo %% 5 == 0])
## [1] 914500

Menggambar Pola

Buatlah pola seperti gambar di atas menggunakan algoritme iteratif!

# Membuat pola dasar
mat <- matrix(NA,6,6)

# Membuat pola sesuai yang diinginkan
for (i in 1:6) {
  for (j in 1:(6-i+1)) {
    mat[i,j] <- 1
  }
}

# Hasil Akhir
mat
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    1    1    1    1    1
## [2,]    1    1    1    1    1   NA
## [3,]    1    1    1    1   NA   NA
## [4,]    1    1    1   NA   NA   NA
## [5,]    1    1   NA   NA   NA   NA
## [6,]    1   NA   NA   NA   NA   NA

Root Finding

Penentuan akar persamaan dari suatu fungsi dapat dilakukan dengan berbagai metode, seperti :

  • Metode Fixed Point
  • Metode Newton Raphson
  • Metode Secant
  • Metode Bisection

Metode Fixed Point

Metode fixed-point adalah metode iteratif untuk menyelesaikan persamaan dengan bentuk f(x) = x. Metode ini membangkitkan sekuen titik-titik x0, x1, x2, … xn yang diharapkan konvergen terhadap suatu titik sehingga f(a)=a. Suatu perkiraan nilai awal x0 diperlukan untuk menentukan dugaan berikutnya, yaitu x1=f(x0), x2=f(x1), … , xn+1=f(xn), terus berulang. Kelemahan dari metode ini adalah proses iterartif yang relatif lambat.

Formula dari Metode Fixed Point adalah :

\[\mathrm{x_{n+1}=g(x_n)}\]

fixedpoint <- function(ftn, x0, tol = 1e-9, max.iter = 100) {
  xold <- x0
  xnew <- ftn(xold)
  iter <- 1
  cat("At iteration 1 value of x is:", xnew, "\n")
  
  while ((abs(xnew-xold) > tol) && (iter < max.iter)) {
    # xold digunakan untuk menyimpan akar-akar persamaan dari iterasi sebelumnya
    xold <- xnew
    # xnew digunakan untuk menghitung akar-akar persamaan pada iterasi yang sedang berjalan
    xnew <- ftn(xold) #ingat ftn bentuknya fungsi
    # iter digunakan untuk menyimpan banyaknya iterasi
    iter <- iter + 1
    cat("At iteration", iter, "value of x is:", xnew, "\n")
  }
  # output bergantung pada kesuksesan algoritma
  if (abs(xnew-xold) > tol) {
    cat("Algorithm failed to converge\n")
    return(NULL)
  } else {
    cat("Algorithm converged\n")
    return(xnew)
  }
}

Metode Newton Raphson

Metode Newton-Raphson adalah pendekatan untuk menemukan akar persamaan nonlinier dan merupakan salah satu algoritma pencarian akar yang paling umum karena kesederhanaan dan kecepatannya. Dalam metode newton-raphson, akar dari suatu fungsi adalah titik di mana f(x) = 0. Newton-Raphson adalah metode berulang yang dimulai dengan tebakan awal akar (x0). Metode ini menggunakan turunan dari fungsi f′(x) serta fungsi asli f(x) sehingga hanya dapat digunakan jika turunannya dapat ditentukan.

Formula dari Metode Newton-Raphson adalah :

\[\mathrm{x_{n+1}=x_n - \frac{f(x_1)}{f'(x_1)}}\]

newtonraphson <- function(ftn, x0, tol = 1e-5, max.iter = 100) {
  x <- x0
  fx <- ftn(x)
  iter <- 0
  while ((abs(fx[1]) > tol) && (iter < max.iter)) {
    x <- x - fx[1]/fx[2]
    fx <- ftn(x)
    iter <- iter + 1
    cat("At iteration", iter, "value of x is:", x, "\n")
  }
  # output bergantung pada kesuksesan algoritma
  if (abs(fx[1]) > tol) {
    cat("Algorithm failed to converge\n")
    return(NULL)
  } else {
    cat("Algorithm converged\n")
    return(x)
  }
}

Metode Secant

Metode secant adalah metode berulang yang membutuhkan dua tebakan awal akar, (sebelumnya metode Newton-Raphson hanya menggunakan satu tebakan awal akar). Karena membutuhkan dua pendekatan akar awal, metode secant lebih lambat menuju nilai konvergen dibandingkan dengan metode Newton-Raphson, namun biasanya hasilnya lebih stabil. Metode secant juga memiliki kelebihan lainnya dibandingkan dengan metode Newton-Raphson karena tidak memerlukan turunan dari fungsi yang akan dicari akar persamaannya. Metode secant menggunakan garis-garis potong (karena itu diperlukan dua nilai awal awal) untuk menemukan akar suatu fungsi sedangkan metode Newton-Raphson mendekati akar dengan garis singgung.

Formula dari Metode Secant adalah :

\[\mathrm{x_{n+1}=x_n - f(x_n)\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})}}\]

secant <- function(ftn, x0, x1, tol = 1e-9, max.iter = 100) {
  xa <- x0
  fxa <- ftn(xa)
  xb <- x1
  fxb <- ftn(xb)
  iter <- 0
  while ((abs(xb-xa) > tol) && (iter < max.iter)) {
    xt <- fxb*((xb-xa)/(fxb-fxa))
    xc <- xb-xt
    # fxc <- ftn(xc)
    iter <- iter + 1
    cat("At iteration", iter, "value of x is:", xc, "\n")
    xa <- xb
    xb <- xc
    fxa <- ftn(xa)
    fxb <- ftn(xb)
  }
  if (abs(xb-xa) > tol) {
    cat("Algorithm failed to converge\n")
    return(NULL)
  } else {
    cat("Algorithm converged\n")
    return(xb)
  }
}

Metode Bisection

Metode Bisection adalah pendekatan lain untuk mencari akar dari fungsi kontinu f (x) pada interval [a, b]. Metode ini memanfaatkan teorema nilai tengah yang disebut teorema Bolzano yang menyatakan bahwa jika nilai f (a) dan f (b) memiliki tanda yang berlawanan, selang/interval harus mengandung setidaknya satu akar. Kelebihan dari metode bisection adalah langkah-langkah iterasi yang relatif lebih mudah, namun kekurangannya adalah konvergensi menuju solusi lebih lambat dibandingkan dengan metode pencarian akar lainnya.

Algoritma - Metode Bisection

Metode bisection dimulai dengan menghitung titik tengah, \(\mathrm{m=\frac {a+b}{2}}\) dari interval. Selanjutnya, fungsi dievaluasi pada titik f(m). Sesuai teorema Bolzano, jika f(a) dan f(m) memiliki tanda yang berlawanan, metode bisection menggantikan nilai b dengan titik tengah yang dihitung (m). Begitu juga jika f (b) dan f(m) memiliki tanda yang berlawanan, a akan diganti dengan titik tengah (m). Langkah ini memastikan masih ada root yang terkandung dalam interval. Prosedur kemudian berlanjut ke iterasi berikutnya. Solusinya dikatakan ditemukan setelah fungsinya sama dengan 0 pada f(m).

Sehingga jika dirangkum cara kerja dari metode bisection tersebut adalah :

Tetapkan dua nilai a dan b sehingga f(a)f(b) < 0.

  1. Jika |b – al| ≤ ε maka stop.
  2. Tetapkan \(\mathrm{m=\frac {a+b}{2}}\) ; jika f(m) = 0 maka stop.
  3. Jika f(a)f(m) < 0 maka tetapkan b = m; selainnya, tetapkan a = m.
  4. Lanjutkan ke langkah 1.
bisection <- function(ftn, xl, xr, tol = 1e-7,max.iter = 100) {
  # Jika tanda dari hasil evaluasi fungsi pada titik xl dan xr berbeda maka program dihentikan
  if ((ftn(xl) > 0) && (ftn(xr) < 0)) {
    stop('signs of ftn(xl) and ftn(xr) differ')
  }
  
  for (iter in 1:max.iter) {
    xm <- (xl + xr) / 2 # menghitung nilai tengah

    cat("At iteration", iter, "value of x is:", xm, "\n")       
    # Jika fungsinya sama dengan 0 di titik tengah atau titik tengah di bawah toleransi
    # yang diinginkan, hentikan dan keluarkan nilai tengah sebagai akar
    if ((ftn(xm) == 0) || abs(xr - xl) < tol) {
      return(xm)
    }
    # Jika diperlukan iterasi lain, periksa tanda-tanda fungsi pada titik xm dan xl
    # lalu tetapkan kembali xl atau xr sebagai titik tengah yang akan digunakan
    # pada iterasi berikutnya.
    ifelse(sign(ftn(xm)) == sign(ftn(xl)), xl <- xm, xr <- xm)
  }
}

Berikut adalah persamaan fungsi dan parameter yang diperlukan untuk masing-masing metode :

Metode Fungsi Parameter
Fixed Point Fungsi awal x0 : nilai awal
Newton Raphson Fungsi awal dan Turunan Fungsi x0 : nilai awal
Secant Fungsi Awal x0, x1 : 2 nilai awal
Bisection Fungsi Awal xl, xr : nilai rentang

Selain parameter tersebut, sebenarnya juga terdapat parameter nilai toleransi dna nilai maksimal iterasi dengan nilai default yaitu tol=1e-9 dan max.iter=100.

Contoh soal Root Finding

Contoh 1

Misalkan diberikan suatu persamaan \(\mathrm{x^3-2x-5=0}\).

  1. Gambarkan fungsi tersebut

Yang pertama harus dilakukan adalah mendefinisikan fungsi

# mendefinikan persamaan dalam bentuk fungsi
fx <- function(x) {
  x^3 - 2*x - 5
}

Selanjutnya adalah menggambarkan fungsi tersebut

curve(fx,xlim=c(-3,3), col='steelblue',lwd=2)
abline(h=0)
abline(v=0)

  1. Carilah akar-akar persamaan dengan menggunakan metode fixed-point, newton-raphson, secant dan bisection!

Sebelum mencari akar persamaan dengan berbagai metode, kita definisikan dahulu fungsi nya ke dalam berbagai bentuk persamaan yang sesuai untuk masing-masing metode yang digunakan

# mendefinisikan persamaan dalam bentuk fungsi untuk fixed-iteration
fx_fi <- function(x){
  5/(x^2-2)
}

# mendefinikan persamaan dalam bentuk fungsi untuk newton-raphson
fx_nr <- function(x) {
  rumus_fungsi <- function(x) x^3-2*x-5
  fungsi <- rumus_fungsi(x)
  # mencari turunan dengan menggunakan fungsi Deriv
  # dari package Deriv
  rumus_turunan <- Deriv::Deriv(fx)
  turunan <- rumus_turunan(x)
  return(c("fungsi"=fungsi,"turunan"=turunan))
}

Metode Fixed Point

# fixed-point
fixedpoint(ftn = fx_fi,x0=0)
## At iteration 1 value of x is: -2.5 
## At iteration 2 value of x is: 1.176471 
## At iteration 3 value of x is: -8.117978 
## At iteration 4 value of x is: 0.07824535 
## At iteration 5 value of x is: -2.507676 
## At iteration 6 value of x is: 1.165925 
## At iteration 7 value of x is: -7.804949 
## At iteration 8 value of x is: 0.08486483 
## At iteration 9 value of x is: -2.509035 
## At iteration 10 value of x is: 1.164075 
## At iteration 11 value of x is: -7.752778 
## At iteration 12 value of x is: 0.08605028 
## At iteration 13 value of x is: -2.50929 
## At iteration 14 value of x is: 1.163728 
## At iteration 15 value of x is: -7.743083 
## At iteration 16 value of x is: 0.08627333 
## At iteration 17 value of x is: -2.509339 
## At iteration 18 value of x is: 1.163662 
## At iteration 19 value of x is: -7.741248 
## At iteration 20 value of x is: 0.08631566 
## At iteration 21 value of x is: -2.509348 
## At iteration 22 value of x is: 1.16365 
## At iteration 23 value of x is: -7.740899 
## At iteration 24 value of x is: 0.08632371 
## At iteration 25 value of x is: -2.50935 
## At iteration 26 value of x is: 1.163647 
## At iteration 27 value of x is: -7.740833 
## At iteration 28 value of x is: 0.08632524 
## At iteration 29 value of x is: -2.50935 
## At iteration 30 value of x is: 1.163647 
## At iteration 31 value of x is: -7.74082 
## At iteration 32 value of x is: 0.08632553 
## At iteration 33 value of x is: -2.50935 
## At iteration 34 value of x is: 1.163647 
## At iteration 35 value of x is: -7.740818 
## At iteration 36 value of x is: 0.08632558 
## At iteration 37 value of x is: -2.50935 
## At iteration 38 value of x is: 1.163647 
## At iteration 39 value of x is: -7.740817 
## At iteration 40 value of x is: 0.08632559 
## At iteration 41 value of x is: -2.50935 
## At iteration 42 value of x is: 1.163647 
## At iteration 43 value of x is: -7.740817 
## At iteration 44 value of x is: 0.0863256 
## At iteration 45 value of x is: -2.50935 
## At iteration 46 value of x is: 1.163647 
## At iteration 47 value of x is: -7.740817 
## At iteration 48 value of x is: 0.0863256 
## At iteration 49 value of x is: -2.50935 
## At iteration 50 value of x is: 1.163647 
## At iteration 51 value of x is: -7.740817 
## At iteration 52 value of x is: 0.0863256 
## At iteration 53 value of x is: -2.50935 
## At iteration 54 value of x is: 1.163647 
## At iteration 55 value of x is: -7.740817 
## At iteration 56 value of x is: 0.0863256 
## At iteration 57 value of x is: -2.50935 
## At iteration 58 value of x is: 1.163647 
## At iteration 59 value of x is: -7.740817 
## At iteration 60 value of x is: 0.0863256 
## At iteration 61 value of x is: -2.50935 
## At iteration 62 value of x is: 1.163647 
## At iteration 63 value of x is: -7.740817 
## At iteration 64 value of x is: 0.0863256 
## At iteration 65 value of x is: -2.50935 
## At iteration 66 value of x is: 1.163647 
## At iteration 67 value of x is: -7.740817 
## At iteration 68 value of x is: 0.0863256 
## At iteration 69 value of x is: -2.50935 
## At iteration 70 value of x is: 1.163647 
## At iteration 71 value of x is: -7.740817 
## At iteration 72 value of x is: 0.0863256 
## At iteration 73 value of x is: -2.50935 
## At iteration 74 value of x is: 1.163647 
## At iteration 75 value of x is: -7.740817 
## At iteration 76 value of x is: 0.0863256 
## At iteration 77 value of x is: -2.50935 
## At iteration 78 value of x is: 1.163647 
## At iteration 79 value of x is: -7.740817 
## At iteration 80 value of x is: 0.0863256 
## At iteration 81 value of x is: -2.50935 
## At iteration 82 value of x is: 1.163647 
## At iteration 83 value of x is: -7.740817 
## At iteration 84 value of x is: 0.0863256 
## At iteration 85 value of x is: -2.50935 
## At iteration 86 value of x is: 1.163647 
## At iteration 87 value of x is: -7.740817 
## At iteration 88 value of x is: 0.0863256 
## At iteration 89 value of x is: -2.50935 
## At iteration 90 value of x is: 1.163647 
## At iteration 91 value of x is: -7.740817 
## At iteration 92 value of x is: 0.0863256 
## At iteration 93 value of x is: -2.50935 
## At iteration 94 value of x is: 1.163647 
## At iteration 95 value of x is: -7.740817 
## At iteration 96 value of x is: 0.0863256 
## At iteration 97 value of x is: -2.50935 
## At iteration 98 value of x is: 1.163647 
## At iteration 99 value of x is: -7.740817 
## At iteration 100 value of x is: 0.0863256 
## Algorithm failed to converge
## NULL

Metode Newton-Raphson

# newton raphson
newtonraphson(ftn=fx_nr,x0=0)
## At iteration 1 value of x is: -2.5 
## At iteration 2 value of x is: -1.567164 
## At iteration 3 value of x is: -0.5025924 
## At iteration 4 value of x is: -3.820706 
## At iteration 5 value of x is: -2.549393 
## At iteration 6 value of x is: -1.608111 
## At iteration 7 value of x is: -0.5761004 
## At iteration 8 value of x is: -4.59771 
## At iteration 9 value of x is: -3.083543 
## At iteration 10 value of x is: -2.022194 
## At iteration 11 value of x is: -1.123764 
## At iteration 12 value of x is: 1.208652 
## At iteration 13 value of x is: 3.58079 
## At iteration 14 value of x is: 2.655233 
## At iteration 15 value of x is: 2.216106 
## At iteration 16 value of x is: 2.102125 
## At iteration 17 value of x is: 2.094584 
## At iteration 18 value of x is: 2.094551 
## Algorithm converged
##   fungsi 
## 2.094551

Metode Secant

# secant
secant(fx,x0=1,x1=3)
## At iteration 1 value of x is: 1.545455 
## At iteration 2 value of x is: 1.859163 
## At iteration 3 value of x is: 2.20035 
## At iteration 4 value of x is: 2.079799 
## At iteration 5 value of x is: 2.093704 
## At iteration 6 value of x is: 2.094559 
## At iteration 7 value of x is: 2.094551 
## At iteration 8 value of x is: 2.094551 
## At iteration 9 value of x is: 2.094551 
## Algorithm converged
## [1] 2.094551

Metode Bisection

# bisection
bisection(fx,xl=-3,xr=3)
## At iteration 1 value of x is: 0 
## At iteration 2 value of x is: 1.5 
## At iteration 3 value of x is: 2.25 
## At iteration 4 value of x is: 1.875 
## At iteration 5 value of x is: 2.0625 
## At iteration 6 value of x is: 2.15625 
## At iteration 7 value of x is: 2.109375 
## At iteration 8 value of x is: 2.085938 
## At iteration 9 value of x is: 2.097656 
## At iteration 10 value of x is: 2.091797 
## At iteration 11 value of x is: 2.094727 
## At iteration 12 value of x is: 2.093262 
## At iteration 13 value of x is: 2.093994 
## At iteration 14 value of x is: 2.09436 
## At iteration 15 value of x is: 2.094543 
## At iteration 16 value of x is: 2.094635 
## At iteration 17 value of x is: 2.094589 
## At iteration 18 value of x is: 2.094566 
## At iteration 19 value of x is: 2.094555 
## At iteration 20 value of x is: 2.094549 
## At iteration 21 value of x is: 2.094552 
## At iteration 22 value of x is: 2.094551 
## At iteration 23 value of x is: 2.094551 
## At iteration 24 value of x is: 2.094552 
## At iteration 25 value of x is: 2.094552 
## At iteration 26 value of x is: 2.094551 
## At iteration 27 value of x is: 2.094551
## [1] 2.094551

Contoh 2

Misalkan diberikan suatu persamaan \(\mathrm{x^3-x-1=0}\)

  1. Gambarkan fungsi tersebut!

  2. Carilah akar-akar persamaan dengan menggunakan metode fixed-point, newton-raphson, secant dan bisection!

Jawaban

a. Grafik Fungsi

# mendefinikan persamaan dalam bentuk fungsi
fx <- function(x) {
  x^3 - x - 1
}

curve(fx,xlim=c(-3,3), col='steelblue',lwd=2)
abline(h=0)
abline(v=0)

b. Akar Persamaan

# mendefinisikan persamaan dalam bentuk fungsi untuk fixed-point
fx_fp <- function(x){
  1/(x^2-1)
}

# mendefinikan persamaan dalam bentuk fungsi untuk newton-raphson
fx_nr <- function(x) {
  rumus_fungsi <- function(x) x^3-x-1
  fungsi <- rumus_fungsi(x)
  # mencari turunan dengan menggunakan fungsi Deriv
  # dari package Deriv
  rumus_turunan <- Deriv::Deriv(fx)
  turunan <- rumus_turunan(x)
  return(c("fungsi"=fungsi,"turunan"=turunan))
}

Metode Fixed Point

# fixed-point
fixedpoint(ftn = fx_fp,x0=0)
## At iteration 1 value of x is: -1 
## At iteration 2 value of x is: Inf 
## At iteration 3 value of x is: 0 
## At iteration 4 value of x is: -1 
## At iteration 5 value of x is: Inf 
## At iteration 6 value of x is: 0 
## At iteration 7 value of x is: -1 
## At iteration 8 value of x is: Inf 
## At iteration 9 value of x is: 0 
## At iteration 10 value of x is: -1 
## At iteration 11 value of x is: Inf 
## At iteration 12 value of x is: 0 
## At iteration 13 value of x is: -1 
## At iteration 14 value of x is: Inf 
## At iteration 15 value of x is: 0 
## At iteration 16 value of x is: -1 
## At iteration 17 value of x is: Inf 
## At iteration 18 value of x is: 0 
## At iteration 19 value of x is: -1 
## At iteration 20 value of x is: Inf 
## At iteration 21 value of x is: 0 
## At iteration 22 value of x is: -1 
## At iteration 23 value of x is: Inf 
## At iteration 24 value of x is: 0 
## At iteration 25 value of x is: -1 
## At iteration 26 value of x is: Inf 
## At iteration 27 value of x is: 0 
## At iteration 28 value of x is: -1 
## At iteration 29 value of x is: Inf 
## At iteration 30 value of x is: 0 
## At iteration 31 value of x is: -1 
## At iteration 32 value of x is: Inf 
## At iteration 33 value of x is: 0 
## At iteration 34 value of x is: -1 
## At iteration 35 value of x is: Inf 
## At iteration 36 value of x is: 0 
## At iteration 37 value of x is: -1 
## At iteration 38 value of x is: Inf 
## At iteration 39 value of x is: 0 
## At iteration 40 value of x is: -1 
## At iteration 41 value of x is: Inf 
## At iteration 42 value of x is: 0 
## At iteration 43 value of x is: -1 
## At iteration 44 value of x is: Inf 
## At iteration 45 value of x is: 0 
## At iteration 46 value of x is: -1 
## At iteration 47 value of x is: Inf 
## At iteration 48 value of x is: 0 
## At iteration 49 value of x is: -1 
## At iteration 50 value of x is: Inf 
## At iteration 51 value of x is: 0 
## At iteration 52 value of x is: -1 
## At iteration 53 value of x is: Inf 
## At iteration 54 value of x is: 0 
## At iteration 55 value of x is: -1 
## At iteration 56 value of x is: Inf 
## At iteration 57 value of x is: 0 
## At iteration 58 value of x is: -1 
## At iteration 59 value of x is: Inf 
## At iteration 60 value of x is: 0 
## At iteration 61 value of x is: -1 
## At iteration 62 value of x is: Inf 
## At iteration 63 value of x is: 0 
## At iteration 64 value of x is: -1 
## At iteration 65 value of x is: Inf 
## At iteration 66 value of x is: 0 
## At iteration 67 value of x is: -1 
## At iteration 68 value of x is: Inf 
## At iteration 69 value of x is: 0 
## At iteration 70 value of x is: -1 
## At iteration 71 value of x is: Inf 
## At iteration 72 value of x is: 0 
## At iteration 73 value of x is: -1 
## At iteration 74 value of x is: Inf 
## At iteration 75 value of x is: 0 
## At iteration 76 value of x is: -1 
## At iteration 77 value of x is: Inf 
## At iteration 78 value of x is: 0 
## At iteration 79 value of x is: -1 
## At iteration 80 value of x is: Inf 
## At iteration 81 value of x is: 0 
## At iteration 82 value of x is: -1 
## At iteration 83 value of x is: Inf 
## At iteration 84 value of x is: 0 
## At iteration 85 value of x is: -1 
## At iteration 86 value of x is: Inf 
## At iteration 87 value of x is: 0 
## At iteration 88 value of x is: -1 
## At iteration 89 value of x is: Inf 
## At iteration 90 value of x is: 0 
## At iteration 91 value of x is: -1 
## At iteration 92 value of x is: Inf 
## At iteration 93 value of x is: 0 
## At iteration 94 value of x is: -1 
## At iteration 95 value of x is: Inf 
## At iteration 96 value of x is: 0 
## At iteration 97 value of x is: -1 
## At iteration 98 value of x is: Inf 
## At iteration 99 value of x is: 0 
## At iteration 100 value of x is: -1 
## Algorithm failed to converge
## NULL

Metode Newton-Raphson

newtonraphson(ftn=fx_nr,x0=0)
## At iteration 1 value of x is: -1 
## At iteration 2 value of x is: -0.5 
## At iteration 3 value of x is: -3 
## At iteration 4 value of x is: -2.038462 
## At iteration 5 value of x is: -1.390282 
## At iteration 6 value of x is: -0.9116119 
## At iteration 7 value of x is: -0.3450285 
## At iteration 8 value of x is: -1.427751 
## At iteration 9 value of x is: -0.9424179 
## At iteration 10 value of x is: -0.4049494 
## At iteration 11 value of x is: -1.706905 
## At iteration 12 value of x is: -1.155756 
## At iteration 13 value of x is: -0.6941918 
## At iteration 14 value of x is: 0.7424943 
## At iteration 15 value of x is: 2.781296 
## At iteration 16 value of x is: 1.982725 
## At iteration 17 value of x is: 1.536927 
## At iteration 18 value of x is: 1.357262 
## At iteration 19 value of x is: 1.325663 
## At iteration 20 value of x is: 1.324719 
## Algorithm converged
##   fungsi 
## 1.324719

Metode Secant

secant(fx,x0=1,x1=3)
## At iteration 1 value of x is: 1.083333 
## At iteration 2 value of x is: 1.148686 
## At iteration 3 value of x is: 1.379925 
## At iteration 4 value of x is: 1.314886 
## At iteration 5 value of x is: 1.324227 
## At iteration 6 value of x is: 1.324722 
## At iteration 7 value of x is: 1.324718 
## At iteration 8 value of x is: 1.324718 
## At iteration 9 value of x is: 1.324718 
## Algorithm converged
## [1] 1.324718

Metode Bisection

bisection(fx,xl=-3,xr=3)
## At iteration 1 value of x is: 0 
## At iteration 2 value of x is: 1.5 
## At iteration 3 value of x is: 0.75 
## At iteration 4 value of x is: 1.125 
## At iteration 5 value of x is: 1.3125 
## At iteration 6 value of x is: 1.40625 
## At iteration 7 value of x is: 1.359375 
## At iteration 8 value of x is: 1.335938 
## At iteration 9 value of x is: 1.324219 
## At iteration 10 value of x is: 1.330078 
## At iteration 11 value of x is: 1.327148 
## At iteration 12 value of x is: 1.325684 
## At iteration 13 value of x is: 1.324951 
## At iteration 14 value of x is: 1.324585 
## At iteration 15 value of x is: 1.324768 
## At iteration 16 value of x is: 1.324677 
## At iteration 17 value of x is: 1.324722 
## At iteration 18 value of x is: 1.324699 
## At iteration 19 value of x is: 1.324711 
## At iteration 20 value of x is: 1.324717 
## At iteration 21 value of x is: 1.324719 
## At iteration 22 value of x is: 1.324718 
## At iteration 23 value of x is: 1.324717 
## At iteration 24 value of x is: 1.324718 
## At iteration 25 value of x is: 1.324718 
## At iteration 26 value of x is: 1.324718 
## At iteration 27 value of x is: 1.324718
## [1] 1.324718

Contoh 3

Lakukan Modifikasi pada fungsi newtonraphson sehingga output dari fungsi tersebut berupa list dengan tiga object:

  1. Object pertama bernama root berisi akar persamaan.

  2. Object kedua bernama x_values berisi semua nilai x yang dicobakan

  3. Object ketiga bernama error berisi nilai dari hasil operasi abs(fx[1]) Gunakan persamaan pada Latihan 1 untuk menujukkan hasil outputnya!

Jawaban

newtonraphson2 <- function(ftn, x0, tol = 1e-5, max.iter = 100) {
    x <- x0
    fx <- ftn(x)
    iter <- 0
    vector <- c()
   while ((abs(fx[1]) > tol) && (iter < max.iter)) {
         x <- x - fx[1]/fx[2]
         fx <- ftn(x)
         iter <- iter + 1
         vector <- c(vector, x)
         }
   # output depends on success of algorithm
   if (abs(fx[1]) > tol) {
       cat("Algorithm failed to converge\n")
       return(NULL)
       } else {
            cat("Algorithm converged\n")
            list.NR <- list("Root" = x, "X_values" = vector, "error" = abs(fx[1])) # output berupa list
            return(list.NR)
            }
     }

Penggunaan dengan contoh soal sebelumnya :

fx_nr <- function(x) {
  rumus_fungsi <- function(x) x^3-x-1
  fungsi <- rumus_fungsi(x)
  # mencari turunan dengan menggunakan fungsi Deriv
  # dari package Deriv
  rumus_turunan <- Deriv::Deriv(fx)
  turunan <- rumus_turunan(x)
  return(c("fungsi"=fungsi,"turunan"=turunan))
}

newtonraphson2(fx_nr, 0)
## Algorithm converged
## $Root
##   fungsi 
## 1.324719 
## 
## $X_values
##     fungsi     fungsi     fungsi     fungsi     fungsi     fungsi     fungsi 
## -1.0000000 -0.5000000 -3.0000000 -2.0384615 -1.3902821 -0.9116119 -0.3450285 
##     fungsi     fungsi     fungsi     fungsi     fungsi     fungsi     fungsi 
## -1.4277507 -0.9424179 -0.4049494 -1.7069046 -1.1557564 -0.6941918  0.7424943 
##     fungsi     fungsi     fungsi     fungsi     fungsi     fungsi 
##  2.7812959  1.9827252  1.5369274  1.3572625  1.3256631  1.3247188 
## 
## $error
## fungsi.fungsi 
##  3.545493e-06

Contoh 4

Selidiki penggunaan fungsi browser() di R. Gunakan fungsi newtonraphson untuk ilustrasi penggunaan. Laporkan hasil penyelidikanmu beserta ilustrasi koding!

Jawaban

Fungsi browser() dalam R berfungsi untuk menginterupsi suatu ekspresi/proses dan memeriksa/menginspeksi lingkungan pada ekspresi/proses tersebut.

Misal jika ilustrasi penggunaan pada fungsi newtonraphson, maka saat fungsi tersebut dijalankan dan fungsi browser dipanggil, maka fungsi browser akan memeriksa lingkungan dari fungsi newtonraphson atau bisa dikatakan melakukan debugging.

Pada syntax newtonraphson, diselipkan fungsi browser() dalam proses iterasi while.

newtonraphson <- function(ftn, x0, tol = 1e-5, max.iter = 100) {
  x <- x0
  fx <- ftn(x)
  iter <- 0
  while ((abs(fx[1]) > tol) && (iter < max.iter)) {
    browser() #Memasukkan fungsi browser()
    x <- x - fx[1]/fx[2]
    fx <- ftn(x)
    iter <- iter + 1
    cat("At iteration", iter, "value of x is:", x, "\n")
  }
  # output bergantung pada kesuksesan algoritma
  if (abs(fx[1]) > tol) {
    cat("Algorithm failed to converge\n")
    return(NULL)
  } else {
    cat("Algorithm converged\n")
    return(x)
  }
}

Contoh 5

Buatlah gambar fungsi \(\mathrm{x^3-x-1}\) dengan menggunakan package ggplot dengan menggunakan fungsi geom_curve!

df.x <- c()
df.y <- c()
for (i in -3:3) {
  x = i
  y = x^3 - x - 1
  df.x <- c(df.x, x)
  df.y <- c(df.y, y)
}
df.5 <- data.frame(df.x, df.y)


library(ggplot2)
ggplot (df.5)+
  geom_curve(aes(x=df.5[1,1], y=df.5[1,2], xend = df.5[4,1], yend = df.5[4,2]), curvature = -0.33, color="red")+
  geom_curve(aes(x=df.5[4,1], y=df.5[4,2], xend = df.5[7,1], yend = df.5[7,2]), curvature = 0.33, color="red")+
  xlab("x") + ylab("y") +
  geom_vline(xintercept = 0, col="grey50")+
  geom_hline(yintercept = 0, col="grey50")

Integral Numerik

Integrasi numerik berarti komputasi suatu integral menggunakan teknik numerik. Komputasi numerik suatu integral peubah tunggal disebut juga kuadratur (quadrature).

Beberapa pendekatan dalam integrasi numerik adalah :

  1. Trapezoidal
  2. Simpson
  3. Gauss (Adaptive) Quadrature
  4. Monte Carlo Integration

Dalam metode Trapezoid atau Simpson, algoritma berdasarkan tahapan :

  1. pembagian interval integrasi menjadi sejumlah sub interval dengan lebar yang sama dan bentuk sederhana berupa segi empat (trapesium)

  2. menghitung luas setiap sub interval

  3. menjumlahkan luas semua sub interval

Pendekatan yang lebih akurat adalah kuadratur Gauss (Gaussian Quadrature) yang membagi integral tidak sama lebarnya berdasarkan bobot (weight) untuk setiap titik potongnya dan menghitung integral sebagai hasil kali bobot dengan nilai fungsi pada setiap titiknya.

Pendekatan tersebut dapat digunakan untuk adaptive quadrature. Nilai integral didekati dengan berbagai sub interval yang lebih kecil; proses ini berlangsung sampai suatu nilai kesalahan (error) tertentu.

R telah menyiapkan fungsi integrasi yaitu dengan fungsi integrate(). Misalkan akan dicari integral dari \(\mathrm{f(x)=e^{-x}cos(x)}\)

f = function(x) exp(-x) * cos(x)

q <- integrate(f, 0, pi)
q
## 0.521607 with absolute error < 7.6e-15

Metode Trapezoidal

Pada metode trapezoidal, dilakukan pengintegralan berbatas \(\mathrm{\int_{a}^{b} f(x)= F(b)-F(a)}\)

Interval [a,b] dibagi n sehingga lebar setiap sub interval adalah \(\mathrm{h=\frac{(b-a)}{n}}\) dan diperoleh sekuen titik-titik a = x0, x1, x2, … , xn = b.

Sehingga \(\mathrm{\int_{a}^{b} f(x) dx = T}\) adalah :

### Program R - Metode Trapezoidal

trapezoid <- function(ftn, a, b, n = 100) {
     h <- (b-a)/n
     x.vec <- seq(a, b, by = h)
     f.vec <- sapply(x.vec, ftn)     # ftn(x.vec)
     T <- h*(f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)
     return(T)
}

Selain menggunakan fungsi tersebut, dapat pula digunakan fungsi dalam R yaitu fungsi integrate() dan fungsi trapz dalam package pracma

Contoh soal - Metode Trapezoidal

  1. Hitunglah \(\mathrm{\int_{0}^{1} 4x^3 dx }\)
fx_trapezoidal <- function(x) return(4*x^3)

trapezoid(fx_trapezoidal, 0, 1, n = 20)
## [1] 1.0025

Jika dilakukan penghitungan dengan fungsi integrate() dan fungsi trapz dalam package pracma hasilnya

integrate(fx_trapezoidal,0,1)
## 1 with absolute error < 1.1e-14
library(pracma)
## Warning: package 'pracma' was built under R version 4.0.5
## 
## Attaching package: 'pracma'
## The following object is masked _by_ '.GlobalEnv':
## 
##     secant
fx_trapezoidal <- function(x) return(4*x^3)
xs = seq(0, 1)
ys = fx_trapezoidal(xs)

trapz(xs, ys)
## [1] 2

Contoh lain penggunaan fungsi trapz untuk menyelesaikan integral \(\mathrm{\int_{0}^{\pi} e^{-x}cos(x) dx }\)

library(pracma)
f = function(x) exp(-x) * cos(x)
xs = seq(0, pi, length.out = 101)
ys = f(xs)

trapz(xs, ys)
## [1] 0.5216928
library(pracma)
f = function(x) exp(-x) * cos(x)
xs = seq(0, pi, length.out = 101)
ys = f(xs)
h = pi/100
ya = (ys[2] - ys[1])
ye = (ys[101] - ys[100])


trapz(xs, ys) - h/12 * (ye - ya)
## [1] 0.521607
  1. Gunakan metode trapezoidal untuk mencari integral \(\mathrm{\int_{0}^{\pi}sin(x)^2dx}\) untuk n=6
# Mendefinisikan fungsi f(x)
f <- function(x){
  sin(x)^2
}

# Menggambar kurva f(x)
curve(f,0,pi)
abline(v=seq(0,pi,length.out = 7),col=2,lty=2)

# Menghitung integral menggunakan trapezoid
trapezoid(f,0,pi,n = 6)
## [1] 1.570796

Metode Simpson

Metode simpson berdasarkan pada konsep jumlah luas trapezium pada setiap sub interval [xi, xi+1]. Interval [a,b] dibagi n dan n harus merupakan bilangan genap.

Fungsi f(x) diduga dengan suatu fungsi polinomial derajat 2 pada setiap sub interval [xi, xi+1]. Pada sub interval h, tetapkan 3 titik yaitu u, v dan w sehingga u < v < w.

Program R - Metode Simpson

# integral fungsi ftn peubah tunggal dengan n sub interval (genap)
# integral dengan batas bawah a dan batas atas b

simpson_n <- function(ftn, a, b, n = 100) {
    n <- max(c(2*(n %/% 2), 4))
    h <- (b-a)/n
    x.vec1 <- seq(a+h, b-h, by = 2*h)       # ganjil
    x.vec2 <- seq(a+2*h, b-2*h, by = 2*h)   # genap
    f.vec1 <- sapply(x.vec1, ftn)   # ganjil
    f.vec2 <- sapply(x.vec2, ftn)   # genap
    S <- h/3*(ftn(a) + ftn(b) + 4*sum(f.vec1) + 2*sum(f.vec2))
    return(S)
}

Contoh Soal - Metode Simpson

  1. Untuk pengerjaan contoh yang sama dengan metode trapezoid, akan dicari integral dari \(\mathrm{\int_{0}^{\pi} e^{-x}cos(x) dx }\) dengan metode simpson dengan batas 0-20 dan 0-100
ftn <- function(x) return(exp(-x) * cos(x))

simpson_n(ftn, 0, pi, 20) 
## [1] 0.5215999
simpson_n(ftn, 0, pi, 100) 
## [1] 0.5216069
  1. Gunakan metode simpson untuk mencari integral \(\mathrm{\int_{0}^{8}\sqrt(x)dx}\) dengan n=4
# Mendefinisikan fungsi f(x)
f2 <- function(x){
  sqrt(x)
}

# Menggambar kurva f(x)
curve(f2,0,8)
abline(v=seq(0,8,length.out = 5),col=2,lty=2)

# Menghitung integral menggunakan metode Simpson
simpson_n(f2,0,8,n = 4)
## [1] 14.85549

Metode Gauss (Adaptive) Quadrature

Pada Gauss kuadratur batasan integral tidak dibagi sama sesuai dengan pembobot (weight; wi) untuk setiap titik potongnya dan menghitung integral sebagai hasil kali pembobot dengan nilai fungsi pada setiap titiknya.

Nilai integral didekati dengan sub interval yang berbeda dan lebih kecil; proses ini berlangsung sampai suatu nilai kesalahan (error) tertentu.

Di dalam R ada fungsi legendre.quadrature.rules untuk integrasi kuadratur ini. Fungsi ini memerlukan library orthopolynom, polynom, dan gaussquad.

\[\mathrm{\int_{a}^{b} f(x) dx \approx \sum_{i=1}^{n}w_if(x_i) }\] ### Contoh soal - Metode Gauss (Adaptive) Quadrature

  1. Hitung dari \(\mathrm{f(x)=x^6}\) pada interval [-1, 1] dengan Legendre node dan pembobot ordo 4 !
library(polynom) 
## 
## Attaching package: 'polynom'
## The following object is masked from 'package:pracma':
## 
##     integral
library(orthopolynom)
## Warning: package 'orthopolynom' was built under R version 4.0.5
library(gaussquad)
## Warning: package 'gaussquad' was built under R version 4.0.5
# Fungsi Legendre
f = function(x) x^6
Lq = legendre.quadrature.rules(4)[[4]]  # Legendre of order 4
xi = Lq$x; wi = Lq$w            # nodes and weights

sum(wi * f(xi))     # quadrature
## [1] 0.2857143
#Fungsi Integrate
f = function(x) x^6
integrate(f,-1,1)
## 0.2857143 with absolute error < 3.2e-15
  1. Gunakan metode four-point Gaussian untuk mengintegralkan \(\mathrm{\int_{0}^{3}exp(-x^2)dx}\) !
# Mendefinisikan fungsi f(x) dan f(t)
f3a <- function(x){
  exp(-x^2)
}
f3b <- function(t){
  (3/2)*exp((-(3/2*(t+1))^2))
}

# Legendre order 4
Lq = legendre.quadrature.rules(4)[[4]]
xi = Lq$x
wi = Lq$w 

# Menghitung integral berdasarkan gaussian quadratur
sum(wi * f3b(xi)) 
## [1] 0.8841359

Metode Monte Carlo Integration

Monte Carlo integration adalah metode statistik berdasarkan penarikan contoh acak

Program R - Metode Monte Carlo

# User-defined function
mc_integral <- function(ftn, a, b,m=1000){
  # Membangkitkan x berdistribusi U(a,b)
  x <- runif(m,a,b)
  # Menghitung rata-rata dari output fungsi
  Gx <- ftn(x)
  Gx_m <- mean(Gx)
  theta.hat <- (b-a)*Gx_m
  return(theta.hat)
}

Monte Carlo dengan batas [0,1]

Misalkan g(x) adalah suatu fungsi dan akan ditentukan \(\mathrm{\int_{a}^{b} g(x) dx}\)

Jika X adalah suatu peubah acak dengan kepekatan f(x), maka nilai harapan peubah acak Y=g(X) adalah

\[\mathrm{E[g(X)]=\int_{-\infty}^{\infty} g(x) dx}\] Sehingga, penduga monte carlo nya

### Monte Caelo dengan batas bukan [0,1] maka,

Sehingga, algoritma Monte Carlo untuk integral \(\mathrm{\theta=\int_{a}^{b}g(x)dx}\) adalah :

  1. Bangkitkan X1, X2, … , Xm dari Uniform(a,b).

  2. Hitung \(\mathrm{\bar{g_m(X)}=\frac{1}{m}\sum_{i=1}^{m}g(X_i)}\)

  3. \(\mathrm{\hat \theta = (b-a)\bar g(X)}\)

Monte Carlo dengan batas tak terhingga

Sehingga dapat dilakukan substitusi

Sehingga penduga CDF Normal Baku

Bagian 1

set.seed(123)
x <- seq(.1, 2.5, length = 10)
m <- 10000
u <- runif(m)
cdf <- numeric(length(x))
for (i in 1:length(x)) {
      g <- x[i] * exp(-(u * x[i])^2 / 2)
      cdf[i] <- mean(g) / sqrt(2 * pi) + 0.5
}


x <- seq(.1, 2.5, length = 10)
Phi <- pnorm(x) # fungsi pnorm


print(round(rbind(x, cdf, Phi), 3))
##     [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
## x   0.10 0.367 0.633 0.900 1.167 1.433 1.700 1.967 2.233 2.500
## cdf 0.54 0.643 0.737 0.816 0.879 0.925 0.957 0.977 0.989 0.996
## Phi 0.54 0.643 0.737 0.816 0.878 0.924 0.955 0.975 0.987 0.994

Bagian 2

set.seed(123)
x <- seq(.1, 2.5, length = 10)
m <- 10000
u <- runif(m)
cdf <- numeric(length(x))
for (i in 1:length(x)) {
  g <- -x[i] * exp(-(u * (-x[i]))^2 / 2)    # Φ(−𝑥)
  cdf[i] <- mean(g) / sqrt(2 * pi) + 0.5
}

round(cdf,3)
##  [1] 0.460 0.357 0.263 0.184 0.121 0.075 0.043 0.023 0.011 0.004

Contoh soal - Metode Monte Carlo

  1. Tentukan Penduga Monte Carlo untuk \(\mathrm{\theta=\int_{0}^{1}e^{-x}dx}\)
m <- 10000
x <- runif(m)   # U(0,1)

theta.hat <- mean(exp(-x))

print(theta.hat)
## [1] 0.6343679
print(1 - exp(-1))
## [1] 0.6321206
m <- 10000
x <- runif(m, min=0, max=1)

theta.hat <- mean(exp(-x))*1

print(theta.hat)
## [1] 0.6319533
print(exp(-0) - exp(-1))
## [1] 0.6321206

Sehingga \(\mathrm{\hat\theta}=0.6355\) dan \(\mathrm{\theta=1-e^{-1}=0.6321}\)

  1. Tentukan penduga MC untuk \(\mathrm{\theta=\int_{2}^{4}e^{-x}}\)
m <- 10000
x <- runif(m, min=2, max=4) # a=2 dan b=4 ; U(a,b)

theta.hat <- mean(exp(-x)) * 2  # (𝑏−𝑎) (𝑔(𝑌))

print(theta.hat)            # penduga
## [1] 0.118435
print(exp(-2) - exp(-4))        # hasil riil
## [1] 0.1170196

Sehingga \(\mathrm{\hat\theta}=0.1160195\) dan \(\mathrm{\theta=1-e^{-1}=0.1170196}\)

  1. Tentukan integral \(\mathrm{\int_{0}^{1/2}arcsin(\sqrt x) dx}\) dan tentukan berapa banyak sampel yang dibutuhkan agar hasil integral ini sesuai dengan hasil exactnya yaitu 1/4 !
# Mendefinisikan fungsi f(x)
f <- function(x){
  asin(sqrt(x))
}

# Menggambar kurva f(x)
curve(f,0,0.5)
abline(v=seq(0,0.5,length.out = 5),col=2,lty=2)

mc_integral(f,0,0.5, m=1000)
## [1] 0.2537197

Sehingga, dibutuhkan m sebanyak 1000 untuk mendapatkan hasil mendekati nilai exact sebesar 0.25.