32 lines
941 B
R
32 lines
941 B
R
rk4 <- function(f, x0, y0, x1, n) {
|
|
vx <- double(n + 1)
|
|
vy <- double(n + 1)
|
|
vx[1] <- x <- x0
|
|
vy[1] <- y <- y0
|
|
h <- (x1 - x0)/n
|
|
for(i in 1:n) {
|
|
k1 <- h*f(x, y)
|
|
k2 <- h*f(x + 0.5*h, y + 0.5*k1)
|
|
k3 <- h*f(x + 0.5*h, y + 0.5*k2)
|
|
k4 <- h*f(x + h, y + k3)
|
|
vx[i + 1] <- x <- x0 + i*h
|
|
vy[i + 1] <- y <- y + (k1 + k2 + k2 + k3 + k3 + k4)/6
|
|
}
|
|
cbind(vx, vy)
|
|
}
|
|
|
|
sol <- rk4(function(x, y) x*sqrt(y), 0, 1, 10, 100)
|
|
cbind(sol, sol[, 2] - (4 + sol[, 1]^2)^2/16)[seq(1, 101, 10), ]
|
|
|
|
vx vy
|
|
[1,] 0 1.000000 0.000000e+00
|
|
[2,] 1 1.562500 -1.457219e-07
|
|
[3,] 2 3.999999 -9.194792e-07
|
|
[4,] 3 10.562497 -2.909562e-06
|
|
[5,] 4 24.999994 -6.234909e-06
|
|
[6,] 5 52.562489 -1.081970e-05
|
|
[7,] 6 99.999983 -1.659460e-05
|
|
[8,] 7 175.562476 -2.351773e-05
|
|
[9,] 8 288.999968 -3.156520e-05
|
|
[10,] 9 451.562459 -4.072316e-05
|
|
[11,] 10 675.999949 -5.098329e-05
|