Tugas Praktikum Keempat STA561 Pemrograman Statistika

Klik disini untuk ke halaman rpubs.

Algoritma Iteratif (Praktikum Pertemuan Kedelapan)

Sorting

Latihan 1

Misalkan saja ada sebuah vektor x={3,1,4,1,5,9,2,6,5}, buat program untuk mengurutkan vektor tersebut dari kecil ke besar, kemudian bandingkan hasilnya dengan fungsi sort! Gunakan algoritme pada slide kuliah ke-16 !

Jawaban:

Ada beberapa algoritme iteratif yang bisa diterapkan untuk melakukan pengurutan:

  • Bubble Sort

  • Selection Sort

  • Insertion Sort

  • Merge Sort

  • Quick Sort

Algoritme yang digunakan kali ini adalah Bubble Sort. 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

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

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

# Mendefinisikan objek baru yang akan digunakan untuk menampilkan hasil sortir
x_sort <- x

# Banyaknya amatan 
nx <- length(x)

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

x_sort
## [1] 1 1 2 3 4 5 5 6 9
sort(x)
## [1] 1 1 2 3 4 5 5 6 9

Latihan 2

Berdasarkan jawaban Latihan 1, buatlah fungsi yang bernama urut yang membutuhkan input vektor x dan outputnya merupakan vektor yang sudah terurut!

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

# output fungsi urut
x <- c(3,1,4,1,5,9,2,6,5)
urut(x)
## [1] 1 1 2 3 4 5 5 6 9

Menggambar Pola

Latihan 3

Buatlah pola seperti gambar di atas menggunakan algoritme iteratif!

# Membuat pola dasar
mat <- matrix(NA,6,6)
mat
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]   NA   NA   NA   NA   NA   NA
## [2,]   NA   NA   NA   NA   NA   NA
## [3,]   NA   NA   NA   NA   NA   NA
## [4,]   NA   NA   NA   NA   NA   NA
## [5,]   NA   NA   NA   NA   NA   NA
## [6,]   NA   NA   NA   NA   NA   NA
# Membuat pola sesuai yang diinginkan
for (i in 1:6) { #baris
  for (j in 1:(6-i+1)) { #kolom
    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

Deret Bilangan

Latihan 4

Susunlah syntax R untuk menghitung penjumlahan deret tak hingga berikut!

z=2−6/5+8/10−10/17+12/26−…

Hint: Cari formula/pola dari deret tersebut terlebih dahulu

Jawaban:

Formula dari deret tersebut adalah

# 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
# Banyaknya Iterasi
n
## [1] 200002

Latihan 5

Susunlah syntax R untuk menghitung penjumlahan deret tak hingga berikut!

z=2/5+3/40+6/135+11/320+18/625+27/1080+…

Hint: Cari formula/pola dari deret tersebut terlebih dahulu

Jawaban:

Formula dari deret tersebut adalah

# 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

Latihan 6

Berdasarkan jawaban Latihan 4 dan 5, buatlah 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.

deret <- function(x,n){
  
  if(missing(n)) {
  #Input Fungsi Deret Bilangan (1 input)
      # 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
}
output <- matrix(c(z,n),nrow =2, ncol =1)
rownames(output)<- c("jumlah deret","banyaknya iterasi")
return(output)
}
else{
  output <- matrix(c(n,fz2(n)),nrow =2, ncol =1)
  rownames(output)<- c("n","suku ke-n")
  return(output)
}
}
#Fungsi Latihan 4
# Menuliskan rumus fungsi dalam bentuk function
fz2 <- function(n){
  (-1)^(n+1)*((2*n+2)/((n^2)+1))
}

deret(fz2)
##                           [,1]
## jumlah deret      1.267192e+00
## banyaknya iterasi 2.000002e+06
deret(fz2,4)
##                 [,1]
## n          4.0000000
## suku ke-n -0.5882353
#Fungsi Latihan 5
fz2 <- function(n){
  ((n^2) - 2*n + 3) / (5*n^3)
}

deret(fz2)
##                          [,1]
## jumlah deret      2.61992e+00
## banyaknya iterasi 2.00000e+05
deret(fz2,4)
##               [,1]
## n         4.000000
## suku ke-n 0.034375

Bilangan Fibonacci

Latihan 7

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

Metode Roof Finding (Praktikum Pertemuan Kesembilan)

Berikut adalah user-defined function dari metode-metode Root Finding dengan mengikuti algoritme masing-masing metode.

Metode Fixed Point

Fungsi fixedpoint memiliki beberapa argumen yaitu ftn yang berupa function di R, x0 merupakan nilai awal, tol merupakan nilai tolerasi dan max.iter adalah iterasi maksimum. tol dan max.iter masing-masing memiliki nilai default 1e-9 dan 100.

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

Fungsi newtonraphson memiliki beberapa argumen yaitu ftn yang berupa function di R, x0 merupakan nilai awal, tol merupakan nilai tolerasi dan max.iter adalah iterasi maksimum. tol dan max.iter masing-masing memiliki nilai default 1e-9 dan 100.

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

Fungsi secant memiliki beberapa argumen yaitu ftn yang berupa function di R, x0 dan x1 merupakan nilai awal, tol merupakan nilai tolerasi dan max.iter adalah iterasi maksimum. tol dan max.iter masing-masing memiliki nilai default 1e-9 dan 100.

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

Fungsi bisection memiliki beberapa argumen yaitu ftn yang berupa function di R, xl dan xr merupakan nilai awal, tol merupakan nilai tolerasi dan max.iter adalah iterasi maksimum. tol dan max.iter masing-masing memiliki nilai default 1e-9 dan 100.

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

Misalkan diberikan suatu persamaan x^3−2x−5=0.

Gambarkan fungsi tersebut!

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

Jawaban

# 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))
}
  1. Gambarkan fungsi tersebut!
curve(fx,xlim=c(-3,3), col='steelblue',lwd=2)
abline(h=0)
abline(v=0)

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

# 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

Latihan 2

Misalkan diberikan suatu persamaan x^3−x−1=0. a. Gambarkan fungsi tersebut!

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

Jawaban

# 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))
}
  1. Gambarkan fungsi tersebut!
curve(fx,xlim=c(-3,3), col='steelblue',lwd=2)
abline(h=0)
abline(v=0)

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

# 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
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
secant(fx,x0=1,x1=3)
## At iteration 1 value of x is: 1.083333 
## 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: 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

Latihan 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!

newtonraphson.modif <- 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.newtonraphson.modif <- list("Root" = x, "X_values" = vector, "error" = abs(fx[1])) # output berupa list
            return(list.newtonraphson.modif)
            }
     }
#Menentukan akar ciri dari persamaan  1

ftn <- function(x) {
# returns function value and its derivative at x
     fx <- x^3 - 2*x - 5  # fungsi f(x)
     dfx <- 3*x^2 - 2    # turunan pertama f(x)
     return(c(fx, dfx))
}
newtonraphson.modif(ftn, 0)
## Algorithm converged
## $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] 6.472653e-09

Latihan 4

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

Jawaban:

Fungsi \(browser()\) dapat digunakan seperti breakpoint yaitu eksekusi kode akan berhenti pada titik yang dipanggil. Kemudian pengguna dapat memeriksa nilai variabel, mengeksekusi kode R arbitrer, dan menelusuri kode baris demi baris. Jadi, fungsi browser() ini digunakan untuk melakukan debugging.

Setelah fungsi \(browser()\) masuk ke dalam kode, interpreter interaktif akan mulai. Kode R apa pun dapat dijalankan seperti biasa, seperti contoh akan digunakan fungsi browser() untuk melihat program fungsi newtonraphson.

Persamaan yang digunakan untuk melakukan debugging adalah:

dengan nilai x0 yang dimasukkan adalah 3 dan nilai tol diisi dengan 1e-04.

newtonraphson <- function(ftn, x0, tol = 1e-05, max.iter = 100) {
  x <- x0
  fx <- ftn(x)
  iter <- 0
  while ((abs(fx[1]) > tol) && (iter < max.iter)) {
    browser() #ditambah 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)
  }
}
ftn <- function(x) {
    fx <- log(x) - exp(-x) # fungsi f(x)
    dfx <- 1/x + exp(-x) # turunan pertama f(x)
return(c(fx, dfx))
}
newtonraphson(ftn,x0=3,1e-04)
## Called from: newtonraphson(ftn, x0 = 3, 1e-04)
## 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: 0.2624136 
## debug at <text>#5: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#6: browser()
## 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: 0.7224658 
## debug at <text>#5: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#6: browser()
## 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: 1.156032 
## debug at <text>#5: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#6: browser()
## 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: 1.299908 
## debug at <text>#5: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## debug at <text>#6: browser()
## 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: 1.309759 
## debug at <text>#5: (while) (abs(fx[1]) > tol) && (iter < max.iter)
## Algorithm converged
## [1] 1.309759

Dengan menambahkan fungsi browser() dalam syntax maka dapat dilihat jalannya program tersebut baris per baris.

Ketika program dijalankan dengan argumen yang dimasukkan x0 = 3, didapati nilai x sudah terisi dengan nilai x0 yaitu 3, nilai f(x) juga sudah terisi sesuai dengan nilai f(X) pada saat x=3.

Selanjutnya akan dijalankan syntax pada baris ke-7 yang akan mengganti nilai x sebelumnya.

Maka nilai x yang sebelumnya berisi 3 berganti menjadi 0.26

Proses dapat terus dijalankan sampai debugging selesai dilakukan, sehingga dapat dilakukan pemeriksaan apakah syntax yang kita buat pada setiap baris sudah sesuai hasilnya dengan yang diinginkan.

Latihan 5

Buatlah gambar fungsi pada Latihan 2 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 (Praktikum Pertemuan Kesepuluh)

Packages

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

Metode Trapezoidal

Metode trapezoidal bisa menggunakan user-defined function dan fungsi trapzfun dari package pracma.

Memerlukan 2 titik.

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

Metode Simpson

Memerlukan 3 titik.

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

Metode Gauss (Adaptive) Quadrature

Metode Gauss (Adaptive) Quadrature bisa menggunakan fungsi integrate dan fungsi legendre.quadrature.rules dari package gaussquad.

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

Latihan 1

Gunakan metode trapezoidal untuk mengintegralkan fungsi berikut:

untuk n=6 dengan menggunakan Cara Manual dan Program di R

Jawaban

Cara Manual

Diketahui

Program di R

# 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
# Jika menggunakan n lainnya
trapezoid(f,0,pi,n = 1)
## [1] 7.066985e-32
# Menghitung integral menggunakan trapzfun
trapzfun(f,0,pi)
## $value
## [1] 1.570796
## 
## $iter
## [1] 2
## 
## $rel.err
## [1] 0

Latihan 2

Gunakan metode Simpson untuk mengintegralkan fungsi berikut

  1. Cara Manual

Panjang setiap sub interval adalah sebagai berikut:

  1. Program di R
# 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
# Jika menggunakan n lainnya
simpson_n(f2,0,8,n = 15)
## [1] 15.04988

Latihan 3

  1. Program di R
# 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

Latihan 4

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.

ftn <- function(x) {
  asin(sqrt(x))
}

curve(ftn,0,0.5)

round(mc_integral(ftn, 0, 0.5,m=1000),4) #banyak sampel 1.000
## [1] 0.255
round(mc_integral(ftn, 0, 0.5,m=100000),4) #banyak sampel 100.000
## [1] 0.2503
round(mc_integral(ftn, 0, 0.5,m=1000000),4) #banyak sampel 1.000.000
## [1] 0.25
hasil.exact<- (1/4) 
round(hasil.exact,4)
## [1] 0.25

Banyak sampel yang dibutuhkan agar hasil integral ini sesuai dengan hasil exactnya yaitu 1/4 adalah 1.000.000

TERIMA KASIH

Referensi

Dito, Gerry Alfa dan Raharjo, Mulianto. 2021. Algoritma Iteratif. Retrieved From https://newlms.ipb.ac.id

Dito, Gerry Alfa dan Raharjo, Mulianto. 2021. Root Finding. Retrieved From https://newlms.ipb.ac.id

Dito, Gerry Alfa dan Raharjo, Mulianto. 2021. Integral Numerik. Retrieved From https://newlms.ipb.ac.id

Gilat, Amos and Subramaniam, Vish. 2013. Numerical Methods for Engineers and Scientists: An Introduction with Applications Using MATLAB 3rd Ed. Wiley

Wigena, Aji Hamim. 2021. STA561 Pemrograman Statistika: Algoritma Iteratif. Retrieved from https://newlms.ipb.ac.id/

Wigena, Aji Hamim. 2021. STA561 Pemrograman Statistika: Penentuan Akar Persamaan (Root Finding). Retrieved from https://newlms.ipb.ac.id/

Wigena, Aji Hamim. 2021. STA561 Pemrograman Statistika: Integrasi Numerik (Numerical Integration). Retrieved from https://newlms.ipb.ac.id/

  1. Tugas Praktikum 4 STA561 Pemrograman Statistika, Reni Amelia, Mahasiswa Pascasarjana Statistika dan Sains Data, IPB University, ↩︎