diff --git a/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs b/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs index 2ba52584..bdeb210d 100644 --- a/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs +++ b/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs @@ -39,7 +39,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests { readonly double[] _t = { -2.0, -1.0, 0.0, 1.0, 2.0 }; readonly double[] _y = { 1.0, 2.0, -1.0, 0.0, 1.0 }; - + //test data for min max values + readonly double[] _x = { -4, -3, -2, -1, 0, 1, 2, 3, 4 }; + readonly double[] _z = { -7, 2, 5, 0, -3, -1, -4, 0, 6 }; + /// /// Verifies that the interpolation matches the given value at all the provided sample points. /// @@ -202,5 +205,37 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests CollectionAssert.DoesNotContain(yipol, Double.NaN); } + + /// + /// Tests that the cubic spline returns the correct t values where the derivative is 0 + /// + [Test] + public void NaturalSplineGetHorizontalDerivativeTValues() + { + CubicSpline it = CubicSpline.InterpolateBoundaries(_x, _z, SplineBoundaryCondition.SecondDerivative, -4.0, SplineBoundaryCondition.SecondDerivative, 4.0); + var horizontalDerivatives = it.GetHorizontalDerivativeTValues(); + //readonly double[] _x = { -4, -3, -2, -1, 0, 1, 2, 3, 4 }; + //readonly double[] _z = { -7, 2, 5, 0, -3, -1, -4, 0, 6 }; + Assert.AreEqual(4, horizontalDerivatives.Length, "Incorrect number of points with derivative value equal to 0"); + Assert.IsTrue(horizontalDerivatives[0]>=-3 && horizontalDerivatives[0] <= -2,"Spline returns wrong t value: "+horizontalDerivatives[0]+" for first point"); + Assert.IsTrue(horizontalDerivatives[1] >= -1 && horizontalDerivatives[1] <= 0, "Spline returns wrong t value: " + horizontalDerivatives[1] + " for second point"); + Assert.IsTrue(horizontalDerivatives[2] >= 0 && horizontalDerivatives[2] <= 1, "Spline returns wrong t value: " + horizontalDerivatives[2] + " for third point"); + Assert.IsTrue(horizontalDerivatives[3] >= 2 && horizontalDerivatives[3] <= 3, "Spline returns wrong t value: " + horizontalDerivatives[3] + " for fourth point"); + Console.WriteLine("GetHorizontalDerivativeTValues checked out ok for cubic spline."); + } + + /// + /// Tests that the min and max values for the natural spline are correct + /// + [Test] + public void NaturalSplineGetMinMaxTvalues() + { + CubicSpline it = CubicSpline.InterpolateBoundaries(_x, _z, SplineBoundaryCondition.SecondDerivative, -4.0, SplineBoundaryCondition.SecondDerivative, 4.0); + var minMax = it.GetMinMaxTValues(); + Assert.AreEqual(-4, minMax.Item1, "Spline returns wrong t value for global minimum"); + Assert.AreEqual(4, minMax.Item2, "Spline returns wrong t value for global maximum"); + Console.WriteLine("GetMinMaxTValues checked out ok for cubic spline."); + + } } } diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index 1830da86..2d00bce9 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics @@ -623,5 +623,87 @@ namespace MathNet.Numerics.Interpolation return Math.Min(Math.Max(index, 0), _x.Length - 2); } + + /// + /// Gets all the t values where the derivative is 0 + /// + /// An array of t values (in the domain of teh function) where the derivative of the spline is 0 + public double[] GetHorizontalDerivativeTValues() + { + List points = new List(); + 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(); + } + + /// + /// Returns the t values in the domain of the spline for which it takes the minimum and maximum value. + /// + /// A tuple containing the t value for which the spline is minimum in the first component and maximum in the second component + public Tuple GetMinMaxTValues() + { + double max = double.MinValue; + double min = double.MaxValue; + double minT = 0; + double maxT = 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; + maxT = p; + } + if (y < min) + { + min = y; + minT = p; + } + } + //go through the inflexion, local minimums and local maximums + foreach (double p in GetHorizontalDerivativeTValues()) + { + double y = Interpolate(p); + if (y > max) + { + max = y; + maxT = p; + } + if (y < min) + { + min = y; + minT = p; + } + } + return new Tuple(minT, maxT); + } } }