# Using a simple exponential function as an example, this script demonstrates the simple Euler 
# method, Runge-Kutta 2, and Runge-Kutta 4

# Consider an example of exponential population growth where the growth rate is 0.23 (23%), with an
# initial population of 100. The integrated rate equation is P=100*exp(0.23*t). We can evaluate this
# equation at any time we like and plot it's exact solution. But in many cases we can't derive such 
# an equation we therefore have to estimate the value of the population by taking repeated time 
# steps and evaluating an expression that approximates the value of the population after each step.

# Setup the integrated growth equation and evaluate it at a number of time points for comparison to 
# our numerical results
tp = seq(0, 10, 0.005)	# vector of time points at which the equation will be evaluated
P0 = 100				# initial population
k = 0.10				# growth rate
Pa = P0*exp(k*tp)		# Population analytical solution (vector of the same length as tp)
plot(tp, Pa, type="l")	# Plot our exact solution


# Simple Euler method - predict dependent variable at time (t), by using the value at time (t-dt) and 
# the rate of change at time (t-dt)
plot(tp, Pa, type="l", lwd=3)	# plot of the analytical solution for population growth
points(0, 100, pch=19, cex=2)	# point at initial population at time zero 
abline(P0, (k*P0), col="blue")	# tangent line at point above
dx = 5							# time at which we will evaluate population using Euler method
dP1 = k*P0*dx					# estimated change in population over interval dt 
P1 = P0 + dP1
lines(c(0,dx), c(P0, P0), lty="dashed")	# dashed line showing tim einterval (dt)
lines(c(dx,dx), c(P0, P1), lty="dashed")	# dashed line showing cahnge in population
points(dx, P1, pch=19, cex=2, col="blue")	# point that is Euler estimate of pop at (t+dt)
# redraw the plot with a line segment and point showing our estimate so far
plot(tp, Pa, type="l", lwd=3)	# analytical solution
lines(c(0,dx), c(P0, P1))		# line segment showing estimated change in population
points(dx, P1, pch=19, cex=2)	# estimaed population after dt
b = P1 - (k*P1*dx)				# y-intercept for tangent line through first estimated population
abline(b, (k*P1), col="blue")	# tangent line at point above
dP2 = k*P1*dx					# estimated change in population over interval dt 
P2 = P1 + dP2
lines(c(dx,2*dx), c(P1, P1), lty="dashed")	# dashed line showing tim einterval (dt)
lines(c(2*dx,2*dx), c(P1, P2), lty="dashed")	# dashed line showing cahnge in population
points(2*dx, P2, pch=19, cex=2, col="blue")	# point that is Euler estimate of pop at (t+dt)

plot(tp, Pa, type="l", lwd=3)	# analytical solution
lines(c(0,dx), c(P0, P1))		# line segment showing estimated change in population
points(dx, P1, pch=19, cex=2)
lines(c(dx, 2*dx), c(P1, P2))		# line segment showing estimated change in population
points(2*dx, P2, pch=19, cex=2)


y = function(yinit,x){
	y = c(yinit)			# Initial value passed to function
	
	# Divide x interval into n steps of length dx
	# Doing this here, inside the function, will allow us to evaluate y easily over different 
	# intervals by passing different vectors of time points to the function
	n = length(x)			# calculate how many times "y" will need to be evaluated
	dx=(x[n]-x[1])/(n-1)	# calculate the size of the time step, this assumes that all time-steps
							# are the same size
	
	# Set-up a loop to calculate "y" at each time point. We don't need to calculate it at the first 
	# time point, t=0, because we provided an initial value, so we start at the second time point
	for (i in 2:n) {
		dy = f(x[i-1])*dx	# the function "f" is defined below for the differential rate equation
		y[i] = y[i-1]+dy	# calculate "y" as HAVE = HAD + CHANGE
		}
	return(y)
}

f = function(x) {x^2}	# This is the differntial rate equation (the slope, or first derivative,  
						# at time xi), that describes the rate of change in y given a value of x

# Create a vector of time points
xinit = 0 	# start time 
xmax = 2	# end time
dx = .1		# time-step
x = seq(xinit,xmax,dx)	# vector of time points

yinit = 0	# initial value of the dependent variable
x1 = x		# I don't know why this assignment was necessary

y1 = y(yinit,x)	# call our function "y" with it's needed parameters, and assign the result to "y1"

plot(x1,y1, type="l") # Plot the result as a line plot, (again, I don't know why "x1" is necessary)
lines(spline(x, x^3/3, n=201), col="red")