# I use the function knn.reg for simple comparisons below:
#
if( !require('FNN') ){
    install.packages("FNN", dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(FNN)

source('euler.R')


# Problem 1:
#
x_grid = seq( 0, 0.4, by=0.1 ) # where do we want to evaluate our solution

dy.dx = function(x, y){ 3 + x - y }

h=0.1
res = euler(dy.dx, h=h, start=0, y0=1, end=0.4)
yhat_h1 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

h=0.05
res = euler(dy.dx, h=h, start=0, y0=1, end=0.4)
yhat_h2 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

h=0.025
res = euler(dy.dx, h=h, start=0, y0=1, end=0.4)
yhat_h3 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

# compute with the exact solution:
#
h = 0.0
y_exact = x_grid + 2 - exp(-x_grid)

A = rbind( yhat_h1$pred, yhat_h2$pred, yhat_h3$pred, y_exact )
rownames(A) = NULL
R = cbind( c( 0.1, 0.05, 0.025, 0 ), A )
colnames(R) = c('h', 'y(0)', 'y(0.1)', 'y(0.2)', 'y(0.3)', 'y(0.4)')
print('Problem 1:')
print(R)


# Problem 2:
#
x_grid = seq( 0, 0.4, by=0.1 ) # where do we want to evaluate our solution

dy.dx = function(x, y){ 2*y - 1 }

h=0.1
res = euler(dy.dx, h=h, start=0, y0=1, end=0.4)
yhat_h1 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

h=0.05
res = euler(dy.dx, h=h, start=0, y0=1, end=0.4)
yhat_h2 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

h=0.025
res = euler(dy.dx, h=h, start=0, y0=1, end=0.4)
yhat_h3 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

# compute with the exact solution:
#
h = 0.0
y_exact = (1.+exp(2*x_grid))/2

A = rbind( yhat_h1$pred, yhat_h2$pred, yhat_h3$pred, y_exact )
rownames(A) = NULL
R = cbind( c( 0.1, 0.05, 0.025, 0 ), A )
colnames(R) = c('h', 'y(0)', 'y(0.1)', 'y(0.2)', 'y(0.3)', 'y(0.4)')
print('Problem 2:')
print(R)



# Problem 3:
#
x_grid = seq( 0, 0.4, by=0.1 ) # where do we want to evaluate our solution

dy.dx = function(x, y){ 0.5 - x + 2*y }

h=0.1
res = euler(dy.dx, h=h, start=0, y0=1, end=0.4)
yhat_h1 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

h=0.05
res = euler(dy.dx, h=h, start=0, y0=1, end=0.4)
yhat_h2 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

h=0.025
res = euler(dy.dx, h=h, start=0, y0=1, end=0.4)
yhat_h3 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

# compute with the exact solution:
#
h = 0.0
y_exact = x_grid/2 + exp(2*x_grid)

A = rbind( yhat_h1$pred, yhat_h2$pred, yhat_h3$pred, y_exact )
rownames(A) = NULL
R = cbind( c( 0.1, 0.05, 0.025, 0 ), A )
colnames(R) = c('h', 'y(0)', 'y(0.1)', 'y(0.2)', 'y(0.3)', 'y(0.4)')
print('Problem 3:')
print(R)



# Problem 4:
#
x_grid = seq( 0, 0.4, by=0.1 ) # where do we want to evaluate our solution

dy.dx = function(x, y){ 3*cos(x) - 2*y }

h=0.1
res = euler(dy.dx, h=h, start=0, y0=0, end=0.4)
yhat_h1 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

h=0.05
res = euler(dy.dx, h=h, start=0, y0=0, end=0.4)
yhat_h2 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

h=0.025
res = euler(dy.dx, h=h, start=0, y0=0, end=0.4)
yhat_h3 = knn.reg( as.matrix( res$xs ), y=res$ys, test=as.matrix( x_grid ), k=1, algorithm='brute' )

# compute with the exact solution:
#
h = 0.0
y_exact = (3/5) * ( sin(x_grid) + 2*cos(x_grid) ) - (6/5) * exp(-2*x_grid)

A = rbind( yhat_h1$pred, yhat_h2$pred, yhat_h3$pred, y_exact )
rownames(A) = NULL
R = cbind( c( 0.1, 0.05, 0.025, 0 ), A )
colnames(R) = c('h', 'y(0)', 'y(0.1)', 'y(0.2)', 'y(0.3)', 'y(0.4)')
print('Problem 4:')
print(R)