Laporan 4 Praktikum Pemrograman Statistik - Tugas Satria June Adwendi

Praktikum Pertemuan 8 - Algoritma Iteratif

Latihan 1

Sorting

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 pada slide kuliah adalah Bubble Sort. Untuk melakukan pengurutan, iterasi yang akan digunakan terdiri dari dua tingkat:

Iterasi untuk memastikan semua elemen diperiksa

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

Latihan 2

Berdasarkan jawaban Latihan 1, buatlah fungsi yang bernama urut dengan 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

Latihan 3

Menggambar Pola

##      [,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

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

Latihan 4

Deret Bilangan

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("hasil akhir","banyaknya iterasi")
return(output)
}
else{
  output <- matrix(c(n,fz2(n)),nrow =2, ncol =1)
  rownames(output)<- c("n","suku ke-n")
  return(output)
}
}


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

deret(fz2)
##                           [,1]
## hasil akhir       1.267192e+00
## banyaknya iterasi 2.000002e+06
deret(fz2,4)
##                 [,1]
## n          4.0000000
## suku ke-n -0.5882353
# Menuliskan rumus fungsi latihan 5 dalam bentuk function
fz2 <- function(n){
  ((n^2) - 2*n + 3) / (5*n^3)
}

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

Latihan 7

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

Praktikum Pertemuan 9 - Root Finding

Metode Root Finding

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

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

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

3. 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-9,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 .

  1. Gambarkan fungsi tersebut!

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

Jawab:

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

# mendefinisikan persamaan dalam bentuk fungsi untuk fixed-point
fx_fp <- 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
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!
# fixed-point
fixedpoint(ftn = fx_fp,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 
## 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
# 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 
## At iteration 28 value of x is: 2.094551 
## At iteration 29 value of x is: 2.094551 
## At iteration 30 value of x is: 2.094551 
## At iteration 31 value of x is: 2.094551 
## At iteration 32 value of x is: 2.094551 
## At iteration 33 value of x is: 2.094551 
## At iteration 34 value of x is: 2.094551
## [1] 2.094551

Latihan 2

Misalkan diberikan suatu persamaan 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

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

# 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))
}
  1. Gambarkan Fungsi
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!
# 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
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 
## 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
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 
## At iteration 28 value of x is: 1.324718 
## At iteration 29 value of x is: 1.324718 
## At iteration 30 value of x is: 1.324718 
## At iteration 31 value of x is: 1.324718 
## At iteration 32 value of x is: 1.324718 
## At iteration 33 value of x is: 1.324718 
## At iteration 34 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 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.

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

Latihan 5

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

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

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.0.5
ggplot (df.lat5)+
  geom_curve(aes(x=df.lat5[1,1], y=df.lat5[1,2], xend = df.lat5[4,1], yend = df.lat5[4,2]), curvature = -0.33, color="blue")+
  geom_curve(aes(x=df.lat5[4,1], y=df.lat5[4,2], xend = df.lat5[7,1], yend = df.lat5[7,2]), curvature = 0.33, color="blue")+
  xlab("x") + ylab("y") +
  geom_vline(xintercept = 0, col="grey50")+
  geom_hline(yintercept = 0, col="grey50")

Praktikum Pertemuan 10 - Integral Numerik

Package

library(pracma)
## Warning: package 'pracma' was built under R version 4.0.5
## 
## 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 Integral Numerik

1. Metode Trapezoidal

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

T=(h/2) (f(x0)+2f(x1)+2f(x2)+…+2f(xn−1)+f(xn)) dengan x0=a,xn=b

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

2. Metode Simpson

S=(h/3) (f(x0)+4f(x1)+2f(x2)+…+4f(xn−1)+f(xn)) dengan x0=a,xn=b

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

3. Metode Gauss (Adaptive) Quadrature

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

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

∫(batas 0,π) sin(x)^2 dx

untuk n=6 dengan menggunakan

  1. Cara Manual

  2. Program di R

Jawaban:

a. Cara Manual

Diketahui

Panjang setiap sub interval adalah sebagai berikut:

Menghitung fungsi f(x) pada setiap xi

Kemudian berdasarkan metode trapezoidal berikut adalah hasil akhirnya Menghitung fungsi f(x) pada setiap xi

Jika kita bandingkan dengan hasil yang didapatkan secara exact adalah sebagai berikut:

b. Program 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

∫(batas 0, 8) (√x) dx

untuk n=4 dengan menggunakan

  1. Cara Manual

  2. Program di R

Jawaban

a. Cara Manual

Panjang setiap sub interval adalah sebagai berikut:

sehinga

Menghitung fungsi untuk setiap titik xi:

Menggantikan semua nilai-nilai diatas ke dalam rumus Simpson

Solusi exact dari integral ini adalah sebagai berikut

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

Gunakan metode four-point gaussian quadratur untuk mengintegralkan fungsi berikut

∫(batas 0, 3) exp(−x^2) dx

a. Cara Manual

Rumus Umum

Langkah pengerjaanya adalah sebagai berikut

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

Gunakan metode Monte Carlo untuk mengintegralkan fungsi berikut

∫(batas 0, 1/2) arcsin(√x) dx

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.2499
round(mc_integral(ftn, 0, 0.5,m=100000),4) #banyak sampel 100.000
## [1] 0.2501
round(mc_integral(ftn, 0, 0.5,m=1000000),4) #banyak sampel 1.000.000
## [1] 0.2499
hasil.exact<- (1/4) 
round(hasil.exact,4)
## [1] 0.25
#Banyak sampel yang dibutuhkan agar hasil integral sesuai dengan hasil exactnya yaitu 1/4 adalah sebanya 1.000.000 sampel

Demikian, Terima Kasih

  1. Satria June Adwendi, IPB University, ↩︎