如何求出插值表达式

一、原因

在实现拉格朗日插值时遇到一个问题，即：

当我们手算时用拉格朗日很容易得到具体表达式，但是当我们用机器实现时，却很难得到，因为机器不能表达出未知数x（虽然可以用TensorFlow的占位符，但杀鸡焉用牛刀），因此我们能很容易的得到具体的近似解，却无法得到表达式。

---

但是，python的scipy库中的lagrange实现了这一功能，如图：

---

所以我不禁思考如何实现这一功能，后来我换了一个思路，将系数当成方程组解，用高斯消元法求解。但是其机器的截断误差时不能避免的，所以应该还有其他更好的方法。

---

如果其他人能够有更好的思路，不吝赐教。

二、实现代码（python）

import numpy as npimport matplotlib.pyplot as pltfrom scipy.interpolate import lagrange'''    列主元高斯消元法    A:系数增广矩阵    n:未知数个数'''def main_element_gauss(A,n):    for i in range(0,n-1):        if(np.max(A[i:,i])!=A[i,i]):    #如果当前系数不是最大值，则列主元            temp_i=int(np.where(A==np.max(A[i:,i]))[0]) #temp_i为最大值所在的行索引            A[[i,temp_i],:]=A[[temp_i,i],:] #交换        for j in range(1,(n-i)):            A[i+j,:]=A[i+j,:]*A[i,i]-A[i,:]*(A[i+j,i]) #消元        #print("第%d次消元系数矩阵为：\n"%(i+1),A)    #回代得解    X=np.zeros((n,1))    for i in range(n-1,-1,-1):        temp=0        for j in range(0,n):            temp+=A[i,j]*X[j,0]        X[i,0]=(A[i,n]-temp)/A[i,i]    return X   #print("方程组的解为：\n",X)def main():    X=[float(i) for i in (input("请输入X的对应值：").split())]    Y=[float(i) for i in (input("请输入Y的对应值：").split())]    #scipy    #lagrange(X,Y)    X=np.array(X)    Y=np.array(Y).reshape(len(Y),1)    X1=np.zeros((len(X),len(X)))    #求出表达式的增广矩阵    for i in range(len(X)):        for j in range(len(X)):            X1[i,j]=X[i]**(len(X)-j-1)              C=np.hstack((X1,Y))    #得到系数解    R=main_element_gauss(C,len(X))    #在最小值至最大值区间取1000点    temp=np.linspace(np.min(X),np.max(X),1000,endpoint=True)    #将1000个点对应的y值求出，绘图；另把表达式用字符串表达    expression=np.zeros(temp.shape)    count=0    title='y='    tag=True #控制+号    for i in range(len(R)-1,-1,-1):        expression+=R[count]*(temp**i)        if(float(R[count])>0 ):            if(tag):                tag=False            else:                title=title+'+'        if(float(R[count])!=0):            if(i==0):                title=title+str(float(R[count]))            elif(i==1):                title=title+str(float(R[count]))+'x'            else:                title=title+str(float(R[count]))+'x**'+str(i)        count=count+1    plt.title(title)    plt.plot(X,Y,'or') #红点表示插的值    plt.plot(temp,expression)    plt.show()if __name__=='__main__':    main()

三、结果展示