f <- function(t, y) -2*t*y^2
y.exact <- function(t) 1/(1+t^2)
euler <- function(t, y, h, f) f(t, y)
heun <- function(t, y, h, f) 0.5*(f(t, y) + f(t+h, y + h*f(t, y)))
rk4 <- function(t, y, h, f) {
K1 <- f(t, y)
K2 <- f(t + 0.5*h, y + 0.5*h*K1)
K3 <- f(t + 0.5*h, y + 0.5*h*K2)
K4 <- f(t + h, y + h*K3)
1/6*(K1 + 2*(K2+K3) + K4)
}
odesolve <- function(f, t0, T, y0, n, Phi, y.exact) {
t <- seq(t0, T, length.out=n+1)
h <- (T-t0)/n
y <- rep(0, n+1)
y[1] <- y0
for (j in 1:n) {
y[j+1] <- y[j] + h*Phi(t[j], y[j], h, f)
}
cbind(t, y, y - y.exact(t))
}
error <- function(n, Phi, y.exact) {
err <- rep(0, length(n))
for (k in seq_along(h)) {
y <- odesolve(f, 0, 1, 1, n[k], Phi, y.exact)
err[k] <- max(abs(y[,3]))
}
err
}
n <- 2^(0:12)
h <- 1/n
svg(point=16)
oldpar <- par(mar=c(5.1,4.1,1.1,1.1))
plot(c(1e-4, 1), c(1e-16, 1), xlab="h", ylab="Fehler", log="xy", type="n")
points(h, error(n, euler, y.exact), col="magenta")
lines(h, error(n, euler, y.exact), lty=2, col="magenta")
points(h, error(n, heun, y.exact), col="darkgreen")
lines(h, error(n, heun, y.exact), lty=2, col="darkgreen")
points(h, error(n, rk4, y.exact), col="gold4")
lines(h, error(n, rk4, y.exact), lty=2, col="gold4")
legend("bottomright", legend= c("expl. Euler", "Heun", "klass. Runge-Kutta"), pch=1, lty=2, col=c("magenta","darkgreen","gold4"))
par(oldpar)
dev.off()