我知道这可能是一个有点“请帮我做作业”的问题,但我正在尝试按照这里所解释的De Boor算法来实现:http://en.wikipedia.org/wiki/De_Boor's_algorithm
我已经有了一个通用类,可以接受任何类型的数组,并使用两个函数为每个点生成x和y。
当我向GetPoint(x)提供非常接近_d中最低x值的x时,我得到的插值非常不准确。
有人能看出我哪里做错了吗?非常感谢您的帮助。
我正在使用的数据是:
代码
当我向GetPoint(x)提供非常接近_d中最低x值的x时,我得到的插值非常不准确。
有人能看出我哪里做错了吗?非常感谢您的帮助。
我正在使用的数据是:
x y
0.002739726 0.999988548
0.019178082 0.999917365
0.038356164 0.999835128
0.057534247 0.999753381
0.082191781 0.999648832
0.167123288 0.999286638
0.249315068 0.998939303
0.334246575 0.998583164
0.419178082 0.998222171
0.495890411 0.997899634
0.747945205 0.99680615
1 0.99559971
1.249315068 0.994187653
1.495890411 0.9926708
1.747945205 0.99066351
2 0.988439344
2.249315068 0.985377424
2.498630137 0.981972695
2.750684932 0.978188863
3.002739726 0.974069647
4.002739726 0.952415995
5.002739726 0.926226504
6.002739726 0.897271422
7.005479452 0.866241091
8.005479452 0.834348616
9.005479452 0.802369725
代码
public class Spline<T>
{
/// <summary>
/// Coordinate
/// </summary>
private struct TwoDPoint
{
public double x;
public double y;
};
/// <summary>
/// Unused apart from debugging.
/// </summary>
private T [] _originalDataPoints;
/// <summary>
/// Stores the points we will perform the algorithm on.
/// </summary>
private TwoDPoint [] _d;
private int _n; // Value for the "order" used in the De Boor algorithm
public Spline(T [] dataPoints, Func<T, double> xFunc, Func<T, double> yFunc, int n)
{
SortAndAssignPoints(dataPoints, xFunc, yFunc);
_n = n;
}
/// <summary>
/// Populates points in d sorted ascending by their x value
/// </summary>
private void SortAndAssignPoints(T [] dataPoints, Func<T, double> xFunc, Func<T, double> yFunc)
{
_originalDataPoints = dataPoints;
_d = dataPoints
.Select(point => new TwoDPoint()
{
x = xFunc(point),
y = yFunc(point)
})
.OrderBy(point => point.x).ToArray();
}
/// <summary>
/// Function which gets an interpolated value
/// </summary>
public double GetPoint(double x)
{
double sum = 0;
for (int i = 0; i < _d.Length; i++)
{
sum += _d[i].y * N_General(_n, i, x);
}
return sum;
}
/// <summary>
/// General N function as given in the algorithm
/// </summary>
private double N_General(int n, int i, double x)
{
if (n == 0)
{
return N_Base(i, x);
}
double firstCoefficient;
if ((i + n) > (_d.Length - 1))
{
firstCoefficient = 0D;
}
else
{
firstCoefficient =
(x - _d[i].x) / (_d[i + n].x - _d[i].x);
}
double secondCoefficient;
if ((i + n + 1) > (_d.Length - 1))
{
secondCoefficient = 0D;
}
else
{
secondCoefficient =
(_d[i + n + 1].x - x) / (_d[i + n + 1].x - _d[i + 1].x);
}
return firstCoefficient * N_General(n - 1, i, x) +
secondCoefficient * N_General(n - 1, i + 1, x);
}
/// <summary>
/// N_i^0 as given in the algorithm
/// </summary>
private double N_Base(int i, double x)
{
//Bound check
if (
(i >= (_d.Length - 1)) ||
(i < 0)
)
{
return 0;
}
if ((x >= _d[i].x) &&
(x < _d[i+1].x))
{
return 1;
}
else
{
return 0;
}
}
}