Tugas Praktikum Metode Numerik (Turunan & Integral Numerik)

rm(list = ls())

# Pertemuan 12: Contoh Metode Pias

cat("Pertemuan 12: Contoh Metode Pias\n")

## Pertemuan 12: Contoh Metode Pias

# Definisi fungsi dan parameter
f <- function(x) { x^2 + x }
a <- 0
b <- 1
n <- 4
h <- (b - a) / n
x <- a

# Nilai sejati
nilai_sejati <- 5/6
cat("Nilai Sejati: 5/6 =", nilai_sejati, "\n\n")

## Nilai Sejati: 5/6 = 0.8333333

# 1.1 Kaidah Segiempat Kiri (Left Riemann Sum)
IP.SEK <- function(f, x, h, n){
  X <- c()
  X[1] <- x
  for (i in 2:(n+1)){
    X[i] <- X[i-1] + h
  }
  h * sum(f(X[1:n]))
}

# 1.2 Kaidah Segiempat Kanan (Right Riemann Sum)
IP.SEka <- function(f, x, h, n){
  X <- c()
  X[1] <- x
  for (i in 2:(n+1)){
    X[i] <- X[i-1] + h
  }
  h * sum(f(X[2:(n+1)]))
}

# 1.3 Kaidah Titik Tengah (Midpoint Rule)
IP.SEC <- function(f, x, h, n){
  n2 <- n * 2
  X <- c()
  X[1] <- x
  for (i in 2:(n2+1)){
    X[i] <- X[i-1] + h/2
  }
  X <- X[c(1:n2) %% 2 == 0]
  h * sum(f(X))
}

# 1.4 Kaidah Trapesium (Trapezoidal Rule) - DIPERBAIKI
IP.T <- function(f, a, b, x, h, n){
  fa <- f(a)
  fb <- f(b)
  X <- seq(from = a, to = b, by = h)
  fc <- 2 * sum(f(X[2:n]))
  (h/2) * (fa + fb + fc)
}

# Hitung nilai-nilai
hasil_segiempat_kiri <- IP.SEK(f, x, h, n)
hasil_segiempat_kanan <- IP.SEka(f, x, h, n)
hasil_titik_tengah <- IP.SEC(f, x, h, n)
hasil_trapesium <- IP.T(f, a, b, x, h, n)

cat("Kaidah Segiempat Kiri  :", hasil_segiempat_kiri, " | Error:", abs(nilai_sejati - hasil_segiempat_kiri), "\n")

## Kaidah Segiempat Kiri  : 0.59375  | Error: 0.2395833

cat("Kaidah Segiempat Kanan :", hasil_segiempat_kanan, " | Error:", abs(nilai_sejati - hasil_segiempat_kanan), "\n")

## Kaidah Segiempat Kanan : 1.09375  | Error: 0.2604167

cat("Kaidah Titik Tengah    :", hasil_titik_tengah, " | Error:", abs(nilai_sejati - hasil_titik_tengah), "\n")

## Kaidah Titik Tengah    : 0.828125  | Error: 0.005208333

cat("Kaidah Trapesium       :", hasil_trapesium, " | Error:", abs(nilai_sejati - hasil_trapesium), "\n")

## Kaidah Trapesium       : 0.84375  | Error: 0.01041667

# Pertemuan 13: Metode Romberg

cat("Pertemuan 13: Metode Romberg\n")

## Pertemuan 13: Metode Romberg

# Fungsi trapesium untuk n pias 
trapezoid <- function(f, a, b, n) {
  h <- (b - a) / n
  x <- seq(a, b, by = h)
  y <- f(x)
  (h/2) * (y[1] + 2 * sum(y[2:n]) + y[n+1])
}

# Definisi fungsi dan parameter
f <- function(x) { 1 / (1 + x) }
a <- 0
b <- 1

# Hitung A0, A1, A2, A3 (kaidah trapesium dengan n=1,2,4,8)
A0 <- trapezoid(f, a, b, 1)   # n=1
A1 <- trapezoid(f, a, b, 2)   # n=2
A2 <- trapezoid(f, a, b, 4)   # n=4
A3 <- trapezoid(f, a, b, 8)   # n=8

cat("A0 (n=1):", A0, "\n")

## A0 (n=1): 2.25

cat("A1 (n=2):", A1, "\n")

## A1 (n=2): 0.7083333

cat("A2 (n=4):", A2, "\n")

## A2 (n=4): 0.6970238

cat("A3 (n=8):", A3, "\n")

## A3 (n=8): 0.6941219

# Ekstrapolasi Richardson
# A (orde h^2) -> B (orde h^4): B_k = A_k + (A_k - A_{k-1})/(2^2 - 1)
B1 <- A1 + (A1 - A0) / (2^2 - 1)
B2 <- A2 + (A2 - A1) / (2^2 - 1)
B3 <- A3 + (A3 - A2) / (2^2 - 1)

cat("\nB1:", B1, "\n")

## 
## B1: 0.1944444

cat("B2:", B2, "\n")

## B2: 0.693254

cat("B3:", B3, "\n")

## B3: 0.6931545

# B (orde h^4) -> C (orde h^6): C_k = B_k + (B_k - B_{k-1})/(2^4 - 1)
C2 <- B2 + (B2 - B1) / (2^4 - 1)
C3 <- B3 + (B3 - B2) / (2^4 - 1)

cat("\nC2:", C2, "\n")

## 
## C2: 0.7265079

cat("C3:", C3, "\n")

## C3: 0.6931479

# C (orde h^6) -> D (orde h^8): D_k = C_k + (C_k - C_{k-1})/(2^6 - 1)
D3 <- C3 + (C3 - C2) / (2^6 - 1)

cat("\nD3:", D3, "\n")

## 
## D3: 0.6926184

# Nilai eksak untuk perbandingan
eksak <- log(2)  # ∫_0^1 1/(1+x) dx = ln(2)
cat("\nNilai eksak (ln(2)):", eksak, "\n")

## 
## Nilai eksak (ln(2)): 0.6931472

cat("Error D3:", abs(eksak - D3), "\n")

## Error D3: 0.0005288034

# Tabel Romberg
cat("\n--- TABEL ROMBERG ---\n")

## 
## --- TABEL ROMBERG ---

tabel <- matrix(NA, 4, 4)
tabel[1, 1] <- A0
tabel[2, 1] <- A1
tabel[2, 2] <- B1
tabel[3, 1] <- A2
tabel[3, 2] <- B2
tabel[3, 3] <- C2
tabel[4, 1] <- A3
tabel[4, 2] <- B3
tabel[4, 3] <- C3
tabel[4, 4] <- D3

colnames(tabel) <- c("k=0", "k=1", "k=2", "k=3")
rownames(tabel) <- c("n=1", "n=2", "n=4", "n=8")
print(round(tabel, 8))

##           k=0       k=1       k=2       k=3
## n=1 2.2500000        NA        NA        NA
## n=2 0.7083333 0.1944444        NA        NA
## n=4 0.6970238 0.6932540 0.7265079        NA
## n=8 0.6941219 0.6931545 0.6931479 0.6926184

# Pertemuan 13: Kuadratur Gauss-Legendre 2-titik

cat("Pertemuan 13: Kuadratur Gauss-Legendre 2-titik\n")

## Pertemuan 13: Kuadratur Gauss-Legendre 2-titik

# Fungsi Gauss-Legendre 2-titik
gauss_legendre_2 <- function(f, a, b) {
  t1 <- 1/sqrt(3)
  t2 <- -1/sqrt(3)
  x1 <- (a+b)/2 + (b-a)/2 * t1
  x2 <- (a+b)/2 + (b-a)/2 * t2
  
  cat("  t1 =", t1, ", t2 =", t2, "\n")
  cat("  x1 =", x1, ", x2 =", x2, "\n")
  cat("  f(x1) =", f(x1), ", f(x2) =", f(x2), "\n")
  
  (b-a)/2 * (f(x1) + f(x2))
}

# Definisi fungsi dan parameter
f <- function(x) { x^2 + 1 }
a <- 1
b <- 2

hasil_gauss <- gauss_legendre_2(f, a, b)

##   t1 = 0.5773503 , t2 = -0.5773503 
##   x1 = 1.788675 , x2 = 1.211325 
##   f(x1) = 4.199359 , f(x2) = 2.467308

cat("Hasil Kuadratur Gauss-Legendre 2-titik:", hasil_gauss, "\n")

## Hasil Kuadratur Gauss-Legendre 2-titik: 3.333333

# Nilai eksak
eksak_gauss <- (b^3/3 + b) - (a^3/3 + a)
cat("Nilai eksak:", eksak_gauss, "\n")

## Nilai eksak: 3.333333

cat("Error:", abs(eksak_gauss - hasil_gauss), "\n")

## Error: 0

Tugas Praktikum Metode Numerik (Turunan & Integral Numerik)

Nabila Putri Amandah

2025-12-01