Matlab实验指导书(实验六)

MATLAB实验指导书

编著：李新平

二零零八年三月十四日

实验六、数据插值和数据拟合

6.1 实验目的

1）掌握用 MA TLAB 计算拉格朗日、分段线性、三次样条三种插值的方法，改变节点 的数目，对三种插值结果进行初步分析。

2）掌握用 MA TLAB 进行多项式最小二乘拟合，会选择合适的函数及转化为线性函数。 3）通过实例学习用数据插值和数据拟合解决实际问题。

6.2 分段线性插值

设给定一元未知函数 ) (x f y = 的 1 + n 个结点的数据 b x x a n = < < = L 0 对应的函数 值 n y y , , 0 L ，根据这些结点数据求其余 ) ( i j x j 1 点的函数值 j y ，可将相邻两个节点之间用 直线连接起来，如此形成的一条折线(见右图)构成的分段线性函数 ) (x I n 来近似表示未知函 数 ) (x f ，从而解决该插值问题的方法就称为分段线性插值。可用如下公式表示：

)

( ) ( ) ( 0

x f x l y x I n

j j j n ? = ? = 其余

, 0 , , ) ( 1

1

1

1 1 1

+ + + - - - ￡ ￡ - - ￡ ￡ - - ? ? ?

? ? í

ì

= j j j j j j

j j j j j x x x x x x x x x x x x x x x l 可用 MA TLAB 命令 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 拉格朗日插值

设未知函数 ) (x g y = 是n 次多项式，给定该n 次多项式 1 + n 个结点的数据 ),

, {( i i y x ， } , , 0 n i L = 根据这些结点数据求其余 ) ( i j x j 1 点的函数值 j y ，可考虑如下构造：

) ( ) (

) ( 0 x g x l y x L n

i i i n ? = ? = )

( ) )( ( ) ( )

( ) )( ( ) ( ) (

1 1 0 1 1 0 n i i i i i i n i i i x x x x x x x x x x x x x x x x x l - - - - - - - - = + - + - L L L L 来近似表示n 次多项式 ) (x g ，从而解决该插值问题的方法就称为拉格朗日插值。

用 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 三次样条插值

设给定一元未知函数 ) (x h y = 的 1 + n 个结点的数据 } , , 0 ), , {( n i y x i i L = ，根据这些结 点数据求其余 ) ( i j x j 1 点的函数值 j y ，可考虑做如下近似构造 ) ( ) ( x h x S ? ：

(1) 在每个小区间 ] , [ 1 i i x x - 上都是 3 次多项式, 即 i i i i i d x c x b x a x s + + + = 2

3

) ( ? (2) i i y x S = ) ( , n i , ,

0L = ? (3) ) (x S 在 n x x x ￡ ￡ 0 上的二阶导数连续, 即 ) ( ) ( 1 i i i i x s x s + = , ) ( ' ) ( ' 1 i i i i x s x s + = ,

)

( " ) ( " 1 i i i i x s x s + = ? 在 MA TLAB 中可以采用函数 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 多项式最小二乘数拟合

为了估计指定函数(如 x a a y 1 0 + = , 2

2 1 0 x a x a a y + + = , x

a a y 1 0 +

= , x

a e a y 1 0 = 等) 函数的参数， 从而确定该函数， 即作参数估计求得该函数， 可用先将以上函数转化了以 1 0 ,a a 为自变量的线性函数，再将给定的 1 + n 个结点数据 } , , 0 ),

, {( n i y x i i L = 代入线性函数，求

解非齐次线性方程组 Y XA = ，使距离 i d 的平方和最小的解 T

a a A ) , ( 1 0 = 即为所求参数。该 方法就称为线性最小二乘拟合。

在 MA TLAB 中可用函数 aa=ployfit(x,y,n)来求得参数a ， 从而实现线性最小二乘拟合(当 然也可以实现多项式函数拟合)，其中参数 x 为给定结点数据的横坐标向量，y 为对应的纵 坐标向量，n 为多项式的次数(如线性拟合则 n 为 1，二次多项式拟合为 2 等)，函数返回值 aa 为拟合的多项式系数。

在求得多项式系数后，为了求得多项式的值，可用 MA TLAB 函数 y=polyval(aa,x)求得 系数为 aa 的多项式在指定点 x 的函数值 y 。

【例】用多项式最小二乘拟合求电阻R 与温度t 之间的关系 b at R + = ： >>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) 1 1 , 1 2

￡ ￡ - - = x x y ?

(2) 2

2 , 2

￡ ￡ - = - x e

y x 3．用 x y = 在 x =0, 1, 4, 9, 16 产生 5个节点函数值 5 2 1 , , , y y y L 。用不同的节点构造插值 公式来计算 5 = x 处的插值（如用 5 2 4 1 5 1 ~ ? ~ ? ~ y y y y y y 等），与精确值比较并分析。 4．下表给出的 x , y 数据位于机翼断面的轮廓线上，y 1 和 y 2 分别对应轮廓的上下线。假设 需要得到x 坐标每改变 0.1时的 y 坐标。试完成加工所需数据，并画出曲线。 x 0 3 5 7 9 11 12 13 14 15 y 1 0 1.8 2.2 2.7 3.0 3.1 2.9 2.5 2.0 1.6 y 2

1.2

1.7

2.0

2.1

2.0

1.8

1.2

1.0

1.6

5．对以下数据分别作二次和三次多项式拟合，求得多项式系数，并画出图形。

x 2 4 5 6 7 8 10 12 14 16 y

6.4

8.4

9.28

9.5

9.7

9.86

10.2

10.4

10.5

10.6

6．根据以下数据，确定1790~1900 年美国人口的指数增长模型 rt

e x t x 0

) ( = 中的人口年增长 率r 和初始人口数 0 x 。 t 1790 1800 1810 1820 1830 1840 1850 1860 1870 1880 1890 1900 x

3.9

5.3

7.2

9.6

12.9

17.1

23.2

31.4

38.6

50.2

62.9

76.0

附 lagr1.m 文件内容如下, 请仔细读懂各语句含义