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