import numpy as np

 from sympy import *

 import matplotlib.pyplot as plt

 def f(x):

     return  x ** 4

 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(0, 2, 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(0, 2, 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 =   n ** 4

     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插值多项式')

     lossCal(P)