2.3 函數(shù)插值
看其他篇章到目錄選擇
在Commons Math中的analysis.interpolation包中有所有的與函數(shù)插值相關(guān)的類和接口定義。這一篇主要從這個包分析，來研究一下函數(shù)插值的應(yīng)用。在2.1的api doc中添加了很多新的接口和類實現(xiàn)，但是2.0的source code里還是只有少量的實現(xiàn)。這里以2.0的source code為標準，輔助以2.1的api doc（其實這都是不影響的）。
插值是數(shù)學領(lǐng)域數(shù)值分析中的通過已知的離散數(shù)據(jù)求未知數(shù)據(jù)的過程或方法。給定n個離散數(shù)據(jù)點（稱為節(jié)點）(xk,yk)，k= 1,2,...,n。對于，求x所對應(yīng)的y的值稱為內(nèi)插。f(x)為定義在區(qū)間[a,b]上的函數(shù)。x1,x2,x3...xn為[a,b]上n個互不相同的點，G為給定的某意函數(shù)類。若G上有函數(shù)g(x)滿足： g(xi) = f(xi),k = 1,2,...n
則稱g(x)為f(x)關(guān)于節(jié)點x1,x2,x3...xn在G上的插值函數(shù)
在正式開始之前，不得不提在上一篇文章中講多項式函數(shù)時提到的PolynomialFunctionLagrangeForm和PolynomialFunctionNewtonForm，雖然把它們都放到了polynomials包里，但是其實它們完成的還是插值的工作。上一節(jié)沒有細講，這節(jié)來展開說一下。它們都是多項式插值，多項式插值用多項式對一組給定數(shù)據(jù)進行插值的過程。換句話說就是，對于一組給定的數(shù)據(jù)（如來自于采樣的數(shù)據(jù)），其目的就是尋找一個恰好通過這些數(shù)據(jù)點的多項式。
PolynomialFunctionLagrangeForm和PolynomialFunctionNewtonForm均實現(xiàn)UnivariateRealFunction接口，它們內(nèi)部都定義了一個double coefficients[]用來表示多項式的系數(shù)。以拉格朗日形式為例來說明具體用法，它定義了double x[]和y[]表示采樣的數(shù)據(jù)點。構(gòu)造方法支持PolynomialFunctionLagrangeForm(double[] x, double[] y)，在實現(xiàn)接口方法value()的時候調(diào)用了evaluate()方法，這一方法的原型是static double evaluate(double[] x, double[] y, double z)。其中z就是未知變量，返回值則是插值。這一方法運用了Neville's Algorithm算法實現(xiàn)，具體見源代碼：
1
public static double evaluate(double x[], double y[], double z) throws
2
DuplicateSampleAbscissaException, IllegalArgumentException 
{
3
4 int i, j, n, nearest = 0;
5 double value, c[], d[], tc, td, divider, w, dist, min_dist;
6
7 verifyInterpolationArray(x, y);
8
9 n = x.length;
10 c = new double[n];
11 d = new double[n];
12 min_dist = Double.POSITIVE_INFINITY;
13
for (i = 0; i < n; i++) {
14 // initialize the difference arrays
15 c[i] = y[i];
16 d[i] = y[i];
17 // find out the abscissa closest to z
18 dist = Math.abs(z - x[i]);
19 if (dist < min_dist) {
20 nearest = i;
21 min_dist = dist;
22
}
23 }
24
25 // initial approximation to the function value at z
26 value = y[nearest];
27
28 for (i = 1; i < n; i++) {
29 for (j = 0; j < n-i; j++) {
30 tc = x[j] - z;
31 td = x[i+j] - z;
32 divider = x[j] - x[i+j];
33 if (divider == 0.0) {
34 // This happens only when two abscissas are identical.
35 throw new DuplicateSampleAbscissaException(x[i], i, i+j);
36 }
37 // update the difference arrays
38 w = (c[j+1] - d[j]) / divider;
39 c[j] = tc * w;
40 d[j] = td * w;
41 }
42 // sum up the difference terms to get the final value
43 if (nearest < 0.5*(n-i+1)) {
44 value += c[nearest]; // fork down
45 } else {
46 nearest--;
47 value += d[nearest]; // fork up
48 }
49 }
50
51 return value;
52
}
53
54
55
測試代碼示例如下：
1/** *//**
2 *
3 */
4package algorithm.math;
5
6import org.apache.commons.math.FunctionEvaluationException;
7import org.apache.commons.math.MathException;
8import org.apache.commons.math.analysis.UnivariateRealFunction;
9import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
10import org.apache.commons.math.analysis.interpolation.UnivariateRealInterpolator;
11import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
12import org.apache.commons.math.analysis.polynomials.PolynomialFunctionLagrangeForm;
13import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
14
15/** *//**
16 * @author Jia Yu
17 * @date 2010-11-21
18 */
19public class InterpolationTest {
20
21 /** *//**
22 * @param args
23 */
24 public static void main(String[] args) {
25 // TODO Auto-generated method stub
26 polynomialsInterpolation();
27 System.out.println("-------------------------------------------");
28 interpolatioin();
29 }
30
31 private static void interpolatioin() {
32 // TODO Auto-generated method stub
33 // double x[] = { 0.0, 0.5, 1.0 };
34 // double y[] = { 0.0, 0.5, 1.0 };
35 double x[] = { 0.0, Math.PI / 6d, Math.PI / 2d, 5d * Math.PI / 6d,
36 Math.PI, 7d * Math.PI / 6d, 3d * Math.PI / 2d,
37 11d * Math.PI / 6d, 2.d * Math.PI };
38 double y[] = { 0d, 0.5d, 1d, 0.5d, 0d, -0.5d, -1d, -0.5d, 0d };
39 UnivariateRealInterpolator i = new SplineInterpolator();
40 UnivariateRealFunction f = null;
41 // interpolate y when x = 0.5
42 try {
43 f = i.interpolate(x, y);
44 System.out.println("when x = 0.5, y = " + f.value(0.5));
45 } catch (MathException e) {
46 // TODO Auto-generated catch block
47 e.printStackTrace();
48 }
49
50 //check polynomials functions
51 PolynomialFunction polynomials[] = ((PolynomialSplineFunction) f).getPolynomials();
52 for(int j=0;j<polynomials.length;j++){
53 System.out.println("cubic spline:f"+j+"(x) = "+polynomials[j]);
54 }
55 }
56
57 private static void polynomialsInterpolation() {
58 // TODO Auto-generated method stub
59 double x[] = { 0.0, -1.0, 0.5 };
60 double y[] = { -3.0, -6.0, 0.0 };
61 PolynomialFunctionLagrangeForm p = new PolynomialFunctionLagrangeForm(
62 x, y);
63 // output directly
64 System.out.println("ugly output is " + p);
65 // interpolate y when x = 1.0
66 try {
67 System.out.println("when x = 1.0, y = " + p.value(1.0));
68 } catch (FunctionEvaluationException e) {
69 // TODO Auto-generated catch block
70 e.printStackTrace();
71 }
72 // degree
73 System.out.println("polynomial degree is " + p.degree());
74 // coefficients
75 for (int i = 0; i < p.getCoefficients().length; i++) {
76 System.out.println("coeff[" + i + "] is " + p.getCoefficients()[i]);
77 }
78 //
79 }
80
81}
82
83
輸出結(jié)果：
ugly output is org.apache.commons.math.analysis.polynomials.PolynomialFunctionLagrangeForm@1b67f74
when x = 1.0, y = 4.0
polynomial degree is 2
coeff[0] is -3.0
coeff[1] is 5.0
coeff[2] is 2.0
-------------------------------------------
when x = 0.5, y = 0.47956828499706067
cubic spline:f0(x) = 1.0026761414789407 x - 0.17415828593927732 x^3
cubic spline:f1(x) = 0.5 + 0.8594366926962349 x - 0.2735671958343121 x^2 - 0.08707914296963848 x^3
cubic spline:f2(x) = 1.0 + 5.551115123125783E-17 x - 0.5471343916686237 x^2 + 0.08707914296963837 x^3
cubic spline:f3(x) = 0.5 - 0.859436692696235 x - 0.27356719583431244 x^2 + 0.17415828593927762 x^3
cubic spline:f4(x) = -1.002676141478941 x + 6.938893903907228E-17 x^2 + 0.1741582859392774 x^3
cubic spline:f5(x) = -0.5 - 0.8594366926962351 x + 0.2735671958343122 x^2 + 0.08707914296963865 x^3
cubic spline:f6(x) = -1.0 + 3.3306690738754696E-16 x + 0.5471343916686243 x^2 - 0.08707914296963899 x^3
cubic spline:f7(x) = -0.5 + 0.8594366926962345 x + 0.27356719583431127 x^2 - 0.1741582859392767 x^3
PolynomialFunctionLagrangeForm類沒有重寫toString方法。這個類的兩個亮點是計算coefficients的算法（這里不貼了，有興趣的讀源碼）和Neville算法。
下面該正式的interpolation包了，這個包在2.0的source code里只有一個接口定義UnivariateRealInterpolator，用來表示一個單元實函數(shù)的插值。接口只定義了一個接口方法public UnivariateRealFunction interpolate(double xval[], double yval[])，用來構(gòu)建插值函數(shù)，可以看到，這個方法返回一個單元實函數(shù)。這個接口下面有四個類來實現(xiàn)，表達四種插值方法：DividedDifferenceInterpolator, LoessInterpolator, NevilleInterpolator, SplineInterpolator。其中的NevilleInterpolator插值其實就是前面提到的拉格朗日插值，DividedDifferenceInterpolator調(diào)用的其實是polynomials包中的牛頓插值。我們略過不講，以樣條插值為例來講解interpolation包中的插值類。
樣條插值是使用一種名為樣條的特殊分段多項式進行插值的形式。由于樣條插值可以使用低階多項式樣條實現(xiàn)較小的插值誤差，這樣就避免了使用高階多項式所出現(xiàn)的龍格現(xiàn)象，所以樣條插值得到了流行。
SplineInterpolator類實現(xiàn)了UnivariateRealInterpolator接口，它實現(xiàn)的是一種三次樣條插值，它的interpolate(double []x, double []y)方法返回一個PolynomialSplineFunction對象，該樣條函數(shù)包含n個三次多項式，分段區(qū)分的標準就是輸入?yún)?shù)x數(shù)組就是樣條函數(shù)PolynomialSplineFunction中的knots數(shù)組（忘記了這一概念的同學請回到上一節(jié)）。實現(xiàn)源碼這里不貼了，因為這個類也就這一個方法，非常容易閱讀。它的實現(xiàn)來源是1989年的第四版《數(shù)值分析》，作者是R.L. Burden和 J.D. Faire。
插值代碼是保護性的，在插值前有檢驗參數(shù)的操作，如果參數(shù)正確，那么插值會拋出異常。
當然在2.1版本中，Math包又加入了新的插值方法，主要定義了新的多元實函數(shù)插值接口，實現(xiàn)類里加入了雙三次樣條插值等新方法，使得插值應(yīng)用更加方便。
相關(guān)資料：
插值：http://zh.wikipedia.org/zh-cn/%E6%8F%92%E5%80%BC
多項式插值：http://zh.wikipedia.org/zh-cn/%E5%A4%9A%E9%A1%B9%E5%BC%8F%E6%8F%92%E5%80%BC
拉格朗日插值法：http://zh.wikipedia.org/zh-cn/%E6%8B%89%E6%A0%BC%E6%9C%97%E6%97%A5%E5%A4%9A%E9%A1%B9%E5%BC%8F
Neville's Algorithm：http://mathworld.wolfram.com/NevillesAlgorithm.html
樣條插值：http://zh.wikipedia.org/zh-cn/%E6%A0%B7%E6%9D%A1%E6%8F%92%E5%80%BC
龍格現(xiàn)象：http://zh.wikipedia.org/zh-cn/%E9%BE%99%E6%A0%BC%E7%8E%B0%E8%B1%A1
Commons math包：http://commons.apache.org/math/index.html