从拉格朗日插值谈到牛顿均差插值、切比雪夫插值

      在数值分析中，多项式插值用多项式对一组给定的数据进行插值。本文将从拉格朗日插值

讲到牛顿均差插值、切比雪夫插值。

    一、拉格朗日插值

      1.1、什么是拉格朗日插值？

                  假设有n个数据点 (x1, y1),...,(xn,yn) ，对于每一个1<=k<=n,定义 (n-1)阶多项式

                                   

                   则定义 n-1阶多项式

                                     

                  这个多项式就是拉格朗日插值多项式。而给定n个数据点，求出一个n-1阶(准确地说，不超过n-1阶)

            多项式，使得这个多项式对应的函数经过这n个数据点的过程就叫拉格朗日插值。

                  将任何一个数据点(xk,yk)(1<=k<=n)带入到拉格朗日多项式Pn-1(x)中，很明显Pn-1(xk)=0。

                  举个例子，对于三个数据点(0,1), (2, 2), (3, 4), 对应的拉格朗日插值多项式为

                          

                 即 2阶函数 f(x)=经过点(0, 1), (2, 2), (3, 4), 也就是f(x)插值这三个点。

                 这里需要注意的是，对于n个数据点，对应的拉格朗日插值多项式不一定是n-1阶，可能小于n-1阶。

                 例如，经过点(0, 2), (1, 1), (2, 0) 的拉格朗日插值多项式是 P(x) = -x + 2 ,是一阶的。

                 

     二、牛顿均差插值

       2.1 牛顿均差插值介绍

                 正如上面所示，拉格朗日插值法能够很方便写出经过n个点的插值多项式。但是在实际中，

         拉格朗日插值法很少用于计算，因为有多种可供选择的更加便于管理和计算复杂度更低的形式。

         牛顿均差插值法就是其中一种。

                给定n个点，定义n-1阶多项式

                         

           这个多项式就是牛顿均值插值多项式。

                  其中, 各项的系数可以有下面对应公式计算得到

                                   

            其它各项系数，以此类推。

                  上面计算各项系数的公式有点不好理解，说白了，可以将所有的系数存放在一个二维数组a里,

            其中系数从x的0次项到x的(n-1)次项的系数依次等于a[0][0], a[1][1], ...., a[n-1][n-1]。

            对于j=0, a[i][j]=f(xi)(0<=i<n); 当i>=1且j>=1时, a[i][j] = (a[i][j-1] - a[i-1][j-1])  /  (x(i) - x(i-j)) ;

                  这样可能还是有点抽象，下面举个例子来加深理解。

                  对于点(0, 1), (2, 2), (3, 4), 多项式各项系数对应的二维数组a，

                             a[0][0] = f(x0) = f(0) = 1,

                             a[1][0] = f(x1)  = f(2) = 2,     a[1][1] = (a[1][0] - a[0][0])  /  (x1 - x(1-0)) = (2-1)/(2-0) = 1/2;

                             a[2][0] = f(x2)  = f(3) = 4,     a[2][1] = (a[2][0] - a[1][0]) / (x2 - x(2-1)) = (4-2)/(3-2) = 2,

                                                                         a[2][2] = (a[2][1] - a[1][1]) / (x2 - x(2-2)) = (2-1/2) / (3-0) = 1/2;

                  则经过这三个点的牛顿均差插值多项式是

                                     P(x) =  1 + (x-0)/2 + (x-0)(x-2)/2

                  由上面的叙述，我们很容易就能用程序实现牛顿均差插值。

         2.2 牛顿均差插值程序实现

                  牛顿均差插值程序实现主要代码如下:

public class NewtonMean {
	private double[][] points;  //points to get theirs Newton Interpolation
	
	public NewtonMean(double[][] points) {
		this.points = points;
	}
	
	/* Get Newton Interpolation of points given */
	public double[] newtonMean() {
		if(points[0].length != 2) {
			throw new IllegalArgumentException("Each row of points can only " +
				"have two elements, the first for x, the other for y.");
		}		
		return computeCoefficients();		
	}

	/* Compute coefficient for each item of Polynomial */
	private double[] computeCoefficients() {
		/* Compute mean difference */
		double[][] meanDifference = new double[points.length][points.length];
		for(int i=0; i<meanDifference.length; i++) {
			meanDifference[i][0] = points[i][1];
		}						
		for(int i=0; i<meanDifference.length; i++) {
			for(int j=1; j<=i; j++) {
				meanDifference[i][j] = (meanDifference[i][j-1] - meanDifference[i-1][j-1]) / 
						(points[i][0] - points[i-j][0]);			    
			}		   
		}
		
		/* Get coefficient of each item */
		double[] coefficients = new double[points.length];
		for(int i=0; i<coefficients.length; i++) {
			coefficients[i] = meanDifference[i][i];
		}
		return coefficients;
    }	
}

      牛顿均差插值法相对拉格朗日插值法另一个优点是，增加一个节点时，牛顿均差插值法只需在当前得到的

插值多项式的基础上添加一项即可，而Lagrange 插值所有的基函数都需要重新计算。

    

  三、切比雪夫插值

      切比雪夫插值的提出是为了改善插值在差值区间上的最大误差控制

                      

      3.1 切比雪夫插值节点的选取

               对于一个插值区间 [a, b], 如果要在这个插值区间上选取n个点作为插值基点,使得上面的最大误差

         最小，则基点的选法如下:

                                               

               对于 i = 1, 2, ...., n. 下面的不等式在区间 [a, b] 上满足

                                                   

    3.2 切比雪夫插值程序实现                        

              上面已经实现了牛顿均差插值,切比雪夫插值只要在插值区间上选取出插值基点后，再调用

     牛顿均值差值法即可实现。

                    其代码实现如下(以函数 f(x)=sinx 的插值为例):                                 

public class ChebyshevInterpolation {
	
	private double rangeStart;  //Start value of range
	private double rangeEnd;    //End value of range
	private int numberOfBasicPoints;  
			
	public ChebyshevInterpolation(double rangeStart, double rangeEnd,
			int numberOfBasicPoints) {
		this.rangeStart = rangeStart;
		this.rangeEnd = rangeEnd;
		this.numberOfBasicPoints = numberOfBasicPoints;
	}

        /* Chebyshev Interpolation  */
	public double[] chebyshevInterpolation() {
		double[][] generatedBasicPoints = generateBasicPoints();
		NewtonMean newton = new NewtonMean(generatedBasicPoints);
		return newton.newtonMean();
	}		
	
	/* Generate basic points */	 
	private double[][] generateBasicPoints() {
		double[][] points = new double[numberOfBasicPoints][2];
		double scaling = (rangeEnd - rangeStart)/2;
		double displacement = (rangeStart + rangeEnd)/2;
		double factor = Math.PI/(2*numberOfBasicPoints);
		for(int i=0; i<points.length; i++) {
			points[i][0] = scaling * Math.cos((2*i + 1)*factor) + displacement ;
			points[i][1] = func(points[i][0]);
		}
		return points;
	} 
	
	/* Get function value of x */
	private double func(double x) {
		return Math.sin(x);
	}	
}