Added computational methods, including root finding methods and cubic spline interpolation.

Author: adam-rumpf

.gitignore

@@ -1,3 +1,4 @@
+_resources/
 datafiles/
 objects/
 options/

GameMakerScripts.yyp

@@ -23,6 +23,8 @@ Common array functions (such as searching and counting).
 
 Numerical algorithms for computational mathematics.
 
+* [`_cubic_spline_coefficients`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_cubic_spline_coefficients/_cubic_spline_coefficients.gml): Calculates vectors of coefficients for the natural cubic polynomial interpolating a given data set.
+* [`_cubic_spline_evaluate`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_cubic_spline_evaluate/_cubic_spline_evaluate.gml): Evaluates a piecewise cubic polynomial at a specific value given a given set of coefficient vectors and intervals.
 * [`_root_bisection`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_root_bisection/_root_bisection.gml): Finds a function root on a specified interval using the bisection method.
 * [`_root_newton`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_root_newton/_root_newton.gml): Finds a function root given a derivative and an initial guess using Newton's method.
 

README.md

@@ -0,0 +1,79 @@
+/// @func _cubic_spline_coefficients(x, y)
+/// @desc Calculates the coefficients for a natural cubic spline interpolating a set of points.
+/// @requires _tridiagonal_solve
+/// @param {real[]} x Sample inputs, as a 1D array in ascending order.
+/// @param {real[]} y Sample outputs, as a 1D array the same length as the input array.
+/// @return {real[][]} Set of coefficients to define each segment of the natural cubic spline (see documentation for details), undefined in case of error.
+/**
+ * For a data set of the form [x0, x1, ..., xn], [y0, y1, ..., yn], the natural cubic spline
+ * consists of n-1 piecewise polynomials of the form:
+ *     si(x) = ai (x-xi)^3 + bi (x-xi)^2 + ci (x-xi) + di
+ * on interval (xi, x(i+1)). The spline is twice differentiable with a second derivative
+ * that vanishes at x0 and xn.
+ *
+ * The return value of this function is a list of 4 lists, each of length n-1 and consisting
+ * of the coefficients of the form:
+ *     [
+ *         [a0, a1, ..., a(n-1)],
+ *         [b0, b1, ..., b(n-1)],
+ *         [c0, c1, ..., c(n-1)],
+ *         [d0, d1, ..., d(n-1)]
+ *     ]
+ */
+
+function _cubic_spline_coefficients(xx, yy)
+{
+	// Verify that argument lists have equal length
+	if (array_length(xx) != array_length(yy))
+		return undefined;
+	
+	// Get length of data lists
+	var n = array_length(xx);
+	
+	// Calculate interval differences in x and y
+	var dx, dy;
+	dx = array_create(n-1);
+	dy = array_create(n-1);
+	for (var i = 0; i < n-1; i++)
+	{
+		dx[i] = xx[i+1] - xx[i];
+		dy[i] = yy[i+1] - yy[i];
+	}
+	
+	// Define a tridiagonal system to calculate the quadratic coefficients
+	var upper, lower, diag, rhs;
+	upper = array_create(n-1, 0); // upper diagonal
+	lower = array_create(n-1, 0); // lower diagonal
+	diag = array_create(n, 1); // main diagonal
+	rhs = array_create(n, 0); // constant vector
+	for (var i = 1; i < n-1; i++)
+	{
+		lower[i-1] = dx[i-1];
+		upper[i] = dx[i];
+		diag[i] = 2*(dx[i-1] + dx[i]);
+		rhs[i] = 3*((dy[i]/dx[i]) - (dy[i-1]/dx[i-1]));
+	}
+	
+	// Solve for quadratic coefficients
+	var bb = _tridiagonal_solve(lower, diag, upper, rhs);
+	
+	// Substitute to solve for remaining coefficients, and trim long coefficient arrays
+	var a, b, c, d;
+	a = array_create(n-1);
+	b = array_create(n-1);
+	c = array_create(n-1);
+	d = array_create(n-1);
+	for (var i = 0; i < n-1; i++)
+	{
+		// Calculate cubic and linear coefficients
+		a[i] = (bb[i+1] - bb[i])/(3*dx[i]);
+		c[i] = (dy[i]/dx[i]) - (dx[i]/3)*(2*bb[i] + bb[i+1]);
+		
+		// Trim quadratic and constant coefficients (constants come directly from data)
+		b[i] = bb[i];
+		d[i] = yy[i];
+	}
+	
+	// Return all coefficient arrays
+	return [a, b, c, d];
+}

scripts/_cubic_spline_coefficients/_cubic_spline_coefficients.gml

@@ -0,0 +1,31 @@
+/// @func _cubic_spline_evaluate(a, b, c, d, x, t)
+/// @desc Evaluate a natural cubic spline at point t for coefficient vectors a, b, c, and d and input data vector x.
+/// @param {real[]} a Vector of cubic coefficients [a0, a1, ..., a(n-1)].
+/// @param {real[]} b Vector of quadratic coefficients [b0, b1, ..., b(n-1)].
+/// @param {real[]} c Vector of linear coefficients [c0, c1, ..., c(n-1)].
+/// @param {real[]} d Vector of contant coefficients [d0, d1, ..., d(n-1)].
+/// @param {real[]} x Vector of sample input values [x0, x1, ..., xn], in ascending order.
+/// @param {real} t Value at which to evaluate the spline.
+/// @return {real} Value of spline function evaluated at point t, depending on which interval if falls into, undefined in case of array size mismatch.
+/**
+ * For a data set of the form [x0, x1, ..., xn], [y0, y1, ..., yn], the cubic spline consists of a
+ * consists of n-1 piecewise polynomials of the form:
+ *     si(x) = ai (x-xi)^3 + bi (x-xi)^2 + ci (x-xi) + di
+ * on interval (xi, x(i+1)).
+ */
+
+function _cubic_spline_evaluate(a, b, c, d, xx, t)
+{
+	// Verify that all arguments have the appropriate length
+	var n = array_length(xx);
+	if ((array_length(a) != n-1) || (array_length(b) != n-1) || (array_length(c) != n-1) || (array_length(d) != n-1))
+		return undefined;
+	
+	// Find which subinterval the input value belongs to
+	var ii = 0; // index of subinterval
+	while ((t >= xx[ii+1]) && (ii < n-2))
+		ii++;
+	
+	// Evaluate the function using the corresponding coefficients
+	return a[ii]*power((t-xx[ii]), 3) + b[ii]*sqr(t-xx[ii]) + c[ii]*(t-xx[ii]) + d[ii];
+}

scripts/_cubic_spline_evaluate/_cubic_spline_evaluate.gml

@@ -2,7 +2,7 @@
 /// @desc Generates an array of equally-spaced values over a specified range with a specified step size.
 /// @param {real} start First value of array.
 /// @param {real} stop Final value of array (if the step size evenly divides the array; otherwise gives an upper bound to the final value).
-/// @param {real} [step=1] Step size. May be positive if start <= stop and negative if start >= stop.
+/// @param {real} [step=1] Step size. May be positive if start <= stop and negative otherwise.
 /// @return {real[]} Array of values beginning with start, incrementing by step, and ending with the last value that does not extend beyond stop.
 
 function _range(start, stop)

scripts/_range/_range.gml

@@ -1,5 +1,5 @@
 /// @func _root_bisection(f, a, b[, err[, cutoff]])
-/// @func Finds the root of a real-valued function on a specified interval using the bisection method.
+/// @func Finds the root of a real-valued, continuous function on a specified interval using the bisection method.
 /// @param {function} f Function (which maps reals to reals) to find the root of. The signs of f(a) and f(b) must be opposite.
 /// @param {real} a Lower bound of search interval.
 /// @param {real} b Upper bound of search interval.

scripts/_root_bisection/_root_bisection.gml

@@ -1,5 +1,5 @@
 /// @func _root_newton(f, fp, x0[, err[, cutoff]])
-/// @func Finds the root of a real-valued, differentiable function for an initial guess using Newton's method.
+/// @func Finds the root of a real-valued, differentiable function from an initial guess using Newton's method.
 /// @param {function} f Differentiable function (which maps reals to reals) to find the root of.
 /// @param {function} fp First derivative of function.
 /// @param {real} x0 Initial guess.