Web 技术研究所

我一直坚信着,Web 将会成为未来应用程序的主流

多项式插值的计算

  给定数列的前n个项,我们都可以找到一个满足所有给定项的n次多项式。也就是说,在任何有限项的数列中都可以找出多项式作为它的通式。要找出这个多项式也不困难,只需要使用线性代数的知识就可以找出我们需要的多项式,本文使用JavaScript来实现它。
  首先,我们可以把给定的数列当做一个n阶等差数列,得到的n阶差是0。这就意味着存在一个函数,它的n阶导数是0。那么根据微积分的基本原理,这个函数就是一个n次多项式。但是我们并不知道这个多项式每一项的系数,只知道给定的n个x与y的对应值。这就可以把它的各项系数当做未知数,再带入已知的各组x和y,构成一个多元一次方程组。解出这个多元一次方程组就可以得到多项式的各次项系数。有了这些系数,我们最终要求的函数就出来了。
  计算的关键就是解这个多元一次方程组。如果人工去解这玩意儿肯定蛋疼死,我们让程序来完成,思路依然按照人工的解法。先把这个方程组写成矩阵,然后使用高斯消除法化简。下面是代码 function findFunction(s,l){
  var i,j,k,m=[];
  //构造矩阵
  for(i in s){
    tmp=[];
    for(j=0;j<l;j++)tmp.push(Math.pow(i,j));
    m.push(tmp.concat(s[i]));
  };
  //高斯消除
  for(i=0;i<l;i++){ //项
    //下三角
    for(j=i;j<l;j++){ //行
      tmp=m[j][i];
      for(k=i;k<=l;k++){ //列
        if(tmp)m[j][k]/=tmp;
        if(j!=i)m[j][k]-=m[i][k];
      };
    };
    //上三角
    for(j=0;j<i;j++){ //行
      tmp=m[j][i];
      for(k=i;k<=l;k++){ //列
        m[j][k]-=m[i][k]*tmp;
      };
    };
  };
  //生成函数
  k=[];
  for(i=0;i<l;i++)k.push(m[i][l]);
  return function(x){
    var i,r=0;
    for(i=0;i<l;i++)r+=Math.pow(x,i)*k[i];
    return r;
  };
};
  这个函数有两个参数,s是存放已知项的JavaScript对象,l是已知项的个数。构造矩阵要做的其实就是通过多项式的原型来算出每一个未知项目在给定x时的系数。我的测试数据是这样的 s={0:100,50:30,120:60,180:70,230:130,280:150,310:180};
f=findFunction(s,7);
  s中的key是x值,value是对应x的y值。把这个s通过上面定义的findFunction函数来处理。在执行完构造矩阵的步骤时会得到这样一个矩阵。

  这就是把方程组写成矩阵的形式,注意这里的定义是“0的0次方等于1”,JavaScript的pow函数本身也是这个定义。左边的7列是多项式的每一项带入x后的东西,由于每一项的系数不同,所以每一项写在单独的一列上。最右边的一列是相应的y值。我们接下来要做的是把左边的这7行7列变换成单位矩阵,所有变换过程中产生的副作用都堆积到最右边的那一列上,这就是高斯消除法。当然也有其他方法,不过这个方法是最容易理解的了。我使用的算法也是刚刚才写出来的,没做过啥优化,也许还存在一些问题。这个算法应该自己理解了写一遍,要不然对整个概念都清晰不起来。做完这个变换后可以得到另一个矩阵。

  这个矩阵已经把多项式的所有系数都变得互不关联了,唯一对应的就是右边的值。所以我们需要的多项式的每一项系数就是右边的值。现在只要按照我们生产这个矩阵的顺序把右边的值取出来作为多项式的各项系数即可。我的findFunction函数返回的是一个函数,就是这个多项式函数了。我的例子使用了CANVAS把这个函数的图像画出来了,可以发现这条曲线是平滑的,并且经过了所有给定的点。

  不过毕竟是多项式插值,当次数高时它在数值上很不稳定。也就是说,当次数高时,我只要让一个点改变很小的值,曲线就会发生很大的变化。而且你也可以看到我们上面矩阵中的数值非常大,因为它的次数增大时数值是指数增长的。如果我们使用数值法计算,随着次数的增加,误差也会指数增长,这就是传说中的龙格现象。所以多项式插值不太实用,不过有了多项式插值的理念,理解起其他东西就很容易。这篇文章就说到这儿吧,最后是文中例子的完整代码。

  多项式插值的计算
网名:
52.91.185.*
电子邮箱:
仅用于接收通知
提交 悄悄的告诉你,Ctrl+Enter 可以提交哦
神奇海螺
[查看全部主题]
各类Web技术问题讨论区
发起新主题
本模块采用即时聊天邮件通知的模式
让钛合金F5成为历史吧!
次碳酸钴的技术博客,文章原创,转载请保留原文链接 ^_^