Browse Source

Added function go get points where cubic spline derivative is 0 (GetHorizontalPoints) and also to get the maximum value of the cubic spline function within its domain (GetMaxPoints)

pull/888/head
Eduardo Rojas 5 years ago
parent
commit
a89fcd09ac
  1. 70
      src/Numerics/Interpolation/CubicSpline.cs

70
src/Numerics/Interpolation/CubicSpline.cs

@ -623,5 +623,75 @@ namespace MathNet.Numerics.Interpolation
return Math.Min(Math.Max(index, 0), _x.Length - 2);
}
/// <summary>
/// Gets all the points (x values) where the derivative is 0
/// </summary>
/// <returns>An array of X values where the derivative of the spline is 0</returns>
public double[] GetHorizontalPoints()
{
List<double> points = new List<double>();
for (int index = 0; index < _x.Length-1; index++)
{
double a = 6 * _c3[index]; //derive ax^3 and multiply by 2
double b = 2 * _c2[index]; //derive bx^2
double c = _c1[index];//derive cx
double d = b * b - 2 * a * c;
//first check if a is 0, if so its a linear function, this happens with quadratic condition
if (a.AlmostEqual(0))
{
double x = _x[index] - c / b;
//check if the result is in the domain
if (_x[index] <= x && x <= _x[index + 1]) points.Add(x);
}
else if (d.AlmostEqual(0))//its a quadratic with a single solution
{
double x = _x[index] - b / a;
if (_x[index] <= x && x <= _x[index + 1]) points.Add(x);
}
else if (d > 0)//only has a solution if d is greater than 0
{
d = (double)System.Math.Sqrt(d);
//apply quadratic equation
double x1 = _x[index]+ (-b + d) / a;
double x2 = _x[index]+ (-b - d) / a;
//Add any solution points that fall within the domain to the list
if ((_x[index] <= x1) && (x1 <= _x[index + 1])) points.Add(x1);
if ((_x[index] <= x2) && (x2 <= _x[index + 1])) points.Add(x2);
}
}
return points.ToArray();
}
/// <summary>
/// Returns the minimum value for the spline within its domain.
/// </summary>
/// <returns></returns>
public double GetMaxPoint()
{
double max = double.MinValue;
double x=0;
//go through the functions points to check if one of them has a higher value
foreach (double p in _x)
{
double y = Interpolate(p);
if (y > max)
{
max = y;
x = p;
}
}
//go through the inflexion, local minimums and local maximums
foreach (double p in GetHorizontalPoints())
{
double y = Interpolate(p);
if (y > max)
{
max = y;
x = p;
}
}
return x;
}
}
}

Loading…
Cancel
Save