song.yz@foxmail.com wechat: math-box

计算数学



牛顿插值

拉格朗日插值多项式每增加一个插值基点，原先计算的插值多项式pn(x) 对 pn+1(x) 没有用， 这样必然增加计算工作量，尤其在没有计算机的情况下这个问题更加突出。 我们期望增加插值基点时原先的计算结果对后面的计算仍然有用. 在介绍牛顿插值之前先要了解一下差商的基本概念。

f[xk]=f(xk)f[xk−1,xk]=f[xk]−f[xk−1]xk−1−xkf[xk−2,xk−1,xk]=f[xk−1,xk]−f[xk−2,xk−1]xk−2−xk⋮f[xk−j,xk−j+1,⋯,xk]=f[xk−j+1,⋯,xk]−f[xk−j,⋯,xk−1]xk−xk−j 有了上面的准备，便可以给出牛顿插值多项式为： Nn(x)=f[x0]+f[x0,x1](x−x0)+f[x0,x1,x2](x−x0)(x−x1)+⋯+f[x0,x1,⋯,xn](x−x0)(x−x1)⋯(x−xn−1) 如果记wk=(x−x0)(x−x1)⋯(x−xk−1),k=1,2,⋯,n 则牛顿插值可以表达为： Nn(x)=f[x0]+f[x0,x1]w1+f[x0,x1,x2]w2+⋯+f[x0,x1,⋯,xn]wn

下面给出牛顿插值的C#源代码。 


using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
class cnewton
/*--------------------------------------class comment
 Version   :  V1.0
 Coded by  :  syz
 Date      :  2011-07-16 21:33:04          *星期六*
----------------------------------------------------
Desciption :  牛顿插值法
 *
parameters :
 *   fout--------输出文件
 * 
Methods    :
 *  newton-------牛顿插值法（对多数据一次性完成插值）
 *  newton0------对单个点插值法
 *  drive--------驱动程序
--------------------------------------------------*/
{
    public StreamWriter fout = new StreamWriter("fout.txt");
    public VEC newton(VEC x, VEC y, VEC xx)
    /*------------------------------------------ comment
    Author      : syz 
    Date        : 2011-07-16 21:34:30         
    ---------------------------------------------------
    Desciption  : 牛顿插值法
     *
    Post Script :
     *   详细算法参看《Fortran科学计算与工程》
    paramtegers :
     *  x--------样本自变量
     *  y--------样本函数值
     *  xx-------要插值的向量
     *  返回值-------插值结果（维数同xx）
    -------------------------------------------------*/
    {
        int N = x.dim;
        int M = xx.dim;
        VEC yy=new VEC(M);
        for (int i = 0; i < M; i++)
            yy[i] = newton0(x, y, xx[i]);
        return yy;

    }
    public double newton0(VEC x, VEC y, double t)
    /*------------------------------------------ comment
    Author      : syz 
    Date        : 2011-07-16 21:36:41         
    ---------------------------------------------------
    Desciption  :  对单点牛顿插值方法
     *
    Post Script :
     *     编程时需要注意向量与矩阵的指标
    paramtegers :
     * x,y-------样本自变量与应变量
     * t---------要插值的点
    -------------------------------------------------*/
    {
        //获取样本维数
        int N = x.dim;        
        //声明矩阵
        MAT Q = new MAT(N, N);
        VEC b = new VEC(N);
        //把样本函数值复制给Q第一列
        for (int i = 0; i < N; i++)
            Q[i, 0] = y[i];
        for (int i = 1; i < N; i++)
            for (int j = 1; j <= i; j++)
                Q[i, j] = (Q[i, j - 1] - Q[i - 1, j - 1]) / (x[i] - x[i - j]);
        b[N - 1] = Q[N - 1, N - 1];
        for (int k = N - 1; k >= 1; k--)
            b[k - 1] = Q[k - 1, k - 1] + b[k] * (t - x[k - 1]);
        //返回运算结果
        return b[0];

    }
}