Added Differentiate3 to Interpolation

mathnet · Sep 10, 2024 · 3cc5b65 · 3cc5b65
1 parent f196418
commit 3cc5b65
Show file tree

Hide file tree

Showing 10 changed files with 109 additions and 1 deletion.
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
@@ -576,7 +576,18 @@ public double Differentiate2(double t)
         {
             int k = LeftSegmentIndex(t);
             var x = t - _x[k];
-            return 2*_c2[k] + x*6*_c3[k];
+            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>

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>