Komstat Aproksimasi Integral

#APROKSIMASI INTEGRAL

FUNGSI POLINOM-> INTEGRAL X^N DX- 1/(N+1) x^N+1 + c

METODE APROKSIMASI YANG DAPAT DIGUNAKAN UNTUK MEMBERIKAN NILAI HAMPIRAN DARI INTEGRAL TERTENTU. METODE-> RIEMANN, NEWTON-COTES, DAN POMBERG

Secara umum untuk metode aproksimasi integral menggunakan cara mempartisi selang [a,b] menjadi n sub-selang.

𝑎=𝑥_0<𝑥_1<…<𝑥_(𝑛−1)<𝑥_𝑛=𝑏

SUB-SELANG-> BATAS ATAS DAN BATAS BAWAH TERPENUHI

#RIEMANN -> MENGHITUNG LUAS PERSEGI PANJANG

Jumlahan Riemann didefinisikan sebagai: 𝑅=∑(𝑖=0)^(𝑛−1)▒〖𝑓(𝑐_𝑖)(𝑥(𝑖+1)−𝑥_𝑖)〗

Untuk Riemann Kiri ->𝑐_𝑖=𝑥_𝑖-> YANG PALING KIRI DARI SUB-SELANG Untuk Riemann Kanan -> 𝑐_𝑖=𝑥_(𝑖+1)-> YANG PALING KANAN DARI SUB-SELANG Untuk Riemann Tengah -> 𝑐_𝑖=(𝑥_𝑖+𝑥_(𝑖+1))/2 -> MERATA2KANNYA

Dipilih pembagian sub-selang sehingga setiap sub-selang memiliki jarak yang sama: 𝑥_(𝑖+1)−𝑥_𝑖=(𝑏−𝑎)/𝑛

*SELANG-> BATAS SELANG/ JUMLAH PARTISI

riemann<-function(f,a,b,n){
  d<-(b-a)/n
  sum=0
  for(i in 1:n){
    sum=sum+f(a+(i-1)*d)*d
  }
  return(sum)
}

f<-function(x) 1/x^2
riemann(f,1,5,100)

## [1] 0.8194644

integrate(f,1,5)

## 0.8 with absolute error < 1.1e-07

riemannkanan<-function(f,a,b,n){
  d<-((b-a)/n)
  sum=0
  for(i in 1:n){
    sum=sum(f(seq(a, b - d, length = n)) * d)
  }
  return(sum)
}

riemannkanan(f,1,5,100)

## [1] 0.8194644

integrate(f,1,5)

## 0.8 with absolute error < 1.1e-07

#NEWTON COTES

trapezoid<-function(f,a,b){
  return (((b-a)/2)*(f(a)+f(b)))
}
f<-function(x) exp(x)
trapezoid(f,0,3)

## [1] 31.62831

integrate(f,0,3)

## 19.08554 with absolute error < 2.1e-13

Latihan: Buatlah fungsi aproksimasi dari metode Simpson 1/3, Simpson 3/8, dan Boole, kemudian bandingkan hasilnya dengan menggunakan fungsi dasar integrate.(yang di praktikum latihan!)

#simpson 1/3

simpson13=function(f,a,b){
  return(((b-a)/6)*f(a)+4*f((a+b)/2)+f(b))
}
simpson13(f,0,3)

## [1] 38.51229

integrate(f,0,3)

## 19.08554 with absolute error < 2.1e-13

#Simpson 3/8

simpson38=function(f,a,b){
  h=(b-a)/3
  return (((b-a)/8)*(f(a)+3*f(a+h)+3*f(b-h)+f(b)))
}
simpson38(f,0,3)

## [1] 19.27783

integrate(f,0,3)

## 19.08554 with absolute error < 2.1e-13

#Boole

boole <- function(f,a,b){
return ((b-a)/90*(7*f(a)+32*f((3*a+b)/4)+12*f((a+b)/2)+32*f((a+3*b)/4)+7*f(b)))
}
boole(f,0,3)

## [1] 19.09102

integrate(f,0,3)

## 19.08554 with absolute error < 2.1e-13

Keempat metode Newton-Cotes dapat diaplikasikan ke dalam bentuk kelipatan (Komposit). Diterapkan untuk sub-interval hasil dari partisi.

Multiple Trapezoid* ∑(𝑖=0)^(𝑛−1)▒〖(𝑥(𝑖+1)−𝑥_𝑖)/2 (𝑓(𝑥_𝑖 )+𝑓(𝑥_(𝑖+1) )) 〗 dijabarkan menjadi (𝑏−𝑎)/2𝑛 (𝑓(𝑥_0 )+2𝑓(𝑥_1 )+2𝑓(𝑥_2 )+…+2𝑓(𝑥_(𝑛−1) )+𝑓(𝑥_𝑛 ))

MulTrap<-function(f,a,b,n){
  d<-(b-a)/n
  sum=f(a)
  for(i in 1:n-1){
    sum=sum+2*f(a+i*d)
  }
  sum=sum+f(a+n*d)
  return (d*sum/2)
}
f<-function(x) exp(x)
MulTrap(f,0,3,1000)

## [1] 19.08855

integrate(f,0,3)

## 19.08554 with absolute error < 2.1e-13

Untuk metode Simpson-1/3 juga dapat diaplikasikan dalam bentuk kelipatan Syaratnya, jumlah partisi haruslah merupakan kelipatan 2. Jumlahkan bentuk Simpson-1/3 dari indeks genap 2ℎ/6 [𝑓(𝑥_0 )+4𝑓(𝑥_1 )+2𝑓(𝑥_2 )+4𝑓(𝑥_3 )+…+2𝑓(𝑥_(𝑛−1) )+𝑓(𝑥_𝑛 )]

MulSimp1<-function(f,a,b,n){
  d<-(b-a)/n
  sum=f(a)
  for(i in 1:n-1){
    #jika index ganjil
    if(i %% 2 == 1){
      sum=sum+4*f(a+i*d)
    } else {
      #jika indexgenap
      sum=sum+2*f(a+i*d)
    }
  }
    sum=sum+f(a+n*d)
    return (d*sum/3)
}
f<-function(x) 2^x
MulSimp1(f,0,3,1000)

## [1] 10.10087

integrate(f,0,3)

## 10.09887 with absolute error < 1.1e-13

Untuk metode Simpson-3/8 Syaratnya partisinya harus kelipatan 3 Nilai yang digunakan adalah 𝑥_(𝑖+3)−𝑥_𝑖=3ℎ

MulSimp3_8 <- function(f, a, b, n) {
  h <- (b - a) / n
  sum <- f(a) + f(b)

  for (i in 1:(n-1)) {
    xi <- a + i * h

    # Jika i kelipatan 3, gunakan koefisien 2
    if (i %% 3 == 0) {
      sum <- sum + 2 * f(xi)
    } else {
      # Jika i bukan kelipatan 3, gunakan koefisien 3
      sum <- sum + 3 * f(xi)
    }
  }

  result <- (3 * h / 8) * sum
  return(result)
}

# Fungsi uji
f <- function(x) 2^x
MulSimp3_8(f, 0, 3, 1000)

## [1] 10.09287

integrate(f,0,3)

## 10.09887 with absolute error < 1.1e-13

Metode Boole: Syarat partisinya harus kelipatan 4 Nilai yang digunakan adalah 𝑥_(𝑖+4)−𝑥_𝑖=4ℎ

MulBoole <- function(f, a, b, n) {
  h <- (b - a) / n
  sum <- f(a) + f(b)

  for (i in 1:(n-1)) {
    xi <- a + i * h

    # Jika i kelipatan 4, gunakan koefisien 14
    if (i %% 4 == 0) {
      sum <- sum + 14 * f(xi)
    } else {
      # Jika i bukan kelipatan 4 dan i kelipatan 2, gunakan koefisien 32
      if (i %% 2 == 0) {
        sum <- sum + 32 * f(xi)
      } else {
        # Jika i bukan kelipatan 4 dan i bukan kelipatan 2, gunakan koefisien 12
        sum <- sum + 12 * f(xi)
      }
    }
  }

  result <- (2 * h / 45) * sum
  return(result)
}

# Fungsi uji
f <- function(x) 2^x
MulBoole(f, 0, 3, 1000)

## [1] 7.84747

integrate(f,0,3)

## 10.09887 with absolute error < 1.1e-13

#ROMBERG

-> Improvisasi dari metode Newton-Cotes -> Gabungan antara Multiple Trapezoid dan Ekstrapolasi Richardson -> Notasi: 𝑅_(𝑖,𝑗)->diperoleh dari integrasi menggunakan metode Multiple Trapezoid dengan banyaknya partisi adalah 2^(𝑖−1) ->Pada halaman sebelumnya sudah dibuat fungsi MulTrap, jalankan fungsi untuk setiap iterasinya yaitu MulTrap(f,a,b,2^i-1) ->Kemudian, untuk iterasi pada indeks j, digunakan formula sebagai berikut: 𝑅_(𝑖,𝑗)=(4^(𝑗−1) 𝑅_(𝑖,𝑗−1)−𝑅_(𝑖−1,𝑗−1))/(4^(𝑗−1)−1)

data tpype: matrix

romberg<-function(f,a,b,n){
  R<-matrix(NA,nrow=n,ncol=n)
  R[1,1]<-MulTrap(f,a,b,1)
  for (i in 2:n){
    R[i,1]<-MulTrap(f,a,b,2^(i-1))
    for (j in 2:i){
      R[i,j]<-(4^(j-1)*R[i,j-1]-R[i-1,j-1])/(4^(j-1)-1)
    }
  }
  return(R)
}

Misalkan fungsi f(x)=exp(x) Deklarasikan dengan f<-function(x) exp(x) Bagaimanakah nilai dari Romberg(f,0,3,5) Bandingkan dengan jika memakai fungsi dasar integrate.

MulTrap <- function(f, a, b, n) {
  h <- (b - a) / n
  sum <- 0.5 * (f(a) + f(b))

  for (i in 1:(n-1)) {
    xi <- a + i * h
    sum <- sum + f(xi)
  }

  result <- h * sum
  return(result)
}

romberg <- function(f, a, b, n) {
  R <- matrix(NA, nrow = n, ncol = n)
  R[1, 1] <- MulTrap(f, a, b, 1)

  for (i in 2:n) {
    R[i, 1] <- MulTrap(f, a, b, 2^(i-1))

    for (j in 2:i) {
      R[i, j] <- (4^(j-1) * R[i, j-1] - R[i-1, j-1]) / (4^(j-1) - 1)
    }
  }

  return(R)
}

# Fungsi uji
f <- function(x) exp(x)
result_romberg <- romberg(f, 0, 3, 5)
result_integrate <- integrate(f, 0, 3)

result_romberg[nrow(result_romberg), ncol(result_romberg)]

## [1] 19.08562

result_integrate$value

## [1] 19.08554

##LATIHAN

Gunakan semua teknik yang sudah dipelajari untuk memberikan aproksimasi dari integral di bawah ini: 1. ∫0^{(𝜋/4)▒〖𝑒}3𝑥 sin⁡2𝑥 𝑑𝑥〗 2. ∫(−1)^1▒(cos⁡𝑥 )^2 𝑑𝑥 3. ∫_3^{3.5▒𝑥/√(𝑥}2−4) 𝑑𝑥 4. ∫_𝑒^2𝑒▒1/(𝑥 ln⁡𝑥 ) 𝑑𝑥

Jawab:

∫_0^{(𝜋/4)▒〖𝑒}3𝑥 sin⁡2𝑥 𝑑𝑥〗

Teknik Aproksimasi

# Fungsi yang diintegralkan
f <- function(x) exp(3 * x) * sin(2 * x)

# Metode trapesium untuk menghitung integral
MulTrap <- function(f, a, b, n) {
  h <- (b - a) / n
  sum <- 0.5 * (f(a) + f(b))

  for (i in 1:(n-1)) {
    xi <- a + i * h
    sum <- sum + f(xi)
  }

  result <- h * sum
  return(result)
}

# Batas integrasi
a <- 0
b <- pi/4

# Jumlah partisi
n <- 1000

# Menggunakan metode trapesium untuk mengaproksimasi integral
approximation <- MulTrap(f, a, b, n)

# Menampilkan hasil aproksimasi
print(approximation)

## [1] 2.58863

Teknik Riemann

# Fungsi yang diintegralkan
f <- function(x) exp(3 * x) * sin(2 * x)

# Batas integrasi
a <- 0
b <- pi/4

# Jumlah subinterval
n <- 1000

# Lebar subinterval
delta_x <- (b - a) / n

# Menggunakan teknik Riemann untuk mengaproksimasi integral
riemann_approximation <- sum(f(seq(a + 0.5 * delta_x, b - 0.5 * delta_x, delta_x)) * delta_x)

# Menampilkan hasil aproksimasi
print(riemann_approximation)

## [1] 2.588628

Teknik Newton cotes

# Fungsi yang diintegralkan
f <- function(x) exp(3 * x) * sin(2 * x)

# Batas integrasi
a <- 0
b <- pi/4

# Jumlah partisi
n <- 1000

# Lebar subinterval
h <- (b - a) / n

# Metode trapesium untuk mengaproksimasi integral
newton_cotes_approximation <- h * (0.5 * (f(a) + f(b)) + sum(f(seq(a + h, b - h, h))))

# Menampilkan hasil aproksimasi
print(newton_cotes_approximation)

## [1] 2.58863

Teknik Romberg

# Fungsi yang diintegralkan
f <- function(x) exp(3 * x) * sin(2 * x)

# Metode trapesium untuk menghitung integral
MulTrap <- function(f, a, b, n) {
  h <- (b - a) / n
  sum <- 0.5 * (f(a) + f(b))

  for (i in 1:(n-1)) {
    xi <- a + i * h
    sum <- sum + f(xi)
  }

  result <- h * sum
  return(result)
}

# Metode Romberg untuk mengaproksimasi integral
romberg <- function(f, a, b, n) {
  R <- matrix(NA, nrow = n, ncol = n)
  R[1, 1] <- MulTrap(f, a, b, 1)

  for (i in 2:n) {
    R[i, 1] <- MulTrap(f, a, b, 2^(i-1))

    for (j in 2:i) {
      R[i, j] <- (4^(j-1) * R[i, j-1] - R[i-1, j-1]) / (4^(j-1) - 1)
    }
  }

  return(R)
}

# Batas integrasi
a <- 0
b <- pi/4

# Jumlah iterasi Romberg
n_iter <- 5

# Menggunakan metode Romberg untuk mengaproksimasi integral
result_romberg <- romberg(f, a, b, n_iter)

# Menampilkan hasil aproksimasi dari elemen terakhir matriks Romberg
print(result_romberg[n_iter, n_iter])

## [1] 2.58864

Note: Dari hasil aproksimasi integral menggunakan 4 metode diatas, yang paling akurat adalah pada teknik keempat yaitu Teknik Romberg dengan hasil aproksimasi sebesar 2.58864.

∫_(−1)^1▒(cos⁡𝑥 )^2 𝑑𝑥

Teknik Aproksimasi

# Fungsi yang diintegralkan
f <- function(x) (cos(x))^2

# Batas integrasi
a <- -1
b <- 1

# Jumlah partisi
n <- 1000

# Lebar subinterval
h <- (b - a) / n

# Metode trapesium untuk mengaproksimasi integral
trapezoidal_approximation <- h * (0.5 * (f(a) + f(b)) + sum(f(seq(a + h, b - h, h))))

# Menampilkan hasil aproksimasi
print(trapezoidal_approximation)

## [1] 1.454648

Teknik Riemann

# Fungsi yang diintegralkan
f <- function(x) (cos(x))^2

# Batas integrasi
a <- -1
b <- 1

# Jumlah subinterval
n <- 1000

# Lebar subinterval
delta_x <- (b - a) / n

# Menggunakan teknik Riemann (titik tengah) untuk mengaproksimasi integral
riemann_approximation <- sum(f(seq(a + 0.5 * delta_x, b - 0.5 * delta_x, delta_x)) * delta_x)

# Menampilkan hasil aproksimasi
print(riemann_approximation)

## [1] 1.454649

Teknik Newton cotes

# Fungsi yang diintegralkan
f <- function(x) (cos(x))^2

# Batas integrasi
a <- -1
b <- 1

# Jumlah partisi
n <- 1000

# Lebar subinterval
h <- (b - a) / n

# Metode trapesium (Newton-Cotes) untuk mengaproksimasi integral
newton_cotes_approximation <- h * (0.5 * (f(a) + f(b)) + sum(f(seq(a + h, b - h, h))))

# Menampilkan hasil aproksimasi
print(newton_cotes_approximation)

## [1] 1.454648

Teknik Romberg

# Fungsi yang diintegralkan
f <- function(x) (cos(x))^2

# Metode trapesium untuk menghitung integral
MulTrap <- function(f, a, b, n) {
  h <- (b - a) / n
  sum <- 0.5 * (f(a) + f(b))

  for (i in 1:(n-1)) {
    xi <- a + i * h
    sum <- sum + f(xi)
  }

  result <- h * sum
  return(result)
}

# Metode Romberg untuk mengaproksimasi integral
romberg <- function(f, a, b, n) {
  R <- matrix(NA, nrow = n, ncol = n)
  R[1, 1] <- MulTrap(f, a, b, 1)

  for (i in 2:n) {
    R[i, 1] <- MulTrap(f, a, b, 2^(i-1))

    for (j in 2:i) {
      R[i, j] <- (4^(j-1) * R[i, j-1] - R[i-1, j-1]) / (4^(j-1) - 1)
    }
  }

  return(R)
}

# Batas integrasi
a <- -1
b <- 1

# Jumlah iterasi Romberg
n_iter <- 5

# Menggunakan metode Romberg untuk mengaproksimasi integral
result_romberg <- romberg(f, a, b, n_iter)

# Menampilkan hasil aproksimasi dari elemen terakhir matriks Romberg
print(result_romberg[n_iter, n_iter])

## [1] 1.45465

Note: Dari hasil aproksimasi integral menggunakan 4 metode diatas, yang paling akurat adalah pada teknik kedua yaitu Teknik Riemann dengan hasil aproksimasi sebesar 1.454649.

∫_3^{3.5▒𝑥/√(𝑥}2−4) 𝑑𝑥

Teknik Aproksimasi

# Fungsi yang diintegralkan
f <- function(x) x / sqrt(x^2 - 4)

# Batas integrasi
a <- 3
b <- 3.5

# Jumlah partisi
n <- 1000

# Lebar subinterval
h <- (b - a) / n

# Metode trapesium untuk mengaproksimasi integral
trapezoidal_approximation <- h * (0.5 * (f(a) + f(b)) + sum(f(seq(a + h, b - h, h))))

# Menampilkan hasil aproksimasi
print(trapezoidal_approximation)

## [1] 0.6362133

Teknik Riemann

# Fungsi yang diintegralkan
f <- function(x) x / sqrt(x^2 - 4)

# Batas integrasi
a <- 3
b <- 3.5

# Jumlah subinterval
n <- 1000

# Lebar subinterval
delta_x <- (b - a) / n

# Menggunakan teknik Riemann (titik tengah) untuk mengaproksimasi integral
riemann_approximation <- sum(f(seq(a + 0.5 * delta_x, b - 0.5 * delta_x, delta_x)) * delta_x)

# Menampilkan hasil aproksimasi
print(riemann_approximation)

## [1] 0.6362133

Teknik Newton cotes

# Fungsi yang diintegralkan
f <- function(x) x / sqrt(x^2 - 4)

# Batas integrasi
a <- 3
b <- 3.5

# Jumlah partisi
n <- 1000

# Lebar subinterval
h <- (b - a) / n

# Metode trapesium (Newton-Cotes) untuk mengaproksimasi integral
newton_cotes_approximation <- h * (0.5 * (f(a) + f(b)) + sum(f(seq(a + h, b - h, h))))

# Menampilkan hasil aproksimasi
print(newton_cotes_approximation)

## [1] 0.6362133

Teknik Romberg

# Fungsi yang diintegralkan
f <- function(x) x / sqrt(x^2 - 4)

# Metode trapesium untuk menghitung integral
MulTrap <- function(f, a, b, n) {
  h <- (b - a) / n
  sum <- 0.5 * (f(a) + f(b))

  for (i in 1:(n-1)) {
    xi <- a + i * h
    sum <- sum + f(xi)
  }

  result <- h * sum
  return(result)
}

# Metode Romberg untuk mengaproksimasi integral
romberg <- function(f, a, b, n) {
  R <- matrix(NA, nrow = n, ncol = n)
  R[1, 1] <- MulTrap(f, a, b, 1)

  for (i in 2:n) {
    R[i, 1] <- MulTrap(f, a, b, 2^(i-1))

    for (j in 2:i) {
      R[i, j] <- (4^(j-1) * R[i, j-1] - R[i-1, j-1]) / (4^(j-1) - 1)
    }
  }

  return(R)
}

# Batas integrasi
a <- 3
b <- 3.5

# Jumlah iterasi Romberg
n_iter <- 5

# Menggunakan metode Romberg untuk mengaproksimasi integral
result_romberg <- romberg(f, a, b, n_iter)

# Menampilkan hasil aproksimasi dari elemen terakhir matriks Romberg
print(result_romberg[n_iter, n_iter])

## [1] 0.6362151

Note: Dari hasil aproksimasi integral menggunakan 4 metode diatas, yang paling akurat adalah pada teknik keempat yaitu Teknik Romberg dengan hasil aproksimasi sebesar 0.6362151.

∫_𝑒^2𝑒▒1/(𝑥 ln⁡𝑥 ) 𝑑𝑥

Teknik Aproksimasi

# Fungsi yang diintegralkan
f <- function(x) 1 / (x * log(x))

# Batas integrasi
a <- exp(1)  # e
b <- exp(2)  # e^2

# Jumlah partisi
n <- 1000

# Lebar subinterval
h <- (b - a) / n

# Metode trapesium untuk mengaproksimasi integral
trapezoidal_approximation <- h * (0.5 * (f(a) + f(b)) + sum(f(seq(a + h, b - h, h))))

# Menampilkan hasil aproksimasi
print(trapezoidal_approximation)

## [1] 0.6931476

Teknik Riemann

# Fungsi yang diintegralkan
f <- function(x) 1 / (x * log(x))

# Batas integrasi
a <- exp(1)  # e
b <- exp(2)  # e^2

# Jumlah subinterval
n <- 1000

# Lebar subinterval
delta_x <- (b - a) / n

# Menggunakan teknik Riemann (titik tengah) untuk mengaproksimasi integral
riemann_approximation <- sum(f(seq(a + 0.5 * delta_x, b - 0.5 * delta_x, delta_x)) * delta_x)

# Menampilkan hasil aproksimasi
print(riemann_approximation)

## [1] 0.6931469

Teknik Newton cotes

# Fungsi yang diintegralkan
f <- function(x) 1 / (x * log(x))

# Batas integrasi
a <- exp(1)  # e
b <- exp(2)  # e^2

# Jumlah partisi
n <- 1000

# Lebar subinterval
h <- (b - a) / n

# Metode trapesium (Newton-Cotes) untuk mengaproksimasi integral
newton_cotes_approximation <- h * (0.5 * (f(a) + f(b)) + sum(f(seq(a + h, b - h, h))))

# Menampilkan hasil aproksimasi
print(newton_cotes_approximation)

## [1] 0.6931476

Teknik Romberg

# Fungsi yang diintegralkan
f <- function(x) 1 / (x * log(x))

# Metode trapesium untuk menghitung integral
MulTrap <- function(f, a, b, n) {
  h <- (b - a) / n
  sum <- 0.5 * (f(a) + f(b))

  for (i in 1:(n-1)) {
    xi <- a + i * h
    sum <- sum + f(xi)
  }

  result <- h * sum
  return(result)
}

# Metode Romberg untuk mengaproksimasi integral
romberg <- function(f, a, b, n) {
  R <- matrix(NA, nrow = n, ncol = n)
  R[1, 1] <- MulTrap(f, a, b, 1)

  for (i in 2:n) {
    R[i, 1] <- MulTrap(f, a, b, 2^(i-1))

    for (j in 2:i) {
      R[i, j] <- (4^(j-1) * R[i, j-1] - R[i-1, j-1]) / (4^(j-1) - 1)
    }
  }

  return(R)
}

# Batas integrasi
a <- exp(1)  # e
b <- exp(2)  # e^2

# Jumlah iterasi Romberg
n_iter <- 5

# Menggunakan metode Romberg untuk mengaproksimasi integral
result_romberg <- romberg(f, a, b, n_iter)

# Menampilkan hasil aproksimasi dari elemen terakhir matriks Romberg
print(result_romberg[n_iter, n_iter])

## [1] 0.6931516

Note: Dari hasil aproksimasi integral menggunakan 4 metode diatas, yang paling akurat adalah pada teknik keempat yaitu Teknik Romberg dengan hasil aproksimasi sebesar 0.6931516.

Komstat Aproksimasi Integral

Jasmine Georgina Sekartaji_121450159

2023-11-13