Change Dir

Commons Math學習筆記——函數插值

2.3 函數插值

看其他篇章到目錄選擇
在Commons Math中的analysis.interpolation包中有所有的與函數插值相關的類和接口定義。這一篇主要從這個包分析，來研究一下函數插值的應用。在2.1的api doc中添加了很多新的接口和類實現，但是2.0的source code里還是只有少量的實現。這里以2.0的source code為標準，輔助以2.1的api doc（其實這都是不影響的）。

插值是數學領域數值分析中的通過已知的離散數據求未知數據的過程或方法。給定n個離散數據點（稱為節點）(x_k,y_k)，k= 1,2,...,n。對于，求x所對應的y的值稱為內插。f(x)為定義在區間[a,b]上的函數。x₁,x₂,x₃...x_n為[a,b]上n個互不相同的點，G為給定的某意函數類。若G上有函數g(x)滿足： g(x_i) = f(x_i),k = 1,2,...n

則稱g(x)為f(x)關于節點x₁,x₂,x₃...x_n在G上的插值函數

在正式開始之前，不得不提在上一篇文章中講多項式函數時提到的PolynomialFunctionLagrangeForm和PolynomialFunctionNewtonForm，雖然把它們都放到了polynomials包里，但是其實它們完成的還是插值的工作。上一節沒有細講，這節來展開說一下。它們都是多項式插值，多項式插值用多項式對一組給定數據進行插值的過程。換句話說就是，對于一組給定的數據（如來自于采樣的數據），其目的就是尋找一個恰好通過這些數據點的多項式。

 1public 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

輸出結果：
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，用來表示一個單元實函數的插值。接口只定義了一個接口方法public UnivariateRealFunction interpolate(double xval[], double yval[])，用來構建插值函數，可以看到，這個方法返回一個單元實函數。這個接口下面有四個類來實現，表達四種插值方法：DividedDifferenceInterpolator, LoessInterpolator, NevilleInterpolator, SplineInterpolator。其中的NevilleInterpolator插值其實就是前面提到的拉格朗日插值，DividedDifferenceInterpolator調用的其實是polynomials包中的牛頓插值。我們略過不講，以樣條插值為例來講解interpolation包中的插值類。

樣條插值是使用一種名為樣條的特殊分段多項式進行插值的形式。由于樣條插值可以使用低階多項式樣條實現較小的插值誤差，這樣就避免了使用高階多項式所出現的龍格現象，所以樣條插值得到了流行。

SplineInterpolator類實現了UnivariateRealInterpolator接口，它實現的是一種三次樣條插值，它的interpolate(double []x, double []y)方法返回一個PolynomialSplineFunction對象，該樣條函數包含n個三次多項式，分段區分的標準就是輸入參數x數組就是樣條函數PolynomialSplineFunction中的knots數組（忘記了這一概念的同學請回到上一節）。實現源碼這里不貼了，因為這個類也就這一個方法，非常容易閱讀。它的實現來源是1989年的第四版《數值分析》，作者是R.L. Burden和 J.D. Faire。

插值代碼是保護性的，在插值前有檢驗參數的操作，如果參數正確，那么插值會拋出異常。

當然在2.1版本中，Math包又加入了新的插值方法，主要定義了新的多元實函數插值接口，實現類里加入了雙三次樣條插值等新方法，使得插值應用更加方便。

相關資料：

插值：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

龍格現象：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

posted on 2010-12-16 22:30 changedi 閱讀(4507) 評論(0)  編輯  收藏 所屬分類: 數學