# I use the function knn.reg to evaluate the numerical solution at a fixed grid of points # if( !require('FNN') ){ install.packages('FNN', dependencies=TRUE, repos='http://cran.rstudio.com/') } library(FNN) # Use the ode function in the deSolve library to integrate the differential equations: # if( !require('deSolve') ){ install.packages('deSolve') } library(deSolve) my_yprime = function(t, y, parameters){ x = y[1] y = y[2] dy = numeric(2) dy[1] = x*(1 - 0.5*x - 0.5*y) dy[2] = y*(-0.25 + 0.5*x) list(dy) } t0 = 0.0 t1 = 1.0 y0 = c(4, 1) # Where do we want to evaluate the solution: # t_grid = seq( t0, t1, by=0.2 ) # Euler's method: # h = 0.1 compute_grid = seq( t0, t1, by=h ) res = ode(y0, compute_grid, my_yprime, parms=c(), method='euler') res = as.data.frame(res) colnames(res) = c('time', 'x', 'y') e_xhat = knn.reg( as.matrix(res$time), y=res$x, test=as.matrix(t_grid), k=1, algorithm='brute' ) e_yhat = knn.reg( as.matrix(res$time), y=res$y, test=as.matrix(t_grid), k=1, algorithm='brute' ) R = cbind( t_grid, e_xhat$pred, e_yhat$pred ) colnames(R) = c('time', 'x_hat', 'y_hat') print(sprintf('Eulers method (h=%f)', h)) print(R) # RK4 method: # h = 0.2 compute_grid = seq( t0, t1, by=h ) res = ode(y0, compute_grid, my_yprime, parms=c(), method='rk4') res = as.data.frame(res) colnames(res) = c('time', 'x', 'y') e_xhat = knn.reg( as.matrix(res$time), y=res$x, test=as.matrix(t_grid), k=1, algorithm='brute' ) e_yhat = knn.reg( as.matrix(res$time), y=res$y, test=as.matrix(t_grid), k=1, algorithm='brute' ) R = cbind( t_grid, e_xhat$pred, e_yhat$pred ) colnames(R) = c('time', 'x_hat', 'y_hat') print(sprintf('RK4 method (h=%.2f)', h)) print(R) # RK4 method smaller h: # h = 0.1 compute_grid = seq( t0, t1, by=h ) res = ode(y0, compute_grid, my_yprime, parms=c(), method='rk4') res = as.data.frame(res) colnames(res) = c('time', 'x', 'y') e_xhat = knn.reg( as.matrix(res$time), y=res$x, test=as.matrix(t_grid), k=1, algorithm='brute' ) e_yhat = knn.reg( as.matrix(res$time), y=res$y, test=as.matrix(t_grid), k=1, algorithm='brute' ) R = cbind( t_grid, e_xhat$pred, e_yhat$pred ) colnames(R) = c('time', 'x_hat', 'y_hat') print(sprintf('RK4 method (h=%.2f)', h)) print(R)