R/plot_Simpson.R

In alexjweil/integrateIt: Integrates Functions Numerically

#' Computes the quadratic interpolation of the provided
#' values.
#'
#' @param u  The first value to interpolate
#' @param fu The y-value corresponding to the first
#' @param v  The midpoint to interpolate on
#' @param fv The y-value corresponding to the midpoint value
#' @param w  The last value to interpolate
#' @param fw The y-value corresponding to the last value
#' @export
quadratic.interpolation = function(u, fu, v, fv, w, fw) {
	curve(fu * (x-v) * (x-w) / ((u-v) * (u-w)) +
		  fv * (x-u) * (x-w) / ((v-u) * (v-w)) +
		  fw * (x-u) * (x-v) / ((w-u) * (w-v)), u, w, add = TRUE, col="red")
}


#' Plots the parabolas created in a Simpson's rule
#' numerical integration.
#'
#' @param Simpson The Simpson object to plot.
#' @export
setMethod(f="plot",
		  signature="Simpson",
		  definition=function(x=NULL, y=x, ...){
		  	# The range of valid values.
		  	range = which(x@x>=x@start & x@x<=x@end)
		  	trimmed.x = x@x[range]
		  	trimmed.y = x@y[range]

		  	# Create the plot
		  	plot(trimmed.x, trimmed.y, pch=16, main="Numerical Integration (Simpson's Rule)", ylim=range(trimmed.y), ylab="y", xlab="x")

		  	# Map each index to its interpolation
		  	sapply(seq(1,length(x@y),2), function(i) {quadratic.interpolation(x@x[i], x@y[i], x@x[i+1], x@y[i+1], x@x[i+2], x@y[i+2])})

		  }
)