diff --git a/src/Managed.UnitTests/InterpolationTests/InterpolationTest.cs b/src/Managed.UnitTests/InterpolationTests/InterpolationTest.cs index 04624134..9be797d8 100644 --- a/src/Managed.UnitTests/InterpolationTests/InterpolationTest.cs +++ b/src/Managed.UnitTests/InterpolationTests/InterpolationTest.cs @@ -50,5 +50,12 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests Factory = Interpolation.CreateRationalPoleFree, Order = new[] { 1, 2, 6 } }; + + [VerifyContract] + public readonly IContract NevillePolynomialContractTests = new InterpolationContract() + { + Factory = (t, x) => new NevillePolynomialInterpolation(t, x), + Order = new[] { 1, 2, 6 } + }; } } diff --git a/src/Managed/Interpolation/Algorithms/NevillePolynomialInterpolation.cs b/src/Managed/Interpolation/Algorithms/NevillePolynomialInterpolation.cs new file mode 100644 index 00000000..58ee1542 --- /dev/null +++ b/src/Managed/Interpolation/Algorithms/NevillePolynomialInterpolation.cs @@ -0,0 +1,225 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.Interpolation.Algorithms +{ + using System; + using System.Collections.Generic; + + /// + /// Lagrange Polynomial Interpolation using Neville's Algorithm. + /// + /// + /// + /// This algorithm supports differentiation, but doesn't support integration. + /// + /// + /// When working with equidistant or Chebyshev sample points it is + /// recommended to use the barycentric algorithms specialized for + /// these cases instead of this arbitrary Neville algorithm. + /// + /// + public class NevillePolynomialInterpolation : IInterpolation + { + /// + /// Sample Points t. + /// + private IList points; + + /// + /// Spline Values x(t). + /// + private IList values; + + /// + /// Initializes a new instance of the NevillePolynomialInterpolation class. + /// + public NevillePolynomialInterpolation() + { + } + + /// + /// Initializes a new instance of the NevillePolynomialInterpolation class. + /// + /// Sample Points t + /// Sample Values x(t) + public NevillePolynomialInterpolation( + IList samplePoints, + IList sampleValues) + { + this.Initialize(samplePoints, sampleValues); + } + + /// + /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative). + /// + /// + /// + bool IInterpolation.SupportsDifferentiation + { + get { return true; } + } + + /// + /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature). + /// + /// + bool IInterpolation.SupportsIntegration + { + get { return false; } + } + + /// + /// Initialize the interpolation method with the given spline coefficients. + /// + /// Sample Points t + /// Sample Values x(t) + public void Initialize( + IList samplePoints, + IList sampleValues) + { + if (null == samplePoints) + { + throw new ArgumentNullException("samplePoints"); + } + + if (null == sampleValues) + { + throw new ArgumentNullException("sampleValues"); + } + + if (samplePoints.Count != sampleValues.Count) + { + throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLengths); + } + + this.points = samplePoints; + this.values = sampleValues; + } + + /// + /// Interpolate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t). + public double Interpolate(double t) + { + double[] x = new double[this.values.Count]; + this.values.CopyTo(x, 0); + + for (int level = 1; level < x.Length; level++) + { + for (int i = 0; i < x.Length - level; i++) + { + double hp = t - this.points[i + level]; + double ho = this.points[i] - t; + double den = this.points[i] - this.points[i + level]; + x[i] = ((hp * x[i]) + (ho * x[i + 1])) / den; + } + } + + return x[0]; + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated first derivative at point t. + /// + /// + public double Differentiate(double t) + { + double[] x = new double[this.values.Count]; + double[] dx = new double[this.values.Count]; + this.values.CopyTo(x, 0); + + for (int level = 1; level < x.Length; level++) + { + for (int i = 0; i < x.Length - level; i++) + { + double hp = t - this.points[i + level]; + double ho = this.points[i] - t; + double den = this.points[i] - this.points[i + level]; + 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 dx[0]; + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t) + /// Interpolated second derivative at point t. + /// Interpolated first derivative at point t. + /// + /// + public double Differentiate( + double t, + out double interpolatedValue, + out double secondDerivative) + { + double[] x = new double[this.values.Count]; + double[] dx = new double[this.values.Count]; + double[] d2x = new double[this.values.Count]; + this.values.CopyTo(x, 0); + + for (int level = 1; level < x.Length; level++) + { + for (int i = 0; i < x.Length - level; i++) + { + double hp = t - this.points[i + level]; + double ho = this.points[i] - t; + double den = this.points[i] - this.points[i + level]; + d2x[i] = ((hp * d2x[i]) + (ho * d2x[i + 1]) + (2 * dx[i]) - (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; + } + } + + interpolatedValue = x[0]; + secondDerivative = d2x[0]; + return dx[0]; + } + + /// + /// Integrate up to point t. + /// + /// Right bound of the integration interval [a,t]. + /// Interpolated definite integral over the interval [a,t]. + /// + double IInterpolation.Integrate(double t) + { + throw new NotSupportedException(); + } + } +} diff --git a/src/Managed/Managed.csproj b/src/Managed/Managed.csproj index d076ff04..58e0c140 100644 --- a/src/Managed/Managed.csproj +++ b/src/Managed/Managed.csproj @@ -56,6 +56,7 @@ + diff --git a/src/Native/Native.csproj b/src/Native/Native.csproj index 5b250a92..5a6e697f 100644 --- a/src/Native/Native.csproj +++ b/src/Native/Native.csproj @@ -77,6 +77,9 @@ Interpolation\Algorithms\LinearSplineInterpolation.cs + + Interpolation\Algorithms\NevillePolynomialInterpolation.cs + Interpolation\Algorithms\RationalPoleFreeInterpolation.cs