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包里，但是其实它们完成的还是插值的工作。上一节没有细讲，这节来展开说一下。它们都是多项式插值，多项式插值用多项式对一组给定数据进行插值的过程。换句话说就是，对于一组给定的数据（如来自于采样的数据），其目的就是寻找一个恰好通过这些数据点的多项式。

PolynomialFunctionLagrangeForm和PolynomialFunctionNewtonForm均实现UnivariateRealFunction接口，它们内部都定义了一个double coefficients[]用来表示多项式的系数。以拉格朗日形式为例来说明具体用法，它定义了double x[]和y[]表示采样的数据点。构造方法支持PolynomialFunctionLagrangeForm(double[] x, double[] y)，在实现接口方法value()的时候调用了evaluate()方法，这一方法的原型是static double evaluate(double[] x, double[] y, double z)。其中z就是未知变量，返回值则是插值。这一方法运用了Neville's Algorithm算法实现，具体见源代码：

public static double evaluate(double x[], double y[], double z) throws
2

DuplicateSampleAbscissaException, IllegalArgumentException {
3

int i, j, n, nearest = 0;
5

double value, c[], d[], tc, td, divider, w, dist, min_dist;
6

verifyInterpolationArray(x, y);
8

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

// initial approximation to the function value at z
26

value = y[nearest];
27

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

return value;
52

}
53

测试代码示例如下：

/**
2

*
3

*/
4

package algorithm.math;
5

import org.apache.commons.math.FunctionEvaluationException;
7

import org.apache.commons.math.MathException;
8

import org.apache.commons.math.analysis.UnivariateRealFunction;
9

import org.apache.commons.math.analysis.interpolation.SplineInterpolator;
10

import org.apache.commons.math.analysis.interpolation.UnivariateRealInterpolator;
11

import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
12

import org.apache.commons.math.analysis.polynomials.PolynomialFunctionLagrangeForm;
13

import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
14

/**
16

* @author Jia Yu
17

* @date 2010-11-21
18

*/
19

public class InterpolationTest {
20

/**
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

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

//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

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

}
82

输出结果：
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包又加入了新的插值方法，主要定义了新的多元实函数插值接口，实现类里加入了双三次样条插值等新方法，使得插值应用更加方便。

Change Dir

导航

公告

随笔分类(125)

随笔档案(123)

统计

留言簿(18)

积分与排名

“牛”们的博客

各个公司技术

我的链接

淘宝技术

阅读排行榜

评论排行榜

Commons Math学习笔记——函数插值

2.3 函数插值


只有注册用户登录后才能发表评论。




网站导航: 博客园 IT新闻 Chat2DB C++博客博问管理
相关文章: Logistic Regression——用线解决问题基数估计 Commons Math学习笔记——聚类和回归 Commons Math学习笔记——随机生成和统计初步 Commons Math学习笔记——分数和复数 Commons Math学习笔记——分布 Commons Math学习笔记——函数方程求解 Commons Math学习笔记——函数积分 Commons Math学习笔记——函数插值 Commons Math学习笔记——多项式函数