# 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('beuler.R') source('runge_kutta.R') t0 = 0.0 y0 = 1.0 t1 = 0.1 # Where do we want to evaluate the solution: # x_grid = c( t1/2, t1 ) y_exact = function(t){ exp(-100*t) + t } y_exact_grid = y_exact( x_grid ) # Eulers method: # dy.dt = function(t, y){ -100*y + 100*t + 1 } h=0.1 while( TRUE ){ res = euler(dy.dt, h=h, start=t0, y0=y0, end=t1) e_yhat = knn.reg( as.matrix( res\$xs ), y=res\$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' ) error = e_yhat\$pred - y_exact_grid small_enough = abs(error)<0.0005 if( all(small_enough) ){ break } h = h/2 } print(sprintf('Euler h= %f', h)) # Backward Eulers method: # ynp1_res_fn = function(ynp1, yn, tn, h){ ynp1 - yn - h * ( -100 * ynp1 + 100*(tn + h) + 1 ) } ynp1_res_prime_fn = function(ynp1, yn, tn, h){ 1 + 100 * h } h = 0.1 while( TRUE ){ be_res = beuler(ynp1_res_fn, ynp1_res_prime_fn, t0, y0, t1, h) be_yhat = knn.reg( as.matrix( be_res\$t ), y=res\$y, test=as.matrix( x_grid ), k=1, algorithm='brute' ) error = be_yhat\$pred - y_exact_grid small_enough = abs(error)<0.0005 if( all(small_enough) ){ break } h = h/2 } print(sprintf('Backwards Euler h= %f', h)) # Runge-Kutta method: # h=0.1 while( TRUE ){ res = runge_kutta(dy.dt, h=h, start=t0, y0=y0, end=t1) e_yhat = knn.reg( as.matrix( res\$xs ), y=res\$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' ) error = e_yhat\$pred - y_exact_grid small_enough = abs(error)<0.0005 if( all(small_enough) ){ break } h = h/2 } print(sprintf('RK4 h= %f', h))