Web 技术研究所

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

三次hermite插值

  三次hermite插值所做的事就是通过平面上的两个点和它们的切线斜率来确定一个三次函数。具体的实现方式有很多种,书本上介绍的基本都是借助lagrange插值的方法,但是lagrange插值不容易直接得到多项式函数,所以这里我依然使用了矩阵的方法。
  与之前的矩阵法做多项式插值类似,同样是把多项式的各次项系数拿出来作为需要求解的目标,利用给定的条件构造矩阵来解决问题。但与普通的多项式插值不同的是,hermite插值的条件是带导数的,甚至可以带高阶导数。但是现在咱只讨论给定两个点和一阶导数来求三次函数的情况,搞明白这个的做法就自然能搞明白高阶高次做法。我就讨厌教科书里不该说的说一大堆,重点的就是提不到。我手上的《数值分析》就是蛋疼的书,说lagrange插值完全没说道重点,一堆式子都不知道哪儿就突然出现了,然后无聊的例子一大堆,明明都是可以自己举一反三的东西。嘛,不喷它了。
  其实,带不带导数并没有太大影响,因为我们的解决方法是构造等式,然后把一堆等式做成方程组来求解多项式的各项系数。无论是几阶导数,它最终都可以写成等式,只要给定的条件足够就可以求出足够高次的函数。比如我要说的三次hermite插值就是要求三次函数,三次函数有四个系数吧?于是我们需要四个或四个以上条件才可以求出它。这里我们有平面上两个点的坐标,这就是两个条件,还有这个两个坐标的切线斜率,也就是一阶导数在该点上的值,这也是两个条件。现在我们只要把这些条件写成等式做成方程组,用矩阵求解就好了。
  比如现在有两个点p₁=(x₁,y₁)和p₂=(x₂,y₂),然后这两点上的切线斜率分别是k₁和k₂。我们要求的是过p₁和p₂两点并满足给定切线斜率的一个三次函数。于是我们先来找等式,点的等式很容易找到吧?由于我们要求的是三次函数,因此只要把点的x和y值带入三次函数即可。而三次函数是过这个点的,所以可以划等号。比如p₁的等式就是y₁=ax₁³+bx₁²+cx₁+d,这里的a,b,c,d是三次函数的系数,也是我们最终要求得的东西。p₂和p₁是一样的我就不写了,这样我们就得到了两个等式。然后我们还知道这些点上的一阶导数值,要把它写成等式也不难,把三次函数求导了就行。比如k₁是p₁点的切线斜率,那么也就是这个三次函数的导数在p₁点的值,于是可以写成k₁=3ax₁²+2bx₁+c。这个二次函数就是三次函数的导数,它的导数在x₁上的值就是p₁的切线斜率,也就是k₁,于是可以划上等号。k₂和p₂也是同样的情况,这样我们就有了四个等式。 y₁=ax₁³+bx₁²+cx₁+d
y₂=ax₂³+bx₂²+cx₂+d
k₁=3ax₁²+2bx₁+c
k₂=3ax₂²+2bx₂+c
  把这四个等式视为方程组,就可以用矩阵法求出a,b,c,d这些系数了,这也味着三次函数被求出来了。不过要算矩阵也有多种方法,之前的多项式插值就直接用了高斯消除法把矩阵一点点的化成单位矩阵最终直接得出结果。其实还可以使用逆矩阵的方式来做,之前的文章我也有说过行列式的计算与矩阵求逆,求得逆矩阵后乘上y₁、y₂、k₁、k₂,这些值组成的列向量就可以得到a,b,c,d这些系数了。这回的例子我就是用逆矩阵的方式来做的。
  基本思路就是这些,下面是实现的核心代码,至于绘图代码我就不在文章贴出来了,可以在末尾的链接中看到。 function hermite(x1,y1,k1,x2,y2,k2){
  /**************************************************
  Author:次碳酸钴(admin@web-tinker.com)
  参数:
    x1,y1,k1 第一个点的x和y坐标与切线斜率
    x2,y2,k2 第二个点的x和y坐标与切线斜率
  返回:
    一个三次函数
  **************************************************/
  var e,i,j,m,n,r,o,s,f;
  e=[ //构造矩阵
    1,x1,x1*x1,x1*x1*x1,1,x2,x2*x2,x2*x2*x2,
    0,1,2*x1,3*x1*x1,0,1,2*x2,3*x2*x2
  ],f=function(e){ //三阶行列式计算
    return e[0]*e[4]*e[8]+e[1]*e[5]*e[6]+e[2]*e[3]*e[7]-
           e[2]*e[4]*e[6]-e[1]*e[3]*e[8]-e[0]*e[5]*e[7]
  };
  //计算矩阵的四阶行列式
  for(i=s=0;i<4;s+=f(o)*e[i]*(i%2?-1:1),i++)if(e[i])
    for(j=4,o=[];j<e.length;j++)if(j%4!=i)o.push(e[j]);
  //计算大矩阵部分,并除以上面的四阶行列式结果
  for(r=[],i=0;i<4;i++)
    for(j=0;j<4;r[i+j*4]=((i+j)%2?-1:1)*f(o)/s,j++)
      for(o=[],m=0;m<3;m++)for(n=0;n<3;n++)
        o.push(e[((i+m+1)%4)*4+(j+n+1)%4]);
  //计算三次函数的系数(矩阵乘向量)
  var a,b,c,d;
  d=y1*r[0]+y2*r[1]+k1*r[2]+k2*r[3];
  c=y1*r[4]+y2*r[5]+k1*r[6]+k2*r[7];
  b=y1*r[8]+y2*r[9]+k1*r[10]+k2*r[11];
  a=y1*r[12]+y2*r[13]+k1*r[14]+k2*r[15];
  //返回一个三次函数
  return function(x){return a*x*x*x+b*x*x+c*x+d};
};

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