javascript 如何使用拉格朗日插值计算多项式的系数
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/9860937/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
How to calculate coefficients of polynomial using Lagrange interpolation
提问by jcubic
I need to calculate coefficients of polynomial using Lagrange interpolation polynomial, as my homework, I decide to do this in Javascript.
我需要使用拉格朗日插值多项式计算多项式的系数，作为我的作业，我决定在 Javascript 中执行此操作。
here is definition of Lagrange polynomial (L(x))
这是拉格朗日多项式 (L(x)) 的定义
Lagrange basis polynomials are defined as follows
拉格朗日基多项式定义如下
Calculate y value for specific X (W(x) function) is simple but I need to calculate coefficients of polynomial (array of [a0, a1, ..., an]) I need to do this to n<=10 but it will be nice to have arbitrary n, then I can put that function into horner function and draw that polynomial.
计算特定 X 的 y 值（W(x) 函数）很简单，但我需要计算多项式的系数（[a0, a1, ..., an] 数组）我需要这样做 n<=10 但它拥有任意 n 会很好，然后我可以将该函数放入 horner 函数并绘制该多项式。
I have function that calculate denominator in first equation
我有在第一个方程中计算分母的函数
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;
}
}
return result;
}
and function that return y using horner method (I also have drawing function using canvas)
和使用 horner 方法返回 y 的函数（我也有使用画布的绘图功能）
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;
};
}
anybody know algorithm to do this, or idea how to calculate those coefficients
任何人都知道这样做的算法，或者知道如何计算这些系数
采纳答案by Daniel Fischer
Well, you can do it the naive way. Represent a polynomial by the array of its coefficients, the array
好吧，你可以以天真的方式做到这一点。用系数数组表示多项式，数组
[a_0,a_1,...,a_n]
corresponding to a_0 + a_1*X + ... + a_n*X^n. I'm no good with JavaScript, so pseudocode will have to do:
对应于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
Start with the constant polynomial 1/((x_1-x_0)* ... *(x_i-x_{i-1})*(x_i-x_{i+1})*...*(x_i-x_n))and multiply with X - x_kfor all k != i. So that gives the coefficients for Li, then you just multiply them with yi(you could do that by initialising coefficientsto y_i/denominator(i,points)if you pass the y-values as parameters) and add all the coefficients together finally.
从常数多项式开始1/((x_1-x_0)* ... *(x_i-x_{i-1})*(x_i-x_{i+1})*...*(x_i-x_n))并乘以X - x_kfor all k != i。所以，让L用系数我，那么你只要将它们相乘，其中y我（你可以做到这一点通过初始化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]
Calculating each Liis O(n2), so the total calculation is O(n3).
计算每个L i是O(n2)，所以总计算是O(n3)。
Update:In your jsFiddle, you had an error in the polynomial multiplication loop in addition to (the now corrected) mistake with the start index I made, it should be
更新：在你的 jsFiddle 中，除了我做的开始索引的错误（现在更正）之外，你在多项式乘法循环中还有一个错误，它应该是
for (var j= (k < i) ? (k+1) : k; j--;) {
new_coefficients[j+1] += coefficients[j];
new_coefficients[j] -= points[k].x*coefficients[j];
}
Since you decrement jwhen testing, it needs to start one higher.
由于您j在测试时递减，因此需要从更高的开始。
That doesn't produce a correct interpolation yet, but it's at least more sensible than before.
这还没有产生正确的插值，但至少比以前更明智。
Also, in your hornerfunction,
此外，在您的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;
};
}
you multiply the highest coefficient twice with x, it should be
你乘以最高系数两次x，它应该是
if (i == 0) {
return array[0];
}
instead. Still no good result, though.
反而。尽管如此，仍然没有好结果。
Update2:Final typo fixes, the following works:
更新 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;
}
回答by guest
You can find coefficients of Lagrange interpolation polynomial relatively easy if you use a matrix form of Lagrange interpolation presented in "Beginner's guide to mapping simplexes affinely ", section "Lagrange interpolation". I'm afraid, I don't know JavaScript to provide you with appropriate code, but I worked with Python a little bit, maybe the following can help (sorry for bad codestyle -- I'm mathematician, not programmer)
如果您使用“仿射单纯形映射初学者指南 ”，“拉格朗日插值”部分中介绍的拉格朗日插值矩阵形式，您可以相对容易地找到拉格朗日插值多项式的系数。恐怕，我不知道 JavaScript 可以为您提供适当的代码，但是我使用过 Python 一点点，也许以下内容可以提供帮助（抱歉代码风格不好——我是数学家，不是程序员）
import numpy as np
# input
x = [0, 2, 4, 5] # <- x's
y = [2, 5, 7, 7] # <- y's
# calculating coefficients
M = [[_x**i*(-1)**(i*len(x)) for _x in x] for i in range(len(x))]
C = [np.linalg.det((M+[y]+M)[d:d+len(x)]) for d in range(len(x)+1)]
C = (C / C[0] * (-1)**(len(x)+1) )[1:]
# polynomial lambda-function
poly = lambda _x: sum([C[i] * _x**i for i in range(len(x))])
# output and tests
print("Coefficients:\n", C)
print("TESTING:")
for _x, _y in zip(x, y):
result = "[OK]" if np.allclose(_y, poly(_x)) else "[ERROR]"
print(_x, " mapped to: ", poly(_x), " ; expected: ", _y, result)
This code calculates coefficients of the Lagrange interpolation polynomial, prints them, and tests that given x's are mapped into expected y's. You can test this code with Google colab, so you don't have to install anything. Probably, you can translate it to JS.
此代码计算拉格朗日插值多项式的系数，打印它们，并测试给定的 x 是否映射到预期的 y。您可以使用Google colab测试此代码，因此您无需安装任何东西。也许，您可以将其翻译为 JS。
回答by Kevin Bertman
This code will determine the coefficients starting with the constant term.
此代码将确定从常数项开始的系数。
var xPoints=[2,4,3,6,7,10]; //example coordinates
var yPoints=[2,5,-2,0,2,8];
var coefficients=[];
for (var m=0; m<xPoints.length; m++) coefficients[m]=0;
for (var m=0; m<xPoints.length; m++) {
var newCoefficients=[];
for (var nc=0; nc<xPoints.length; nc++) newCoefficients[nc]=0;
if (m>0) {
newCoefficients[0]=-xPoints[0]/(xPoints[m]-xPoints[0]);
newCoefficients[1]=1/(xPoints[m]-xPoints[0]);
} else {
newCoefficients[0]=-xPoints[1]/(xPoints[m]-xPoints[1]);
newCoefficients[1]=1/(xPoints[m]-xPoints[1]);
}
var startIndex=1;
if (m==0) startIndex=2;
for (var n=startIndex; n<xPoints.length; n++) {
if (m==n) continue;
for (var nc=xPoints.length-1; nc>=1; nc--) {
newCoefficients[nc]=newCoefficients[nc]*(-xPoints[n]/(xPoints[m]-xPoints[n]))+newCoefficients[nc-1]/(xPoints[m]-xPoints[n]);
}
newCoefficients[0]=newCoefficients[0]*(-xPoints[n]/(xPoints[m]-xPoints[n]));
}
for (var nc=0; nc<xPoints.length; nc++) coefficients[nc]+=yPoints[m]*newCoefficients[nc];
}
