Persamaan Diferensial (Turunan) pada RStudio

Lembaga : Universitas Islam Negeri Maulana Malik Ibrahim Malang

Jurusan : Teknik Informatika

Fakultas : Sains dan Teknologi # Chapter 10 Persamaan Diferensial

Persamaan diferensial merupakan persoalan matematis yang sering dijumpai dalam bidang teknik lingkungan. Sering kali suatu persamaan diferensial tidak dapat diselesaikan secara analitik sehingga diperlukan metode numerik untuk menyelesaikannya. Pada Chapter 10, kita akan membahas masalah-masalah dalam persamaan diferensial dan metode penyelesaiannya. Adapun yang akan dibahas pada Chapter 10 kali ini antara lain:

Initial value problems
Sistem persamaan diferensial
Persamaan diferensial parsial

0.1 10.1 Initial value problems

Initial value problems merupakan permasalahan yang sering ditemukan pada proses dekomposisi zat kimia atau polutan dalam reaktor. Penyelesaiaan persamaan diferensial biasanya dipersulit dengan tidak tersedianya informasi yang cukup untuk menyelesaikannya.Sebuah persamaan diferensial f′(x,…)f′(x,…) merupakan hasil diferensiasi beberapa fungsi f(x,…)f(x,…). Proses penyelesaian persamaan diferensial, dan menemukan nilai f(x,…)f(x,…) untuk beberapa nilai x,… tidak dimungkinkan karena integral dari f′(x,…)f′(x,…) hanya digambarkan bentuk umum. Pergeseran vertikal, atas atau bawah, tidak diketahui. Pergeseran vertikal ini menghasilkan konstanta integrasi.

Selama proses diferensiasi, nilai apapun dari proses pergeseran vertikal (integrasi) akan hilang sebagai akibat dari eliminasi konstanta yang memiliki turunan 0. Kita biasa melakukannya ketika mengintegrasikan fungsi dengan menambahkan konstanta +C+C pada proses integrasi ke integral yang tidak terbatas. Hal ini terkadang bukan menjadi permasalahan sebab jika menemukan nilai integrasi pada suatu batas tertentu syarat +C+C dibatalkan dan konstanta integrasi tidak diperlukan.

Untuk persamaan diferensial biasa, tidak ada pembatalan yang nyaman, yang mengarah ke initial value problems. Initial value problems memberikan nilai f(x0,…)f(x0,…), di mana x0x0 biasanya bernilai 0, meski tidak diharuskan. Nilai awal ini memberikan informasi yang cukup untuk menyelesaikan persamaan dan menemukan nilai aktual dari f(x,…)f(x,…) untuk sejumlah nilai xx. Terdapat beberapa metode yang akan dibahas pada Chapter 10.1, antara lain:

Metode Euler
Metode Heun
Metode Titik Tengah
Metode Runge-Kutta Orde 4
Metode multistep linier

0.1.1 10.1.1 Metode Euler

Metode Euler merupakan metode paling sederhana yang diturunkan dari deret Taylor. Penyelesaian initial value problems menggunakan metode Euler dilakukan melalui Persamaan (10.1).

yi+1=yi+f(xi,yi)h(10.1)(10.1)yi+1=yi+f(xi,yi)h

dimana ii merupakan tahapan iterasi.

Algoritma Metode Euler

Tentukan titik awal integrasi x0x0 dan y0y0.
Tentukan jumlah iterasi nn dan step size hh yang digunakan.
Lakukan integrasi menggunakan Persamaan (10.1).

Algoritma tersebut, selanjutnya dapat disusun ke dalam sebuah fungsi R. Fungsi tersebut adalah sebagai berikut:

euler <- function(f, x0, y0, h, n){
  x <- x0
  y <- y0
  
  for(i in 1:n){
    y0 <- y0 + h*f(x0, y0)
    x0 <- x0 + h
    x <- c(x,x0)
    y <- c(y, y0)
  }
  
  return(data.frame(x=x, y=y))
}

Contoh 10.1 Selesaikan persamaan diferensial di bawah ini, jika diketahui f(0)=1 menggunakan h=0,05 dan n=100!

f′(x,y)=y2x+1f′(x,y)=y2x+1

Jawab:

Penyelesaian secara analitik persamaan tersebut untuk nilai f(0)=1f(0)=1 sebagai berikut:

f(x)=√2x+1f(x)=2x+1

Secara numerik persamaan tersebut dapat diselesaikan sebagai berikut:

iterasi 1

y(0,5)=1+0,05×1(2⋅0)+1=1,05y(0,5)=1+0,05×1(2⋅0)+1=1,05

iterasi 2

y(0,1)=1,05+0,05×1(2⋅0,05)+1=1.097727y(0,1)=1,05+0,05×1(2⋅0,05)+1=1.097727

Kita dapat juga menggunakan fungsi euler() untuk menyelesaikan persamaan tersebut secara numerik. Hasil yang diperoleh selanjutnya diplotkan dengan hasil yang diperoleh menggunakan metode analitik. Berikut adalah sintaks yang digunakan:

# metode numerik
f1 <- function(x,y){y/(2*x+1)}
num <- euler(f1, x0=0, y0=1, h=0.05, n=100)

# metode analitik
f2 <- function(x){sqrt(2*x+1)}
x0 <- 0
y0 <- 1
x <- x0
y <- y0

for(i in 1:100){
  y0 <- f2(x0+0.05)
  x0 <- x0+0.05
  x <- c(x, x0)
  y <- c(y, y0)
}
true <- data.frame(x=x, y=y)

Visualisasi integrasi numerik dengan metode Euler dan metode analitik

Gambar 10.1: Visualisasi integrasi numerik dengan metode Euler dan metode analitik

Berdasarkan hasil visualisasi dapat dilihat bahwa metode Euler dapat dengan baik memberikan pendekatan nilai integrasi persamaan. Pembaca dapat mencoba untuk melakukan simulasi kembali dengan nilai hh yang lebih kecil.

0.1.2 10.1.2 Metode Heun

Metode Heun merupakan salah satu peningkatan dari metode Euler. Metode ini melibatkan 2 buah persamaan. Persamaan pertama disebut sebagai persamaan prediktor yang digunakan untuk memprediksi nilai integrasi awal (Persamaan (10.2)). Persamaan kedua disebut sebagai persamaan korektor yang mengoreksi hasil integrasi awal (Persamaan (10.3)). Metode Heun pada Chapter ini merupakan metode prediktor-korektor satu tahapan. Akurasi integrasi dapat ditingkatkan dengan melakukan koreksi ulang terhadap nilai koreksi semula menggunakan persamaan kedua.

y0i+1=yi+f(xi,yi)h(10.2)(10.2)yi+10=yi+f(xi,yi)h

yi+1=yi+f(xi,yi)+f(xi+1,y0i+1)2h(10.3)(10.3)yi+1=yi+f(xi,yi)+f(xi+1,yi+10)2h

Algoritma Metode Heun

Tentukan titik awal integrasi x0x0 dan y0y0.
Tentukan jumlah iterasi nn dan step size hh yang digunakan.
Lakukan prediksi nilai awal dengan Persamaan (10.2).
Lakukan koreksi nilai awal menggunakan Persamaan (10.3).
Lakukan koreksi terhadap nilai koreksi yang dihasilkan sebelumnya menggunakan Persamaan (10.3).

Kita dapat membangun sebuah fungsi yang dapat melakukan proses integrasi menggunakan metode Heun. Berikut adalah sintaks yang digunakan:

heun <- function(f, x0, y0, h, n, iter=1){
  x <- x0
  y <- y0
  
  for(i in 1:n){
    ypred0 <- f(x0,y0)
    ypred1 <- y0 + h*ypred0
    ypred2 <- f(x0+h,ypred1)
    ykor <- y0 + h*(ypred0+ypred2)/2
    if(iter!=1){
      for(i in 1:iter){
        ykor <- y0 + h*(ypred0+f(x0+h,ykor))/2
      }
    }
    y0 <- ykor
    x0 <- x0 + h
    x <- c(x, x0)
    y <- c(y, y0)
  }
  
  return(data.frame(x=x,y=y))
}

Contoh 10.2 Selesaikan kembali persamaan yang ditampilkan pada Contoh 10.1 menggunakan metode Heun!

Jawab:

Contoh perhitungan secara manual menggunakan metode Heun untuk sekali iterasi adalah sebagai berikut:

f′(0;1)=1(2⋅0)+1=1f′(0;1)=1(2⋅0)+1=1

y01=1+0,05⋅1=1,05y10=1+0,05⋅1=1,05

y′1=f′(0,05;1,05)=1,05(2⋅0,05)+1=0,9545455y1′=f′(0,05;1,05)=1,05(2⋅0,05)+1=0,9545455

y1=1+0,05⋅1+0,95454552=1,047727y1=1+0,05⋅1+0,95454552=1,047727

Penyelesaian persamaan tersebut menggunakan fungsi heun() dengan iterasi pada nilai koreksi sebanyak 1 kali disajikan pada sintaks berikut:

num <- heun(f1, x0=0, y0=1, h=0.05, n=100)

Visualisasi integrasi numerik dengan metode Heun dan metode analitik

Gambar 10.2: Visualisasi integrasi numerik dengan metode Heun dan metode analitik

0.1.3 10.1.3 Metode Titik Tengah

Metode titik tengah menggunakan setengah step size pada metode Euler untuk melakukan estimasi terhadap integral suatu persamaan diferensial. Metode ini melakukan perhitungan melalui dua tahapan yaitu: menghitung nilai estimasi integral pada setengah step size(Persamaan (10.4)) dan menghitung nilai integral menggunkan hasil perhitungan setengah step size sebelumnya (Persamaan (10.5)).

yi+12=yi+f(xi,yi)h2(10.4)(10.4)yi+12=yi+f(xi,yi)h2

yi+1=yi+f(xi+12,yi,12)h(10.5)(10.5)yi+1=yi+f(xi+12,yi,12)h

Algoritma Metode Tengah

Tentukan titik awal integrasi x0x0 dan y0y0.
Tentukan jumlah iterasi nn dan step size hh yang digunakan.
Lakukan integrasi pada setengah tahapan iterasi menggunakan Persamaan (10.4).
Lakukan iterasi pada setengah tahapan selanjutnya menggunakan Persamaan (10.5).

Berdasarkan algoritma tersebut, kita dapat membangun sebuah fungsi pada R yang adapat digunakan untuk menyelesaikan persamaan diferensial menggunakan metode titik tengah. Berikut adalah sintaks yang digunakan:

midpt <- function(f, x0, y0, h, n){
  x <- x0
  y <- y0
  
  for(i in 1:n){
    s1 <- y0 + f(x0,y0) * h/2
    s2 <- h * f(x0+h/2,s1)
    y0 <- y0 + s2
    x0 <- x0 + h
    x <- c(x, x0)
    y <- c(y, y0)
  }
  
  return(data.frame(x=x,y=y))
}

Contoh 10.3 Selesaikan kembali persamaan yang ditampilkan pada Contoh 10.1 menggunakan metode titik tengah!

Jawab:

Contoh perhitungan secara manual menggunakan metode titik tengah untuk sekali iterasi adalah sebagai berikut:

y12=1+1(2⋅0)+1⋅0,052=1,025y12=1+1(2⋅0)+1⋅0,052=1,025

y1=1+1,025(2⋅0,025)+1⋅0,05=1,0488y1=1+1,025(2⋅0,025)+1⋅0,05=1,0488

Penyelesaian persamaan tersebut menggunakan fungsi midpt() dengan iterasi pada nilai koreksi sebanyak 1 kali disajikan pada sintaks berikut:

num <- midpt(f1, x0=0, y0=1, h=0.05, n=100)

Visualisasi integrasi numerik dengan metode titik tengah dan metode analitik

Gambar 10.3: Visualisasi integrasi numerik dengan metode titik tengah dan metode analitik

0.1.4 10.1.4 Metode Runge-Kutta Orde 4

Runge-Kutta orde 4 merupakan metode yang paling populer dalam penyelesaian persamaan diferensial. Metode ini dapat memperoleh akurasi deret Taylor tanpa memerlukan diferensiasi orde yang lebih tinggi. Metode Runge-Kutta orde 4 dituliskan ke dalam Persamaan (10.6).

yi+1=yi+16(k1+2k2+2k3+k4)h(10.6)(10.6)yi+1=yi+16(k1+2k2+2k3+k4)h

dimana

k1=f(xi,yi)(10.7)(10.7)k1=f(xi,yi)

k2=f(xi+12h,yi+12k1h)(10.8)(10.8)k2=f(xi+12h,yi+12k1h)

k3=f(xi+12h,yi+12k2h)(10.9)(10.9)k3=f(xi+12h,yi+12k2h)

k4=f(xi+h,yi+k3h)(10.10)(10.10)k4=f(xi+h,yi+k3h)

Algoritma Metode Tengah

Tentukan titik awal integrasi x0x0 dan y0y0.
Tentukan jumlah iterasi nn dan step size hh yang digunakan.
Lakukan integrasi menggunakan Persamaan (10.6).

Berdasarkan algoritma tersebut, kita dapat membangun sebuah fungsi pada R yang adapat digunakan untuk menyelesaikan persamaan diferensial menggunakan metode Runge-Kutta orde 4. Berikut adalah sintaks yang digunakan:

rk4 <- function(f, x0, y0, h, n){
  x <- x0
  y <- y0
  
  for(i in 1:n){
    k1 <- f(x0,y0)
    k2 <- f(x0+0.5*h,y0+0.5*k1*h)
    k3 <- f(x0+0.5*h,y0+0.5*k2*h)
    k4 <- f(x0+h,y0+k3*h)
    y0 <- y0 + (1/6)*(k1+2*k2+2*k3+k4)*h
    x0 <- x0 + h
    x <- c(x, x0)
    y <- c(y, y0)
  }
  
  return(data.frame(x=x,y=y))
}

Contoh 10.4 Selesaikan kembali persamaan yang ditampilkan pada Contoh 10.1 menggunakan metode Runge-Kutta orde 4!

Jawab:

Contoh perhitungan secara manual menggunakan metode titik tengah untuk sekali iterasi adalah sebagai berikut:

k1=1(2⋅0)+1=1k1=1(2⋅0)+1=1

k2=1+12⋅1⋅0,05(2⋅0⋅12⋅0,05)+1=1,025k2=1+12⋅1⋅0,05(2⋅0⋅12⋅0,05)+1=1,025

k2=1+12⋅1,025⋅0,05(2⋅0⋅12⋅0,05)+1=1,025625k2=1+12⋅1,025⋅0,05(2⋅0⋅12⋅0,05)+1=1,025625

k2=1+1,025625⋅0,05(2⋅0⋅0,05)+1=1,051281k2=1+1,025625⋅0,05(2⋅0⋅0,05)+1=1,051281

y1=1+16(1+2⋅1,025+2⋅1,025625+1,051281)0,05=1.051271y1=1+16(1+2⋅1,025+2⋅1,025625+1,051281)0,05=1.051271

1+(1/6)*(1+2*1.025+2*1.025625+1.051281)*0.05

## [1] 1.051

Iterasi dapat pula dilakukan dengan menggunakan fungsi rk4(). Berikut adalah sintaks yang digunakan:

num <- rk4(f1, x0=0, y0=1, h=0.05, n=100)

Visualisasi integrasi numerik dengan metode Runge-Kutta orde 4 dan metode analitik

Gambar 10.4: Visualisasi integrasi numerik dengan metode Runge-Kutta orde 4 dan metode analitik.

0.1.5 10.1.5 Metode Multistep Linier

Jika metode Runge-Kutta mengalami kesulitan karena terlalu banyak evaluasi fungsi yang digunakan, masuk akal untuk bertanya apakah kita dapat menggunakan kembali beberapa evaluasi fungsi sebelumnya, yang sudah kita buat. Sebagai contoh, kita ingin tahu apakah kita dapat menggunakan kembali estimasi f(0,1)f(0,1) dan f(0,2)f(0,2) untuk memperkirakan nilai f(0,3)f(0,3). Jika kita dapat menggunakan kembali perkiraan sebelumnya, kita dapat memperoleh akurasi tambahan tanpa menimbulkan penalti kinerja yang terkait dengan evaluasi fungsi tambahan. Metode multistep linier dikembangkan untuk mengatasi masalah ini.

Di satu sisi, metode multistep linier dasar untuk persamaan diferensial hanya mencakup satu titik xixi, dalam perhitungan xi+1xi+1. Ini persis bagaimana fungsi metode Euler dan metode Euler merupakan metode multistep linier dasar. Metode selanjutnya menggunakan xi−1xi−1 dan xixi untuk menghitung xi+1xi+1. Metode Adams-Bashforth menggunakan tambahan berbobot, termasuk bobot negatif, dari langkah dan poin untuk sampai pada langkah berikutnya. Seperti metode numerik lainnya, bobot muncul dari interpolasi polinomial titik yang tersedia.

Metode Adam-Bashforth orde 2 didasarkan pada Persamaan (10.11).

yi+2=yi+1+h2(3f(xi+1,yi+1)−f(xi,yi))(10.11)(10.11)yi+2=yi+1+h2(3f(xi+1,yi+1)−f(xi,yi))

Pendekatan ini melakukan interpolasi antara titik sebelumnya untuk memperkirakan titik ketiga dalam grup. Titik ketiga ini menjadi titik tengah dari iterasi berikutnya saat seluruh proses berlanjut. Karena nilai sebelumnya disimpan dan digunakan kembali, evaluasi fungsi tambahan tidak diperlukan. Jika metode Runge-Kutta dapat dibandingkan dengan tip-toeing melalui bidang vektor, maka metode Adams-Bashforth dapat sama dibandingkan dengan menjalankan melalui bidang vektor. Namun, ini tidak berarti metode Adams-Bashforth lebih unggul.

Algoritma Metode Tengah

Tentukan titik awal integrasi x0x0 dan y0y0.
Tentukan jumlah iterasi nn dan step size hh yang digunakan.
Lakukan pendekatan pada iterasi ke-1 menggunakan metode Euler.
Lakukan integrasi ke-2 sampai n menggunakan Persamaan (10.11).

Berdasarkan algoritma tersebut, kita dapat membangun sebuah fungsi pada R yang adapat digunakan untuk menyelesaikan persamaan diferensial menggunakan metode multistep linier. Berikut adalah sintaks yang digunakan:

adambashforth <- function(f, x0, y0, h, n){
  # pendekatan Euler untuk x1 dan y1
  y1 <- y0 + h*f(x0,y0)
  x1 <- x0 + h
  
  x <- c(x0,x1)
  y <- c(y0,y1)
  n <- n-1
  
  for(i in 1:n){
    yn <- y1 + 1.5*h*f(x1,y1) - 0.5*h*f(x0,y0)
    xn <- x1 + h
    
    y0 <- y1
    x0 <- x1
    y1 <- yn
    x1 <- xn
    
    y <- c(y,y1)
    x <- c(x,x1)
  }
  
  return(data.frame(x=x,y=y))
}

Contoh 10.5 Selesaikan kembali persamaan yang ditampilkan pada Contoh 10.1 menggunakan metode Runge-Kutta orde 4!

Jawab:

Untuk melakukan iterasi, kita dapat menggunakan fungsi adambashforth(). Berikut adalah sintaks yang digunakan:

num <- rk4(f1, x0=0, y0=1, h=0.05, n=100)

Visualisasi integrasi numerik dengan metode Adam-Bashfoth orde 2 dan metode analitik

Gambar 10.5: Visualisasi integrasi numerik dengan metode Adam-Bashfoth orde 2 dan metode analitik

0.2 10.2 Sistem Persamaan Diferensial

Sistem persamaan diferensial akan sering pembaca temui dalam pemodelan sistem dinamik. Pada proses pemodelan sistem tersebut akan ditemukan aksi-interasi antar komponennya yang dinyatakan ke dalam suatu sistem persamaan diferensial.

Pada fungsi rk4sys() disusun sebuah fungsi pada R untuk melakukan iterasi terhadap sistem persamaan diferensial menggunakan metode Runge-Kutta orde 4. Berikut adalah sintaks yang digunakan:

rk4sys <- function(f, x0, y0, h, n){
  x <- x0
  y <- y0
  
  values <- data.frame(x=x,t(y0))
  for(i in 1:n){
    k1 <- f(x0,y0)
    k2 <- f(x0+0.5*h,y0+0.5*k1*h)
    k3 <- f(x0+0.5*h,y0+0.5*k2*h)
    k4 <- f(x0+h,y0+k3*h)
    y0 <- y0 + (1/6)*(k1+2*k2+2*k3+k4)*h
    x0 <- x0 + h
    values <- rbind(values, data.frame(x=x0,t(y0)))
  }
  
  return(values)
}

Contoh 10.6 Sebuah model populasi yang dibagi menjadi tiga kelompok umur: anak (0-12 tahun), melahirkan anak (13-40 tahun), dan usia (41 tahun atau lebih). Kelompok 0-12 meningkat berdasarkan kelahiran di kelompok 13-40, dan menurun dengan kematian dan lewat bagian ke kelompok 13-40. Kelompok 13–40 meningkat dengan perolehan dari kelompok 0-12, dan berkurang dengan kematian dan lewat bagian ke dalam kelompok> = 41. Kelompok> = 41 meningkat dengan keuntungan dari kelompok 13-40, dan menurun dengan kematian. Parameter dipilih untuk mewakili tingkat kelahiran dan kematian yang cukup tinggi seperti yang ditemukan di banyak masyarakat berkembang. Model interaksi tersebut dinyatakan ke dalam persamaan interkasi di bawah ini! Simulasikan dinamika populasi pada model tersebut pada 10 tahun ke depan, jika diketahui populasi semula masing-masing kelompok usia secara berurutan adalah 200,400, dan 400, serta nilai masing-masing koefisien kelahiran, kematian kelompok 1, kematian kelompok2, dan kematian kelompok 3 secara berurutan adalah 0,5;0,1;0,1;0,25!

pop1′=b⋅pop2′+1112⋅pop1′(1−d1)pop1′=b⋅pop2′+1112⋅pop1′(1−d1)

pop2′=112⋅pop1′⋅(1−d1)+2627⋅pop2′⋅(1−d2)pop2′=112⋅pop1′⋅(1−d1)+2627⋅pop2′⋅(1−d2)

pop3,=127⋅pop2,⋅(1−d2)+pop3′⋅(1−d3)pop3,=127⋅pop2,⋅(1−d2)+pop3′⋅(1−d3)

Jawab:

Sistem persamaan diferensial perlu ditransformasi ke dalam bentuk sebuah fungsi pada R.

Population = function(x,y) {
   # parameter
   b = 0.5 # Birth rate in 13-40 group
   d1 = 0.1 # Death rate of 0-12 group
   d2 = 0.1 # Death rate of 13-40 group
   d3 = 0.25 # Death rate of 41 and older
   
   y1 = y[1] # 0-12 group population
   y2 = y[2] # 13-40 group population
   y3 = y[3] # 41 and older population
   y1.new = b*y2 + (11/12)*y1*(1 - d1)
   y2.new = (1/12)*y1*(1 - d1) + (26/27)*y2*(1 - d2)
   y3.new = (1/27)*y2*(1 - d2) + y3*(1 - d3)
   return(c(y1.new, y2.new, y3.new))
}

Nilai awal dituliskan seperti berikut:

# Nilai awal
y = c(pop1=200,pop2=400,pop3=400)

Simulasi model ditampilkan pada sintaks berikut:

pop <- rk4sys(Population, x0=0, y0=y, h=0.1, n=100)
head(pop)

##     x  pop1  pop2  pop3
## 1 0.0 200.0 400.0 400.0
## 2 0.1 239.0 437.9 432.6
## 3 0.2 283.4 479.6 467.9
## 4 0.3 334.0 525.4 506.1
## 5 0.4 391.4 575.9 547.4
## 6 0.5 456.4 631.3 592.1

Visualisasi hasil simulasi model dinamika populasi

Gambar 10.6: Visualisasi hasil simulasi model dinamika populasi

0.3 10.3 Penyelesaian Persamaan Diferensial dan Sistem Persamaan Diferensial Menggunakan Fungsi `ode()`

Fungsi ode() pada paket deSolve merupakan salah satu fungsi yang dapat digunakan untuk menyelesaikan persamaan dan sistem persamaan diferensial. Fungsi ini mudah digunakan serta menyediakan berbagai macam metode iterasi numerik untuk menyelesaikan persamaan diferensial. Format umum fungsi ini antara lain:

ode(y, times, func, parms,method , ...)

Catatan:

y : nilai awal (kondisi awal) suatu persamaan diferensial.

times : deret waktu yang terdiri dari kapan simulasi dimulai, kapan simulasi berkahir dan berapa step size yang digunakan.

func : fungsi yang berisi persamaan diferensial. Return value pada fungsi haruslah berupa list.

parms : list parameter yang diinputan kedalam func.

method : sebuah string yang berupa metode intgrasi yang digunakan. Metode yang tersedia dan kegunaan metode tersebut antara lain:

“bdf”: menangani persamaan diferensial menggunakan formula diferensiasi mundur (bacward differensiation) dan cocok untuk menangani kondisi stiff. Metode ini setara dengan method="lsode".

“bdf_d”: menggunakan formula diferensiasi mundur yang memanfaatkan iterasi Jacobi-Newton (mengabaikan elemen diagonal Jacobian). Cocok digunakan persamaan atau sistem persamaan diferensial kondisi stiff. Metode ini setara dengan method="lsode",mf=23.

“adams”: metode Adams yang menggunakan iterasi fungsional (tanpa menggunakan Jacobian). Cocok digunakan untuk persamaan atau sistem persamaan non stiff. Metode ini setara dengan method="lsode", mf=10.

“impAdams”: metode Adams implisit yang menggunakan iterasi Newton-Raphson. Metode ini setara dengan method="lsode",mf=12.

“impAdams_d”: metode Adams implisit yang menggunakan iterasi Jacobi-Newton. Metode ini setara dengan method="lsode, mf=13".

“euler”: metode iterasi Euler

“rk4”: metode iterasi Runge-Kutta orde 4.

metode lain: “lsoda”, “lsode”, “lsodes”, “lsodar”, “vode”, “daspk”, “ode23”, “ode45”, “radau”.

...: argumen tambahan untuk integrator atau method.

Contoh penerapan fungsi tersebut menggunakan Contoh 10.6, sebagai berikut:

library(deSolve)

# sistem persamaan diferensial
Population = function(t,y,param) {
  y1 = y[1] # 0-12 group population
  y2 = y[2] # 13-40 group population
  y3 = y[3] # 41 and older population
  y1.new = b*y2 + 11/12*y1*(1 - d1)
  y2.new = 1/12*y1*(1 - d1) + 26/27*y2*(1 - d2)
  y3.new = 1/27*y2*(1 - d2) + y3*(1 - d3)
  return(list(c(y1.new, y2.new, y3.new)))
}

# parameter
b = 0.5 # Birth rate in 13-40 group
d1 = 0.1 # Death rate of 0-12 group
d2 = 0.1 # Death rate of 13-40 group
d3 = 0.25 # Death rate of 41 and older

# nilai awal
y = c(pop1=200, pop2=400, pop3=400)

# waktu
Tahun = seq(0,10,0.1)

# simulasi
out = ode(func = Population, y = y, times = Tahun, 
          parms = c(b,d1,d2,d3), method = "rk4")

Visualisasi hasil simulasi model dinamika populasi menggunakan paket desolve

Gambar 10.7: Visualisasi hasil simulasi model dinamika populasi menggunakan paket desolve

0.4 10.4 Persamaan Diferensial Parsial

Persamaan diferensial parsial (PDE) banyak dijumpai pada pemodelan transport polutan dalam bidang teknik lingkungan. Persamaan diferensial parsial merupakan persamaan diferensial yang melibatkan lebih dari satu variabel independen, biasanya variabel waktu dan satu atau lebih variabel posisi atau beberapa variabel spasial. PDE diklasifikasikan menjadi 3 jenis: parabolik (time-dependent dan difusif), hiperbolik (time-dependent dan gelombang), dan eliptik (time-independent). Dalam penyelesaian PDE pada umumnya kita menggunakan metode FTCS (forward in time, centered in space). Untuk memahami definisi tersebut, pembaca dapat membaca kembali Chapter 9.1.

0.4.1 10.4.1 Persamaan Difusi

Persamaan difusi merupakan contoh PDE parabolik. Persamaan difusi satu dimensi spasial ditampilkan pada Persamaan (10.12).

∂C∂t=D∂2C∂x2(10.12)(10.12)∂C∂t=D∂2C∂x2

Persamaan tersebut mirip dengan persamaan konduksi panas. Pada persamaan konduksi panas, variabel konsentrasi CC diganti dengan variabel temperatur TT, dan koefisien difusi D diganti dengan koefsien difusi termal K.

Untuk menyelesaikan Persamaan (10.12), persamaan tersebut diubah ke dalam bentuk metode beda hingga, dimana turunan waktu menggunakan pendekatan Euler (metode beda hingga maju) dan turunan variabel spasial diubah ke dalam bentuk pendekatan titik pusat. Proses diskretisasi Persamaan (10.12), ditampilkan pada Persamaan (10.13).

C(i+1,j)−C(i,j)Δt=DC(i,j+1)+C(i,j−1)−2C(i,j)Δx2(10.13)(10.13)C(i+1,j)−C(i,j)Δt=DC(i,j+1)+C(i,j−1)−2C(i,j)Δx2

dimana ii merupakan step untuk variabel waktu tt dan jj merupakan step untuk variabel spasial xx.

Persamaan (10.13) dapat disusun kembali sehingga menjadi Persamaan (10.14) yang menyatakan persamaan konsentrasi CC pada saat i+1i+1 pada posisi jj.

C(i+1,j)=C(i,j)+AC(i,j+1)+C(i,j−1)−2C(i,j)(10.14)C(i+1,j)=C(i,j)+A[C(i,j+1)+C(i,j−1)−2C(i,j)]

dimana

A=DΔtΔx2(10.15)(10.15)A=DΔtΔx2

Untuk stabilitas komputasi, pemilihan peningkatan waktu terhadap jarak yang dinyatakan pada nilai AA harus ≤12≤12.

Pada contoh berikut, kita akan melakukan simulasi menggunakan Persamaan (10.14). Parameter yang digunakan dan nilai awal yang digunakan dinyatakan pada sintaks berikut:

dt    <- 3                # Timestep, s
dx    <- 0.1              # Distance step, cm
D     <- 1e-4             # Diffusion coeff, cm^2/s

# Cek apakah syarat stabilitas terpenuhi
D*dt/dx^2 <= 0.5

## [1] TRUE

# Desain grid points
L     <- 1                # Length from -L/2 to L/2
n     <- L/dx + 1         # Number of grid points
x     <- seq(-L/2,L/2,dx) # Location of grid points
steps <- 30               # Number of iterations
time  <- 0:steps

Langkah selanjutnya adalah inisiasi konsentrasi awal, dimana seluruh konsentrasi awal pada tiap grid adalah nol kecuali pada grid pusat.

C <- matrix(rep(0, (steps+1)*n), nrow = steps+1, ncol = n)
C[1, round(n/2)] <- 1/dx  # initial spike at central point

Simulasi selanjutnya dilakukan dengan melakukan iterasi pada grid points konsentrasi yang telah dibuat.

# Loop over desired number of time steps
for(i in 1:(steps-1)){
  for(j in 2:(n-1)){
    C[i+1,j] <- C[i,j] + (D*dt/dx^2)*(C[i,j+1] + C[i,j-1] - 2*C[i,j])
  }
}

Visualisasi dari hasil simulasi tersebut ditampilkan pada Gambar 10.8.

Visualisasi 3D simulasi difusi partikulat

Gambar 10.8: Visualisasi 3D simulasi difusi partikulat

0.4.2 10.4.2 Persamaan Gelombang

Persamaan gelombang 1 dimensi ditampilkan pada Persamaan (10.16).

∂2W∂t2=c2∂2W∂x2(10.16)(10.16)∂2W∂t2=c2∂2W∂x2

dimana WW merupakan pemindahan dan cc merupakan kecepatan gelombang. Persamaan (10.16) merupakan bentuk PDE hiperbolik. Versi lebih sederhana dari Persamaan (10.16) adalah persamaan adveksi yang ditampilkan pada Persamaan (10.17).

∂y∂t=−c∂y∂x(10.17)(10.17)∂y∂t=−c∂y∂x

Persamaan (10.17) menggambarkan evolusi pada bidang skalar y(x,y)y(x,y) dibawa oleh gelombang dengan kecepatan konstan cc dan bergerak dari kiri ke kanan jika c>0c>0. Persamaan adveksi merupakan contoh palin sederhana dari persamaan konservasi flux.

Diskretisasi Persamaan (10.17) ditampilkan pada Persamaan (10.18).

y(i+1,j)=y(i,j)−cΔt2Δxy(i,j+1)−y(i,j−1)(10.18)y(i+1,j)=y(i,j)−cΔt2Δx[y(i,j+1)−y(i,j−1)]

Pada contoh berikut, kita akan melakukan simulasi menggunakan Persamaan (10.17). Parameter yang digunakan dan nilai awal yang digunakan dinyatakan pada sintaks berikut:

dt    <- 0.002              # Timestep, s
L     <- 1                  # Length from -L/2 to L/2
n     <- 50                 # Number of grid points
v     <- 1                  # Wavespeed, cm/s
dx    <- L/n                # Distance step, cm
x     <- (1:n-0.5)*dx-L/2   # Location of grid points
steps <- L/(v*dt)           # Number of iterations
time  <- 0:steps  
sig   <- 0.1                # Standard deviation of initial Gaussian wave
amp0  <- exp(-x^2/(2*sig^2))# Initial Gaussian amplitude

Kondisi awal dan kondisi batas periodik dinyatakan sebagai berikut:

C <- matrix(rep(0, (steps+1)*n), nrow = steps+1, ncol = n)
C[1, ] <- amp0
jplus1 <- c(2:n,1)
jminus1 <- c(n,1:(n-1))

Langkah selanjutnya melakukan interasi untuk melihat perubahan nilai pada grid points.

for(i in 1:steps){
  for(j in 1:n){
    C[i+1,j] <- C[i,j] + (v*dt/(2*dx))*(C[i,jplus1[j]] - C[i,jminus1[j]])
  }
}

Visualisasi dari hasil simulasi tersebut ditampilkan pada Gambar 10.9.

Visualisasi simulasi adveksi

Gambar 10.9: Visualisasi simulasi adveksi

0.4.3 10.4.3 Persamaan Laplace

Persamaan Laplace dalam dua dimensi disajikan pada Persamaan (10.19).

∂V∂x2+∂V∂y2=0(10.19)(10.19)∂V∂x2+∂V∂y2=0

Persamaan (10.19) merupakan contoh tipe ketiga PDE, persamaan elips. Persoalan ini sering muncul dalam bidang elektrostatik, gravitasi, dan bidang lain di mana potensi VV harus dihitung sebagai fungsi posisi. Jika ada muatan atau massa dalam ruang, dan jika kita menggeneralisasi ke tiga dimensi, persamaannya menjadi persamaan Poisson

∂V∂x2+∂V∂y2+∂V∂z2=0(10.20)(10.20)∂V∂x2+∂V∂y2+∂V∂z2=0

Bergantung pada geometri masalahnya, Persamaan (10.20) juga dapat ditulis dalam koordinat bola, silindris, atau lainnya.

Untuk menyelesaikan persamaan eliptik jenis ini, persamaan harus diberi syarat batas. Biasanya dengan menentukan bahwa titik, garis, atau permukaan tertentu dipertahankan pada nilai konstan potensial. Kemudian potensi di titik lain disesuaikan sampai mencapai perkiraan yang diinginkan. (Dalam kasus yang jarang terjadi, persamaan dengan kondisi batas dapat diselesaikan dengan tepat secara analitis; tetapi biasanya solusi perkiraan harus cukup.)

Ada banyak pendekatan untuk solusi numerik dari persamaan Laplace. Mungkin yang paling sederhana adalah Jacobi, di mana titik-titik interior secara berturut-turut didekati dengan rata-rata titik-titik di sekitarnya, sedangkan titik-batas dibatasi pada nilai-nilai tetap dan yang ditentukan. Kita asumsikan sebagai contoh bidang persegi, dibatasi oleh (0,1)(0,1) dalam arah xx dan yy, di mana ujung pada y=1y=1 dipertahankan pada V=1V=1 dan tiga tepi lainnya dipertahankan pada V=0V=0. Kita buat tebakan awal yang paling sederhana untuk potensi di titik interior, tetapi ini akan disamakan saat solusinya menyatu.

Pada contoh berikut, kita akan melakukan simulasi persamaan Laplace menggunakan metode Jacobi. Parameter yang digunakan dalam simulasi adalah sebagai berikut:

n   <- 30                 # Number of grid points/side
L   <- 1                  # Length of a side
dx  <- L/(n-1)            # Grid spacing
x   <- y    <- 0:(n-1)*dx # x and y coorfinates

Tebakan awal untuk profil voltase adalah sebagai berikut:

V0  <- 1
V   <- matrix(V0/2*sin(2*pi*x/L)*sin(2*pi*y/L),
              nrow = n, ncol = n, byrow = TRUE)

Kondisi batas di tetapkan seperti berikut:

V[1,] <- 0
V[n,] <- 0
V[,1] <- 0
V[,n] <- V0*rep(1,n)

Selanjutnya visualisasikan tebakan awal.

Visualisasi tebakan awal solusi persaman Laplace

Gambar 10.10: Visualisasi tebakan awal solusi persaman Laplace

Proses iterasi menggunakan metode Jacobi ditampilkan pada sintaks berikut:

newV  <- V
itmax <- n^2
tol   <- 1e-4
for(it in 1:itmax){
  dVsum <- 0
  for(i in 2:(n-1)){
    for(j in 2:(n-1)){
      newV[i,j] <- 0.25*(V[i-1,j]+V[i+1,j]+V[i,j-1]+V[i,j+1])
      dVsum <- dVsum + abs(1-V[i,j]/newV[i,j])
    }
  }
  V <- newV
  dV = dVsum/(n-2)^2 # Average deviation from prev. value
  if(dV < tol) break # Desired tolerance achived
}

it # Iterations to achieve convergence to tol

## [1] 419

dV

## [1] 9.908e-05

Hasil simulasi selanjutnya di visualisasikan kembali.

Visualisasi hasil simulasi solusi persaman Lplace

Gambar 10.11: Visualisasi hasil simulasi solusi persaman Lplace

## Penyelesaian Persamaan Diferensial Parsial Menggunakan PaketReacTran`

Package ReacTran memfasilitasi pemodelan transport reaktif dalam 1, 2, dan 3 dimensi. Paket ini "berisi rutinitas yang memungkinkan pengembangan model transportasi reaktif dalam sistem perairan (sungai, danau, lautan), media berpori (agregat flok, sedimen, …) dan bahkan organisme yang diidealkan.

Pada paket ReacTran terdapat sejumlah fasilitas fungsi, antara lain:

Fungsi untuk menyiapkan kisi beda hingga (1D atau 2D)
Fungsi untuk melampirkan parameter dan properti ke kisi ini (1D atau 2D)
Fungsi untuk menghitung jangka waktu transport adveksi-difusi di atas grid (1D, 2D, 3D)
Berbagai fungsi lainnya.

Saat paket ReacTran dimuat, paket ini juga memuat dua paket pendukung, yaitu: rootSolve dan deSolve. Paket rootSolve berguna untuk memecahkan persamaan diferensial untuk kondisi tunak baik persamaan diferensial uniform atau multikomponen (1D, 2D, dan 3D). Sedangkan deSolve berguna untuk menyediakan fungsi yang digunakan untuk memperoleh penyelesaian numerik persamaan diferensial biasa (ODE), persamaan diferensial parsial (PDE), persamaan aljabar diferensial (DAE) dan persamaan delay differential.

0.4.4 10.4.4 `setup.grid.1D()`

Fungsi setup.grid.1D() digunakan untuk membentuk kisi satu dimensi. Secara sederhana fungsi ini membagi ruang satu dimensi LL antara x.upx.up dan x.downx.down menjadi sejumlah NN kisi sebesar dx.1. Format umum fungsi tersebut, adalah sebagai berikut:

setup.grid.1D(x.up=0, x.down=NULL,L=NULL, 
              N=NULL, dx.1=NULL, 
              p.dx.1= rep(1,length(L)), 
              max.dx.1=L, dx.N=NULL, 
              p.dx.N=rep(1,length(L)), 
              max.dx.N=L)

Catatan:

x.up : posisi hilir.

x.down : posisi hulu.

L : x.down-x.up.

N : jumlah kisi = L/dx.1.

Pada situasi yang lebih kompleks, ukuran sel atau kisi dapat bervariasi, atau dapat pula lebih dari satu zona. Kondisi ini dijelaskan lebih jauh pada laman bantuan fungsi.

Nilai yang dihasilkan dari fungsi setup.grid.1D() termasuk x.mid (vektor sepanjang NN yang menyatakan posisi titik tengah) dan x.int (vektor sepanjang N+1N+1 yang menyatakan posisi antar muka antara kisi sel dimana flux diukur).

Fungsi plot untuk grid.1D() memvisualissikan baik posisi sel dan ketebalan kotak, menampilkan x.mid dan x.int. Contoh di halaman bantuan menunjukkan perilaku ini.

setup.grid.1D() berfungsi sebagai titik awal untuk setup.grid.2D, yang membuat kisi di atas domain persegi panjang yang ditentukan oleh dua kisi 1D ortogonal.

0.4.5 10.4.5 `setup.prop.1D()`

Banyak model transportasi akan melibatkan kisi-kisi dengan properti konstan. Tetapi jika beberapa properti yang mempengaruhi difusi atau adveksi bervariasi dengan posisi di grid, variasi dapat digabungkan dengan fungsi setup.prop.1D() (atau setup.prop.2D() dalam dua dimensi).

Diberikan fungsi matematis atau matriks data, fungsi setup.prop.1D() menghitung nilai properti yang diminati di tengah sel kisi dan pada antarmuka antar sel. Format fungsi tersebut adalah sebagai berikut:

setup.prop.1D(func=NULL, value=NULL, xy=NULL,
              interpolate="spline",grid, ...)

Catatan:

func : fungsi yang mengatur ketergantungan spasial properti.

value : nilai konstan yang diberikan ke properti jika tidak ada ketergantungan spasial.

xy : matriks data di mana kolom pertama memberikan posisi, dan kolom kedua memberikan nilai-nilai yang diinterpolasi di atas grid.

iterpolate : metode interpolasi yang digunakan (spline atau linier).

grid: objek yang dibentuk melalui fungsi setup.grid.1D().

...: argumen tambahan func.

0.4.6 10.4.6 `tran.1D()`

Fungsi ini menghitung istilah transportasi — laju perubahan konsentrasi akibat difusi dan adveksi — dalam model cairan 1D (fraksi volume = 1) atau padatan berpori (fraksi volume dapat bervariasi dan <1).tran.1D() juga digunakan untuk masalah dalam geometri bola atau silinder, meskipun dalam kasus ini antarmuka sel jaringan akan memiliki area variabel. Format fungsi ini adalah sebagai berikut:

tran.1D(C, C.up = C[1], C.down = C[length(C)], 
        flux.up = NULL, flux.down = NULL, 
        a.bl.up = NULL, a.bl.down = NULL, D = 0, 
        v = 0, AFDW = 1, VF = 1, A = 1, dx, 
        full.check = FALSE, full.output = FALSE)

Catatan:

C : vektor konsentrasi pada titik tengah kisi sel.

C.up, C.down : konsentrasi pada hulu dan hilir batas.

flux.up,flux.down : flux dari dan keluar sistem pada hulu dan hilir batas.

Jika terdapat tranport konvektif sepanjang hulu dan hilir batas, a.bl.up dan a.bl.down merupakan koefisiennya.

D: Koefisien difusi.

v: koefisien adveksi.

VF: fraksi volume.

A: fraksi area.

dx: ketebalan kisi, baik nilai konstan atau vektor.

full.check, full.output: nilai logik untuk memeriksa konsistensi dan mengatur output dari perhitungan. Keduanya FALSE secara default.

Ketika full.output = FALSE, nilai-nilai yang dikembalikan oleh trans.1D() adalah dC, yaitu: laju perubahan C di pusat setiap sel kisi karena transportasi, dan flux.up dan flux.down, yatu: fluks masuk dan keluar dari model di batas hulu dan hilir.

ReacTran juga memiliki fungsi untuk memperkirakan istilah difusi dan adveksi dalam model dua dan tiga dimensi, dan dalam koordinat silinder dan kutub. Jumlah input tumbuh dengan dimensi, tetapi input pada dasarnya sama seperti pada kasus 1D. Lihat halaman bantuan untuk tran.2D(), tran.3D(), tran.cylindrical(), dan tran.polar().

Namun penyempurnaan lain adalah fungsi tran.volume.1D(), yang memperkirakan istilah transportasi volumetrik dalam model 1D. Berbeda dengan tran.1D(), yang menggunakan fluks (massa per satuan luas per satuan waktu), tran.volume.1D() menggunakan aliran (massa per satuan waktu). Ini berguna untuk memodelkan saluran yang area lintas bagiannya berubah, ketika perubahan area tidak perlu dimodelkan secara eksplisit. Ini juga memungkinkan input lateral dari saluran samping.

0.4.7 10.4.7 `ode.1D()` atau `steady.1D()`

Setelah kisi telah diatur dan properti ditetapkan dan model transport telah diformulasikan dengan tran.1D() (atau analog 2D atau 3D-nya), maka ReacTran memanggil ode.1D() dari paket deSolve jika solusi time-dependent diperlukan, atau stable.1D() dari paket rootSolve jika solusi kondisi-stabil diinginkan. Sistem ODE yang dihasilkan dari metode pendekatan garis biasanya sparse dan non-stiff. Integrator dalam deSolve, seperti “lsoda” (metode default 1D) sangat cocok untuk menangani sistem persamaan tersebut. Jika sistem ODE non-stiff, maka “adams” umumnya merupakan metode yang baik.

0.5 10.5 Contoh Penerapan Paket `ReacTran`

0.5.1 10.5.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.5.2 10.5.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.5.3 10.5.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.5.4 10.5.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

Persamaan Diferensial (Turunan) pada RStudio

Heny Rimadana, Prof. Dr. Suhartono, M.Kom

12/3/2021

0.1 10.1 Initial value problems

0.1.1 10.1.1 Metode Euler

0.1.2 10.1.2 Metode Heun

0.1.3 10.1.3 Metode Titik Tengah

0.1.4 10.1.4 Metode Runge-Kutta Orde 4

0.1.5 10.1.5 Metode Multistep Linier

0.2 10.2 Sistem Persamaan Diferensial

0.3 10.3 Penyelesaian Persamaan Diferensial dan Sistem Persamaan Diferensial Menggunakan Fungsi `ode()`

0.4 10.4 Persamaan Diferensial Parsial

0.4.1 10.4.1 Persamaan Difusi

0.4.2 10.4.2 Persamaan Gelombang

0.4.3 10.4.3 Persamaan Laplace

0.4.4 10.4.4 `setup.grid.1D()`

0.4.5 10.4.5 `setup.prop.1D()`

0.4.6 10.4.6 `tran.1D()`

0.4.7 10.4.7 `ode.1D()` atau `steady.1D()`

0.5 10.5 Contoh Penerapan Paket `ReacTran`

0.5.1 10.5.1 Persamaan Adveksi-Difuksi 1D

0.5.2 10.5.2 Persamaan Gelombang 1D

0.5.3 10.5.3 Persamaan Laplace

0.5.4 10.5.4 Persamaan Poisson untuk Dipol

Persamaan Diferensial (Turunan) pada RStudio

Heny Rimadana, Prof. Dr. Suhartono, M.Kom

12/3/2021

0.1 10.1 Initial value problems

0.1.1 10.1.1 Metode Euler

0.1.2 10.1.2 Metode Heun

0.1.3 10.1.3 Metode Titik Tengah

0.1.4 10.1.4 Metode Runge-Kutta Orde 4

0.1.5 10.1.5 Metode Multistep Linier

0.2 10.2 Sistem Persamaan Diferensial

0.3 10.3 Penyelesaian Persamaan Diferensial dan Sistem Persamaan Diferensial Menggunakan Fungsi ode()

0.4 10.4 Persamaan Diferensial Parsial

0.4.1 10.4.1 Persamaan Difusi

0.4.2 10.4.2 Persamaan Gelombang

0.4.3 10.4.3 Persamaan Laplace

0.4.4 10.4.4 setup.grid.1D()

0.4.5 10.4.5 setup.prop.1D()

0.4.6 10.4.6 tran.1D()

0.4.7 10.4.7 ode.1D() atau steady.1D()

0.5 10.5 Contoh Penerapan Paket ReacTran

0.5.1 10.5.1 Persamaan Adveksi-Difuksi 1D

0.5.2 10.5.2 Persamaan Gelombang 1D

0.5.3 10.5.3 Persamaan Laplace

0.5.4 10.5.4 Persamaan Poisson untuk Dipol

0.3 10.3 Penyelesaian Persamaan Diferensial dan Sistem Persamaan Diferensial Menggunakan Fungsi `ode()`

0.4.4 10.4.4 `setup.grid.1D()`

0.4.5 10.4.5 `setup.prop.1D()`

0.4.6 10.4.6 `tran.1D()`

0.4.7 10.4.7 `ode.1D()` atau `steady.1D()`

0.5 10.5 Contoh Penerapan Paket `ReacTran`