数值计算方法实验之Hermite 多项式插值 (Python 代码)

一、实验目的

在已知f(x),x∈[a,b]的表达式，但函数值不便计算,或不知f(x),x∈[a,b]而又需要给出其在[a,b]上的值时，按插值原则f(x_i)= y_i(i= 0,1…….,n)求出简单函数P(x)(常是多项式),使其在插值基点x_i,处成立P(x_i)= y_i(i=0,1,……,n),而在[a,b]上的其它点处成立f(x)≈P(x).

二、实验原理

三、实验内容

求f(x)=x⁴在[0,2]上按5个等距节点确定的Hermite插值多项式.

四、实验程序

 1 import numpy as np
 2 from sympy import *
 3 import matplotlib.pyplot as plt
 4 
 5 
 6 def f(x):
 7     return  x ** 4
 8 
 9 
10 def ff(x):  # f[x0, x1, ..., xk]
11     ans = 0
12     for i in range(len(x)):
13         temp = 1
14         for j in range(len(x)):
15             if i != j:
16                 temp *= (x[i] - x[j])
17         ans += f(x[i]) / temp
18     return ans
19 
20 
21 def draw(L, newlabel= 'Lagrange插值函数'):
22     plt.rcParams['font.sans-serif'] = ['SimHei']
23     plt.rcParams['axes.unicode_minus'] = False
24     x = np.linspace(0, 2, 100)
25     y = f(x)
26     Ly = []
27     for xx in x:
28         Ly.append(L.subs(n, xx))
29     plt.plot(x, y, label='原函数')
30     plt.plot(x, Ly, label=newlabel)
31     plt.xlabel('x')
32     plt.ylabel('y')
33     plt.legend()
34 
35     plt.savefig('1.png')
36     plt.show()
37 
38 
39 def lossCal(L):
40     x = np.linspace(0, 2, 101)
41     y = f(x)
42     Ly = []
43     for xx in x:
44         Ly.append(L.subs(n, xx))
45     Ly = np.array(Ly)
46     temp = Ly - y
47     temp = abs(temp)
48     print(temp.mean())
49 
50 
51 def calM(P, x):
52     Y =   n ** 4
53     dfP = diff(P, n)
54     return solve(Y.subs(n, x[0]) - dfP.subs(n, x[0]), [m,])[0]
55 
56 
57 if __name__ == '__main__':
58     x = np.array(range(11)) - 5
59     y = f(x)
60 
61     n, m = symbols('n m')
62     init_printing(use_unicode=True)
63 
64     P = f(x[0])
65     for i in range(len(x)):
66         if i != len(x) - 1:
67             temp = ff(x[0:i + 2])
68         else:
69             temp = m
70         for j in x[0:i + 1]:
71             temp *= (n - j)
72         P += temp
73     P = expand(P)
74 
75     P = P.subs(m, calM(P, x))
76     draw(P, newlabel='Hermite插值多项式')
77     lossCal(P)

五、运算结果