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