Metode Numerik

Tugas Kelompok

1 Anggota Kelompok

1. Nikita Indriyani
2. Putri Angelina Windjaya
3. Sherly Taurin Siridion

2 Soal 1

Tunjukkan 10 hasil iterasi metode Euler untuk persamaan diferensial \(f′(x, y) = y\) dimana \(x_0 = 0, y_0 = 2\) dan step size \(h = 0, 1\)!

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

f1 <- function(x,y){y}
num_euler <- euler(f1, x0=0, y0=2, h=c(0,1), n=10)
num_euler

##     x    y
## 1   0    2
## 2   0    2
## 3   1    4
## 4   0    2
## 5   2    8
## 6   0    2
## 7   3   16
## 8   0    2
## 9   4   32
## 10  0    2
## 11  5   64
## 12  0    2
## 13  6  128
## 14  0    2
## 15  7  256
## 16  0    2
## 17  8  512
## 18  0    2
## 19  9 1024
## 20  0    2
## 21 10 2048

3 Soal 2

Kerjakan lagi soal no.1 dengan menggunakan metode lainnya dan bandingkan seluruh metode tersebut menggunakan penyelesaian analitiknya!

3.1 Metode Heun

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

num_heun <- heun(f1, x0=0, y0=2, h=c(0,1), n=10)
num_heun

##     x          y
## 1   0     2.0000
## 2   0     2.0000
## 3   1     5.0000
## 4   0     2.0000
## 5   2    12.5000
## 6   0     2.0000
## 7   3    31.2500
## 8   0     2.0000
## 9   4    78.1250
## 10  0     2.0000
## 11  5   195.3125
## 12  0     2.0000
## 13  6   488.2812
## 14  0     2.0000
## 15  7  1220.7031
## 16  0     2.0000
## 17  8  3051.7578
## 18  0     2.0000
## 19  9  7629.3945
## 20  0     2.0000
## 21 10 19073.4863

3.2 Metode Titik Tengah

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

num_midpt <- midpt(f1, x0=0, y0=2, h=c(0,1), n=10)
num_midpt

##     x          y
## 1   0     2.0000
## 2   0     2.0000
## 3   1     5.0000
## 4   0     2.0000
## 5   2    12.5000
## 6   0     2.0000
## 7   3    31.2500
## 8   0     2.0000
## 9   4    78.1250
## 10  0     2.0000
## 11  5   195.3125
## 12  0     2.0000
## 13  6   488.2812
## 14  0     2.0000
## 15  7  1220.7031
## 16  0     2.0000
## 17  8  3051.7578
## 18  0     2.0000
## 19  9  7629.3945
## 20  0     2.0000
## 21 10 19073.4863

3.3 Metode Runge-Kutta Orde 4

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

num_rk4 <- rk4(f1, x0=0, y0=2, h=c(0,1), n=10)
num_rk4

##     x            y
## 1   0     2.000000
## 2   0     2.000000
## 3   1     5.416667
## 4   0     2.000000
## 5   2    14.670139
## 6   0     2.000000
## 7   3    39.731626
## 8   0     2.000000
## 9   4   107.606488
## 10  0     2.000000
## 11  5   291.434237
## 12  0     2.000000
## 13  6   789.301059
## 14  0     2.000000
## 15  7  2137.690367
## 16  0     2.000000
## 17  8  5789.578077
## 18  0     2.000000
## 19  9 15680.107292
## 20  0     2.000000
## 21 10 42466.957249

3.4 Metode Adams-Bashforth

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

num_adam <- adam(f1, x0=0, y0=2, h=c(0,1), n=10)
num_adam

##     x         y
## 1   0    2.0000
## 2   0    2.0000
## 3   1    4.0000
## 4   0    2.0000
## 5   2    9.0000
## 6   0    2.0000
## 7   3   20.5000
## 8   0    2.0000
## 9   4   46.7500
## 10  0    2.0000
## 11  5  106.6250
## 12  0    2.0000
## 13  6  243.1875
## 14  0    2.0000
## 15  7  554.6562
## 16  0    2.0000
## 17  8 1265.0469
## 18  0    2.0000
## 19  9 2885.2891
## 20  0    2.0000
## 21 10 6580.6992

3.5 Metode Analitik dan Perbandingan

f2 <- function(x){2^(x+1)}
x0 <- 0
y0 <- 2
x <- x0
y <- y0
for(i in 1:10){
y0 <- f2(x0+0.05)
x0 <- x0+0.05
x <- c(x, x0)
y <- c(y, y0)
}
num_analitik <- data.frame(x=x, y=y)
num_analitik

##       x        y
## 1  0.00 2.000000
## 2  0.05 2.070530
## 3  0.10 2.143547
## 4  0.15 2.219139
## 5  0.20 2.297397
## 6  0.25 2.378414
## 7  0.30 2.462289
## 8  0.35 2.549121
## 9  0.40 2.639016
## 10 0.45 2.732081
## 11 0.50 2.828427

library(ggplot2)
library(plotly)

num_adam$group     <- "adam"
num_analitik$group <- "analitik"
num_euler$group    <- "euler"
num_heun$group     <- "heun"
num_midpt$group    <- "midpt"
num_rk4$group      <- "rk4"

data               <- rbind(num_adam,
                            num_analitik,
                            num_euler,
                            num_heun,
                            num_midpt,
                            num_rk4)

theplot            <- ggplot(data,
                             aes(x = x,
                                 y = y,
                             color = group))+
                      geom_line()

ggplotly(theplot, tooltip = "text")

4 Soal 3

Bentuk kembali fungsi rk4sys() menggunakan algoritma metode Adams-Bashforth dan namai fungsi tersebut adamsys()!

adamsys <- function(f, x0, y0, h, n){
y1 <- y0 + h*f(x0,y0)
x1 <- x0 + h
x <- c(x0,x1)
y <- c(y0,y1)
n <- n-1
values <- data.frame(x=x,t(y0))
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
  values <- rbind(values, data.frame(x=x0,t(y0)))
  }
return(values)
}