From 243d5166862167dd2c92f2713333638dc3619095 Mon Sep 17 00:00:00 2001 From: osswaldo Date: Tue, 10 Sep 2024 16:28:39 +0200 Subject: [PATCH] Added Differentiate3 to Interpolation --- .../InterpolationTests/CubicSplineTest.cs | 14 +++++ .../InterpolationTests/LinearSplineTest.cs | 11 ++++ .../InterpolationTests/LogLinearTest.cs | 56 +++++++++++++++++++ .../InterpolationTests/QuadraticSplineTest.cs | 56 +++++++++++++++++++ .../StepInterpolationTest.cs | 11 +++- src/Numerics/Interpolation/Barycentric.cs | 7 +++ .../BulirschStoerRationalInterpolation.cs | 7 +++ src/Numerics/Interpolation/CubicSpline.cs | 11 ++++ src/Numerics/Interpolation/IInterpolation.cs | 7 +++ src/Numerics/Interpolation/LinearSpline.cs | 7 +++ src/Numerics/Interpolation/LogLinear.cs | 18 ++++++ .../NevillePolynomialInterpolation.cs | 30 ++++++++++ src/Numerics/Interpolation/QuadraticSpline.cs | 7 +++ .../Interpolation/StepInterpolation.cs | 7 +++ .../Interpolation/TransformedInterpolation.cs | 7 +++ 15 files changed, 255 insertions(+), 1 deletion(-) create mode 100644 src/Numerics.Tests/InterpolationTests/LogLinearTest.cs create mode 100644 src/Numerics.Tests/InterpolationTests/QuadraticSplineTest.cs diff --git a/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs b/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs index 4eb636366..1499a1fc1 100644 --- a/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs +++ b/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs @@ -269,6 +269,20 @@ public void CheckNaturalSplineMinMaxValuesPerformance() } } + /// + /// Verifies that the 3rd derivative matches the given value at all the provided sample points. + /// + [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)); + } + /// /// Generates a set of points representing an oscilating decaying function /// diff --git a/src/Numerics.Tests/InterpolationTests/LinearSplineTest.cs b/src/Numerics.Tests/InterpolationTests/LinearSplineTest.cs index 79980ee31..576c9fa66 100644 --- a/src/Numerics.Tests/InterpolationTests/LinearSplineTest.cs +++ b/src/Numerics.Tests/InterpolationTests/LinearSplineTest.cs @@ -28,6 +28,7 @@ // 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)); } + /// + /// Verifies that the 3rd derivative matches the given value at all the provided sample points. + /// + 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 new file mode 100644 index 000000000..07f7f5257 --- /dev/null +++ b/src/Numerics.Tests/InterpolationTests/LogLinearTest.cs @@ -0,0 +1,56 @@ +// +// 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. +// + +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 }; + + /// + /// Verifies that the 3rd derivative matches the given value at all the provided sample points. + /// + [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 new file mode 100644 index 000000000..6847a6cd1 --- /dev/null +++ b/src/Numerics.Tests/InterpolationTests/QuadraticSplineTest.cs @@ -0,0 +1,56 @@ +// +// 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. +// + +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 }; + + /// + /// Verifies that the 3rd derivative matches the given value at all the provided sample points. + /// + [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 index e677aac4f..202a03cb1 100644 --- a/src/Numerics.Tests/InterpolationTests/StepInterpolationTest.cs +++ b/src/Numerics.Tests/InterpolationTests/StepInterpolationTest.cs @@ -1,4 +1,4 @@ -// +// // 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 index 90b3418d7..2bb50763e 100644 --- a/src/Numerics/Interpolation/Barycentric.cs +++ b/src/Numerics/Interpolation/Barycentric.cs @@ -334,6 +334,13 @@ public double Interpolate(double t) /// Interpolated second derivative at point t. double IInterpolation.Differentiate2(double t) => throw new NotSupportedException(); + /// + /// Differentiate three times at point t. NOT SUPPORTED. + /// + /// Point t to interpolate at. + /// Interpolated third derivative at point t. + double IInterpolation.Differentiate3(double t) => throw new NotSupportedException(); + /// /// Indefinite integral at point t. NOT SUPPORTED. /// diff --git a/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs b/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs index 9dfeb767a..ec84aca90 100644 --- a/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs +++ b/src/Numerics/Interpolation/BulirschStoerRationalInterpolation.cs @@ -182,6 +182,13 @@ public double Interpolate(double t) /// Interpolated second derivative at point t. double IInterpolation.Differentiate2(double t) => throw new NotSupportedException(); + /// + /// Differentiate three times at point t. NOT SUPPORTED. + /// + /// Point t to interpolate at. + /// Interpolated third derivative at point t. + double IInterpolation.Differentiate3(double t) => throw new NotSupportedException(); + /// /// Indefinite integral at point t. NOT SUPPORTED. /// diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index d1f87e46c..ba45e4abe 100644 --- 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]; } + /// + /// Differentiate three times at point t. + /// + /// Point t to interpolate at. + /// Interpolated third derivative at point t. + public double Differentiate3(double t) + { + int k = LeftSegmentIndex(t); + return 6*_c3[k]; + } + /// /// Indefinite integral at point t. /// diff --git a/src/Numerics/Interpolation/IInterpolation.cs b/src/Numerics/Interpolation/IInterpolation.cs index 03c4e66c3..60659dbf6 100644 --- a/src/Numerics/Interpolation/IInterpolation.cs +++ b/src/Numerics/Interpolation/IInterpolation.cs @@ -65,6 +65,13 @@ public interface IInterpolation /// Interpolated second derivative at point t. double Differentiate2(double t); + /// + /// Differentiate three times at point t. + /// + /// Point t to interpolate at. + /// Interpolated third derivative at point t. + double Differentiate3(double t); + /// /// Indefinite integral at point t. /// diff --git a/src/Numerics/Interpolation/LinearSpline.cs b/src/Numerics/Interpolation/LinearSpline.cs index 0f192c3e3..eecddb43f 100644 --- a/src/Numerics/Interpolation/LinearSpline.cs +++ b/src/Numerics/Interpolation/LinearSpline.cs @@ -152,6 +152,13 @@ public double Differentiate(double t) /// Interpolated second derivative at point t. public double Differentiate2(double t) => 0d; + /// + /// Differentiate three times at point t. + /// + /// Point t to interpolate at. + /// Interpolated third derivative at point t. + public double Differentiate3(double t) => 0d; + /// /// Indefinite integral at point t. /// diff --git a/src/Numerics/Interpolation/LogLinear.cs b/src/Numerics/Interpolation/LogLinear.cs index e40702e08..156923b61 100644 --- a/src/Numerics/Interpolation/LogLinear.cs +++ b/src/Numerics/Interpolation/LogLinear.cs @@ -152,6 +152,24 @@ public double Differentiate2(double t) return secondDerivative; } + /// + /// Differentiate three times at point t. + /// + /// Point t to interpolate at. + /// Interpolated third derivative at point t. + 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; + } + /// /// Indefinite integral at point t. /// diff --git a/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs b/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs index 389ea8456..3be6a10a9 100644 --- a/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs +++ b/src/Numerics/Interpolation/NevillePolynomialInterpolation.cs @@ -197,6 +197,36 @@ public double Differentiate2(double t) return ddx[0]; } + /// + /// Differentiate three times at point t. + /// + /// Point t to interpolate at. + /// Interpolated third derivative at point t. + 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]; + } + /// /// Indefinite integral at point t. NOT SUPPORTED. /// diff --git a/src/Numerics/Interpolation/QuadraticSpline.cs b/src/Numerics/Interpolation/QuadraticSpline.cs index f25edf55b..491931915 100644 --- 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]; } + /// + /// Differentiate three times at point t. + /// + /// Point t to interpolate at. + /// Interpolated third derivative at point t. + public double Differentiate3(double t) => 0d; + /// /// Indefinite integral at point t. /// diff --git a/src/Numerics/Interpolation/StepInterpolation.cs b/src/Numerics/Interpolation/StepInterpolation.cs index c439ae87a..0d33bfbec 100644 --- a/src/Numerics/Interpolation/StepInterpolation.cs +++ b/src/Numerics/Interpolation/StepInterpolation.cs @@ -139,6 +139,13 @@ public double Differentiate(double t) /// Interpolated second derivative at point t. public double Differentiate2(double t) => Differentiate(t); + /// + /// Differentiate three times at point t. + /// + /// Point t to interpolate at. + /// Interpolated third derivative at point t. + public double Differentiate3(double t) => Differentiate2(t); + /// /// Indefinite integral at point t. /// diff --git a/src/Numerics/Interpolation/TransformedInterpolation.cs b/src/Numerics/Interpolation/TransformedInterpolation.cs index 36b58a6f8..d0099c889 100644 --- a/src/Numerics/Interpolation/TransformedInterpolation.cs +++ b/src/Numerics/Interpolation/TransformedInterpolation.cs @@ -146,6 +146,13 @@ public static TransformedInterpolation Interpolate( /// Interpolated second derivative at point t. double IInterpolation.Differentiate2(double t) => throw new NotSupportedException(); + /// + /// Differentiate three times at point t. NOT SUPPORTED. + /// + /// Point t to interpolate at. + /// Interpolated third derivative at point t. + double IInterpolation.Differentiate3(double t) => throw new NotSupportedException(); + /// /// Indefinite integral at point t. NOT SUPPORTED. ///