插值复习(计算方法_python实现)

1. 题目

4、给出x,y数据位于机翼剖面的轮廓线上，y1和y2分别对应轮廓的上下线。假设需要得到x坐标每改变0.1时的y坐标。试完成加工所需数据，画出曲线，并求机翼剖面的面积。

2.算法描述

解题思路:运用学过的插值算法，对区间x所在的区间 [0,
15]内进行插值拟合，为了提高解题的准确度，把步长由0.1的150个点增加到了1000个点，具体做法为:
通过调用numpy库里的linspace函数随机产生了1000个点得等差数列，所求的结果存入另外两个列表里面，最后运用梯形积分求出不同插值情况下的面积,
算法流程图如下所示:

1.拉格朗日插值:

2.牛顿插值:

3.分段线性插值:

4.三次样条插值

3.源代码

#!/usr/bin/env python
# coding=utf-8

import matplotlib.pyplot as plt
import numpy as np

#牛顿插值
def newton_interpolation(X, Y, x):

    sum = Y[0]
    temp = np.zeros((len(X), len(X)))

    for i in range(0, len(X)):
        temp[i, 0] = Y[i]
    temp_sum = 1.0

    for i in range(1, len(X)):
        temp_sum = temp_sum * (x - X[i - 1])
        for j in range(i, len(X)):
            temp[j, i] = (temp[j, i - 1] - temp[j - 1, i - 1]) / (X[j] - X[j - i])
        sum += temp_sum * temp[i, i]
    return sum

#拉格朗日插值
def laGeLangRi(x, y, a) :
    ans = 0.0
    for i in range(len(y)) :
        t = y[i]
        for j in range(len(y)) :
            if i != j :
                t *= (a - x[j]) / (x[i] - x[j])
        ans += t
    return ans

#分段线性插值
def fenDuanXianXingChaZhi(X, Y, a):
    a = [a] if isinstance(a, (float, int)) else a
    x_loc = np.searchsorted(X, a)

    res_lst = []
    for x, i in zip(a, x_loc):
        L_i = Y[i - 1] * (x - X[i]) / (X[i - 1] - X[i]) + \
              Y[i] * (x - X[i - 1]) / (X[i] - X[i - 1])
        #res_lst.append(L_i)
        res_lst = L_i
    return res_lst

X = [0, 3, 5, 7, 9, 11, 12, 13, 14, 15]
Y = [0, 1.8, 2.2, 2.7, 3.0, 3.1, 2.9, 2.5, 2.0, 1.6]
Y1 = [0, 1.2, 1.7, 2.0, 2.1, 1.8, 1.2, 1.1, 1.0, 1.6]

import scipy.interpolate as spi

xs = np.linspace(np.min(X), np.max(X), 1000, endpoint=True)
#xs = np.arange(np.min(X), np.max(X), 0.1)

print(xs)
#三次样条插值
#ipo1 = spi.splrep(X, Y, k = 1)
#iy1 = spi.splev(xs, ipo1)

ipo3 = spi.splrep(X, Y, k = 3)
Y1_ipo3 = spi.splrep(X, Y1, k = 3)
san1 = spi.splev(xs, ipo3)
Y1_ipo3 = spi.splrep(X, Y1, k = 3)
san2 = spi.splev(xs, Y1_ipo3)

#A = get_diff_table(X, Y)
#df = pd.DataFrame(A)

#xs = np.linspace(np.min(X), np.max(X), 1000, endpoint=True)
#print(xs)
#xs = np.arange(np.min(X), np.max(X), 0.1)
niu1 = []
niu2 = []
la1 = []
la2 = []
f1 = []
f2 = []

for x in xs:
    niu1.append(newton_interpolation(X, Y, x))
    niu2.append(newton_interpolation(X, Y1, x))
    
    la1.append(laGeLangRi(X, Y, x))
    la2.append(laGeLangRi(X, Y1, x))
    f1.append(fenDuanXianXingChaZhi(X, Y, x))
    f2.append(fenDuanXianXingChaZhi(X, Y1, x))

#print(len(xs))
#print(f2)

test = []
test1 = []
test2 = []
test3 = []

tmp1 = []
#梯形面积公式
def Area() :
    sum, sum1, sum2, sum3 = 0.0, 0.0, 0.0, 0.0
    #for i in range(0, 27) :
    for i in range(0, len(xs)) :
        test1.append(la2[i] - la1[i])
        test2.append(niu2[i] - niu1[i])
        test.append(f1[i] - f2[i])
        test3.append(san1[i] - san2[i])

    for i in range(1, len(xs)) :
        sum += (xs[i] - xs[i - 1]) / 2.0 * (test[i - 1] + test[i])
        sum1 += (xs[i] - xs[i - 1]) / 2.0 * (test1[i - 1] + test1[i])
        sum2 += (xs[i] - xs[i - 1]) / 2.0 * (test2[i - 1] + test2[i])
        sum3 += (xs[i] - xs[i - 1]) / 2.0 * (test3[i - 1] + test3[i])

    tmp = 0.0
    for i in range(0, len(X)) :
        tmp1.append(Y[i] - Y1[i])
    for i in range(1, len(X)) :
        tmp += (X[i] - X[i - 1]) / 2.0 * (tmp1[i - 1] + tmp1[i])

    print('最初的面积     : %.15f' % tmp)
    print('拉格朗日插值   : %.15f' % sum1)
    print('牛顿插值       : %.15f' % sum2)
    print('分段插值       :', sum)
    print('三次样条插值   :', sum3)

Area()

#拉格朗日插值_图像
plt.title("laGeLangRi")
plt.plot(X, Y, 'r', label="original")
plt.plot(X, Y1, 'r', label = 'original')
plt.plot(xs, la1, 'b', linestyle = '--', label='la1')
plt.plot(xs, la2, 'b', linestyle = '--', label='la2')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.savefig('laGeLangRi')
plt.show()
#plt.clf()

#牛顿插值_图像
plt.title('Newton')
plt.plot(X, Y, 'r', label="original")
plt.plot(X, Y1, 'r', label = 'original')
plt.plot(xs, niu1, 'gray', linestyle = '--', linewidth = 3, label='niu1')
plt.plot(xs, niu2, 'gray', linestyle = '--', linewidth = 3, label='niu2')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.savefig('Newton')
plt.show()
#plt.clf()

#分段线性插值_图像
plt.title('fenDuanXianXing')
plt.plot(X, Y, 'r', label="original")
plt.plot(X, Y1, 'r', label = 'original')
plt.plot(xs, f1, 'gray', linestyle = '--', linewidth = 3, label='f1')
plt.plot(xs, f2, 'gray', linestyle = '--', linewidth = 3, label='f2')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.savefig('fenDuanXianXing')
plt.show()
#plt.clf()

#三次样条插值_图像
plt.title('sanCiYangTiao')
plt.plot(X, Y, 'b', label='original')
plt.plot(X, Y1, 'b', label='original')
##plt.scatter(new_x, iy1, label = 'Y = iy1')
plt.plot(xs, san1, 'r', linestyle = '--', label = 'san1')
plt.plot(xs, san2, 'r', linestyle = '--', label = 'san2')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.savefig('sanCiYangTiao')
plt.show()
plt.clf()

4.实验结果分析

应用梯形积分求得11个点所对应的面积为:11.750000000000000，对比拉格朗日插值、牛顿插值、分段线性插值、三次样条插值的图形和面积后可以得到，拉格朗日插值法和牛顿插值法产生了非常严重的龙格现象，与实际不符，拉格朗日插值法和牛顿插值法的收敛性无法得到保证， 二者求得的面积分别为:36.239333846993190, 36.239333846993254；分段线性插值的大致图形与原图相近，求得的面积为:11.749972945918891，但它得到的曲线不是很光滑，用带棱角的曲线制作飞机机翼是不可取的；为了解决在结点处保持光滑的问题，引入了三次样条插值，观察图形后可以发现，三次样条插值较为光滑，产生的误差较小，求得的面积为:12.240309390279387。