# 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('predictor_corrector.R') source('adams_moulton.R') source('backwards_differentiation.R') t0 = 0.0 y0 = 3.0 t1 = 2.0 # Where do we want to evaluate the solution: # x_grid = c( 0.5, 1.0, 1.5, 2.0 ) dy.dt = function(t, y){ sqrt(t+y) } dy.dt.dy = function(t, y){ (1/2*sqrt(t+y)) } # The predictor corrector method: # h=0.05 res = predictor_corrector(dy.dt, h=h, start=t0, y0=y0, end=t1) pc_yhat = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' ) # The Adams-Moulton method: # h=0.05 res = adams_moulton(dy.dt, dy.dt.dy, h=h, start=t0, y0=y0, end=t1) am_yhat = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' ) # Backwards differentiation: # h=0.05 res = backwards_differentiation(dy.dt, dy.dt.dy, h=h, start=t0, y0=y0, end=t1) bd_yhat = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' ) A = rbind( pc_yhat$pred, am_yhat$pred, bd_yhat$pred ) rownames(A) = NULL R = cbind( c( 'Predictor_Corrector', 'Adams_Moulton', 'Backwards_Difference' ), A ) colnames(R) = c( 'method', sprintf( 'y(%.2f)', x_grid ) ) print(R)