5.2 插值方法

from scipy import interpolate
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif']=['SimHei']    # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False      # 用来正常显示负号

# 代码5-4
x = np.array([200,220,250,270,280])
y = np.array([4,4.5,4.7,4.8,5.2])
la = interpolate.lagrange(x,y)
print('产量为240吨时，费用为：',la(240))

 产量为240吨时，费用为： 4.7142857142827665

# 代码5-5
# 自定义一阶跳跃差分函数
def diff_self (xi,k):
    '''
    xi：接收array。表示自变量x。无默认，不可省略。
    k：接收int。表示差分的次数，无默认，不可省略
    '''
    diffValue = []
    for i in range(len(xi)-k):
        diffValue.append(xi[i+k]-xi[i])
    return diffValue

diff_self(x,1)

[20, 30, 20, 10]

# 自定义求取差商函数
def diff_quot(xi,yi):
    '''
    xi：接收array。表示自变量x。无默认，不可省略。
    yi：接收array。表示因变量y。无默认，不可省略。
    '''
    length = len(xi)
    quot = []
    temp = yi
    for i in range(1,length): 
        #print("f(x)第" + str(i) + "阶差为：" + str(np.diff(temp,1)))
        #print("x第" + str(i) + "阶差为：" + str(diff_self(xi,i)))
        tem = np.diff(temp,1)/diff_self(xi,i)  # 此处需要numpy广播特性支持
        #print("第" + str(i) + "阶差商为：" +str(tem))
        #只有[0]位置的才是该阶差商，其余位置的是中间值
        quot.append(tem[0])
        temp = tem
    return(quot)

diff_quot(x,y)

[0.025, -0.0003666666666666666, 4.761904761904755e-06, 1.9047619047619106e-07]

# 自定义求取(x-x0)*(x-x1).....*(x-x0)
def get_Wi(k = 0, xi = []):
    '''
    xi：接收array。表示自变量x。无默认，不可省略。
    yi：接收array。表示因变量y。无默认，不可省略。
    '''
    def Wi(x):
        '''
        x：接收int，float，ndarray。表示插值节点。无默认。
        '''
        result = 1.0
        for each in range(k):
            if each == 0:
                strResult = str(result) + '*' + '(' + str(x) + '-' + str(xi[each]) + ')'
            else:
                strResult = str(strResult) + '*' + '(' + str(x) + '-' + str(xi[each]) + ')'
            #print(strResult)
            result *= (x - xi[each])
        #print(strResult)
        return result
    return Wi

get_Wi(5,x)(240)

-9600000.0 

# 自定义牛顿插值公式
def get_Newton_inter(xi,yi):
    '''
    xi：接收array。表示自变量x。无默认，不可省略。
    yi：接收array。表示因变量y。无默认，不可省略。
    '''
    diffQuot = diff_quot(xi,yi)
    def Newton_inter(x):
        '''
        x：接收int，float，ndarray。表示插值节点。无默认。
        '''
        result = yi[0]
        for i in range(0, len(xi)-1):
            result += get_Wi(i+1,xi)(x)*diffQuot[i]
        return result
    return Newton_inter

print('产量为240吨时，费用为：',get_Newton_inter(x,y)(240))

 产量为240吨时，费用为： 4.714285714285715

# 代码5-6
x1 = np.array([200,220,250,270,280,230])
y1 = np.array([4,4.5,4.7,4.8,5.2,4.6])
print('产量为240吨时，费用为：',get_Newton_inter(x1,y1)(240))

x0 = np.linspace(0,np.pi/2,4)
x1 = np.linspace(0,np.pi/2,20)
y0 = np.sin(x0)
y1 = np.sin(x1)
newton = [get_Newton_inter(x0,y0)(i) for i in x1]
plt.figure(figsize=(4,3))
plt.plot(x0,y0,'o')
plt.plot(x1,y1,'r-')
plt.plot(x1,newton,'b-')
plt.title('Newton插值sin函数')
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(0,1.7)
plt.ylim(0,1.2)
plt.legend(['插值点','sinx','Newton插值'])
plt.show()

#插值点
#plt.plot(x0,y0)
plt.plot(x0,y0,'D')

newton = [get_Newton_inter(x0,y0)(i) for i in x1]
#######################################################
#    get_Newton_inter(x0,y0)生成牛顿插值函数表达式    #
#    (i)作为新的自变量进入函数表达式                  #
#    生成一系统关于（i）的牛顿插值函数值              #
#    使用牛顿插值函数值拟合正弦函数                   #
#    画出x1和对应牛顿插值函数值                       #
#######################################################
plt.plot(x1,newton,'b-')

#真实曲线
plt.plot(x1,y1,'r-')