算法名称：多项式内插法和外推法
Lagrangian算法描述：
我们知道，通过任意两点有且只有一条直线，通过任意三点只能确定唯一一个平面。通过N个点y0=f(x0),y1=f(x1),...,yN-1=f(xN-1)的N-1维内插多项式，可用拉格朗日（Lagrangian）多项式明确的表示出来：
（1）其中共有N个项，每一项都是N-1维多项式。对每一项除了一个xi使其值为yi外，对其他所有xi值均为零。然而，这个试子不够直观，于是我们改写成如下形式：
（2）但是，拉格朗日公示并不是非常好的方法，原因在于没有给出误差估计，并且不适合编程，原因在于乘法和除法的次数过多，精度损失可能比较严重。所以，后面我们将讨论Neville方法，进行插值求解。
Neville算法描述：
令P0是经过点（x0,y0）的x零阶多项式的值（即常数），故P0=y0。用类似的方法定义P1, P2, ..., PN-1。再令P01为过点（x0,y0）和（x1,y1）的x一阶多项式的值。用类似的方法定义P12, P23, ..., P(N-2)(N-1)。对高阶多项式进行类似定义，直至P012...(N-1)，它是对所有N个点的插值多项式的值，即所求的解。各个不同的P形成一张“表”，左边为“祖先”，一直到右边的一个“后代”。例如，当N=4时，
Neville方法用递推的方法填写表中的数据，每次一列，从左到右。它基于“子女”P和两个“双亲”之间的关系，得
(3)
因为两个双亲值已由点xi+1...xi+m-1得到，故该递推式可求。
对（3）式的一个改进是，记下双亲值与子女值之间的微笑“差值”，即定义（m=1, 2, ..., N-1）
(4)由（3）式很容易推出关系式
(5)对每一层m, 修正后的C和D都使插值高一级，也就是说，算出来的结果又存入了原来的位置。最后的解P0...(N-1)等于yi加上一些C和一些D，也可能只是一些C或者只是一些D。这些C和D相成一条穿过家族树到达右端子女的路径。
运行示例：
Data x:
32.0
22.2
41.6
10.1
50.5
--------------------------
Data f(x):
0.52992
0.37784
0.66393
0.17537
0.63608
--------------------------
P(x=27.5)
0.4575364991917163
--------------------------
示例程序：
package com.nc4nr.chapter03.neville;
public class PolInt ...{
double[] xa = ...{ 32.0, 22.2, 41.6, 10.1, 50.5 };
double[] ya = ...{ 0.52992, 0.37784, 0.66393, 0.17537, 0.63608 };
int count = 5;
private void polint(double x) ...{
int ns = 0, n = count;
double[] c = new double[n],
d = new double[n];
double dif = Math.abs(x-xa[0]);
for (int i = 0; i < n; i++) ...{
double dift = Math.abs(x-xa[i]);
if (dift < dif) ...{
ns = i;
dif = dift;
}
// 初始误差为yi，也就是说，后面产生的差值将累加到y上，最后产生结果
c[i] = ya[i];
d[i] = ya[i];
}
double y = ya[ns--]; // 从最接近的xi开始逼近
for (int m = 1; m < n; m++) ...{
for (int i = 0; i < n-m; i++) ...{
double ho = xa[i] - x,
hp = xa[i+m] - x,
w = c[i+1] - d[i],
den = ho - hp;
if (den == 0.0) System.out.println("polint: error.");
den = w/den;
d[i] = hp * den;
c[i] = ho * den;
}
double dy = (2*(ns+1)<(n-m) ? c[ns+1] : d[ns--]);
y = y + dy;
}
System.out.println(y);
System.out.println("--------------------------");
}
private void output(double[] v, int n) ...{
for (int i = 0; i < n; i++) ...{
System.out.println(v[i]);
}
System.out.println("--------------------------");
}
public PolInt() ...{
System.out.println("Data x:");
output(xa,count);
System.out.println("Data f(x):");
output(ya,count);
double x = 27.5;
System.out.println("P(x=" + x + ")");
polint(x);
}
public static void main(String[] args) ...{
new PolInt();
}
}
发表于 @ 2007年12月12日 13:14:00|评论(loading...)|编辑