埃尔米特插值（等距节点，只用一个点的导数构造n+1阶Hermite多项式）Python实现

柔光的暖阳◎ 2022-05-24 12:40 180阅读 0赞

## 函数 ##

y=11\+x2  y = 1 1 + x 2

## 埃尔米特插值 ##

埃尔米特多项式构造方法有很多种。  
这里只是用最简单的一种，通过均差来进行构造，最后再通过任意一个点的导数来计算出一个待定系数（**这里假设的是m**）。

下面代码中使用的是第一个点的导数相等来作为限制，算出这个多项式。

## 插值效果 ##

![这里写图片描述][70]

## 代码 ##

**下面代码就是根据之前的拉格朗日插值改进得到的。所以那个label那个就是用拉格朗日的做默认值**

有必要解释一下，**ff(x)**这个函数，只能传一个数组类型，然后算的这个数组的均差。

注意到，这里算出的来的跟之前直接用拉格朗日算出的来数值loss的绝对值均值基本相等。

import numpy as np
    from sympy import *
    import matplotlib.pyplot as plt
    
    
    def f(x):
        return 1 / (1 + x ** 2)
    
    
    def ff(x):  # f[x0, x1, ..., xk]
        ans = 0
        for i in range(len(x)):
            temp = 1
            for j in range(len(x)):
                if i != j:
                    temp *= (x[i] - x[j])
            ans += f(x[i]) / temp
        return ans
    
    
    def draw(L, newlabel= 'Lagrange插值函数'):
        plt.rcParams['font.sans-serif'] = ['SimHei']
        plt.rcParams['axes.unicode_minus'] = False
        x = np.linspace(-5, 5, 100)
        y = f(x)
        Ly = []
        for xx in x:
            Ly.append(L.subs(n, xx))
        plt.plot(x, y, label='原函数')
        plt.plot(x, Ly, label=newlabel)
        plt.xlabel('x')
        plt.ylabel('y')
        plt.legend()
    
        plt.savefig('1.png')
        plt.show()
    
    
    def lossCal(L):
        x = np.linspace(-5, 5, 101)
        y = f(x)
        Ly = []
        for xx in x:
            Ly.append(L.subs(n, xx))
        Ly = np.array(Ly)
        temp = Ly - y
        temp = abs(temp)
        print(temp.mean())
    
    
    def calM(P, x):
        Y = 1 / (1 + n ** 2)
        dfP = diff(P, n)
        return solve(Y.subs(n, x[0]) - dfP.subs(n, x[0]), [m,])[0]
    
    
    if __name__ == '__main__':
        x = np.array(range(11)) - 5
        y = f(x)
    
        n, m = symbols('n m')
        init_printing(use_unicode=True)
    
        P = f(x[0])
        for i in range(len(x)):
            if i != len(x) - 1:
                temp = ff(x[0:i + 2])
            else:
                temp = m
            for j in x[0:i + 1]:
                temp *= (n - j)
            P += temp
        P = expand(P)
    
        P = P.subs(m, calM(P, x))
        draw(P, newlabel='Hermite插值多项式（-5点导数相等）')
        lossCal(P)

[70]: /images/20220524/04d349b935924dacb9331575176f5f12.png