vector_predictor_corrector = function(vector_f=function(t,x,y){}, h=1E-7, xy0=c(0, 0), start=0, end=1){ # # xy0 = initial condition [x(t_0), y(t_0)] # start = t_0 = start of time interation # end = t_1 = end of time integration # h = stepsize # nsteps = (end-start)/h ts = rep(NA, nsteps+1) xys = matrix(rep(NA, 2*(nsteps+1)), nrow=2, ncol=nsteps+1) # (2=n_states) x n_timesteps fs = matrix(rep(NA, 2*(nsteps+1)), nrow=2, ncol=nsteps+1) # (2=n_states) x n_timesteps ts[1] = start xys[1, 1] = xy0[1] xys[2, 1] = xy0[2] fs[, 1] = vector_f(ts[1], xys[1, 1], xys[2, 1]) # Get the first 4 values from Euler's method (better methods could be used): # for (i in 2:5){ ts[i] = start + (i-1)*h f = vector_f(ts[i-1], xys[1, i-1], xys[2, i-1]) xys[1, i] = xys[1, i-1] + h * f[1] xys[2, i] = xys[2, i-1] + h * f[2] fs[, i] = vector_f(ts[i], xys[1, i], xys[2, i]) } # Get the next values using a predictive corrector method: # for (i in 5:nsteps){ x = start + (i-1)*h ts[i+1] = x + h f = vector_f(ts[i-1], xys[1, i-1], xys[2, i-1]) # this is f_n # Take a 4th order Adams-Bashforth step: # xy_predict = xys[, 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 = vector_f(ts[i+1], xy_predict[1], xy_predict[2]) xys[, i+1] = xys[, i] + (h/24)*(9*f_predict + 19*fs[, i] - 5*fs[, i-1] + fs[, i-2]) fs[, i+1] = vector_f(ts[i+1], xys[1, i+1], xys[2, i+1]) } list( ts=ts, xys=xys, fs=fs ) }