vanilla RK4 for ODE2, written in R (07Dec2011)

rk4-ode2.r

# vanilla RK4 written in R (07Dec2011)
#	evaluates a 2nd order ODE by splitting into 2- 1st order ODE 
#	and solved by Runge-Kutta Method
# Concepts and algorithm from:
#  http://en.wikipedia.org/wiki/Runge-Kutta_methods
#  http://www.cms.livjm.ac.uk/etchells/notes/ma200/2ndodes.pdf
#
#	f1(t,y) : substituted function from dz(y')/dt = F(z,y,t), z = y'
#	
#	rk4(f, h, tinit, yinit, limit): main RK4 routine
#		f - function defining the differential equation of
#		    the form dy/dt = f(t,y)
#		h - interval step
#		t/y/zinit - boundary values, where y(tinit) = yinit
#		limit - number of iterations the algorithm is done
#		 
#	TODO: fix documentation, instructions; verify using test cases
#         and compare with plots using WolframAlpha


f1 <- function(z,t,y) {
	value <- 3*cos(t)*z - sin(t) + exp(t) 
	value
	}

rk4 <- function(f1 , h, tinit, yinit, zinit, limit) {
	j <- c(0,0,0,0)
	y <- yinit
	t <- tinit
	z <- zinit
	tplot <- c()
	yplot <- c()
	zplot <- c()
	while (limit > 0) {
		
		z = z + (h/6)*(t + 2*(t + 0.5*h) + 2*(t + 0.5*h) + (t+h))
		
		j[1] <- h*f1(z,t,y)
		j[2] <- h*f1(z,t + 0.5*h, y + 0.5*j[1])
		j[3] <- h*f1(z,t + 0.5*h, y + 0.5*j[2])
		j[4] <- h*f1(z,t + h, y + j[3])
		y = y + (1./6)*(j[1] + 2*j[2] + 2*j[3] + j[4])
		t = t + h
		
		limit <- limit - 1
		tplot[limit] = t
		yplot[limit] = y
		zplot[limit] = z
		}
	plot(tplot,yplot)  ## plots y(t)  vs. t
	#plot(yplot,zplot) ## plots y'(t) vs y(t)
		
	}