ufkc7e3d
/
mathnet-numerics
forked from tsai/mathnet-numerics
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Interpolation: migrate cubic spline to CubicSpline class

optimization-3
Christoph Ruegg 13 years ago
parent
fb35625def
commit
94b5c448e6
2 changed files with 173 additions and 34 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 155
      src/Numerics/Interpolation/CubicSpline.cs
  2. 52
      src/UnitTests/InterpolationTests/CubicSplineTest.cs

155
src/Numerics/Interpolation/CubicSpline.cs
View File

@ -86,6 +86,11 @@ namespace MathNet.Numerics.Interpolation
throw new ArgumentException(Resources.ArgumentVectorsSameLength); throw new ArgumentException(Resources.ArgumentVectorsSameLength);
} }
if (xx.Length < 2)
{
throw new ArgumentOutOfRangeException("x");
}
Sorting.Sort(xx, yy, dd); Sorting.Sort(xx, yy, dd);
var c0 = new double[xx.Length - 1]; var c0 = new double[xx.Length - 1];
@ -122,6 +127,11 @@ namespace MathNet.Numerics.Interpolation
throw new ArgumentException(Resources.ArgumentVectorsSameLength); throw new ArgumentException(Resources.ArgumentVectorsSameLength);
} }
if (xx.Length < 5)
{
throw new ArgumentOutOfRangeException("x");
}
Sorting.Sort(xx, yy); Sorting.Sort(xx, yy);
/* Prepare divided differences (diff) and weights (w) */ /* Prepare divided differences (diff) and weights (w) */
@ -160,6 +170,124 @@ namespace MathNet.Numerics.Interpolation
return InterpolateHermite(xx, yy, dd); return InterpolateHermite(xx, yy, dd);
} }
public static CubicSpline Interpolate(IEnumerable<double> x, IEnumerable<double> y)
{
return InterpolateBoundaries(x, y, SplineBoundaryCondition.SecondDerivative, 0.0, SplineBoundaryCondition.SecondDerivative, 0.0);
}
public static CubicSpline InterpolateBoundaries(IEnumerable<double> x, IEnumerable<double> y,
SplineBoundaryCondition leftBoundaryCondition, double leftBoundary,
SplineBoundaryCondition rightBoundaryCondition, double rightBoundary)
{
var xx = (x as double[]) ?? x.ToArray();
var yy = (y as double[]) ?? y.ToArray();
if (xx.Length != yy.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength);
}
if (xx.Length < 2)
{
throw new ArgumentOutOfRangeException("x");
}
Sorting.Sort(xx, yy);
int n = xx.Length;
// normalize special cases
if ((n == 2)
&& (leftBoundaryCondition == SplineBoundaryCondition.ParabolicallyTerminated)
&& (rightBoundaryCondition == SplineBoundaryCondition.ParabolicallyTerminated))
{
leftBoundaryCondition = SplineBoundaryCondition.SecondDerivative;
leftBoundary = 0d;
rightBoundaryCondition = SplineBoundaryCondition.SecondDerivative;
rightBoundary = 0d;
}
if (leftBoundaryCondition == SplineBoundaryCondition.Natural)
{
leftBoundaryCondition = SplineBoundaryCondition.SecondDerivative;
leftBoundary = 0d;
}
if (rightBoundaryCondition == SplineBoundaryCondition.Natural)
{
rightBoundaryCondition = SplineBoundaryCondition.SecondDerivative;
rightBoundary = 0d;
}
var a1 = new double[n];
var a2 = new double[n];
var a3 = new double[n];
var b = new double[n];
// Left Boundary
switch (leftBoundaryCondition)
{
case SplineBoundaryCondition.ParabolicallyTerminated:
a1[0] = 0;
a2[0] = 1;
a3[0] = 1;
b[0] = 2*(yy[1] - yy[0])/(xx[1] - xx[0]);
break;
case SplineBoundaryCondition.FirstDerivative:
a1[0] = 0;
a2[0] = 1;
a3[0] = 0;
b[0] = leftBoundary;
break;
case SplineBoundaryCondition.SecondDerivative:
a1[0] = 0;
a2[0] = 2;
a3[0] = 1;
b[0] = (3*((yy[1] - yy[0])/(xx[1] - xx[0]))) - (0.5*leftBoundary*(xx[1] - xx[0]));
break;
default:
throw new NotSupportedException(Resources.InvalidLeftBoundaryCondition);
}
// Central Conditions
for (int i = 1; i < xx.Length - 1; i++)
{
a1[i] = xx[i + 1] - xx[i];
a2[i] = 2*(xx[i + 1] - xx[i - 1]);
a3[i] = xx[i] - xx[i - 1];
b[i] = (3*(yy[i] - yy[i - 1])/(xx[i] - xx[i - 1])*(xx[i + 1] - xx[i])) + (3*(yy[i + 1] - yy[i])/(xx[i + 1] - xx[i])*(xx[i] - xx[i - 1]));
}
// Right Boundary
switch (rightBoundaryCondition)
{
case SplineBoundaryCondition.ParabolicallyTerminated:
a1[n - 1] = 1;
a2[n - 1] = 1;
a3[n - 1] = 0;
b[n - 1] = 2*(yy[n - 1] - yy[n - 2])/(xx[n - 1] - xx[n - 2]);
break;
case SplineBoundaryCondition.FirstDerivative:
a1[n - 1] = 0;
a2[n - 1] = 1;
a3[n - 1] = 0;
b[n - 1] = rightBoundary;
break;
case SplineBoundaryCondition.SecondDerivative:
a1[n - 1] = 1;
a2[n - 1] = 2;
a3[n - 1] = 0;
b[n - 1] = (3*(yy[n - 1] - yy[n - 2])/(xx[n - 1] - xx[n - 2])) + (0.5*rightBoundary*(xx[n - 1] - xx[n - 2]));
break;
default:
throw new NotSupportedException(Resources.InvalidRightBoundaryCondition);
}
// Build Spline
double[] dd = SolveTridiagonal(a1, a2, a3, b);
return InterpolateHermite(xx, yy, dd);
}
/// <summary> /// <summary>
/// Three-Point Differentiation Helper. /// Three-Point Differentiation Helper.
/// </summary> /// </summary>
@ -185,6 +313,33 @@ namespace MathNet.Numerics.Interpolation
return (2*a*t) + b; return (2*a*t) + b;
} }
/// <summary>
/// Tridiagonal Solve Helper.
/// </summary>
/// <param name="a">The a-vector[n].</param>
/// <param name="b">The b-vector[n], will be modified by this function.</param>
/// <param name="c">The c-vector[n].</param>
/// <param name="d">The d-vector[n], will be modified by this function.</param>
/// <returns>The x-vector[n]</returns>
static double[] SolveTridiagonal(double[] a, double[] b, double[] c, double[] d)
{
for (int k = 1; k < a.Length; k++)
{
double t = a[k]/b[k - 1];
b[k] = b[k] - (t*c[k - 1]);
d[k] = d[k] - (t*d[k - 1]);
}
var x = new double[a.Length];
x[x.Length - 1] = d[d.Length - 1]/b[b.Length - 1];
for (int k = x.Length - 2; k >= 0; k--)
{
x[k] = (d[k] - (c[k]*x[k + 1]))/b[k];
}
return x;
}
/// <summary> /// <summary>
/// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative). /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
/// </summary> /// </summary>

52
src/UnitTests/InterpolationTests/CubicSplineTest.cs
View File

@ -33,21 +33,11 @@ using NUnit.Framework;
namespace MathNet.Numerics.UnitTests.InterpolationTests namespace MathNet.Numerics.UnitTests.InterpolationTests
{ {
/// <summary>
/// CubicSpline Test case.
/// </summary>
[TestFixture, Category("Interpolation")] [TestFixture, Category("Interpolation")]
public class CubicSplineTest public class CubicSplineTest
{ {
/// <summary>
/// Sample points.
/// </summary>
readonly double[] _t = { -2.0, -1.0, 0.0, 1.0, 2.0 }; 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 };
/// <summary>
/// Sample values.
/// </summary>
readonly double[] _x = { 1.0, 2.0, -1.0, 0.0, 1.0 };
/// <summary> /// <summary>
/// Verifies that the interpolation matches the given value at all the provided sample points. /// Verifies that the interpolation matches the given value at all the provided sample points.
@ -55,11 +45,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
[Test] [Test]
public void NaturalFitsAtSamplePoints() public void NaturalFitsAtSamplePoints()
{ {
IInterpolation interpolation = new CubicSplineInterpolation(_t, _x); IInterpolation it = CubicSpline.Interpolate(_t, _y);
for (int i = 0; i < _y.Length; i++)
for (int i = 0; i < _x.Length; i++)
{ {
Assert.AreEqual(_x[i], interpolation.Interpolate(_t[i]), "A Exact Point " + i); Assert.AreEqual(_y[i], it.Interpolate(_t[i]), "A Exact Point " + i);
} }
} }
@ -85,9 +74,8 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
[TestCase(-10.0, 677, 1e-12)] [TestCase(-10.0, 677, 1e-12)]
public void NaturalFitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError) public void NaturalFitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError)
{ {
IInterpolation interpolation = new CubicSplineInterpolation(_t, _x); IInterpolation it = CubicSpline.Interpolate(_t, _y);
Assert.AreEqual(x, it.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t);
Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t);
} }
/// <summary> /// <summary>
@ -96,11 +84,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
[Test] [Test]
public void FixedFirstDerivativeFitsAtSamplePoints() public void FixedFirstDerivativeFitsAtSamplePoints()
{ {
IInterpolation interpolation = new CubicSplineInterpolation(_t, _x, SplineBoundaryCondition.FirstDerivative, 1.0, SplineBoundaryCondition.FirstDerivative, -1.0); IInterpolation it = CubicSpline.InterpolateBoundaries(_t, _y, SplineBoundaryCondition.FirstDerivative, 1.0, SplineBoundaryCondition.FirstDerivative, -1.0);
for (int i = 0; i < _y.Length; i++)
for (int i = 0; i < _x.Length; i++)
{ {
Assert.AreEqual(_x[i], interpolation.Interpolate(_t[i]), "A Exact Point " + i); Assert.AreEqual(_y[i], it.Interpolate(_t[i]), "A Exact Point " + i);
} }
} }
@ -126,9 +113,8 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
[TestCase(-10.0, 1330.1428571428571429, 1e-12)] [TestCase(-10.0, 1330.1428571428571429, 1e-12)]
public void FixedFirstDerivativeFitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError) public void FixedFirstDerivativeFitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError)
{ {
IInterpolation interpolation = new CubicSplineInterpolation(_t, _x, SplineBoundaryCondition.FirstDerivative, 1.0, SplineBoundaryCondition.FirstDerivative, -1.0); IInterpolation it = CubicSpline.InterpolateBoundaries(_t, _y, SplineBoundaryCondition.FirstDerivative, 1.0, SplineBoundaryCondition.FirstDerivative, -1.0);
Assert.AreEqual(x, it.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t);
Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t);
} }
/// <summary> /// <summary>
@ -137,11 +123,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
[Test] [Test]
public void FixedSecondDerivativeFitsAtSamplePoints() public void FixedSecondDerivativeFitsAtSamplePoints()
{ {
IInterpolation interpolation = new CubicSplineInterpolation(_t, _x, SplineBoundaryCondition.SecondDerivative, -5.0, SplineBoundaryCondition.SecondDerivative, -1.0); IInterpolation it = CubicSpline.InterpolateBoundaries(_t, _y, SplineBoundaryCondition.SecondDerivative, -5.0, SplineBoundaryCondition.SecondDerivative, -1.0);
for (int i = 0; i < _y.Length; i++)
for (int i = 0; i < _x.Length; i++)
{ {
Assert.AreEqual(_x[i], interpolation.Interpolate(_t[i]), "A Exact Point " + i); Assert.AreEqual(_y[i], it.Interpolate(_t[i]), "A Exact Point " + i);
} }
} }
@ -167,9 +152,8 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
[TestCase(-10.0, -37, 1e-12)] [TestCase(-10.0, -37, 1e-12)]
public void FixedSecondDerivativeFitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError) public void FixedSecondDerivativeFitsAtArbitraryPointsWithMaple(double t, double x, double maxAbsoluteError)
{ {
IInterpolation interpolation = new CubicSplineInterpolation(_t, _x, SplineBoundaryCondition.SecondDerivative, -5.0, SplineBoundaryCondition.SecondDerivative, -1.0); IInterpolation it = CubicSpline.InterpolateBoundaries(_t, _y, SplineBoundaryCondition.SecondDerivative, -5.0, SplineBoundaryCondition.SecondDerivative, -1.0);
Assert.AreEqual(x, it.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t);
Assert.AreEqual(x, interpolation.Interpolate(t), maxAbsoluteError, "Interpolation at {0}", t);
} }
/// <summary> /// <summary>
@ -183,10 +167,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
{ {
double[] x, y, xtest, ytest; double[] x, y, xtest, ytest;
LinearInterpolationCase.Build(out x, out y, out xtest, out ytest, samples); LinearInterpolationCase.Build(out x, out y, out xtest, out ytest, samples);
IInterpolation interpolation = new CubicSplineInterpolation(x, y); IInterpolation it = CubicSpline.Interpolate(x, y);
for (int i = 0; i < xtest.Length; i++) for (int i = 0; i < xtest.Length; i++)
{ {
Assert.AreEqual(ytest[i], interpolation.Interpolate(xtest[i]), 1e-15, "Linear with {0} samples, sample {1}", samples, i); Assert.AreEqual(ytest[i], it.Interpolate(xtest[i]), 1e-15, "Linear with {0} samples, sample {1}", samples, i);
} }
} }
} }

Write Preview
Loading…
Cancel
Save