Iterative Algorithms

Algoritma Iteratif

Algoritma buble short

Misalkan ada objek baru yang berisi data tidak terurut c(5,9,1,5,4,2,1,4,1) selanjutnya akan dilakukan pengurutan dengan buble short. Untuk melakukan pengurutan, iterasi yang akan digunakan terdiri dari dua tingkat:

  1. Iterasi untuk memastikan semua elemen diperiksa

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

Program R yang digunakan adalah sebagai berikut:

x <-  c(5,9,1,5,4,2,1,4,1)

# Mendefinisikan objek baru yang akan digunakan untuk menampilkan hasil ururt
x_urut <- x

# Banyaknya amatan 
panjangx <- length(x)

# Iterasi tingkat 1
for (i in 1:(panjangx-1)) {
  # Iterasi tingkat 2
  for (j in 1:(panjangx-i)) {
    # Pemeriksaan kondisi apakah elemen ke-j lebih besar dari elemen ke-j+1
    if (x_urut[j] > x_urut[j+1]) {
      # Jika terpenuhi, lakukan algoritme penukaran
      tmp <- x_urut[j]
      x_urut[j] <- x_urut[j+1]
      x_urut[j+1] <- tmp
    }
  }
}

x_urut
## [1] 1 1 1 2 4 4 5 5 9

Agar memudahkan dan dapat menggunakan Program R tersebut kembali, maka dapat disimpan dalam bentuk function

urut=function(x){x_urut <- x

# Banyaknya amatan 
panjangx <- length(x)
for (i in 1:(panjangx-1)) {
  # Iterasi tingkat 2
  for (j in 1:(panjangx-i)) {
    # Pemeriksaan kondisi apakah elemen ke-j lebih besar dari elemen ke-j+1
    if (x_urut[j] > x_urut[j+1]) {
      # Jika terpenuhi, lakukan algoritme penukaran
      tmp <- x_urut[j]
      x_urut[j] <- x_urut[j+1]
      x_urut[j+1] <- tmp
    }
  }
}
return(x_urut)
}
urut(x)
## [1] 1 1 1 2 4 4 5 5 9

Menggambar Pola

Misalkan terdapat matriks berukuran 8x8 diinginkan sisi sebelah kanan atas membentuk segitiga dengan angka 1. Program R yang digunakan:

# Membuat pola dasar
mat <- matrix(NA,8,8)
mat
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]   NA   NA   NA   NA   NA   NA   NA   NA
## [2,]   NA   NA   NA   NA   NA   NA   NA   NA
## [3,]   NA   NA   NA   NA   NA   NA   NA   NA
## [4,]   NA   NA   NA   NA   NA   NA   NA   NA
## [5,]   NA   NA   NA   NA   NA   NA   NA   NA
## [6,]   NA   NA   NA   NA   NA   NA   NA   NA
## [7,]   NA   NA   NA   NA   NA   NA   NA   NA
## [8,]   NA   NA   NA   NA   NA   NA   NA   NA
# Membuat pola sesuai yang diinginkan
for (i in 1:8) {
  for (j in 1:i) {
    mat[j,i] <- 1
  }
}

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

Deret Bilangan

Misalkan terdapat deret dengan x = 2 -6/5 + 8/10 - 10/17 + 12/26 … sehingga program R yang digunakan untuk menghitung deret tsb:

# Menuliskan rumus fungsi dalam bentuk function
fz1 <- function(n){
  (-1)^(n+1)*((2*n+2)/((n^2)+1))
}

# 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>0.00001) {
  z <- z0 + fz1(n)
  error <- abs(z-z0)
  z0 <- z
  n <- n+1
}

# Hasil Akhir
z 
## [1] 1.267197
n
## [1] 200002

Deret Bilangan dengan fungsi bernama deret

fungsi yang bernama deret dengan ketentuan:

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

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

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

# Menuliskan rumus fungsi dalam bentuk function
deret <- function(fz1,n=0){
  # Mendefinisikan nilai-nilai awal (initial value)
  error <- 10 #Nilai awal kriteria berhenti
  z0 <- 0 #Nilai awal deret
  iterasi <- 1 #Nilai awal iterasi
  hasil = NULL
  # Membuat iterasi untuk mencari jumlah deret
  while (error>0.00001) {
    z <- z0 + fz1(iterasi)
    error <- abs(z-z0)
    z0 <- z
    iterasi <- iterasi+1
    hasil = c(hasil,z)
  }
  if (n==0) {
    return(list(`Hasil Akhir`=z,`Banyaknya Iterasi`=iterasi))
  } else{return(list(`Suku ke`=n,`Nilai`=hasil[n]))}
}
deret(fz2,5) #inputnya fungsi deret bilangan dan bilangan n
## $`Suku ke`
## [1] 5
## 
## $Nilai
## [1] 0.5826194
deret(fz2) #inputnya hanya fungsi deret bilangan
## $`Hasil Akhir`
## [1] 2.159406
## 
## $`Banyaknya Iterasi`
## [1] 19999

Bilangan Fibonacci

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

# 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 3
fibo[fibo %% 3 == 0]
## [1]      0      3     21    144    987   6765  46368 317811
# menjumlahkan bilangan Fibonacci kelipatan 3
sum(fibo[fibo %% 3 == 0])
## [1] 372099

Root Finding

Metode Fixed Point

Metode fiexd-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 , untuk menentukan dugaan berikutnya, yaitu x1=f(x0), dan terus berulang. Proses iterartif relatif lambat.

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

Proses konvergen lebih cepat karena nilai awal iterasi dapat diduga untuk memperoleh nilai pendekatan terhadap nilai sebenarnya. Sekuen nilai-nilai x0, x1, x2, x3, . . . diperoleh secara iteratif yang konvergen terhadap nilai akar yang sebenarnya.

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

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

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)
  }
}

Latihan 1

Gambarkan fungsi x^3−2x−5=0!

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

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

# 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))
}

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

# 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
# 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
# secant
secant(fx,x0=1,x1=3)
## At iteration 1 value of x is: 1.545455 
## Algorithm failed to converge
## NULL
# 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

Berdasarkan hasil dari keempat metode tersebut, metode yang tidak dapat menghasilkan solusi konvergen adalah metode fixed-point sedangkan tiga metode lainnya, akar yang diperoleh adalah 2.094551 dengan iterasi terbanyak ada pada metode bisection. Metode bisection membutuhkan waktu paling lama dibandingkan metode newton raphson dan secant sedangkan metode secant merupakan metode yang paling cepat untuk mendapatkan solusi konvergen karena iterasi yang dibutuhkan paling sedikit diantara metode newton raphson dan bisection.

Latihan 2

Misalkan diberikan suatu persamaan x^3−x−1=0

Gambarkan fungsi tersebut!

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

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

# mendefinisikan persamaan dalam bentuk fungsi untuk fixed-iteration
fx_fi <- 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))
}

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

# fixed-point
fixedpoint(ftn = fx_fi,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
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
secant(fx,x0=1,x1=3)
## At iteration 1 value of x is: 1.083333 
## Algorithm failed to converge
## NULL
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

Berdasarkan hasil dari keempat metode tersebut, metode yang tidak dapat menghasilkan solusi konvergen adalah metode fixed-point sedangkan tiga metode lainnya, akar yang diperoleh adalah 1.324718 dengan iterasi terbanyak ada pada metode bisection. Metode bisection membutuhkan waktu paling lama dibandingkan metode newton raphson dan secant sedangkan metode secant merupakan metode yang paling cepat untuk mendapatkan solusi konvergen karena iterasi yang dibutuhkan paling sedikit diantara metode newton raphson dan bisection.

Latihan 3

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

Object pertama bernama root berisi akar persamaan.

Object kedua bernama x_values berisi semua nilai x yang dicobakan

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

Program R yang digunakan untuk menjawab soal tersebut adalah sebagai berikut:

newtonraphson <- function(ftn, x0, tol = 1e-5, max.iter = 100) {
  x <- x0
  x_values <- NULL
  error <- NULL
  fx <- ftn(x)
  iter <- 0
  while ((abs(fx[1]) > tol) && (iter < max.iter)) {
    error <- c(error,abs(fx[1]))
    x <- x - fx[1]/fx[2]
    x_values <- c(x_values,x)
    fx <- ftn(x)
    iter <- iter + 1
  }
  return(list(root = x, x_values = x_values, error = error))
}

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

#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, turunan))
}

newtonraphson(ftn=fx_nr,x0=0)
## $root
## [1] 2.094551
## 
## $x_values
##  [1] -2.5000000 -1.5671642 -0.5025924 -3.8207065 -2.5493934 -1.6081115
##  [7] -0.5761004 -4.5977096 -3.0835431 -2.0221943 -1.1237641  1.2086516
## [13]  3.5807900  2.6552332  2.2161063  2.1021250  2.0945836  2.0945515
## 
## $error
##  [1] 5.000000e+00 1.562500e+01 5.714632e+00 4.121770e+00 5.313249e+01
##  [6] 1.647076e+01 5.942390e+00 4.039002e+00 9.299526e+01 2.815198e+01
## [11] 9.224909e+00 4.171613e+00 5.651658e+00 3.375152e+01 8.409627e+00
## [16] 1.451367e+00 8.489238e-02 3.582354e-04

Berdasarkan hasil tersebut, diperoleh solusi akar 2.094551 dengan 18 iterasi dan masing-masing nilai x yang dicobakan dan eror adalah seperti di atas.

Latihan 4

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

Program R yang digunakan untuk menjawab soal tersebut adalah sebagai berikut:

newtonraphson_b <- function(ftn, x0, tol = 1e-5, max.iter = 100) {
  browser()
  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)
  }
}
#mendefinikan persamaan dalam bentuk fungsi
fx <- function(x) {
  x^3 - 2 * x - 5
}
#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))
}
# newton raphson
newtonraphson_b(ftn=fx_nr,x0=0)
## Called from: newtonraphson_b(ftn = fx_nr, x0 = 0)
## debug at <text>#3: x <- x0
## debug at <text>#4: fx <- ftn(x)
## debug at <text>#5: iter <- 0
## debug at <text>#6: 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")
## }
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 1 value of x is: -2.5 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 2 value of x is: -1.567164 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 3 value of x is: -0.5025924 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 4 value of x is: -3.820706 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 5 value of x is: -2.549393 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 6 value of x is: -1.608111 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 7 value of x is: -0.5761004 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 8 value of x is: -4.59771 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 9 value of x is: -3.083543 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 10 value of x is: -2.022194 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 11 value of x is: -1.123764 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 12 value of x is: 1.208652 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 13 value of x is: 3.58079 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 14 value of x is: 2.655233 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 15 value of x is: 2.216106 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 16 value of x is: 2.102125 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 17 value of x is: 2.094584 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#7: x <- x - fx[1]/fx[2]
## debug at <text>#8: fx <- ftn(x)
## debug at <text>#9: iter <- iter + 1
## debug at <text>#10: cat("At iteration", iter, "value of x is:", x, "\n")
## At iteration 18 value of x is: 2.094551 
## debug at <text>#6: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#13: if (abs(fx[1]) > tol) {
##     cat("Algorithm failed to converge\n")
##     return(NULL)
## } else {
##     cat("Algorithm converged\n")
##     return(x)
## }
## debug at <text>#17: cat("Algorithm converged\n")
## Algorithm converged
## debug at <text>#18: return(x)
##   fungsi 
## 2.094551

Latihan 5

Buatlah gambar fungsi pada Latihan 2 dengan menggunakan package ggplot dengan menggunakan fungsi geom_curve!

Program R yang digunakan untuk menjawab soal tersebut adalah sebagai berikut:

library(ggplot2)
fx <- function(x) {
  x^3 - x - 1
}

x <- -3
xend <- 3
ggplot() + 
  geom_curve(aes(x = x, y = fx(x), xend = xend, yend= fx(xend)), colour = "blue", curvature = 0.4) +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0)

fx <- function(x) {
  x^3 - x - 1
}

ggplot(data.frame(x = c(-3,3)), aes(x)) +
  stat_function(fun = fx, col = "blue") +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0)

Integral Numerik

Integrasi numerik berarti komputasi suatu integral menggunakan teknik numerik. Komputasi numerik suatu integral peubah tunggal disebut juga kuadratur (quadrature). Beberapa pendekatan integrasi numerik, antara lain Trapezoid atau Simpson berdasarkan pembagian interval integrasi menjadi sejumlah sub interval dengan lebar yang sama dan bentuk sederhana berupa segi empat (trapesium), menghitung luas setiap sub interval, dan menjumlahkan luas semua sub interval.

Pertama panggil package yang akan digunakan

library(pracma)
## Warning: package 'pracma' was built under R version 4.1.2
## 
## Attaching package: 'pracma'
## The following object is masked _by_ '.GlobalEnv':
## 
##     secant
library(gaussquad)
## Loading required package: polynom
## 
## Attaching package: 'polynom'
## The following object is masked from 'package:pracma':
## 
##     integral
## Loading required package: orthopolynom

Metode Trapezoidal

# User-defined function
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)
  Trap <- h*(f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)
  return(Trap)
}
# 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)

trapezoid(f,0,pi,n = 6)
## [1] 1.570796
trapezoid(f,0,pi,n = 1)
## [1] 7.067451e-32
# Menghitung integral menggunakan trapzfun
trapzfun(f,0,pi)
## $value
## [1] 1.570796
## 
## $iter
## [1] 2
## 
## $rel.err
## [1] 0

Berdasarkan hasil tersebut, diperoleh luas integral sebesar 1.570796

Metode Simpson

# User-defined function
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)
}

# Mendefinisikan fungsi f(x)
f2 <- function(x){
  sqrt(x)
}

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

simpson_n(f2,0,8,n = 4)
## [1] 14.85549
# Jika menggunakan n lainnya
simpson_n(f2,0,8,n = 15)
## [1] 15.04988

Berdasarkan hasil tersebut, diperoleh luas integral sebesar 14.85549 jika n yang digunakan sebanyak 4 dan jika n=15 luas integral yang diperoleh sebesar 15.04988

Metode Gauss (Adaptive) Quadrature

Pendekatan yang lebih akurat adalah kuadratur Gauss (Gaussian Quadrature) yang membagi integral tidak sama lebarnya sesuai dengan 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.

# 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
# Menghitung integral dengan fungsi integrate
integrate(f3a,0,3)
## 0.8862073 with absolute error < 3.2e-11

Metode Integrasi Mote Carlo (MC

# 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)
}

Latihan

Gunakan metode Monte Carlo untuk mengintegralkan fungsi dan tentukan berapa banyak sampel yang dibutuhkan agar hasil integral ini sesuai dengan hasil exactnya yaitu 1/4! Silahkan kerjakan sebagai latihan. Gunakan user-defined function yang sudah dibuat di atas.

Program R yang digunakan untuk menjawab soal tersebut adalah sebagai berikut:

# Mendefinisikan fungsi f(x)
f <- function(x){
  asin(sqrt(x))
}

mc_integral(f,0,1/2,10000000)
## [1] 0.2500121
f = function(x) exp(-x) * cos(x)
q <- integrate(f, 0, 1/2)

dengan menggunakan sample sebanyak 10.000.000 maka diperoleh solusi sekitar 0.2500, hasil ini sudah dianggap sesuai dengan hasil exactnya yaitu 1/4