# Implements the backwards differentiation method using Newton iterations to solve for y_{n+1} in any the nonlinear function.
# One needs to pass in two functions
#
# dy.dx => to evaluate the right-hand-side of y' = f(t, y)
# dy.dx.dy => to evaluate the deriative of the right-hand-side of y' = f(t, y) with respect to y
#
# An implementation of Newtons method:
#
#source('../../../../Conte/Code/Chapter3/utils.R')
source('http://www.waxworksmath.com/Authors/A_F/Conte/Code/Chapter3/utils.R')

bd_newton_fn = function(ynp1, dy.dx, tn, h, ys, n){
    fnp1 = dy.dx(tn + h, ynp1)
    ynp1 - ( 48*ys[n] - 36*ys[n-1] + 16*ys[n-2] - 3*ys[n-3] + 12*h*fnp1 )/25
}

bd_newton_fn_prime = function(ynp1, dy.dx.dy, tn, h, ys, n){
    1 - (12/25) * h * dy.dx.dy(tn+h, ynp1)
}

backwards_differentiation = function(dy.dx=function(x,y){}, dy.dx.dy=function(x,y){}, h=1E-7, y0=1, start=0, end=1){

    nsteps = (end-start)/h
    xs = numeric(nsteps+1)
    ys = numeric(nsteps+1)
    fs = numeric(nsteps+1)
    xs[1] = start
    ys[1] = y0
    fs[1] = dy.dx(xs[1], ys[1])

    # Get the first 4 values from the Runge-Kutta method
    #
    rk = runge_kutta(dy.dx, h, y0, start, end=start+4*h)
    for (i in 2:5){
        xs[i] = start + (i-1)*h
        ys[i] = rk$ys[i]
        fs[i] = dy.dx(xs[i], ys[i])
    }

    # Get the next values using the Backwards differentiation formulas:
    #
    for (i in 5:nsteps){
        xs[i+1] = start + i*h

        root_fn = function(y){ bd_newton_fn(y, dy.dx, xs[i], h, ys, i) }
        root_fn_prime = function(y){ bd_newton_fn_prime(y, dy.dx.dy, xs[i], h, ys, i) }
        nm_res = newtons_method(ys[i], root_fn, root_fn_prime, xtol=1.e-3, ftol=1.e-3, ntol=10)
        ys[i+1] = nm_res[1]
        
        fs[i+1] = dy.dx(xs[i+1], ys[i+1])
    }
    list( xs=xs, ys=ys )
}
#
# An example of how to use use this function is:
#
# dy.dx = function(x,y) { 3*x - y + 8 }
# dy.dx.dy = function(x,y) { -1 }
#
# res = backwards_differentiation(dy.dx, dy.dx.dy, start=0, end=0.5, h=0.1, y0=3)
#
# res$ys
# 
# [1] 3.000000 3.490325 3.962538 4.418363 4.859359 5.286939
#