Lembaga : Universitas Islam Negeri Maulana Malik Ibrahim Malang

Jurusan : Teknik Informatika

Fakultas : Sains dan Teknologi

0.1 Contoh Penerapan Paket ReacTran

0.1.1 Persamaan Adveksi-Difuksi 1D

Pada contoh kali ini, kita akan memodifikasi persamaan difusi satu dimensi yang telah dilakukan sebelumnya dengan menambahkan bagian adveksi serta diselesaikan menggunakan paket Reactran.

Siapkan kisi menggunakan fungsi setup.grid.1D() dan persiapkan nilai parameter persamaan.

library(ReacTran)
N     <- 100          # Number of grid cells
xgrid <- setup.grid.1D(x.up = 0, x.down = 1, N = N)
x     <- xgrid$x.mid  # Midpoints of grid cells
D     <- 1e-4         # Diffusion coefficient
v     <- 0.1          # Advection velocity

Bentuk fungsi yang menyatakan persamaan adveksi-difusi.

Diffusion <- function(t, Y, parms) {
    tran<-tran.1D(C=Y,C.up=0,C.down=0,D=D,v=v,dx= xgrid)
    list(dY = tran$dC, flux.up = tran$flux.up,
        flux.down = tran $flux.down)
}

Inisiasi konsentrasi awal pada kisi sel.

Yini <- rep(0,N) # Initial concentration = 0
Yini[2] <- 100   # Except in the second cell

Lakukan perhitungandengn menggunakan time step 0,01.

# Calculate for 5 time units
times <- seq(from = 0, to = 5, by = 0.01)
out <- ode.1D(y = Yini, times = times, func = Diffusion,
             parms = NULL,dimens = N)

Visualisasikan hasil perhitungan.

Visualisasi hasil simulasi persamaan adveksi-difusi menggunakan paket ReacTran

Gambar 10.12: Visualisasi hasil simulasi persamaan adveksi-difusi menggunakan paket ReacTran

0.1.2 Persamaan Gelombang 1D

Persamaan (10.16) dapat diselesaikan dengan cara yang sama dengan persamaan difusi dengan mengatur nilai c2=Dc2=D, memisalkan W=uW=u dan ∂u∂t=v∂u∂t=v, dan menyelesaikannya dengan cara yang familiar untuk variabel berpasangan (u,v)(u,v). Di sini kita misalkan persamaan gelombang 1D untuk senar yang dipetik, yang awalnya ditahan pada 00 amplitudo untuk x<−25x<−25 dan x>25x>25, dan direntangkan secara linear hingga maksimum pada x=0x=0. ode.1D() digunakan untuk menyelesaikan set himpunan ODE simultan dengan c=1c=1.

Langkah pertama yang harus dilakukan adalah melakukan penetapan parameter.

dx    <- 0.2          # Spacing of grid cells
# String extends from -100 to +100
xgrid <- setup.grid.1D(x.up = -100, 
                       x.down = 100, dx.1 = dx)
x     <- xgrid$x.mid  # midpoints of grid cells
N     <- xgrid$N      # number of grid cells

Menetapkan kondisi awal seperti ketinggian senar dan kecepatan.

uini  <- rep(0,N) # String height vector before stretching
vini  <- rep(0,N) # Initial string velocity vector
displ <- 10       # Initial displacement at center of string
# Impose initial triangular height profile on string between - 25
for(i in 1:N) {
  if (x[i] > -25 & x[i] <= 0) uini[i] = displ/25*(25 + x[i])else
  if (x[i] > 0 & x[i] < 25) uini[i] = displ/25*(25 - x[i])
}
yini  <- c(uini, vini)

Menetapkan deret waktu yang digunakan dalam simulasi.

times <- seq(from = 0, to = 50, by = 1)

Membangun fungsi yang akan diselesaikan secara numerik.

 wave <- function(t,y,parms) {
    u <- y[1:N] # Separate displacement and velocity vectors
    v <- y[(N+1):(2*N)]
    du<-v
    dv<-tran.1D(C=u,C.up=0,C.down=0,D=1,dx=xgrid)$dC
    return(list(c(du, dv))) 
}

Selesaikan persamaan menggunakan fungsi ode.1D() dengan metode adams.

out <- ode.1D(func = wave, y = yini, times = times,
             parms = NULL, method = "adams",
             dimens = N, names = c("u", "v"))
u <- subset(out, which = "u") # Extract displacement vector

Visualisasikan hasil perhitungan.

Visualisasi hasil simulasi persamaan gelombang menggunakan paket ReacTran

Gambar 10.13: Visualisasi hasil simulasi persamaan gelombang menggunakan paket ReacTran

0.1.3 Persamaan Laplace

Kita akan kembali melakukan simulasi terhadap Persamaan (10.19) menggunakan metode lainnya. Pada contoh kali ini, gradien pada sumbu yy bernilai -1. Gradien hanyalah fluks, D(∂C=∂x)D(∂C=∂x), dengan DD set sama dengan 1. Fungsi penyelesaian yang digunakan adalah steady.2D(), karena tidak ada ketergantungan waktu dalam persamaan. Sebagai kondisi awal yang arbitrer, kami menggunakan nomor acak Nx×NyNx×Ny yang terdistribusi secara seragam. Kita juga harus menentukan nspecnspec, jumlah spesies dalam model (hanya satu, potensi, dalam hal ini), dimensdimens, vektor dengan 2 nilai dengan jumlah sel dalam arah xx dan yy, dan lrw, panjang array. Lihat halaman bantuan untuk steady.2D() untuk informasi lebih detail.

Langkah pertama dalam melakukan simulasi adalah menetapkan parameter model.

Nx    <- 100
Ny    <- 100
xgrid <- setup.grid.1D(x.up = 0, x.down = 1, N = Nx)
ygrid <- setup.grid.1D(x.up = 0, x.down = 1, N = Ny)
x     <- xgrid$x.mid
y     <- ygrid$x.mid

Bentuk fungsi yang akan diselesaikan secara numerik.

laplace <- function(t, U, parms) {
    w = matrix(nrow = Nx, ncol = Ny, data = U)
    dw = tran.2D(C = w, C.x.up = 0, C.y.down = 0,
    flux.y.up = 0,
    flux.y.down = -1,
    D.x = 1, D.y = 1,
    dx = xgrid, dy = ygrid)$dC
    list(dw) 
}

Mulia dengan bilangan acak uniform sebagai kondisi awal, kemudian selesaikan untuk memperoleh nilai pada kondisi tunak dan buat plot kontur.

out = steady.2D(y = runif(Nx*Ny), func = laplace, 
                  parms = NULL,nspec = 1, 
                  dimens = c(Nx, Ny), lrw = 1e7)

Visualisasi hasil simulasi persamaan laplace pada kondisi tunak menggunakan paket ReacTran

Gambar 10.14: Visualisasi hasil simulasi persamaan laplace pada kondisi tunak menggunakan paket ReacTran

0.1.4 Persamaan Poisson untuk Dipol

Pada contoh kali ini, kita akan menyelesaikan persamaan Poisson untuk dipol yang ditampilkan pada Persamaan (10.21).

∂2w∂x2+∂2w∂y2=−ρε0(10.21)(10.21)∂2w∂x2+∂2w∂y2=−ρε0

untuk dipol yang terletak di tengah selembar persegi jika tidak pada 0 potensial. Untuk mempermudah, kita dapat mengatur semua faktor skala sama dengan satu. Dalam definisi fungsi poisson, nilai-nilai dalam matriks Nx×NyNx×Ny w adalah input melalui vektor data UU. Seperti dalam persamaan Laplace di atas, kita menetapkan nilai awal ww pada sel-sel grid sama dengan angka acak uniform.

Langkah pertama dalam melakukan simulasi adalah menetapkan parameter model.

Nx    <- 100
Ny    <- 100
xgrid <- setup.grid.1D(x.up = 0, x.down = 1, N = Nx)
ygrid <- setup.grid.1D(x.up = 0, x.down = 1, N = Ny)
x     <- xgrid$x.mid
y     <- ygrid$x.mid

Cari nilai xx dan yy pada titik kisi yang mendekati (0,4;0,5)(0,4;0,5) untuk muatan positif dan (0,6;0,5)(0,6;0,5) untuk muatan negatif.

# x and y coordinates of positive and negative charges
ipos <- which.min(abs(x - 0.4))
jpos <- which.min(abs(y - 0.50))
ineg <- which.min(abs(x - 0.6))
jneg <- which.min(abs(y - 0.50))

Bentuk fungsi Poisson yang akan diselesaikan secara numerik.

poisson <- function(t, U, parms) {
    w = matrix(nrow = Nx, ncol = Ny, data = U)
    dw = tran.2D(C = w, C.x.up = 0, C.y.down = 0,
    flux.y.up = 0,
    flux.y.down = 0,
    D.x = 1, D.y = 1,
    dx = xgrid, dy = ygrid)$dC
    dw[ipos,jpos] = dw[ipos,jpos] + 1
    dw[ineg,jneg] = dw[ineg,jneg] - 1
    list(dw) 
}

Selesaikan untuk kondisi tunak dan buat visualisasinya.

out <- steady.2D(y = runif(Nx*Ny), 
                 func = poisson, 
                 parms = NULL, 
                 nspec = 1, 
                 dimens = c(Nx, Ny), 
                 lrw = 1e7)

Visualisasi hasil simulasi persamaan Poisson pada kondisi tunak menggunakan paket ReacTran

Gambar 10.15: Visualisasi hasil simulasi persamaan Poisson pada kondisi tunak menggunakan paket ReacTran

0.2 Studi Kasus

0.2.1 Model Lotka-Voltera

Model Lotka-Voltera merupakan model yang menggambarkan interaksi antara predator dan mangsa. Pada studi kasus kali ini model yang akan digunakan merupakan model dengan 3 bentuk interaksi yaitu: tanaman uu, herbivora vv, dan karnivora ww. Sistem persamaan diferensial sistem tersebut ditampilkan pada kumpulan persamaan berikut:

dudt=au−α1f1(u,v)(10.22)(10.22)dudt=au−α1f1(u,v)

dvdt=−bv+α1f1(u,v)−α2f2(v,w)(10.23)(10.23)dvdt=−bv+α1f1(u,v)−α2f2(v,w)

dwdt=−c(w−w∗)+α2f2(v,w)(10.24)(10.24)dwdt=−c(w−w∗)+α2f2(v,w)

dimana w∗w∗ merupakan tingkat predasi minimum untuk menstabilkan populasi ketika populasi mangsa rendah. Interaksi antar komponen digambarkan ke dalam bentuk persamaan logistik.

fi(x,y)=xy1+kix(10.25)(10.25)fi(x,y)=xy1+kix

Langkah pertama untuk menyelesaikan model adalah melakukan penentuan parameter model dan pembentukan fungsi model.

library(deSolve)
f <- function(x,y,k){x*y/(1+k*x)}
model <- function(t, xx, parms) {
    u = xx[1] # plant resource
    v = xx[2] # herbivore
    w = xx[3] # carnivore
    with(as.list(parms),{
    du = a*u - alpha1*f(u, v, k1)
    dv = -b*v + alpha1*f(u, v, k1) - alpha2*f(v, w, k2)
    dw = -c*(w - wstar) + alpha2*f(v, w, k2)
    list(c(du, dv, dw))
})}

times <- seq(0, 200, 0.1)
parms <- c(a=1, b=1, c=10, alpha1=0.2, alpha2=1,
k1<-0.05, k2=0, wstar=0.006)
xstart <- c(u=10, v=5, w=0.1)

Proses iterasi selanjutnya akan menggunakan metode lsoda, dimana metode iterasi ini dapat berganti secara otomatis dalam menangani sistem persamaan diferensial stiff dan non-stiff. Persamaan diferensial dikatakan stiff apabila variabel dependen berubah berdasarkan 2 atau lebih variabel independen yang sangat berbeda skalanya.

out = ode(xstart, times, model, parms, 
          method ="lsoda")

Ketiga variabel selanjutnya divisualisasikan.

Visualisasi simulasi model Lotka-Voltera

Gambar 10.16: Visualisasi simulasi model Lotka-Voltera

0.3 Referensi

  1. Bloomfield, V.A. 2014. Using R for Numerical Analysis in Science and Engineering. CRC Press.

  2. Chapra, S.C. Canale, R.P. 2015. Numerical Methods For Engineers, Seventh Edition. Mc Graw Hill.

  3. Griffiths, G.W. 2016. Numerical analysis using R : solutions to ODEs and PDEs. Cambridge University Press.

  4. Howard, J.P. 2017. Computational Methods for Numerical Analysis with R. CRC Press.

  5. Kreyszig, E. 2011. Advanced Engineering Mathematics, 10th Edition. John Wiley & Sons.

  6. Soetaert,K., Cash J., Mazzia F. 2012. Solving Differential Equations in R. Springer.

  7. Suparno, S. 2008. Komputasi untuk Sains dan Teknik Edisi II. Departemen Fisika-FMIPA Universitas Indonesia.