In this post I will develop a function on the Polynomial class that takes three points, fits a second order polynomial (a parabola) to them, and then attempts to find the minimum value of the polynomial.
This function works by taking three points as arguments, it then calls the previously developed Solve Polynomial Coefficients function to fit a second order polynomial to them, and then tries to find the minimum of that polynomial.
This is done by noting that a second order polynomial can be expressed as
y = a0 + a1x + a2x2
y’ = a1 +2a2x
y’’ = 2a2
At the minimum, the first derivative is equal to 0 and the second derivative is greater than 0.
Solving for the first derivative equaling zero gives
xm = -a1 / 2a2
If the second derivative is greater than 0, then the function returns the xm value from above, otherwise it returns NaN.
The full code for the Polynomial class to this point is given below.
public class Polynomial
{
public static double[] SolvePolyCoefs(
double[] x,
double[] y)
{
Debug.Assert(x.Length > 0);
Debug.Assert(x.Length == y.Length);
Matrix m = new Matrix(x.Length, x.Length);
Parallel.For(0, x.Length, i =>
{
m[i, 0] = 1.0;
for (int j = 1; j < x.Length; j++)
{
m[i, j] = m[i, j - 1] * x[i];
}
}
);
return m.SolveFor(y);
}
public static double ParabolicMin(double[] x, double[] y)
{
double result = double.NaN;
Debug.Assert(x.Length == 3);
Debug.Assert(y.Length == 3);
double[] values = SolvePolyCoefs(x, y);
if (values[2] > 0)
{
result = -values[1] / (2.0 * values[2]);
}
return result;
}
public static double[] Regress(double[,] z, double[] y)
{
//y=a0 z1 + a1 z1 +a2 z2 + a3 z3 +...
//Z is the functional values.
//Z index 0 is a row, the variables go across index 1.
//Y is the summed value.
//returns the coefficients.
Debug.Assert(z != null && y != null);
Debug.Assert(z.GetLength(0) == y.GetLength(0));
Matrix zMatrix = z;
Matrix zTransposeMatrix = zMatrix.Transpose();
Matrix leftHandSide = zTransposeMatrix * zMatrix;
Matrix rightHandSide = zTransposeMatrix * y;
Matrix coefsMatrix = leftHandSide.SolveFor(rightHandSide);
return coefsMatrix;
}
public static double[] RegressPolynomial(double[] x, double[] y, int polynomialOrder)
{
double[,] z = new double[y.Length, polynomialOrder + 1];
for (int i = 0; i < y.Length; i++)
{
z[i, 0] = 1.0;
for (int j = 1; j < polynomialOrder + 1; j++)
{
z[i, j] = Math.Pow(x[i], j);
}
}
double[] coefs = Polynomial.Regress(z, y);
return coefs;
}
public static double EvaluatePolynomial(double[] coefs, double x)
{
double result = coefs[0];
for (int i = 1; i < coefs.Length; i++)
{
result += coefs[i] * Math.Pow(x, i);
}
return result;
}
}
Subscribe