数值分析(曲线插值)

目录
任务描述及要求 (3)
1 算法设计方案 (4)
1.1 总算法概述 (4)
1.2 插值区间规整化 (4)
1.3 多项式牛顿插值 (5)
1.4 三次样条插值 (5)
2 全部源程序 (7)
2.1 全局常量和变量 (7)
2.2 分段线性插值 (7)
2.3 分段二次多项式插值 (8)
2.4 分段三次多项式插值 (9)
2.5 三次样条插值 (9)
2.6 OpenGL环境下显示曲线 (11)
3各种插值方法细化后数据表 (13)
3.1 分段线性插值法 (13)
3.2 分段二次多项式插值 (18)
3.3 分段三次多项式插值 (23)
3.4 三次样条插值 (29)
4 插值机翼轮廓线演示 (35)
4.1 分段线性插值 (35)
4.2 分段二次多项式插值 (35)
4.3 分段三次多项式插值 (35)
4.4 三次样条插值 (36)
5 不同插值方法的优劣 (37)
5.1 分段线性插值 (37)
5.2 分段高次多项式插值 (37)
5.3 三次样条插值 (37)
任务描述及要求
二、在飞机制造业中，机翼的加工是一项关键技术。

由于机翼的尺寸很大，通常在图纸中只
能标出某些关键点的数据。

下表给出的是某型号飞机的机翼上缘轮廓线的部分数据。

但是，在使用数控机床加工机翼时，由于机床走刀只能沿x方向和y方向走非常小的步长，因此机床编程时需要计算出轮廓线上x坐标每改变1个单位时y的相应坐标。

请根据加工要求分别用分段线性插值法、分段二次多项式插值法、分段三次多项式插值法和三次样条插值法，对上表中的数据进行细化。

打印内容：
1.算法的设计方案。

2.全部源程序（要求注明主程序和每个子程序的功能）。

3.用各种插值方法细化后的全部数据表。

4.画出用各种插值方法细化后的机翼轮廓线。

5.讨论不同插值方法的优劣。

1 算法设计方案
1.1 总算法概述
本次任务可以看作是完成两种类型的曲线插值：多项式插值和样条插值。

下面分别介绍这两种插值的算法。

两种方法里需要统一解决的问题是插值步长和插值区段的划分。

由于机床加工时的步长应该为固定的，不能恰好达到某些设计型值点（如x = 4.74在步长为0.5时不会被取到）。

本次任务中对插值区段的规整化进行了算法上的研究（当然，拟合曲线在计算上仍然通过型值点），目的在于使同一个区段上得到的不同插值结果可以快速融合优化。

例如，高次多项式插值根据计算区间的选取不同，在同一区间内可以有不同的插值结果，将这些结果融合也许可以获得更佳的结果。

插值区间规整化后，实现这种融合更加方便了。

本次任务中，步长选为0.5。

本次任务的图像表达利用了OpenGL的glut库进行了实现。

1.2 插值区间规整化
v o i d h e a d_e n d(f l o a t x1,f l o a t x2,f l o a t*h e a d,f l o a t*e n d)//确定本段取值的起末参数
{
i f(x1-f l o o r(x1)>0.5)*h e a d=c e i l(x1);
e l s e
i f(x1==f l o o r(x1))
*h e a d=x1;
e l s e
*h e a d=f l o o r(x1)+0.5;
i f(x2-f l o o r(x2)>0.5)*e n d=c e i l(x2)-0.5;
e l s e
i f(x2==f l o o r(x2))
*e n d=x2-0.5;
e l s e
*e n d=f l o o r(x2);
}
进行规整后，区间划分情况变为：
[0.0, 4.5]，[5.0, 9.0]，[9.5, 18.5]，[19.0, 37.5]，[38.0, 56.5]，[57.0, 75.5]，[76.0,94.5]，
[95.0, 113.5]，[114.0, 132.5]，[133.0, 151.5]，[152.0, 170.5]，[171.0,189.5]
1.3 多项式牛顿插值
任务中的分段线性插值、分段二次多项式插值、分段三次多项式插值均可以用牛顿插值法来实现。

牛顿插值：
[]001001201011()(),()[,,]()()[,
,]()
()
n n n p x f x f x x x x f x x x x x x x f x x x x x x -=+-+--+
+-- （1.1）
当n 分别为1、2、3时，就分为为所需的线性、二次、三次多项式插值公式。

插值涉及区间分别为n 个连续区间。

具体实现时，将连续区间的端点坐标传入函数，函数自动决定读取后面的若干节点并调用牛顿插值公式计算每个步长位置的y 坐标，把相应结果存于一数组中，供输出文件和OpenGL 绘制时调用。

1.4 三次样条插值
由课本P112页，得三次样本插值函数：
33111111()()()66()()6()(),(1,2,,)
6
i i i i i i i i i i i i i i i i i i M M
s x x x x x h h y M
h x x h y M
h x x x x x i n h ------=
-+-+--+--≤≤= （1.2）
其中：1i i i h x x -=-
而M 的解法为解线性方程组：
112(1,2,,1)i i i i i i
M M M i n γαβ-+++==- （1.3）
其中： 11
i i i i h h h α++=
+ 1i i γα=- 11116
()i i i i i i i i i y y y y h h h h β+-++--=
-+ 任务完成选用了第一种边界条件：
00n αγ== 00n ββ== 00n M M == 于是，线性方程组表达为：
11111
1112002
2
200n n n n M M γαβγαβ----⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣
⎦⎣⎦⎣⎦
（1.4） 由（1.3）和（1.4）得到M 后，将位置x 带入（1.2）即可得x 处的y 值。

任务用了函数initM( )来构造出线性方程组，用函数getM( )利用简化的高斯消去法解得各M 值。

2 全部源程序
程序的基本处理单元是四种插值办法，然后是图形显示操作。

以下分别为相关代码2.1 全局常量和变量
#d e f i n e N13//型值点个数
#d e f i n e d e l t a0.5//步长
F I L E*f i l e;//文件指针
f l o a t X[13]={0.00,4.74,9.50,19.00,38.00,57.00,76.00,
95.00,114.00,133.00,152.00,171.00,190.00};//型值点x坐标
f l o a t Y[13]={0.00,5.32,8.10,11.97,16.15,17.10,16.34,
14.63,12.16,9.69,7.03,3.99,0.00};//型值点y坐标
f l o a t R[(i n t)(190/d e l t a)]={0};//插值结果数组（仅存y坐标）
f l o a t M[13][13]={0};//样条插值线性方程组系数矩阵
f l o a t B[13]={0};//样条插值线性方程组解向量
v o i d h e a d_e n d(f l o a t x1,f l o a t x2,f l o a t*h e a d,f l o a t*e n d)//确定本段取值的起末参数
{
i f(x1-f l o o r(x1)>0.5)*h e a d=c e i l(x1);
e l s e
i f(x1==f l o o r(x1))
*h e a d=x1;
e l s e
*h e a d=f l o o r(x1)+0.5;
i f(x2-f l o o r(x2)>0.5)*e n d=c e i l(x2)-0.5;
e l s e
i f(x2==f l o o r(x2))
*e n d=x2-0.5;
e l s e
*e n d=f l o o r(x2);
}
2.2 分段线性插值
i n t l i n e a r(f l o a t*x,f l o a t*y,f l o a t*r,i n t s)//分段线性插值（区间内插值操作）
{
f l o a t h e a d,e n d;
h e a d_e n d(x[0],x[1],&h e a d,&e n d);//区间规范化
f l o a t k=(y[1]-y[0])/(x[1]-x[0]);
f l o a t t=h e a d;
i n t i=0;
w h i l e(t<=e n d){
r[i]=y[0]+k*(t-x[0]);
t+=0.5;i++;
}
r e t u r n s+i;
}
v o i d l i n e(f l o a t*x,f l o a t*y,f l o a t*r)//分段线性插值（总函数）
{
f i l e=f o p e n("l i n e a r.t x t","w");
i n t i;i n t s=0;//s记录R[]的当前应写位置
f o r(i=0;i<N-1;i++)
s=l i n e a r(&x[i],&y[i],&r[s],s);
f o r(i=0;i<(i n t)(190/d e l t a);i++)//输出数据
f p r i n t f(f i l e,"%f\n",r[i]);
f c l o s e(f i l e);
}
2.3 分段二次多项式插值
i n t s q u a r e(f l o a t*x,f l o a t*y,f l o a t*r,i n t s)//分段二次多项式插值（区间内插值操作）
{
f l o a t h e a d,e n d;
h e a d_e n d(x[0],x[2],&h e a d,&e n d);//区间规范化
f l o a t k1=(y[1]-y[0])/(x[1]-x[0]);
f l o a t k2=(y[2]-y[1])/(x[2]-x[1]);
k2=(k2-k1)/(x[2]-x[0]);
f l o a t t=h e a d;
i n t i=0;
f l o a t t e m p;
w h i l e(t<=e n d)
{
r[i]=y[0]+k1*(t-x[0])+k2*(t-x[0])*(t-x[1]);
t+=0.5;i++;
}
r e t u r n s+i;
}
v o i d s q r(f l o a t*x,f l o a t*y,f l o a t*r)//分段二次多项式插值（总函数）
{
f i l e=f o p e n("s q r.t x t","w");
i n t i;i n t s=0;//s记录R[]的当前应写位置
f o r(i=0;i<N-2;i++)
s=s q u a r e(&x[i*2],&y[i*2],&r[s],s);
f o r(i=0;i<(i n t)(190/d e l t a);i++)
f p r i n t f(f i l e,"%f\n",r[i]);//输出数据
f c l o s e(f i l e);
}
2.4 分段三次多项式插值
i n t c u b i c(f l o a t*x,f l o a t*y,f l o a t*r,i n t s)//分段三次多项式插值（区间内插值操作）
{
f l o a t h e a d,e n d;
h e a d_e n d(x[0],x[3],&h e a d,&e n d);
f l o a t k1=(y[1]-y[0])/(x[1]-x[0]);
f l o a t k2=(y[2]-y[1])/(x[2]-x[1]);
f l o a t k3=(y[3]-y[2])/(x[3]-x[2]);
f l o a t k4=(k2-k1)/(x[2]-x[0]);
k3=(k3-k2)/(x[3]-x[1]);
k3=(k3-k4)/(x[3]-x[0]);
f l o a t t=h e a d;
i n t i=0;
w h i l e(t<=e n d)
{
r[i]=y[0]+k1*(t-x[0])+k4*(t-x[0])*(t-x[1])+k3*(t-x[0])*(t-x[1])*(t-x[2]);
t+=0.5;i++;
}
r e t u r n s+i;
}
v o i d c u b(f l o a t*x,f l o a t*y,f l o a t*r)//分段三次多项式插值（总函数）
{
f i l e=f o p e n("c u b.t x t","w");
i n t i;i n t s=0;
f o r(i=0;i<N-3;i++)
s=c u b i c(&x[i*3],&y[i*3],&r[s],s);
f o r(i=0;i<(i n t)(190/d e l t a);i++)
f p r i n t f(f i l e,"%f\n",r[i]);
f c l o s e(f i l e);
}
2.5 三次样条插值
v o i d i n i t M(v o i d)//构造出线性方程组
{
i n t i,j;
M[0][0]=2;M[N-1][N-1]=2;
f o r(i=1;i<N-1;i++)
{
M[i][i]=2;
M[i][i+1]=(X[i+1]-X[i])/(X[i+1]-X[i-1]);
M[i][i-1]=1-M[i][i+1];
B[i]=6/(X[i+1]-X[i-1])*((Y[i+1]-Y[i])/(X[i+1]-X[i])
-(Y[i]-Y[i-1])/(X[i]-X[i-1]));
}
}
v o i d G e t M()//解得各M值
{
i n t i,j;
f l o a t k;
f o r(i=0;i<N-1;i++){//消去
k=M[i+1][i]/M[i][i];
M[i+1][i]=0;
M[i+1][i+1]-=k*M[i][i+1];
B[i+1]-=k*B[i];
}
f o r(i=N-1;i>1;i--){//回代
k=M[i-1][i]/M[i][i];
M[i-1][i]=0;
B[i-1]-=k*B[i];
}
f o r(i=0;i<N;i++){//求解
B[i]=B[i]/M[i][i];
M[i][i]=1;
}
}
i n t S p i n e(f l o a t*x,f l o a t*y,f l o a t*r,i n t x_i,i n t s)//三次样条插值（区间内插值操作）{
f l o a t h e a d,e n d;
h e a d_e n d(x[0],x[1],&h e a d,&e n d);
f l o a t h=x[1]-x[0];
f l o a t k1=B[x_i]/6/h;
f l o a t k2=B[x_i+1]/6/h;
f l o a t k3=y[0]/h-B[x_i]/6*h;
f l o a t k4=y[1]/h-B[x_i+1]/6*h;
f l o a t t=h e a d;
i n t i=0;
w h i l e(t<=e n d)
{
r[i]=k1*(x[1]-t)*(x[1]-t)*(x[1]-t)+
k2*(t-x[0])*(t-x[0])*(t-x[0])+
k3*(x[1]-t)+
k4*(t-x[0]);
t+=0.5;i++;
}
r e t u r n s+i;
}
v o i d S p i(f l o a t*x,f l o a t*y,f l o a t*r)//三次样条插值（总函数）
{
i n i t M();
G e t M();
f i l e=f o p e n("s p i.t x t","w");
i n t i;i n t s=0;
f o r(i=0;i<N;i++)
s=S p i n e(&x[i],&y[i],&r[s],i,s);
f o r(i=0;i<(i n t)(190/d e l t a);i++)
f p r i n t f(f i l e,"%f\n",r[i]);
f c l o s e(f i l e);
}
2.6 OpenGL环境下显示曲线
#i n c l u d e"s t d a f x.h"
#i n c l u d e<s t d l i b.h>
#i n c l u d e"g l/g l u t.h"
#i n c l u d e"g l/g l.h"
#p r a g m a c o m m e n t(l i b,"o p e n g l32.l i b")//包含O p e n G L库#p r a g m a c o m m e n t(l i b,"g l u32.l i b")
#p r a g m a c o m m e n t(l i b,"g l a u x.l i b")
#p r a g m a c o m m e n t(l i b,"g l u t32.l i b")
v o i d i n i t()
{
g l C l e a r C o l o r(1.0,1.0,1.0,0.0);
g l S h a d e M o d e l(G L_S M O O T H);
}
v o i d d i s p l a y()//显示函数
{
g l C l e a r(G L_C O L O R_B U F F E R_B I T|G L_D E P T H_B U F F E R_B I T);
g l C o l o r3f(0.0,0.0,0.0);
i n t i;
g l L i n e W i d t h(2.0);
g l B e g i n(G L_L I N E_S T R I P);//绘制曲线
f o r(i=0;i<(i n t)(190/d e l t a);i++)
g l V e r t e x2f(i*0.5,R[i]);
g l E n d();
g l P o i n t S i z e(8.0);
g l C o l o r3f(1.0,0.0,0.0);
g l B e g i n(G L_P O I N T S);//绘制型值点
f o r(i=0;i<N;i++)
g l V e r t e x3f(X[i],Y[i],0);
g l E n d();
g l u t S w a p B u f f e r s();
}
v o i d r e s h a p e(i n t w,i n t h)//窗口刷新程序
{
g l V i e w p o r t(0,0,(G L s i z e i)w,(G L s i z e i)h);
g l M a t r i x M o d e(G L_P R O J E C T I O N);
g l L o a d I d e n t i t y();
g l O r t h o(-10.0,200.0,-((f l o a t)h/w*105),(f l o a t)h/w*105,-20.0,20.0);
//定义显示坐标域
g l M a t r i x M o d e(G L_M O D E L V I E W);
g l L o a d I d e n t i t y();
}
i n t_t m a i n(i n t a r g c,c h a r**a r g v)
{
g l u t I n i t(&a r g c,a r g v);
g l u t I n i t D i s p l a y M o d e(G L U T_D E P T H|G L U T_D O U B L E|G L U T_R G B);
g l u t I n i t W i n d o w S i z e(800,500);//窗口初始尺寸
g l u t I n i t W i n d o w P o s i t i o n(150,100);//窗口初始位置
g l u t C r e a t e W i n d o w(a r g v[0]);
i n i t();
l i n e(X,Y,R);//分段线性插值
//s q r(X,Y,R);//分段二次多项式插值
//c u b(X,Y,R);//分段三次多项式插值
//S p i(X,Y,R);//三次样条插值
g l u t D i s p l a y F u n c(d i s p l a y);//注册显示程序
g l u t R e s h a p e F u n c(r e s h a p e);//注册窗口刷新程序
g l u t M a i n L o o p();//开始运行绘制
r e t u r n0;
}
3各种插值方法细化后数据表3.1 分段线性插值法
3.2 分段二次多项式插值
3.3 分段三次多项式插值
3.4 三次样条插值
4 插值机翼轮廓线演示4.1 分段线性插值
整体曲线：
曲线局部：
4.2 分段二次多项式插值
整体曲线：
曲线局部：
4.3 分段三次多项式插值
整体曲线：
曲线局部：
4.4 三次样条插值整体曲线：
曲线局部：
5 不同插值方法的优劣
5.1 分段线性插值
分段线性插值的意义相当于将两个型值点之间用直线连接，曲线在型值点处仅能达到位置连续。

这样的插值方法构造相当简单，但没有实用性。

5.2 分段高次多项式插值
分段高次多项式选择个数等于次数的连续区间进行插值，在这些区间包括内型值点处的曲线都是连续且光滑的。

但缺点是相邻的两个高次多项式插值区间间的型值点处仅能达到位置连续。

虽然可以考虑让各次插值有重叠的区间，然后对重叠区间内的不同插值结果进行融合来改善结果，但这种融合较难实现且效果不佳。

另外，对于每次插值，在插值区段上无法控制得到需要的曲线形状，曲线往往难以满足实际需求。

该种插值方式的优点在于计算非常简便（牛顿插值法）。

5.3 三次样条插值
三次样条插值利用三次多项式插值一个区间内的曲线，可以通过规定型值点处的一阶或二阶导数值来使相邻区间的插值曲线光滑连接，在型值点处达到二阶连续。

该方法的优势：设计时可以为曲线插值提供清晰的几何要求，构造连续性良好的曲线。

缺点在于：计算相对复杂，不具有局部性。