14

我需要使用拉格朗日插值多项式计算多项式的系数,作为我的作业,我决定在 Javascript 中执行此操作。

这是拉格朗日多项式 (L(x)) 的定义

在此处输入图像描述

拉格朗日基多项式定义如下

在此处输入图像描述

计算特定 X 的 y 值(W(x) 函数)很简单,但我需要计算多项式的系数([a0,a1,...,an] 的数组)我需要对 n<=10 执行此操作,但它有任意n会很好,然后我可以将该函数放入horner函数并绘制该多项式。

在此处输入图像描述

我有在第一个方程中计算分母的函数

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;
}

和使用 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;
   };
}

任何人都知道这样做的算法,或者知道如何计算这些系数

4

3 回答 3

10

好吧,你可以用天真的方式来做。通过其系数的数组表示多项式,数组

[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_kfor all k != i。这样就给出了 L i的系数,然后你只需将它们与 y i相乘(如果你将 y 值作为参数传递,你可以通过初始化coefficients来做到这一点y_i/denominator(i,points)),最后将所有系数相加。

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,因此需要从高一开始。

这还不能产生正确的插值,但至少比以前更明智。

此外,在您的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;
}
于 2012-03-25T15:55:26.537 回答
1

如果您使用“仿射单纯形的初学者指南 ”,“拉格朗日插值”部分中介绍的拉格朗日插值矩阵形式,您可以相对容易地找到拉格朗日插值多项式的系数。恐怕,我不知道 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)

此代码计算拉格朗日插值多项式的系数,打印它们,并测试给定的 x 映射到预期的 y。您可以使用Google colab测试此代码,因此您无需安装任何东西。大概,你可以把它翻译成JS。

于 2019-05-31T17:12:25.610 回答
1

此代码将确定从常数项开始的系数。

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];
}

于 2020-04-17T06:39:08.927 回答