Added Differentiate3 to Interpolation

osswaldo · osswaldo · commit 243d51668621 · 2024-09-24T14:28:16.000+02:00
diff --git a/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs b/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs
@@ -269,6 +269,20 @@ public void CheckNaturalSplineMinMaxValuesPerformance()
             }
         }
 
+        /// <summary>
+        /// Verifies that the 3rd derivative matches the given value at all the provided sample points.
+        /// </summary>
+        [TestCase(-10, -8.1428571428571423)]
+        [TestCase(-5, -8.1428571428571423)]
+        [TestCase(0, -10.714285714285715)]
+        [TestCase(5, 2.1428571428571423)]
+        [TestCase(10, 2.1428571428571423)]
+        public void ThirdDerivative(double x, double expected)
+        {
+            IInterpolation it = CubicSpline.InterpolateNatural(_t, _y);
+            Assert.AreEqual(expected, it.Differentiate3(x));
+        }
+
         /// <summary>
         /// Generates a set of points representing an oscilating decaying function
         /// </summary>
diff --git a/src/Numerics.Tests/InterpolationTests/LinearSplineTest.cs b/src/Numerics.Tests/InterpolationTests/LinearSplineTest.cs
@@ -28,6 +28,7 @@
 // </copyright>
 
 using MathNet.Numerics.Interpolation;
+using MathNet.Numerics.Random;
 using NUnit.Framework;
 
 namespace MathNet.Numerics.Tests.InterpolationTests
@@ -54,6 +55,16 @@ public void FirstDerivative()
             Assert.That(ip.Differentiate(3.0), Is.EqualTo(1.0));
         }
 
+        /// <summary>
+        /// Verifies that the 3rd derivative matches the given value at all the provided sample points.
+        /// </summary>
+        public void ThirdDerivative(double x, double expected)
+        {
+            var rnd = new SystemRandomSource(10);
+            IInterpolation it = LinearSpline.Interpolate(_t, _y);
+            Assert.AreEqual(0, it.Differentiate3(rnd.NextDouble()));
+        }
+
         [Test]
         public void DefiniteIntegral()
         {
diff --git a/src/Numerics.Tests/InterpolationTests/LogLinearTest.cs b/src/Numerics.Tests/InterpolationTests/LogLinearTest.cs
@@ -0,0 +1,56 @@
+// <copyright file="CubicSplineTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using MathNet.Numerics.Interpolation;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.Tests.InterpolationTests
+{
+    [TestFixture, Category("Interpolation")]
+    public class LogLinearTest
+    {
+        readonly double[] _t = { 1.0, 2.0, 3.0, 4.0, 5.0 };
+        readonly double[] _y = { 1.0, 4.0, 9.0, 16.0, 25.0 };
+
+        /// <summary>
+        /// Verifies that the 3rd derivative matches the given value at all the provided sample points.
+        /// </summary>
+        [TestCase(0, 0.66604930397785889)]
+        [TestCase(1, 2.6641972159114355)]
+        [TestCase(2, 2.1330961922715468)]
+        [TestCase(3, 1.7142371101090121)]
+        [TestCase(4, 1.4222075877044038d)]
+        [TestCase(5, 2.2221993557881312d)]
+        public void ThirdDerivative(double x, double expected)
+        {
+            IInterpolation it = LogLinear.Interpolate(_t, _y);
+            Assert.AreEqual(expected, it.Differentiate3(x));
+        }
+    }
+}
diff --git a/src/Numerics.Tests/InterpolationTests/QuadraticSplineTest.cs b/src/Numerics.Tests/InterpolationTests/QuadraticSplineTest.cs
@@ -0,0 +1,56 @@
+// <copyright file="QuadraticSplineTest.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using MathNet.Numerics.Interpolation;
+using MathNet.Numerics.Random;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.Tests.InterpolationTests
+{
+    [TestFixture, Category("Interpolation")]
+    public class QuadraticSplineTest
+    {
+        readonly double[] _x = { -20.0, -10.0, -5.0, 0, 5.0, 10.0, 20.0 };
+        readonly double[] _c0 = { -3.0, -2.0, 1.0, 0.0, 1.0, 2.0 };
+        readonly double[] _c1 = { -2.0, -1.0, 0.0, 1.0, 2.0, 3.0 };
+        readonly double[] _c2 = { -1.0, 0.0, 1.0, 2.0, 3.0, 4.0 };
+
+        /// <summary>
+        /// Verifies that the 3rd derivative matches the given value at all the provided sample points.
+        /// </summary>
+        [Test]
+        public void ThirdDerivative()
+        {
+            var rnd = new SystemRandomSource(10);
+            IInterpolation it = new QuadraticSpline(_x, _c0, _c1, _c2);
+
+            Assert.AreEqual(0, it.Differentiate3(rnd.NextDouble()));
+        }
+    }
+}
diff --git a/src/Numerics.Tests/InterpolationTests/StepInterpolationTest.cs b/src/Numerics.Tests/InterpolationTests/StepInterpolationTest.cs
@@ -1,4 +1,4 @@
-﻿// <copyright file="LinearSplineTest.cs" company="Math.NET">
+﻿// <copyright file="StepInterpolationTest.cs" company="Math.NET">
 // Math.NET Numerics, part of the Math.NET Project
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
@@ -29,6 +29,7 @@
 
 using System;
 using MathNet.Numerics.Interpolation;
+using MathNet.Numerics.Random;
 using NUnit.Framework;
 
 namespace MathNet.Numerics.Tests.InterpolationTests
@@ -56,6 +57,14 @@ public void FirstDerivative()
             Assert.That(ip.Differentiate(4.0), Is.EqualTo(0.0));
         }
 
+        public void ThirdDerivative()
+        {
+            var rnd = new SystemRandomSource(10);
+            var x = rnd.NextDouble();
+            IInterpolation ip = new StepInterpolation(_t, _y);
+            Assert.That(ip.Differentiate(x), Is.EqualTo(ip.Differentiate(x)));
+        }
+
         [Test]
         public void DefiniteIntegral()
         {
diff --git a/src/Numerics/Interpolation/Barycentric.cs b/src/Numerics/Interpolation/Barycentric.cs
@@ -334,6 +334,13 @@ public double Interpolate(double t)
         /// <returns>Interpolated second derivative at point t.</returns>
         double IInterpolation.Differentiate2(double t) => throw new NotSupportedException();
 
+        /// <summary>
+        /// Differentiate three times at point t. NOT SUPPORTED.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated third derivative at point t.</returns>
+        double IInterpolation.Differentiate3(double t) => throw new NotSupportedException();
+
         /// <summary>
         /// Indefinite integral at point t. NOT SUPPORTED.
         /// </summary>
diff --git a/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs b/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs
@@ -182,6 +182,13 @@ public double Interpolate(double t)
         /// <returns>Interpolated second derivative at point t.</returns>
         double IInterpolation.Differentiate2(double t) => throw new NotSupportedException();
 
+        /// <summary>
+        /// Differentiate three times at point t. NOT SUPPORTED.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated third derivative at point t.</returns>
+        double IInterpolation.Differentiate3(double t) => throw new NotSupportedException();
+
         /// <summary>
         /// Indefinite integral at point t. NOT SUPPORTED.
         /// </summary>
diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs
@@ -579,6 +579,17 @@ public double Differentiate2(double t)
             return 2*_c2[k] + x*6*_c3[k];
         }
 
+        /// <summary>
+        /// Differentiate three times at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated third derivative at point t.</returns>
+        public double Differentiate3(double t)
+        {
+            int k = LeftSegmentIndex(t);
+            return 6*_c3[k];
+        }
+
         /// <summary>
         /// Indefinite integral at point t.
         /// </summary>
diff --git a/src/Numerics/Interpolation/IInterpolation.cs b/src/Numerics/Interpolation/IInterpolation.cs
@@ -65,6 +65,13 @@ public interface IInterpolation
         /// <returns>Interpolated second derivative at point t.</returns>
         double Differentiate2(double t);
 
+        /// <summary>
+        /// Differentiate three times at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated third derivative at point t.</returns>
+        double Differentiate3(double t);
+
         /// <summary>
         /// Indefinite integral at point t.
         /// </summary>
diff --git a/src/Numerics/Interpolation/LinearSpline.cs b/src/Numerics/Interpolation/LinearSpline.cs
@@ -152,6 +152,13 @@ public double Differentiate(double t)
         /// <returns>Interpolated second derivative at point t.</returns>
         public double Differentiate2(double t) => 0d;
 
+        /// <summary>
+        /// Differentiate three times at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated third derivative at point t.</returns>
+        public double Differentiate3(double t) => 0d;
+
         /// <summary>
         /// Indefinite integral at point t.
         /// </summary>
diff --git a/src/Numerics/Interpolation/LogLinear.cs b/src/Numerics/Interpolation/LogLinear.cs
@@ -152,6 +152,24 @@ public double Differentiate2(double t)
             return secondDerivative;
         }
 
+        /// <summary>
+        /// Differentiate three times at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated third derivative at point t.</returns>
+        public double Differentiate3(double t)
+        {
+            var linearFirstDerivative = _spline.Differentiate(t);
+            var linearSecondDerivative = _spline.Differentiate2(t);
+            var linearThirdDerivative = _spline.Differentiate3(t);
+
+            var thirdDerivative = Differentiate2(t) * linearFirstDerivative +
+                                  2 * Differentiate(t) * linearSecondDerivative +
+                                  Interpolate(t) * linearThirdDerivative;
+
+            return thirdDerivative;
+        }
+
         /// <summary>
         /// Indefinite integral at point t.
         /// </summary>
diff --git a/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs b/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs
@@ -197,6 +197,36 @@ public double Differentiate2(double t)
             return ddx[0];
         }
 
+        /// <summary>
+        /// Differentiate three times at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated third derivative at point t.</returns>
+        public double Differentiate3(double t)
+        {
+            var x = new double[_y.Length];
+            var dx = new double[_y.Length];
+            var ddx = new double[_y.Length];
+            var dddx = new double[_y.Length];
+            _y.CopyTo(x, 0);
+
+            for (int level = 1; level < x.Length; level++)
+            {
+                for (int i = 0; i < x.Length - level; i++)
+                {
+                    double hp = t - _x[i + level];
+                    double ho = _x[i] - t;
+                    double den = _x[i] - _x[i + level];
+                    dddx[i] = ((hp * dddx[i]) + (3 * ddx[i]) + (ho * dddx[i + 1]) - (3 * ddx[i + 1])) / den;
+                    ddx[i] = ((hp * ddx[i]) +  (2 * dx[i])  + (ho * ddx[i + 1]) -  (2 * dx[i + 1])) / den;
+                    dx[i] = ((hp * dx[i]) + x[i]  + (ho * dx[i + 1]) - x[i + 1]) / den;
+                    x[i] = ((hp * x[i]) + (ho * x[i + 1])) / den;
+                }
+            }
+
+            return dddx[0];
+        }
+
         /// <summary>
         /// Indefinite integral at point t. NOT SUPPORTED.
         /// </summary>
diff --git a/src/Numerics/Interpolation/QuadraticSpline.cs b/src/Numerics/Interpolation/QuadraticSpline.cs
@@ -110,6 +110,13 @@ public double Differentiate2(double t)
             return 2*_c2[k];
         }
 
+        /// <summary>
+        /// Differentiate three times at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated third derivative at point t.</returns>
+        public double Differentiate3(double t) => 0d;
+
         /// <summary>
         /// Indefinite integral at point t.
         /// </summary>
diff --git a/src/Numerics/Interpolation/StepInterpolation.cs b/src/Numerics/Interpolation/StepInterpolation.cs
@@ -139,6 +139,13 @@ public double Differentiate(double t)
         /// <returns>Interpolated second derivative at point t.</returns>
         public double Differentiate2(double t) => Differentiate(t);
 
+        /// <summary>
+        /// Differentiate three times at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated third derivative at point t.</returns>
+        public double Differentiate3(double t) => Differentiate2(t);
+
         /// <summary>
         /// Indefinite integral at point t.
         /// </summary>
diff --git a/src/Numerics/Interpolation/TransformedInterpolation.cs b/src/Numerics/Interpolation/TransformedInterpolation.cs
@@ -146,6 +146,13 @@ public static TransformedInterpolation Interpolate(
         /// <returns>Interpolated second derivative at point t.</returns>
         double IInterpolation.Differentiate2(double t) => throw new NotSupportedException();
 
+        /// <summary>
+        /// Differentiate three times at point t. NOT SUPPORTED.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated third derivative at point t.</returns>
+        double IInterpolation.Differentiate3(double t) => throw new NotSupportedException();
+
         /// <summary>
         /// Indefinite integral at point t. NOT SUPPORTED.
         /// </summary>
