编著：李新平
二零零八年三月十四日
湖南理工学院 MATLAB 实验指导书
实验六、数据插值和数据拟合
6.1 实验目的
1）掌握用 MATLAB 计算拉格朗日、分段线性、三次样条三种插值的方法，改变节点 的数目，对三种插值结果进行初步分析。
2）掌握用 MATLAB 进行多项式最小二乘拟合，会选择合适的函数及转化为线性函数。 3）通过实例学习用数据插值和数据拟合解决实际问题。
6.2 分段线性插值
设给定一元未知函数 y = f (x ) 的 n + 1 个结点的数据 a = x 0 < L < x n = b 对应的函数 值 y L , y 根据这些结点数据求其余 x j ( j ¹ i ) 点的函数值 y j ，可将相邻两个节点之间用 0 , n ，直线连接起来，如此形成的一条折线(见右图)构成的分段线性函数 I ) 来近似表示未知函 n (x 数 f (x ) ，从而解决该插值问题的方法就称为分段线性插值。可用如下公式表示：
I n ( x ) = å y j l j ( x ) » f ( x )
j = 0
n
ì
ï ï l j ( x ) = í
ï ï î
x - x j - 1
, x £ x £ x j
x j - x j - 1 j - 1 x - x j + 1
, x j £ x £ x j + 1
x j - x j + 1
0 , 其余
可用 MATLAB 命令 y=interp1(x0,y0,x)来实现， 其中参数 x0 为给定结点数据的横坐标向 量，参数 y0 为 x0 对应的函数值，参数 x 为要未知结点的横坐标向量，函数返回值 y为参数 x 根据分段线性插值得到的函数值。
【例】插值求在[0,15]区间内步长为 0.1 的机床加工数据：
>>x0=[0 3 5 7 9 11 12 13 14 15]; y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
>>x=0:0.1:15; % 插值点
>>y=inpert1(x0, y0, x) % 插值求得函数值
6.3 拉格朗日插值
设未知函数 y = g (x ) 是n次多项式，给定该n次多项式 n + 1 个结点的数据 {( x ， i , y i ),
i = 0 , L , n } 根据这些结点数据求其余 x j ( j ¹ i ) 点的函数值 y j ，可考虑如下构造：
·2·
湖南理工学院
L x ) = å y l x ) » g ( x ) n ( i i (
i = 0 n
MATLAB 实验指导书
l x ) = i (
( x - x L ( x - x x - x L ( x - x 0 ) i - 1 )( i + 1 ) n )
( x L ( x x L ( x i - x 0 ) i - x i - 1 )( i - x i + 1 ) i - x n )
来近似表示n次多项式 g (x ) ，从而解决该插值问题的方法就称为拉格朗日插值。
用 MATLAB 进行拉格朗日插值，要编写拉格朗日插值函数 m 文件，如 lagr1.m，其函 数调用形式为 y=lagr1(x0,y0,x)，其参数形式和函数返回值与分段线性插值相同。
【例】插值求在[0,15]区间内步长为 0.1 的机床加工数据：
>>x0=[0 3 5 7 9 11 12 13 14 15]; y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
>>x=0:0.1:15; % 插值点
>>y=lagr1(x0, y0, x) % 插值求得函数值
6.4 三次样条插值
设给定一元未知函数 y = h (x ) 的 n + 1 个结点的数据 {( x i = 0 , L , n } ，根据这些结 i , y i ), 点数据求其余 x j ( j ¹ i ) 点的函数值 y j ，可考虑做如下近似构造 S ( x ) » h ( x ) ：
(1) 在每个小区间 [ x 上都是 3 次多项式, 即 s x ) = a i - 1 , x i ] i ( i x + b i x + c i x + d i ; (2) S ( x L , n ; i ) = y i , i = 0,
(3) S (x ) 在 x x x ( x ( x 0 £ x £ x n 上的二阶导数连续, 即 s i ( i ) = s i + 1 ( i ) , s i ' i ) = s i + 1 ' i ) ,
3
2
s ( x x ; i \" i ) = s i + 1\" ( i )
在 MATLAB 中可以采用函数 y=interp1(x0,y0,x,’spline’)或 y=spline(x0,y0,x)来实现， 参数 形式和函数返回值跟分段线性插值和拉格朗日插值相同。
【例】插值求在[0,15]区间内步长为 0.1 的机床加工数据：
>>x0=[0 3 5 7 9 11 12 13 14 15]; y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6]; >>x=0:0.1:15; % 插值点
>>y=interp1(x0, y0, x,’spline’) % 或 y=spline(x0, y0, x)
6.5 多项式最小二乘数拟合
2
为了估计指定函数(如 y = a x + a 0 + a 1 x , y = a 0 + a 1 2 x , y = a 0 +
a a 1 1 x
, y = a e 等) 0 x
函数的参数， 从而确定该函数， 即作参数估计求得该函数， 可用先将以上函数转化了以 a 0 ,a 1 为自变量的线性函数，再将给定的 n + 1 个结点数据 {( x i = 0 , L , n } 代入线性函数，求i , y i ),
·3·
湖南理工学院 MATLAB 实验指导书
T
解非齐次线性方程组 XA= Y ，使距离 d a a i 的平方和最小的解 A = ( 0 , 1 ) 即为所求参数。该 方法就称为线性最小二乘拟合。
在 MATLAB 中可用函数 aa=ployfit(x,y,n)来求得参数a， 从而实现线性最小二乘拟合(当 然也可以实现多项式函数拟合)，其中参数 x 为给定结点数据的横坐标向量，y 为对应的纵 坐标向量，n 为多项式的次数(如线性拟合则 n 为 1，二次多项式拟合为 2 等)，函数返回值 aa 为拟合的多项式系数。
在求得多项式系数后，为了求得多项式的值，可用 MATLAB 函数 y=polyval(aa,x)求得 系数为 aa 的多项式在指定点 x的函数值 y。
【例】用多项式最小二乘拟合求电阻R 与温度t之间的关系 R = at + b ： >>t=[20.5 32.5 51 73 95.7]; r=[765 826 873 942 1032]; >>aa=polyfit(t,r,1); % 插值点
>>a=aa(1) % 插值求得一次项系数 a=3.3940 >>b=aa(2) % 插值求得常数项系数 b=702.4918 >>y=polyval(aa,t) >>plot(t,r,’k+’, t,y,’r’)
【上机练习】
1．编制计算拉格朗日插值的函数 m文件 lagr1.m。
2．选择一些函数，在 n 个节点上(n 不要太大，如 5~11)用拉格朗日、分段线性、三次样条 三种插值方法，计算 m个插值点的函数值(m要适中，如 50~100)。通过数值和图形输出， 将三种插值结果与精确值进行比较。适当增加 n，再作比较，由此作初步分析。下列函 数供选择参考：
(1) y = 1 - x , - 1 £ x £ 1 ;
2
(2) y = e
- x 2
, - 2 £ x £ 2
3．用 y = x 在 x =0, 1, 4, 9, 16 产生 5个节点函数值 y L , y 1 , y 2 , 5 。用不同的节点构造插值 公式来计算 x = 5 处的插值（如用 y ，与精确值比较并分析。 1 ~ y 5 ; y 1 ~ y 4 ; y 2 ~ y 5 等）4．下表给出的 x, y 数据位于机翼断面的轮廓线上，y1 和 y2 分别对应轮廓的上下线。假设 需要得到x 坐标每改变 0.1时的 y坐标。试完成加工所需数据，并画出曲线。 x y1 y2 0 0 0 3 1.8 1.2 5 2.2 1.7 7 2.7 2.0 9 3.0 2.1 11 3.1 2.0 12 2.9 1.8 13 2.5 1.2 14 2.0 1.0 15 1.6 1.6 5．对以下数据分别作二次和三次多项式拟合，求得多项式系数，并画出图形。
x y 2 6.4 4 8.4 5 9.28 6 9.5 7 9.7 8 9.86 10 10.2 12 10.4 14 10.5 16 10.6·4·
湖南理工学院 MATLAB 实验指导书
rt
6．根据以下数据，确定1790~1900 年美国人口的指数增长模型 x ( t ) = x e 中的人口年增长 0
率r 和初始人口数 x 0 。 t x 1790 3.9 1800 5.3 1810 7.2 1820 9.6 1830 12.9 1840 17.1 1850 23.2 1860 31.4 1870 38.6 1880 50.2 1890 62.9 1900 76.0 附 lagr1.m文件内容如下, 请仔细读懂各语句含义
·5·
因篇幅问题不能全部显示，请点此查看更多更全内容