source('runge_kutta.R') predictor_corrector = function(dy.dx=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 a predictive corrector method: # for (i in 5:nsteps){ x = start + (i-1)*h xs[i+1] = x + h # Take a 4th order Adams-Bashforth step: # y_predict = ys[i] + (h/24)*(55*fs[i] - 59*fs[i-1] + 37*fs[i-2] - 9*fs[i-3]) # Correct this with a 4th order Adams-Moulton step: # f_predict = dy.dx(x+h, y_predict) ys[i+1] = ys[i] + (h/24)*(9*f_predict + 19*fs[i] - 5*fs[i-1] + fs[i-2]) fs[i+1] = dy.dx(x + h, 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 } # # res = predictor_corrector(dy.dx, 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 #