# 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') source('runge_kutta.R') source('beuler.R') t0 = 0 y0 = 4 t1 = 5 # Where do we want to evaluate the solution: # t_grid = seq( t0, t1, by=0.05 ) # The true solution: # y_exact = function(t){ 0.25 * t^2 + 4 * exp(-10*t) } y_grid = y_exact(t_grid) # Eulers method: # dy.dt = function(t, y){ -10*y + 2.5*t^2 + 0.5*t } h_e=0.201 # Euler's method is unstable for this h #h_e=0.01 # Euler's method is stable for this h res = euler(dy.dt, h=h_e, start=t0, y0=y0, end=t1) e_yhat = knn.reg( as.matrix( res\$xs ), y=res\$ys, test=as.matrix( t_grid ), k=1, algorithm='brute') # Runge-Kutta: # h_rk=0.3 # RK is not stable for this h h_rk=0.1 # RK seems stable for this h res = runge_kutta(dy.dt, h=h_rk, start=t0, y0=y0, end=t1) rk_yhat = knn.reg( as.matrix( res\$xs ), y=res\$ys, test=as.matrix( t_grid ), k=1, algorithm='brute') # Backward Eulers method: # ynp1_res_fn = function(ynp1, yn, tn, h){ ynp1 - yn - h * ( -10*ynp1 + 2.5*(tn + h)^2 + 0.5*(tn + h) ) } ynp1_res_prime_fn = function(ynp1, yn, tn, h){ 1 + 10*h } h_be = 0.3 be_res = beuler(ynp1_res_fn, ynp1_res_prime_fn, t0, y0, t1, h_be) be_yhat = knn.reg( as.matrix( be_res\$t ), y=be_res\$y, test=as.matrix( t_grid ), k=1, algorithm='brute') #postscript("../../WriteUp/Graphics/Chapter8/chap_8_sect_6_prob_4_plots.eps", onefile=FALSE, horizontal=FALSE) all_ys = c( y_grid, e_yhat\$pred, rk_yhat\$pred, be_yhat\$pred ) y_min = min(all_ys) y_max = max(all_ys) plot(t_grid, y_grid,'l', col='green', xlab='t', ylab='y(t)', ylim=c(y_min, y_max)) lines(t_grid, e_yhat\$pred,'p', pch=20, col='red') lines(t_grid, rk_yhat\$pred,'l', col='black') lines(t_grid, be_yhat\$pred,'l', col='blue') legend('topleft', c('Exact Solution', sprintf('Euler (h=%.3f)', h_e), sprintf('RK (h= %.2f)', h_rk), sprintf('BEuler (h= %.2f)', h_be)), lty=c(1, 1, 1, 1), col=c('green','red','black','blue')) grid() #dev.off()