Я использовал один из предложенных алгоритмов, но результаты очень плохие.
Я реализовал алгоритм вики
на Java (код ниже). x(0)
- это points.get(0)
, y(0)
- это values[points.get(0)]
, α
- это alfa
, а μ
- это mi
. В остальном то же, что и в псевдокоде вики.
public void createSpline(double[] values, ArrayList<Integer> points){
a = new double[points.size()+1];
for (int i=0; i <points.size();i++)
{
a[i] = values[points.get(i)];
}
b = new double[points.size()];
d = new double[points.size()];
h = new double[points.size()];
for (int i=0; i<points.size()-1; i++){
h[i] = points.get(i+1) - points.get(i);
}
alfa = new double[points.size()];
for (int i=1; i <points.size()-1; i++){
alfa[i] = (double)3 / h[i] * (a[i+1] - a[i]) - (double)3 / h[i-1] *(a[i+1] - a[i]);
}
c = new double[points.size()+1];
l = new double[points.size()+1];
mi = new double[points.size()+1];
z = new double[points.size()+1];
l[0] = 1; mi[0] = z[0] = 0;
for (int i =1; i<points.size()-1;i++){
l[i] = 2 * (points.get(i+1) - points.get(i-1)) - h[i-1]*mi[i-1];
mi[i] = h[i]/l[i];
z[i] = (alfa[i] - h[i-1]*z[i-1])/l[i];
}
l[points.size()] = 1;
z[points.size()] = c[points.size()] = 0;
for (int j=points.size()-1; j >0; j--)
{
c[j] = z[j] - mi[j]*c[j-1];
b[j] = (a[j+1]-a[j]) - (h[j] * (c[j+1] + 2*c[j])/(double)3) ;
d[j] = (c[j+1]-c[j])/((double)3*h[j]);
}
for (int i=0; i<points.size()-1;i++){
for (int j = points.get(i); j<points.get(i+1);j++){
// fk[j] = values[points.get(i)];
functionResult[j] = a[i] + b[i] * (j - points.get(i))
+ c[i] * Math.pow((j - points.get(i)),2)
+ d[i] * Math.pow((j - points.get(i)),3);
}
}
}
Результат, который я получаю, следующий:
но это должно быть похоже на это:
Я также пытаюсь реализовать алгоритм по-другому в соответствии с: http://www.geos.ed.ac.uk/~yliu23/docs/lect_spline.pdf
Сначала они показывают, как делать линейные сплайны, и это довольно просто. Я создаю функции, которые вычисляют коэффициенты A
и B
. Затем они расширяют линейный сплайн, добавляя вторую производную. Коэффициенты C
и D
тоже легко вычислить.
Но проблемы начинаются, когда я пытаюсь вычислить вторую производную. Я не понимаю, как они их рассчитывают.
Поэтому я реализовал только линейную интерполяцию. Результат:
Кто-нибудь знает, как исправить первый алгоритм или объяснить мне, как вычислить вторую производную во втором алгоритме?