# 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 y0 = 0.5 t1 = 0.5 # Where do we want to evaluate the solution: # x_grid = c( 0.4, 0.5 ) dy.dt = function(t, y){ (y^2 + 2*t*y)/(3 + t^2) } dy.dt.dy = function(t, y){ (2*y + 2*t)/(3 + t^2) } # The predictor corrector method: # h=0.1 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.1 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.1 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)