实数插值方法及其 python 实现

地球被支点撬走啦

已于 2022-01-21 16:17:54 修改

阅读量1.4k

点赞数 1

分类专栏：数学/通信 Python 算法文章标签： python 机器学习

于 2022-01-21 15:38:23 首次发布

本文链接：https://blog.csdn.net/Flag_ing/article/details/122621908

版权

Python 同时被 3 个专栏收录

26 篇文章 9 订阅

订阅专栏

数学/通信

7 篇文章 3 订阅

订阅专栏

算法

7 篇文章 3 订阅

订阅专栏

本文主要介绍两种常用的实数插值方法：拉格朗日（Lagrange）插值 以及 牛顿（Newton）插值 及其python实现。运行效果如下：

1、拉格朗日插值

2、牛顿插值

3、python 实现

1、拉格朗日插值

拉格朗日（Lagrange）插值基函数：

N次插值多项式：

拉格朗日插值多项式的推到如下：

2、牛顿插值

Lagrange插值多项式计算简单便于编程计算，但是如果要增加一个节点，Lagrange插值多项式除增加的一项外，原来的每一项都需要改变。而Newton插值多项式在增加一个节点时只需要在原来的插值多项式中增加一项，而不需要改变原来的项。引出牛顿插值多项式首先要了解差商的概念。

Newton插值多项式为：

3、python 实现

import matplotlib.pyplot as plt
import numpy as np
import warnings


"""
 @功能：Lagrange 插值法 两点一次插值（线性插值）
 @参数：
    L1：需要进行插值的数据列表，待插值地方留空或置零
        如[0, 0, 1, 0, 0, 4, 0, 5], 代表第0, 1, 3、4、6个位置需要进行插值，其余位置为插值节点
        如果插值节点为0时，通过设置 nodeIndex 参数即可进行标定
    nodeIndex ：list2inter 中代表插值节点的位置索引，接上例，nodeIndex=[2,5,7]
"""
def Lagrange1(L1=[], nodeIndex=[]):
    assert (nodeIndex[0] >= 0) and (nodeIndex[-1] < len(L1)), '输入节点索引越界'
    assert len(nodeIndex) >= 2, '输入节点数必须不小于2'
    # 两端位置直接赋边缘值
    if nodeIndex[0] != 0:
        for i in range(nodeIndex[0]):
            L1[i] = L1[nodeIndex[0]]
    if nodeIndex[-1] != len(L1)-1:
        for i in range(nodeIndex[-1], len(L1)):
            L1[i] = L1[nodeIndex[-1]]
    # 需要插值的缝隙数量
    n_inter_space = len(nodeIndex)-1
    for i in range(n_inter_space):
        for j in range(nodeIndex[i]+1, nodeIndex[i+1]):
            # 点斜式  L(x) = y0 + [(y1-y0)/(x1-x0)](x-x0)
            L1[j] = L1[nodeIndex[i]] + \
                            (L1[nodeIndex[i+1]] - L1[nodeIndex[i]]) / (nodeIndex[i+1] - nodeIndex[i]) * \
                            (j - nodeIndex[i])
    return L1


"""
 @功能：Lagrange 三点二次平滑插值，对两端不进行插值
 @参数：
    L2：需要进行插值的数据列表，待插值地方留空或置零
        如[0, 0, 1, 0, 0, 4, 0, 5], 代表第0, 1, 3、4、6个位置需要进行插值，其余位置为插值节点
        如果插值节点为0时，通过设置 nodeIndex 参数即可进行标定
    nodeIndex ：list2inter 中代表插值节点的位置索引，接上例，nodeIndex=[2,5,7]
"""
def Lagrange2_smooth(L2=[], nodeIndex=[]):
    assert (nodeIndex[0] >= 0) and (nodeIndex[-1] < len(L2)), '输入节点索引越界'
    assert len(nodeIndex) >= 3, '输入节点数必须不小于3'
    # 两端位置直接赋边缘值
    if nodeIndex[0] != 0:
        for i in range(nodeIndex[0]):
            L2[i] = L2[nodeIndex[0]]
    if nodeIndex[-1] != len(L2)-1:
        for i in range(nodeIndex[-1], len(L2)):
            L2[i] = L2[nodeIndex[-1]]
    # 需要插值的缝隙数量
    n_inter_space = len(nodeIndex)-1
    for i in range(n_inter_space-1):
        for j in range(nodeIndex[i]+1, nodeIndex[i+2]):
            # 代入公式
            if i != 0 and j < nodeIndex[i+1]:
                L2[j] = 0.5 * \
                        (L2[j] + \
                        (j - nodeIndex[i+1]) * (j - nodeIndex[i+2]) / ((nodeIndex[i  ] - nodeIndex[i+1]) * (nodeIndex[i  ] - nodeIndex[i+2])) * L2[nodeIndex[i  ]] + \
                        (j - nodeIndex[i  ]) * (j - nodeIndex[i+2]) / ((nodeIndex[i+1] - nodeIndex[i  ]) * (nodeIndex[i+1] - nodeIndex[i+2])) * L2[nodeIndex[i+1]] + \
                        (j - nodeIndex[i  ]) * (j - nodeIndex[i+1]) / ((nodeIndex[i+2] - nodeIndex[i  ]) * (nodeIndex[i+2] - nodeIndex[i+1])) * L2[nodeIndex[i+2]] )
            else:
                L2[j] = (j - nodeIndex[i+1]) * (j - nodeIndex[i+2]) / ((nodeIndex[i  ] - nodeIndex[i+1]) * (nodeIndex[i  ] - nodeIndex[i+2])) * L2[nodeIndex[i  ]] + \
                        (j - nodeIndex[i  ]) * (j - nodeIndex[i+2]) / ((nodeIndex[i+1] - nodeIndex[i  ]) * (nodeIndex[i+1] - nodeIndex[i+2])) * L2[nodeIndex[i+1]] + \
                        (j - nodeIndex[i  ]) * (j - nodeIndex[i+1]) / ((nodeIndex[i+2] - nodeIndex[i  ]) * (nodeIndex[i+2] - nodeIndex[i+1])) * L2[nodeIndex[i+2]] 
    return L2


"""
 @功能：Lagrange 三点二次平滑插值,对两端也进行插值
 @参数：
    L2：需要进行插值的数据列表，待插值地方留空或置零
        如[0, 0, 1, 0, 0, 4, 0, 5], 代表第0, 1, 3、4、6个位置需要进行插值，其余位置为插值节点
        如果插值节点为0时，通过设置 nodeIndex 参数即可进行标定
    nodeIndex ：list2inter 中代表插值节点的位置索引，接上例，nodeIndex=[2,5,7]
"""
def Lagrange2_smooth_padding(L2=[], nodeIndex=[]):
    assert (nodeIndex[0] >= 0) and (nodeIndex[-1] < len(L2)), '输入节点索引越界'
    assert len(nodeIndex) >= 3, '输入节点数必须不小于3'
    # 需要插值的缝隙数量
    n_inter_space = len(nodeIndex)-1
    for i in range(n_inter_space-1):
        # 处理开头空白处的数据
        if nodeIndex[0] != 0 and i == 0:
            for j in range(nodeIndex[0]):
                L2[j] = (j - nodeIndex[1]) * (j - nodeIndex[2]) / ((nodeIndex[0] - nodeIndex[1]) * (nodeIndex[0] - nodeIndex[2])) * L2[nodeIndex[0]] + \
                        (j - nodeIndex[0]) * (j - nodeIndex[2]) / ((nodeIndex[1] - nodeIndex[0]) * (nodeIndex[1] - nodeIndex[2])) * L2[nodeIndex[1]] + \
                        (j - nodeIndex[0]) * (j - nodeIndex[1]) / ((nodeIndex[2] - nodeIndex[0]) * (nodeIndex[2] - nodeIndex[1])) * L2[nodeIndex[2]] 
        # 处理结尾空白处的数据
        if nodeIndex[-1] != len(L2)-1 and i == n_inter_space-2:
            for j in range(nodeIndex[-1], len(L2)):
                L2[j] = (j - nodeIndex[n_inter_space-1]) * (j - nodeIndex[n_inter_space  ]) / ((nodeIndex[n_inter_space-2] - nodeIndex[n_inter_space-1]) * (nodeIndex[n_inter_space-2  ] - nodeIndex[n_inter_space  ])) * L2[nodeIndex[n_inter_space-2]] + \
                        (j - nodeIndex[n_inter_space-2]) * (j - nodeIndex[n_inter_space  ]) / ((nodeIndex[n_inter_space-1] - nodeIndex[n_inter_space-2]) * (nodeIndex[n_inter_space-2+1] - nodeIndex[n_inter_space  ])) * L2[nodeIndex[n_inter_space-1]] + \
                        (j - nodeIndex[n_inter_space-2]) * (j - nodeIndex[n_inter_space-1]) / ((nodeIndex[n_inter_space  ] - nodeIndex[n_inter_space-2]) * (nodeIndex[n_inter_space-2+2] - nodeIndex[n_inter_space-1])) * L2[nodeIndex[n_inter_space  ]] 
        # 处理一般位置的数据
        for j in range(nodeIndex[i]+1, nodeIndex[i+2]):
            # 代入公式
            if i != 0 and j < nodeIndex[i+1]:
                L2[j] = 0.5 * \
                        (L2[j] + \
                        (j - nodeIndex[i+1]) * (j - nodeIndex[i+2]) / ((nodeIndex[i  ] - nodeIndex[i+1]) * (nodeIndex[i  ] - nodeIndex[i+2])) * L2[nodeIndex[i  ]] + \
                        (j - nodeIndex[i  ]) * (j - nodeIndex[i+2]) / ((nodeIndex[i+1] - nodeIndex[i  ]) * (nodeIndex[i+1] - nodeIndex[i+2])) * L2[nodeIndex[i+1]] + \
                        (j - nodeIndex[i  ]) * (j - nodeIndex[i+1]) / ((nodeIndex[i+2] - nodeIndex[i  ]) * (nodeIndex[i+2] - nodeIndex[i+1])) * L2[nodeIndex[i+2]] )
            else:
                L2[j] = (j - nodeIndex[i+1]) * (j - nodeIndex[i+2]) / ((nodeIndex[i  ] - nodeIndex[i+1]) * (nodeIndex[i  ] - nodeIndex[i+2])) * L2[nodeIndex[i  ]] + \
                        (j - nodeIndex[i  ]) * (j - nodeIndex[i+2]) / ((nodeIndex[i+1] - nodeIndex[i  ]) * (nodeIndex[i+1] - nodeIndex[i+2])) * L2[nodeIndex[i+1]] + \
                        (j - nodeIndex[i  ]) * (j - nodeIndex[i+1]) / ((nodeIndex[i+2] - nodeIndex[i  ]) * (nodeIndex[i+2] - nodeIndex[i+1])) * L2[nodeIndex[i+2]] 
    return L2

"""
 @功能：Lagrange 三点二次插值, 对两侧也进行插值
 @参数：
    L2：需要进行插值的数据列表，待插值地方留空或置零
        如[0, 0, 1, 0, 0, 4, 0, 5], 代表第0, 1, 3、4、6个位置需要进行插值，其余位置为插值节点
        如果插值节点为0时，通过设置 nodeIndex 参数即可进行标定
    nodeIndex ：list2inter 中代表插值节点的位置索引，接上例，nodeIndex=[2,5,7]
"""
def Lagrange2_padding(L2=[], nodeIndex=[]):
    assert (nodeIndex[0] >= 0) and (nodeIndex[-1] < len(L2)), '输入节点索引越界'
    assert len(nodeIndex) >= 3, '输入节点数必须不小于3'
    # 需要插值的缝隙数量
    n_inter_space = len(nodeIndex)-1
    for i in range(n_inter_space-1):
        # 处理开头空白处的数据
        if nodeIndex[0] != 0 and i == 0:
            for j in range(nodeIndex[0]):
                L2[j] = (j - nodeIndex[1]) * (j - nodeIndex[2]) / ((nodeIndex[0] - nodeIndex[1]) * (nodeIndex[0] - nodeIndex[2])) * L2[nodeIndex[0]] + \
                        (j - nodeIndex[0]) * (j - nodeIndex[2]) / ((nodeIndex[1] - nodeIndex[0]) * (nodeIndex[1] - nodeIndex[2])) * L2[nodeIndex[1]] + \
                        (j - nodeIndex[0]) * (j - nodeIndex[1]) / ((nodeIndex[2] - nodeIndex[0]) * (nodeIndex[2] - nodeIndex[1])) * L2[nodeIndex[2]] 
        # 处理结尾空白处的数据
        if nodeIndex[-1] != len(L2)-1 and i == n_inter_space-2:
            for j in range(nodeIndex[-1], len(L2)):
                L2[j] = (j - nodeIndex[n_inter_space-1]) * (j - nodeIndex[n_inter_space  ]) / ((nodeIndex[n_inter_space-2] - nodeIndex[n_inter_space-1]) * (nodeIndex[n_inter_space-2  ] - nodeIndex[n_inter_space  ])) * L2[nodeIndex[n_inter_space-2]] + \
                        (j - nodeIndex[n_inter_space-2]) * (j - nodeIndex[n_inter_space  ]) / ((nodeIndex[n_inter_space-1] - nodeIndex[n_inter_space-2]) * (nodeIndex[n_inter_space-2+1] - nodeIndex[n_inter_space  ])) * L2[nodeIndex[n_inter_space-1]] + \
                        (j - nodeIndex[n_inter_space-2]) * (j - nodeIndex[n_inter_space-1]) / ((nodeIndex[n_inter_space  ] - nodeIndex[n_inter_space-2]) * (nodeIndex[n_inter_space-2+2] - nodeIndex[n_inter_space-1])) * L2[nodeIndex[n_inter_space  ]] 

        for j in range(nodeIndex[i]+1, nodeIndex[i+2]):
            # 代入公式
            L2[j] = (j - nodeIndex[i+1]) * (j - nodeIndex[i+2]) / ((nodeIndex[i  ] - nodeIndex[i+1]) * (nodeIndex[i  ] - nodeIndex[i+2])) * L2[nodeIndex[i  ]] + \
                            (j - nodeIndex[i  ]) * (j - nodeIndex[i+2]) / ((nodeIndex[i+1] - nodeIndex[i  ]) * (nodeIndex[i+1] - nodeIndex[i+2])) * L2[nodeIndex[i+1]] + \
                            (j - nodeIndex[i  ]) * (j - nodeIndex[i+1]) / ((nodeIndex[i+2] - nodeIndex[i  ]) * (nodeIndex[i+2] - nodeIndex[i+1])) * L2[nodeIndex[i+2]] 
    return L2



"""
 @功能：Lagrange 三点二次插值
 @参数：
    L2：需要进行插值的数据列表，待插值地方留空或置零
        如[1, 0, 0, 4, 0, 5], 代表第2、3、5个位置需要进行插值，其余位置为插值节点
        如果插值节点为0时，通过设置 nodeIndex 参数即可进行标定
    nodeIndex ：list2inter 中代表插值节点的位置索引，接上例，nodeIndex=[0,3,5]
"""
def Lagrange2(L2=[], nodeIndex=[]):
    assert (nodeIndex[0] >= 0) and (nodeIndex[-1] < len(L2)), '输入节点索引越界'
    assert len(nodeIndex) >= 3, '输入节点数必须不小于3'
    # 两端位置直接赋边缘值
    if nodeIndex[0] != 0:
        for i in range(nodeIndex[0]):
            L2[i] = L2[nodeIndex[0]]
    if nodeIndex[-1] != len(L2)-1:
        for i in range(nodeIndex[-1], len(L2)):
            L2[i] = L2[nodeIndex[-1]]
    # 需要插值的缝隙数量
    n_inter_space = len(nodeIndex)-1
    for i in range(n_inter_space-1):
        for j in range(nodeIndex[i]+1, nodeIndex[i+2]):
            # 代入公式
            L2[j] = (j - nodeIndex[i+1]) * (j - nodeIndex[i+2]) / ((nodeIndex[i  ] - nodeIndex[i+1]) * (nodeIndex[i  ] - nodeIndex[i+2])) * L2[nodeIndex[i  ]] + \
                            (j - nodeIndex[i  ]) * (j - nodeIndex[i+2]) / ((nodeIndex[i+1] - nodeIndex[i  ]) * (nodeIndex[i+1] - nodeIndex[i+2])) * L2[nodeIndex[i+1]] + \
                            (j - nodeIndex[i  ]) * (j - nodeIndex[i+1]) / ((nodeIndex[i+2] - nodeIndex[i  ]) * (nodeIndex[i+2] - nodeIndex[i+1])) * L2[nodeIndex[i+2]] 
    return L2


"""
 @功能：计算 n 次差商
 @参数：
    list2inter：需要进行插值的数据列表，待插值地方留空或置零
                如[1, 0, 0, 4, 0, 5], 代表第2、3、5个位置需要进行插值，其余位置为插值节点
                如果插值节点为0时，通过设置 nodeIndex 参数即可进行标定
    nodeIndex ：list2inter 中代表插值节点的位置索引，接上例，nodeIndex=[0,3,5]
           n  ：差商的次数，需要 nodeIndex 中至少包含 n+1 个节点，多余的节点会自动忽略掉
                只计算前 n+1 个节点位置的差商值
"""
def difference_quotient(list2inter=[], nodeIndex=[], n=1):
    if n >= 7:
        warnings.warn("高次插值可能出现 Runge 现象")
    assert (nodeIndex[0] >= 0) and (nodeIndex[-1] < len(list2inter)), '输入节点索引越界'
    assert len(nodeIndex) != 0 or len(list2inter) != 0, '输入节点不能为空'
    # 只截取 n+1 个参数进行计算
    nodeIndex = nodeIndex[:n+1]
    def iter_body(nodeIndex=[]):
        if len(nodeIndex) == 1:
            return list2inter[nodeIndex[0]]
        left = nodeIndex[:-1]  # 取左边 n-1 个值
        right = nodeIndex[1:]  # 取右边 n-1 个值
        return (iter_body(left) - iter_body(right)) / (left[0] - right[-1])

    return iter_body(nodeIndex)


"""
 @功能：计算(x-x0)(x-x1)(x-x2)...(x-xn-1)
 @参数：
    nodeIndex ：list2inter 中代表插值节点的位置索引，接上例，nodeIndex=[0,3,5]
            x ：插值节点的索引值（变量值）
            n ：乘法的次数 如：n=3，则结果为 (x-x0)(x-x1)(x-x2)
"""
def mul_n(nodeIndex=[], x=0, n=0):
    assert len(nodeIndex) >= n, "输入节点数量达不到目标乘法次数的要求"
    result = 1
    for i in range(n):
        result = result * (x - nodeIndex[i])
    return result


"""
 @功能：计算 Newton 插值 建议最高使用 5 次差商，次数太高可能会出现 Runge 现象
 @参数：
    list2inter：需要进行插值的数据列表，待插值地方留空或置零
                如[1, 0, 0, 4, 0, 5], 代表第2、3、5个位置需要进行插值，其余位置为插值节点
                如果插值节点为0时，通过设置 nodeIndex 参数即可进行标定
    nodeIndex ：list2inter 中代表插值节点的位置索引，接上例，nodeIndex=[0,3,5]
           n  ：插值的次数，需要 nodeIndex 中至少包含 n+1 个节点，多余的节点会自动忽略掉
                只计算前 n+1 个节点位置的差商值
"""
def Newton_inter(list2inter=[], nodeIndex=[], n=1):
    assert len(nodeIndex) >= n+1, "输入节点数量达不到目标插值次数的要求"
    # 如果节点数大于插值所需，则取插值到第一个多余节点的前一个数值
    if len(nodeIndex) > n+1:
        num = nodeIndex[n+1]
    else:
        num = len(list2inter)

    dq = [0 for i in range(n+1)]
    for i in range(n+1):
        dq[i] = difference_quotient(list2inter, nodeIndex, i)

    for x in range(num):
        if x in nodeIndex:
            continue
        for j in range(n+1):
            list2inter[x] += dq[j] * mul_n(nodeIndex, x, j)

    return list2inter


L =         [0, 1, 0, 0, 0, 0, 8, 0, 0, 0, 0, 6, 0, 0, 18, 0, 0, 0, 12, 0, 0, 0]
nodeIndex = [   1,             6,            11,       14,          18]

temp = L.copy()
temp1 = L.copy()
temp2 = L.copy()
temp3 = L.copy()
temp4 = L.copy()

plt.figure()
plt.plot(nodeIndex, np.array(L)[nodeIndex], 'mo', label='L0')

inter1 = Lagrange1(L, nodeIndex=nodeIndex)
inter2 = Lagrange2(temp, nodeIndex=nodeIndex)
inter3 = Lagrange2_padding(temp4, nodeIndex=nodeIndex)
inter4 = Lagrange2_smooth(temp1, nodeIndex=nodeIndex)
inter5 = Lagrange2_smooth_padding(temp3, nodeIndex=nodeIndex)
inter6 = Newton_inter(temp2, nodeIndex=nodeIndex, n=4)

plt.plot(range(len(L)), inter1, 'r-', label='L1')
plt.plot(range(len(L)), inter2, 'c:', label='L2')
plt.plot(range(len(L)), inter3, 'g<', label='L2_padding')
plt.plot(range(len(L)), inter4, 'k-.', label='L2_smooth')
plt.plot(range(len(L)), inter5, 'r:', label='L2_smooth_padding')
plt.plot(range(len(L)), inter6, 'b--', label='Newton')

plt.legend(loc='best')
plt.show()

运行结果在文章开头已展示。

地球被支点撬走啦

关注

1
点赞
踩
2

收藏

觉得还不错? 一键收藏
打赏
0
评论
实数插值方法及其 python 实现

本文主要介绍两种常用的实数插值方法：拉格朗日（Lagrange）插值以及牛顿（Newton）插值及其python实现。运行效果如下：1、拉格朗日插值拉格朗日（Lagrange）插值基函数：N次插值多项式：拉格朗日插值多项式的推到如下：2、牛顿插值Lagrange插值多项式计算简单便于编程计算，但是如果要增加一个节点，Lagrange插值多项式除增加的一项外，原来的每一项都需要改变。而Newton插值多项式在增加一个节点时只需要在原来的插值多项式.
复制链接

扫一扫