# 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('../Chapter2/euler.R') t0 = 0 y0 = 1.0 t1 = 2.0 # Where do we want to evaluate the solution: # x_grid = seq( t0, t1, by=0.5 ) # Eulers method: # dy.dt = function(t, y){ 2*t + exp(-t*y) } h=0.025 res = euler(dy.dt, h=h, start=t0, y0=y0, end=t1) e_yhat_h1 = knn.reg( as.matrix( res\$xs ), y=res\$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' ) h=0.0125 res = euler(dy.dt, h=h, start=t0, y0=y0, end=t1) e_yhat_h2 = knn.reg( as.matrix( res\$xs ), y=res\$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' ) A = rbind( e_yhat_h1\$pred, e_yhat_h2\$pred ) rownames(A) = NULL R = cbind( c( 0.025, 0.0125 ), A ) colnames(R) = c( 'h', sprintf('y(%.2f)', x_grid ) ) print('Problem 10 (Euler):') print(R) # Backward Eulers method: # source('beuler.R') ynp1_res_fn = function(ynp1, yn, tn, h){ ynp1 - yn - h * ( 2 * (tn + h) + exp(-(tn+h)*ynp1) ) } ynp1_res_prime_fn = function(ynp1, yn, tn, h){ 1 + h * (tn+h) * exp(-(tn+h)*ynp1) } h = 0.025 be_res = beuler(ynp1_res_fn, ynp1_res_prime_fn, t0, y0, t1, h) be_yhat_h1 = knn.reg( as.matrix( be_res\$t ), y=res\$y, test=as.matrix( x_grid ), k=1, algorithm='brute' ) h = 0.0125 be_res= beuler(ynp1_res_fn, ynp1_res_prime_fn, t0, y0, t1, h) be_yhat_h2 = knn.reg( as.matrix( be_res\$t ), y=res\$y, test=as.matrix( x_grid ), k=1, algorithm='brute' ) A = rbind( be_yhat_h1\$pred, be_yhat_h2\$pred ) rownames(A) = NULL R = cbind( c( 0.025, 0.0125 ), A ) colnames(R) = c( 'h', sprintf( 'y(%.2f)', x_grid ) ) print('Problem 10 (Backwards Euler):') print(R) ymin = min( c(e_yhat_h2\$pred, be_yhat_h2\$pred) ) ymax = max( c(e_yhat_h2\$pred, be_yhat_h2\$pred) ) plot(x_grid, e_yhat_h2\$pred, type='b', col='black', xlab='t', ylab='y(t) approximate', ylim=c(ymin, ymax)) lines(x_grid, be_yhat_h2\$pred, type='b', col='blue') grid() legend( 'topleft', c('Euler', 'Backwards Euler'), lty=c(1, 1), col=c('black', 'blue') )