# 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, u, parameters){ with(as.list(c(u, parameters)),{ t = u[1] y = u[2] du = numeric(2) du[1] = 1 du[2] = lambda * y + 2 * t - lambda * t^2 list(du) }) } t0 = 0.0 t1 = 1.0 y0 = c(0, 0) # Where do we want to evaluate the solution: # t_grid = seq( t0, t1, by=0.1 ) # Where will we compute the solution: # h = 0.01 compute_grid = seq( t0, t1, by=h ) # Integrate the ODE for each value of lambda: # lambdas = c(1, 10, 20, 50) YHAT = c() for( lam in lambdas ){ res = ode(y0, compute_grid, my_yprime, parms=c(lambda=lam), method='rk4') res = as.data.frame(res) colnames(res) = c('time','x','y') yhat = knn.reg(as.matrix(res\$time), y=res\$y, test=as.matrix(t_grid), k=1, algorithm='brute') YHAT = cbind( YHAT, yhat\$pred ) } YHAT = cbind(t_grid, YHAT) colnames(YHAT) = c('t', sprintf('lam= %.0f', lambdas) ) print(YHAT)