# 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) source('runge_kutta.R') t0 = 0.0 y0 = 1 t1 = 2.0 # Where do we want to evaluate the solution: # x_grid = c( seq( t0, 0.5, by=0.1 ), c( 1.0, 1.5, 2.0 ) ) dy.dt = function(t, y){ 1 - t + 4*y } h=0.2 res = runge_kutta(dy.dt, h=h, start=t0, y0=y0, end=t1) rk_yhat_h1 = knn.reg( as.matrix( res\$xs ), y=res\$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' ) h=0.1 res = runge_kutta(dy.dt, h=h, start=t0, y0=y0, end=t1) rk_yhat_h2 = knn.reg( as.matrix( res\$xs ), y=res\$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' ) h=0.05 res = runge_kutta(dy.dt, h=h, start=t0, y0=y0, end=t1) rk_yhat_h3 = knn.reg( as.matrix( res\$xs ), y=res\$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' ) A = rbind( rk_yhat_h1\$pred, rk_yhat_h2\$pred, rk_yhat_h3\$pred ) rownames(A) = NULL R = cbind( c( 0.2, 0.1, 0.05 ), A ) colnames(R) = c( 'h', sprintf( 'y(%.2f)', x_grid ) ) print(R)