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