2021-04-19 2021-04-19 597 words 3 minutes

Contents

插值方法

插值属于数值分析领域中的一种方法，是一种通过已知的离散的数据点，来拟合原函数根据给定的自变量估算因变量的方法。

常用的插值方法有很多，本文章给出三种常见的插值方法的实现。

使用语言: Python

使用下面的数据，预测函数在x=1处的值

1
2


x = np.array([0.5, 0.6, 0.4, 0.7])
y = np.array([-0.6931, -0.5108, -0.9163, -0.3567])

线性插值

线性插值及求一次多项式$p(x)$，满足$p(x_0), p(x_1) = y_1$ 可以根据点斜式方程求解即

$p(x) =y_0 \frac{y_{1}-y_{0}}{x_1-x_0}(x-x_0)$ 还可以将公式整理成如下形式

$p(x) = y_0 \frac{x - x_1}{x_0 - x_1} + y_1 \frac{x- x_0}{x_1 - x_0}$

我们令这里的$l_0(x)=\frac{x-x_1}{x_0 - x_1}, l_1(x)=\frac{x-x_0}{x_1-x_0}$

将其线性组合之后即为

$p(x)=y_0 l_0(x)+y_1 l_1(x)$

是Lagrange插值的特殊形式

此处给出线性插值的代码实现：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17


import numpy as np

def linear(x, y, x0):
    '''
    x : 为横坐标数组
    y : 为纵坐标数组
    x0: 为需要预测点的横坐标
    
    res: 预测结果
    '''
    res = y[1] + (y[1] - y[0])*(x0 - x[1])/(x[1] - x[0]) # 根据公式计算结果
    return res

x = np.array([0.5, 0.6, 0.4, 0.7])
y = np.array([-0.6931, -0.5108, -0.9163, -0.3567])
answer = linear(x, y, 1)
print(np.around(answer, 4)) # 保留四位小数

Largrange插值

根据Lagrange插值基函数$l_k(x)$，其满足如下性质当$i=k,l_k(x_i)=1$

当$i\not ={k}, l_k(x_i) = 0$

其中$l_k(x) = \prod_{i=0,i \not ={k}}^{n} \frac{x-x_i}{x_k-x_i}$

可以得到$p(x)=y_0 l_0(x) + y_1 l_1(x) + … + y_n l_n(x)$

$p(x)$满足$p(x_i)=y_i, i = 0,1,…,n$

即根据插值结点确定的方程，可以使得$p(x_i)=y_i$，是一种可行的插值方法，极大的提高了插值精度。

并且当只有两个插值结点时，Lagrange插值就退化成了线性插值，当有三个结点时，退化成抛物线插值。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25


import numpy as np
from matplotlib import pyplot as plt

def largrange(x, y, x0):
    '''
    x : x数组
    y : y数组
    x0: 预测点横坐标

    res: 预测结果
    '''
    n = len(x)
    res = 0
    for i in range(n):
        param = np.append(x[:i], x[(i+1):]) # 获得除去xi的所有横坐标
        numerator = (x0 - param).prod()     # 分子，使用prod求出矩阵中所有元素的积
        denominator = (x[i] - param).prod() # 计算分母
        res += y[i]*numerator/denominator   # 根据公式计算当前步骤对答案的贡献
    return res

x = np.array([0.5, 0.6, 0.4, 0.7])
y = np.array([-0.6931, -0.5108, -0.9163, -0.3567])
answer = largrange(x, y, 1)

print(np.around(answer, 4))

运行结果： 0.0634

Newton插值

我们引入差商的概念，设有函数$f(x), x_0,…,x_n$

$f[x_i, x_j] = \frac{f(x_j)-f(x_i)}{x_j - x_i}$

称为$f(x)$关于点$x_i, x_j$的一阶差商

$f[x_i, x_j, x_k] = \frac{f[x_j, x_k] - f[x_i, x_j]}{x_k - x_i}$

称为$f(x)$关于点$x_i， x_j， x_k$二阶差商

我们可以得到差商的一般定义，对于k阶差商

$f[x_0, x_1, …, x_k] = \frac{f[x_1, …, x_k] - f[x_0, …, x_{k-1}]}{x_k - x_0}$ 计算差商可以通过差商表来计算

$x_i$	$f(x_i)$	一阶差商	二阶差商	三阶差商	…	n阶差商
$x_0$	$f(x_0)$
$x_1$	$f(x_1)$	$f[x_0, x_1]$
$x_2$	$f(x_2)$	$f[x_1, x_2]$	$f[x_0, x_1, x_2]$
$x_3$	$f(x_3)$	$f[x_2, x_3]$	$f[x_1, x_2, x_3]$	$f[x_0, x_1, x_2, x_3]$
…	…	…	…	…	…
$x_n$	$f(x_n)$	$f[x_{n-1}, x_{n}]$	$f[x_{n-2}, x_{n-1}, x_{n}]$	$f[x_{n-3}, x_{n-2}, x_{n-1}, x_{n}]$		$f[x_0, x_1, …, x_n]$

将对角线上的差商值用来构造插值函数

$f(x) = f(x_0) + f[x_0, x_1](x- x_0) + f[x_0, x_1, x_2](x - x_0)(x - x_1) + … + f[x_0, …, x_n](x - x_0)…(x - x_{n-1}) + f[x,x_0,x_1,x_2,…,x_n](x - x_0)…(x - x_{n})$

其中$R(x) = f[x,x_0,x_1,x_2,…,x_n](x - x_0)…(x - x_{n})$为插值余项

综上，牛顿插值的最终表达式为

$R(x) = N(x)+R(x)$

牛顿插值的插值基函数有可继承的特点，可以优化计算方法

直接利用上一次的计算结果增量更新

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39


import numpy as np


# 计算差商表第j列的元素
def GetColumRes(arr, j):
    n = len(arr[:,0])
    colum = np.zeros(n)
    for i in range(j-1, n):
        colum[i] = (arr[i, j-1] - arr[i-1, j-1])/(arr[i, 0] - arr[i-j+1, 0])
    return colum

# 计算差商
def CalculateDividedDiffernces(x, y):
    n = len(x)
    arr = np.zeros((n, n+1))
    arr[0:n, 0] = x
    arr[0:n, 1] = y
    for i in range(2, n+1):
        # 获得差商表第i列的元素
        arr[:, i] = GetColumRes(arr, i)
    print(arr)
    return arr

def newton(x, y, x0):
    res = y[0]
    n = len(x)
    factor = 1
    DD = CalculateDividedDiffernces(x, y)
    for i in range(1,n):
        # factor因数迭乘
        factor *= (x0-x[i-1])
        # 取出对角线上的元素进行计算
        res += DD[i, i+1]*factor
    return res

x = np.array([0.5, 0.6, 0.4, 0.7])
y = np.array([-0.6931, -0.5108, -0.9163, -0.3567])
answer = newton(x, y, 1)
print(np.around(answer, 4))

运行结果：

差商表

[ 0.5 -0.6931 0. 0. 0. ]

[ 0.6 -0.5108 1.823 0. 0. ]

[ 0.4 -0.9163 2.0275 -2.045 0. ]

[ 0.7 -0.3567 1.86533333 -1.62166667 2.11666667]

运行结果 0.0634

三种插值方法及实现

插值方法

线性插值

Largrange插值

Newton插值