# 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))