Ну, вы можете сделать это наивным способом. Представляем полином массивом его коэффициентов, массивом
[a_0,a_1,...,a_n]
соответствует a_0 + a_1*X + ... + a_n*X^n. Я плохо разбираюсь в JavaScript, поэтому псевдокод должен будет сделать:
interpolation_polynomial(i,points)
coefficients = [1/denominator(i,points)]
for k = 0 to points.length-1
if k == i
next k
new_coefficients = [0,0,...,0] // length k+2 if k < i, k+1 if k > i
if k < i
m = k
else
m = k-1
for j = m downto 0
new_coefficients[j+1] += coefficients[j]
new_coefficients[j] -= points[k]*coefficients[j]
coefficients = new_coefficients
return coefficients
Начните с постоянного многочлена 1/((x_1-x_0)* ... *(x_i-x_{i-1})*(x_i-x_{i+1})*...*(x_i-x_n)) и умножьте на X - x_k для всех k != i. Таким образом, это дает коэффициенты для L i , тогда вы просто умножаете их на y i (вы можете сделать это, инициализируя coefficients в y_i/denominator(i,points), если вы передадите значения y в качестве параметров) и добавьте все коэффициенты вместе, наконец.
polynomial = [0,0,...,0] // points.length entries
for i = 0 to points.length-1
coefficients = interpolation_polynomial(i,points)
for k = 0 to points.length-1
polynomial[k] += y[i]*coefficients[k]
Расчет каждого L i равен O (n²), поэтому общий расчет равен O (n³).
Обновление: В вашем jsFiddle у вас была ошибка в цикле умножения полиномов, в дополнение к (теперь исправленной) ошибке с начальным индексом, который я сделал, должно быть
for (var j= (k < i) ? (k+1) : k; j--;) {
new_coefficients[j+1] += coefficients[j];
new_coefficients[j] -= points[k].x*coefficients[j];
}
Так как при тестировании вы уменьшаете j, нужно начинать на 1 выше.
Это пока не дает правильной интерполяции, но, по крайней мере, более разумно, чем раньше.
Кроме того, в вашей функции horner,
function horner(array, x_scale, y_scale) {
function recur(x, i, array) {
if (i == 0) {
return x*array[0];
} else {
return array[i] + x*recur(x, --i, array);
}
}
return function(x) {
return recur(x*x_scale, array.length-1, array)*y_scale;
};
}
Вы умножаете самый высокий коэффициент в два раза на x, это должно быть
if (i == 0) {
return array[0];
}
вместо этого. Тем не менее, хорошего результата нет.
Обновление2: Окончательные исправления опечаток, следующие работы:
function horner(array, x_scale, y_scale) {
function recur(x, i, array) {
if (i == 0) {
return array[0];
} else {
return array[i] + x*recur(x, --i, array);
}
}
return function(x) {
return recur(x*x_scale, array.length-1, array)*y_scale;
};
}
// initialize array
function zeros(n) {
var array = new Array(n);
for (var i=n; i--;) {
array[i] = 0;
}
return array;
}
function denominator(i, points) {
var result = 1;
var x_i = points[i].x;
for (var j=points.length; j--;) {
if (i != j) {
result *= x_i - points[j].x;
}
}
console.log(result);
return result;
}
// calculate coefficients for Li polynomial
function interpolation_polynomial(i, points) {
var coefficients = zeros(points.length);
// alert("Denominator " + i + ": " + denominator(i,points));
coefficients[0] = 1/denominator(i,points);
console.log(coefficients[0]);
//new Array(points.length);
/*for (var s=points.length; s--;) {
coefficients[s] = 1/denominator(i,points);
}*/
var new_coefficients;
for (var k = 0; k<points.length; k++) {
if (k == i) {
continue;
}
new_coefficients = zeros(points.length);
for (var j= (k < i) ? k+1 : k; j--;) {
new_coefficients[j+1] += coefficients[j];
new_coefficients[j] -= points[k].x*coefficients[j];
}
coefficients = new_coefficients;
}
console.log(coefficients);
return coefficients;
}
// calculate coefficients of polynomial
function Lagrange(points) {
var polynomial = zeros(points.length);
var coefficients;
for (var i=0; i<points.length; ++i) {
coefficients = interpolation_polynomial(i, points);
//console.log(coefficients);
for (var k=0; k<points.length; ++k) {
// console.log(points[k].y*coefficients[k]);
polynomial[k] += points[i].y*coefficients[k];
}
}
return polynomial;
}