Skip to content

Commit 587b4c5

Browse files
Eduardo Rojascdrnet
authored andcommitted
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)
1 parent 5c077f7 commit 587b4c5

1 file changed

Lines changed: 70 additions & 0 deletions

File tree

src/Numerics/Interpolation/CubicSpline.cs

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -623,5 +623,75 @@ int LeftSegmentIndex(double t)
623623

624624
return Math.Min(Math.Max(index, 0), _x.Length - 2);
625625
}
626+
627+
/// <summary>
628+
/// Gets all the points (x values) where the derivative is 0
629+
/// </summary>
630+
/// <returns>An array of X values where the derivative of the spline is 0</returns>
631+
public double[] GetHorizontalPoints()
632+
{
633+
List<double> points = new List<double>();
634+
for (int index = 0; index < _x.Length-1; index++)
635+
{
636+
double a = 6 * _c3[index]; //derive ax^3 and multiply by 2
637+
double b = 2 * _c2[index]; //derive bx^2
638+
double c = _c1[index];//derive cx
639+
double d = b * b - 2 * a * c;
640+
//first check if a is 0, if so its a linear function, this happens with quadratic condition
641+
if (a.AlmostEqual(0))
642+
{
643+
double x = _x[index] - c / b;
644+
//check if the result is in the domain
645+
if (_x[index] <= x && x <= _x[index + 1]) points.Add(x);
646+
}
647+
else if (d.AlmostEqual(0))//its a quadratic with a single solution
648+
{
649+
double x = _x[index] - b / a;
650+
if (_x[index] <= x && x <= _x[index + 1]) points.Add(x);
651+
}
652+
else if (d > 0)//only has a solution if d is greater than 0
653+
{
654+
d = (double)System.Math.Sqrt(d);
655+
//apply quadratic equation
656+
double x1 = _x[index]+ (-b + d) / a;
657+
double x2 = _x[index]+ (-b - d) / a;
658+
//Add any solution points that fall within the domain to the list
659+
if ((_x[index] <= x1) && (x1 <= _x[index + 1])) points.Add(x1);
660+
if ((_x[index] <= x2) && (x2 <= _x[index + 1])) points.Add(x2);
661+
}
662+
}
663+
return points.ToArray();
664+
}
665+
666+
/// <summary>
667+
/// Returns the minimum value for the spline within its domain.
668+
/// </summary>
669+
/// <returns></returns>
670+
public double GetMaxPoint()
671+
{
672+
double max = double.MinValue;
673+
double x=0;
674+
//go through the functions points to check if one of them has a higher value
675+
foreach (double p in _x)
676+
{
677+
double y = Interpolate(p);
678+
if (y > max)
679+
{
680+
max = y;
681+
x = p;
682+
}
683+
}
684+
//go through the inflexion, local minimums and local maximums
685+
foreach (double p in GetHorizontalPoints())
686+
{
687+
double y = Interpolate(p);
688+
if (y > max)
689+
{
690+
max = y;
691+
x = p;
692+
}
693+
}
694+
return x;
695+
}
626696
}
627697
}

0 commit comments

Comments
 (0)