python插值法——牛顿法和拉格朗日法的代码实现

直接给出代码，可以看看。

import sympy as sy
from sympy import *
import matplotlib.pyplot as plt
import numpy as np
def f(a,b,n):#拉格朗日插值法
    global x,y#定义全局变量
    ty=ones(1,n+1);rt=0#给出初始空间
    x=symbols('x')#定义函数变量名
    for i in range(n+1):#确定第i个参数
        for j in range(n+1):
            if i != j:#分母不为零，即排除i
                ty[i]=ty[i]*(x-a[j])/(a[i]-a[j])#求积公式
            else:
                continue
    for p in range(n+1):#求最终函数表达式
        qw=b[p]*ty[p]#确定每个部分
        rt+=qw#求和
    return rt
def g(t1,b,n):#牛顿插值法
    global x,y#定义全局变量
    ty1=ones(1,n+1);ty2=ones(1,n+1);op=0;ty3=ones(1,n+1)
    a=1;rt=0;gt=0
    x=symbols('x')
    for m in range(n):
        a=a*(x-t1[m])
        ty2[m+1]=a
    ty1[0] = b[0]
    for j in range(n - gt):
        ty3[j]=(b[j+1]-b[j])/(t1[j+gt+1]-t1[j])
    gt = gt + 1
    for i in range(1,n+1):
        ty1[i]=ty3[0]
        if n-gt+1==0:
            break
        for j in range(n-gt):
            ty3[j]=(ty3[j+1]-ty3[j])/(t1[j+gt+1]-t1[j])
        gt=gt+1
    for p in range(n+1):
        qw=ty1[p]*ty2[p]
        rt+=qw#求和
    return rt
a=[0.4,0.5,0.6,0.7,0.8];t1=a#自变量
b=[-0.916291,-0.693147,-0.510826,-0.357765,-0.223144]#y值
n=4#几次插值
llp=0.54#需要估计的X值
print('拉格朗日插值：{}次插值法所求的表达式：\nf(x)={}'.format(n,f(a,b,n)))
print('x={}时，函数的估计值：{}\nx={}时，函数的准确值：{}'.format(llp,f(a,b,n).subs(x,llp),llp,float(sy.log(llp))))
print('x={}时，插值法的估计误差：{}'.format(llp,f(a,b,n).subs(x,0.54)-float(sy.log(0.54))))
print('拉格朗日插值：取（0.4,-0.916291）,(0.5,-0.693147),(0.6,-0.510826)时\n{}次插值法所求的表达式：\nf(x)={}'.format(2,f(a,b,2)))
print('x={}时，函数的估计值：{}\nx={}时，函数的准确值：{}'.format(llp,f(a,b,2).subs(x,llp),llp,float(sy.log(llp))))
print('x={}时，插值法的估计误差：{}'.format(llp,f(a,b,2).subs(x,0.54)-float(sy.log(0.54))))
print('拉格朗日插值：取（0.4,-0.916291）,(0.5,-0.693147)时\n{}次插值法所求的表达式：\nf(x)={}'.format(1,f(a,b,1)))
print('x={}时，函数的估计值：{}\nx={}时，函数的准确值：{}'.format(llp,f(a,b,1).subs(x,llp),llp,float(sy.log(llp))))
print('x={}时，插值法的估计误差：{}'.format(llp,f(a,b,1).subs(x,0.54)-float(sy.log(0.54))))
#牛顿
print('牛顿插值：{}次插值法所求的表达式：\nf(x)={}'.format(n,g(a,b,n)))
print('x={}时，函数的估计值：{}\nx={}时，函数的准确值：{}'.format(llp,g(a,b,n).subs(x,llp),llp,float(sy.log(llp))))
print('x={}时，插值法的估计误差：{}'.format(llp,g(a,b,n).subs(x,0.54)-float(sy.log(0.54))))
print('牛顿插值：取（0.4,-0.916291）,(0.5,-0.693147),(0.6,-0.510826)时\n{}次插值法所求的表达式：\nf(x)={}'.format(2,f(a,b,2)))
print('x={}时，函数的估计值：{}\nx={}时，函数的准确值：{}'.format(llp,g(a,b,2).subs(x,llp),llp,float(sy.log(llp))))
print('x={}时，插值法的估计误差：{}'.format(llp,g(a,b,2).subs(x,0.54)-float(sy.log(0.54))))
print('牛顿插值：取（0.4,-0.916291）,(0.5,-0.693147)时\n{}次插值法所求的表达式：\nf(x)={}'.format(1,f(a,b,1)))
print('x={}时，函数的估计值：{}\nx={}时，函数的准确值：{}'.format(llp,g(a,b,1).subs(x,llp),llp,float(sy.log(llp))))
print('x={}时，插值法的估计误差：{}'.format(llp,g(a,b,1).subs(x,0.54)-float(sy.log(0.54))))
x1=np.arange(0.01,4,0.01)
y1=np.log(x1)#真实值
x2=np.arange(0.01,1.6,0.01)
y2=[]#估计值
for i in x2:
    mio=f(a,b,n).subs(x,i)#估计函数值
    y2.append(mio)#加入到列表中
ax=plt.gca()#调用坐标轴
ax.spines['right'].set_color('none')#去掉右边框线
ax.spines['top'].set_color('none')#去掉顶部框线
ax.spines['bottom'].set_position(('data', 0))#原点对齐
ax.spines['left'].set_position(('data', 0))#原点对齐
plt.plot(x1,y1,linestyle='-',c='r',label='True')#真实值图像
plt.plot(x2,y2,linestyle='--',c='b',label='4Approximate')#估计值图像
y2=[]
for i in x2:
    mio=f(a,b,2).subs(x,i)#估计函数值
    y2.append(mio)#加入到列表中
plt.plot(x2,y2,linestyle='--',c='k',label='2Approximate')#估计值图像
y2=[]
for i in x2:
    mio=f(a,b,1).subs(x,i)#估计函数值
    y2.append(mio)#加入到列表中
plt.plot(x2,y2,linestyle='--',c='g',label='1Approximate')#估计值图像
plt.xlabel('X');plt.ylabel('f(X)')#坐标轴
plt.ylim(-3.5,2)#限制y轴的范围
plt.legend(bbox_to_anchor=(1,0.4))#显示标签
plt.show()#展示图像

---

运行的结果 

---

 

拉格朗日插值：4次插值法所求的表达式：
f(x)=38.1787916666667*(5.0 - 10.0*x)*(x - 0.8)*(x - 0.7)*(x - 0.6) + 115.5245*(x - 0.8)*(x - 0.7)*(x - 0.6)*(10.0*x - 4.0) - 255.413*(x - 0.8)*(x - 0.7)*(x - 0.5)*(5.0*x - 2.0) + 178.8825*(x - 0.8)*(x - 0.6)*(x - 0.5)*(3.33333333333333*x - 1.33333333333333) - 37.1906666666666*(x - 0.7)*(x - 0.6)*(x - 0.5)*(2.5*x - 1.0)
x=0.54时，函数的估计值：-0.615984011200000
x=0.54时，函数的准确值：-0.616186139423817
x=0.54时，插值法的估计误差：0.000202128223817155
拉格朗日插值：取（0.4,-0.916291）,(0.5,-0.693147),(0.6,-0.510826)时
2次插值法所求的表达式：
f(x)=4.581455*(5.0 - 10.0*x)*(x - 0.6) + 6.93147*(x - 0.6)*(10.0*x - 4.0) - 5.10826*(x - 0.5)*(5.0*x - 2.0)
x=0.54时，函数的估计值：-0.615319840000000
x=0.54时，函数的准确值：-0.616186139423817
x=0.54时，插值法的估计误差：0.000866299423817218
拉格朗日插值：取（0.4,-0.916291）,(0.5,-0.693147)时
1次插值法所求的表达式：
f(x)=2.23144*x - 1.808867
x=0.54时，函数的估计值：-0.603889400000000
x=0.54时，函数的准确值：-0.616186139423817
x=0.54时，插值法的估计误差：0.0122967394238168
牛顿插值：4次插值法所求的表达式：
f(x)=2.23144*x - 0.309583333333464*(x - 0.7)*(x - 0.6)*(x - 0.5)*(x - 0.4) + 1.92716666666668*(x - 0.6)*(x - 0.5)*(x - 0.4) - 2.04115*(x - 0.5)*(x - 0.4) - 1.808867
x=0.54时，函数的估计值：-0.615984011200000
x=0.54时，函数的准确值：-0.616186139423817
x=0.54时，插值法的估计误差：0.000202128223816933
牛顿插值：取（0.4,-0.916291）,(0.5,-0.693147),(0.6,-0.510826)时
2次插值法所求的表达式：
f(x)=4.581455*(5.0 - 10.0*x)*(x - 0.6) + 6.93147*(x - 0.6)*(10.0*x - 4.0) - 5.10826*(x - 0.5)*(5.0*x - 2.0)
x=0.54时，函数的估计值：-0.615319840000000
x=0.54时，函数的准确值：-0.616186139423817
x=0.54时，插值法的估计误差：0.000866299423816996
牛顿插值：取（0.4,-0.916291）,(0.5,-0.693147)时
1次插值法所求的表达式：
f(x)=2.23144*x - 1.808867
x=0.54时，函数的估计值：-0.603889400000000
x=0.54时，函数的准确值：-0.616186139423817
x=0.54时，插值法的估计误差：0.0122967394238170