|
|
@ -4,7 +4,7 @@ |
|
|
// http://github.com/mathnet/mathnet-numerics
|
|
|
// http://github.com/mathnet/mathnet-numerics
|
|
|
// http://mathnetnumerics.codeplex.com
|
|
|
// http://mathnetnumerics.codeplex.com
|
|
|
//
|
|
|
//
|
|
|
// Copyright (c) 2009-2010 Math.NET
|
|
|
// Copyright (c) 2009-2013 Math.NET
|
|
|
//
|
|
|
//
|
|
|
// Permission is hereby granted, free of charge, to any person
|
|
|
// Permission is hereby granted, free of charge, to any person
|
|
|
// obtaining a copy of this software and associated documentation
|
|
|
// obtaining a copy of this software and associated documentation
|
|
|
@ -77,8 +77,6 @@ namespace MathNet.Numerics.Interpolation |
|
|
/// <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>
|
|
|
/// <seealso cref="Differentiate(double)"/>
|
|
|
|
|
|
/// <seealso cref="DifferentiateAll(double)"/>
|
|
|
|
|
|
bool IInterpolation.SupportsDifferentiation |
|
|
bool IInterpolation.SupportsDifferentiation |
|
|
{ |
|
|
{ |
|
|
get { return true; } |
|
|
get { return true; } |
|
|
@ -87,7 +85,6 @@ namespace MathNet.Numerics.Interpolation |
|
|
/// <summary>
|
|
|
/// <summary>
|
|
|
/// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
|
|
|
/// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
|
|
|
/// </summary>
|
|
|
/// </summary>
|
|
|
/// <seealso cref="Integrate"/>
|
|
|
|
|
|
bool IInterpolation.SupportsIntegration |
|
|
bool IInterpolation.SupportsIntegration |
|
|
{ |
|
|
{ |
|
|
get { return true; } |
|
|
get { return true; } |
|
|
@ -136,16 +133,16 @@ namespace MathNet.Numerics.Interpolation |
|
|
/// <returns>Interpolated value x(t).</returns>
|
|
|
/// <returns>Interpolated value x(t).</returns>
|
|
|
public double Interpolate(double t) |
|
|
public double Interpolate(double t) |
|
|
{ |
|
|
{ |
|
|
int closestLeftIndex = IndexOfClosestPointLeftOf(t); |
|
|
int closestLeftIndex = LeftBracketIndex(t); |
|
|
|
|
|
|
|
|
// Interpolation
|
|
|
// Interpolation
|
|
|
double offset = t - _points[closestLeftIndex]; |
|
|
double offset = t - _points[closestLeftIndex]; |
|
|
int k = closestLeftIndex << 2; |
|
|
int k = closestLeftIndex << 2; |
|
|
|
|
|
|
|
|
return _coefficients[k] |
|
|
return _coefficients[k] |
|
|
+ (offset*(_coefficients[k + 1] |
|
|
+ (offset*(_coefficients[k + 1] |
|
|
+ (offset*(_coefficients[k + 2] |
|
|
+ (offset*(_coefficients[k + 2] |
|
|
+ (offset*_coefficients[k + 3]))))); |
|
|
+ (offset*_coefficients[k + 3]))))); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
/// <summary>
|
|
|
/// <summary>
|
|
|
@ -153,49 +150,38 @@ namespace MathNet.Numerics.Interpolation |
|
|
/// </summary>
|
|
|
/// </summary>
|
|
|
/// <param name="t">Point t to interpolate at.</param>
|
|
|
/// <param name="t">Point t to interpolate at.</param>
|
|
|
/// <returns>Interpolated first derivative at point t.</returns>
|
|
|
/// <returns>Interpolated first derivative at point t.</returns>
|
|
|
/// <seealso cref="IInterpolation.SupportsDifferentiation"/>
|
|
|
|
|
|
/// <seealso cref="DifferentiateAll(double)"/>
|
|
|
|
|
|
public double Differentiate(double t) |
|
|
public double Differentiate(double t) |
|
|
{ |
|
|
{ |
|
|
int closestLeftIndex = IndexOfClosestPointLeftOf(t); |
|
|
int closestLeftIndex = LeftBracketIndex(t); |
|
|
|
|
|
|
|
|
// Differentiation
|
|
|
|
|
|
double offset = t - _points[closestLeftIndex]; |
|
|
double offset = t - _points[closestLeftIndex]; |
|
|
int k = closestLeftIndex << 2; |
|
|
int k = closestLeftIndex << 2; |
|
|
|
|
|
|
|
|
return _coefficients[k + 1] |
|
|
return _coefficients[k + 1] |
|
|
+ (2*offset*_coefficients[k + 2]) |
|
|
+ (2*offset*_coefficients[k + 2]) |
|
|
+ (3*offset*offset*_coefficients[k + 3]); |
|
|
+ (3*offset*offset*_coefficients[k + 3]); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
/// <summary>
|
|
|
/// <summary>
|
|
|
/// Interpolate, differentiate and 2nd differentiate at point t.
|
|
|
/// Differentiate twice at point t.
|
|
|
/// </summary>
|
|
|
/// </summary>
|
|
|
/// <param name="t">Point t to interpolate at.</param>
|
|
|
/// <param name="t">Point t to interpolate at.</param>
|
|
|
/// <returns>Interpolated first derivative at point t.</returns>
|
|
|
/// <returns>Interpolated second derivative at point t.</returns>
|
|
|
/// <seealso cref="IInterpolation.SupportsDifferentiation"/>
|
|
|
public double Differentiate2(double t) |
|
|
/// <seealso cref="Differentiate(double)"/>
|
|
|
|
|
|
public Tuple<double, double, double> DifferentiateAll(double t) |
|
|
|
|
|
{ |
|
|
{ |
|
|
int closestLeftIndex = IndexOfClosestPointLeftOf(t); |
|
|
int closestLeftIndex = LeftBracketIndex(t); |
|
|
double offset = t - _points[closestLeftIndex]; |
|
|
double offset = t - _points[closestLeftIndex]; |
|
|
int k = closestLeftIndex << 2; |
|
|
int k = closestLeftIndex << 2; |
|
|
|
|
|
|
|
|
return new Tuple<double, double, double>( |
|
|
return (2*_coefficients[k + 2]) + (6*offset*_coefficients[k + 3]); |
|
|
_coefficients[k] + (offset*(_coefficients[k + 1] + (offset*(_coefficients[k + 2] + (offset*_coefficients[k + 3]))))), |
|
|
|
|
|
_coefficients[k + 1] + (2*offset*_coefficients[k + 2]) + (3*offset*offset*_coefficients[k + 3]), |
|
|
|
|
|
(2*_coefficients[k + 2]) + (6*offset*_coefficients[k + 3])); |
|
|
|
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
/// <summary>
|
|
|
/// <summary>
|
|
|
/// Integrate up to point t.
|
|
|
/// Indefinite integral at point t.
|
|
|
/// </summary>
|
|
|
/// </summary>
|
|
|
/// <param name="t">Right bound of the integration interval [a,t].</param>
|
|
|
/// <param name="t">Point t to integrate at.</param>
|
|
|
/// <returns>Interpolated definite integral over the interval [a,t].</returns>
|
|
|
|
|
|
/// <seealso cref="IInterpolation.SupportsIntegration"/>
|
|
|
|
|
|
public double Integrate(double t) |
|
|
public double Integrate(double t) |
|
|
{ |
|
|
{ |
|
|
int closestLeftIndex = IndexOfClosestPointLeftOf(t); |
|
|
int closestLeftIndex = LeftBracketIndex(t); |
|
|
|
|
|
|
|
|
// Integration
|
|
|
// Integration
|
|
|
double result = 0; |
|
|
double result = 0; |
|
|
@ -203,18 +189,28 @@ namespace MathNet.Numerics.Interpolation |
|
|
{ |
|
|
{ |
|
|
double w = _points[i + 1] - _points[i]; |
|
|
double w = _points[i + 1] - _points[i]; |
|
|
result += w*(_coefficients[j] |
|
|
result += w*(_coefficients[j] |
|
|
+ ((w*_coefficients[j + 1]*0.5) |
|
|
+ ((w*_coefficients[j + 1]*0.5) |
|
|
+ (w*((_coefficients[j + 2]/3) |
|
|
+ (w*((_coefficients[j + 2]/3) |
|
|
+ (w*_coefficients[j + 3]*0.25))))); |
|
|
+ (w*_coefficients[j + 3]*0.25))))); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
double offset = t - _points[closestLeftIndex]; |
|
|
double offset = t - _points[closestLeftIndex]; |
|
|
int k = closestLeftIndex << 2; |
|
|
int k = closestLeftIndex << 2; |
|
|
|
|
|
|
|
|
return result + (offset*(_coefficients[k] |
|
|
return result + (offset*(_coefficients[k] |
|
|
+ (offset*_coefficients[k + 1]*0.5) |
|
|
+ (offset*_coefficients[k + 1]*0.5) |
|
|
+ (offset*_coefficients[k + 2]/3) |
|
|
+ (offset*_coefficients[k + 2]/3) |
|
|
+ (offset*_coefficients[k + 3]*0.25))); |
|
|
+ (offset*_coefficients[k + 3]*0.25))); |
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
|
|
/// Definite integral between points a and b.
|
|
|
|
|
|
/// </summary>
|
|
|
|
|
|
/// <param name="a">Left bound of the integration interval [a,b].</param>
|
|
|
|
|
|
/// <param name="b">Right bound of the integration interval [a,b].</param>
|
|
|
|
|
|
public double Integrate(double a, double b) |
|
|
|
|
|
{ |
|
|
|
|
|
return Integrate(b) - Integrate(a); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
/// <summary>
|
|
|
/// <summary>
|
|
|
@ -222,7 +218,7 @@ namespace MathNet.Numerics.Interpolation |
|
|
/// </summary>
|
|
|
/// </summary>
|
|
|
/// <param name="t">The value to look for.</param>
|
|
|
/// <param name="t">The value to look for.</param>
|
|
|
/// <returns>The sample point index.</returns>
|
|
|
/// <returns>The sample point index.</returns>
|
|
|
int IndexOfClosestPointLeftOf(double t) |
|
|
int LeftBracketIndex(double t) |
|
|
{ |
|
|
{ |
|
|
// Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included)
|
|
|
// Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included)
|
|
|
int low = 0; |
|
|
int low = 0; |
|
|
|