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。