曲线拟合线性最小二乘法及其MATLAB程序.pdf

《曲线拟合线性最小二乘法及其MATLAB程序.pdf》由会员分享，可在线阅读，更多相关《曲线拟合线性最小二乘法及其MATLAB程序.pdf（11页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、1 1曲线拟合的线性最小二乘法及其曲线拟合的线性最小二乘法及其 MATLABMATLAB程序程序例例 7.2.17.2.1给出一组数据点列入表(xi, yi)7 2 中， 试用线性最小二乘法求拟合曲线，并用（7.2） ， （7.3）和（7.4）式估计其误差，作出拟合曲线.表 72例 7.2.1 的一组数据(xi, yi)xiyi-2.5 -1.7 -1.1 -0.8 0 0.1 1.52.7 3.6-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04解解（1 1）在 MATLAB 工作窗口输入程序 x=-2.5 -1.7 -1.1

2、-0.8 0 0.1 1.5 2.7 3.6; y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.5068.04;plot(x,y,r*),legend(实验数据(xi,yi)xlabel(x), ylabel(y),title(例7.2.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图（略）.（3 3）编写下列 MATLAB程序计算在f (x)(xi, yi)处的函数值，即输入程序 syms a1 a2 a3 a4x=-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6;fi=a1.*x.3+ a2.*x.

3、2+ a3.*x+ a4运行后屏幕显示关于a1,a2, a3和a4的线性方程组fi = -125/8*a1+25/4*a2-5/2*a3+a4,-4913/1000*a1+289/100*a2-17/10*a3+a4,-1331/1000*a1+121/100*a2-11/10*a3+a4,-64/125*a1+16/25*a2-4/5*a3+a4,a4, 1/1000*a1+1/100*a2+1/10*a3+a4,27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4,5832/125*a1+324/25*a2+18/5*

4、a3+a4编写构造误差平方和的MATLAB 程序 y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.5068.04;fi=-125/8*a1+25/4*a2-5/2*a3+a4,-4913/1000*a1+289/100*a2-17/10*a3+a4,-1331/1000*a1+121/100*a2-11/10*a3+a4,-64/125*a1+16/25*a2-4/5*a3+a4,a4,1/1000*a1+1/100*a2+1/10*a3+a4,27/8*a1+9/4*a2+3/2*a3+a4,19683/1000*a1+729/100*a

5、2+27/10*a3+a4,5832/125*a1+324/25*a2+18/5*a3+a4;fy=fi-y; fy2=fy.2; J=sum(fy.2)运行后屏幕显示误差平方和如下J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913/1000*a1+289/100*a2-17/10*a3+a4+171/2)2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/1

6、0*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)2为求使达到a1,a2,a3,a4J最小， 只需利用极值的必要条件J 0(k 1,2,3,4)，ak得到关于的a1,a2,a3,a4线性方程组，这可以由下面的MATLAB 程序完成，即输入程序 syms a1 a2 a3 a4J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)2+(-4913

7、/1000*a1+289/100*a2-17/10*a3+a4.+171/2)2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)2+(a4+91/10)2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/

8、25)2;Ja1=diff(J,a1);Ja2=diff(J,a2);Ja3=diff(J,a3);Ja4=diff(J,a4);Ja11=simple(Ja1),Ja21=simple(Ja2),Ja31=simple(Ja3),Ja41=simple(Ja4),运行后屏幕显示 J 分别对 a1, a2 ,a3 ,a4的偏导数如下Ja11=56918107/10000*a1+32097579/25000*a2+1377283/2500*a3+23667/250*a4-8442429/625Ja21 =32097579/25000*a1+1377283/2500*a2+23667/250*a3

9、+67*a4+767319/625Ja31 =1377283/2500*a1+23667/250*a2+67*a3+18/5*a4-232638/125Ja41 =23667/250*a1+67*a2+18/5*a3+18*a4+14859/25解线性方程组 Ja11 =0，Ja21 =0，Ja31 =0，Ja41=0，输入下列程序A=56918107/10000,32097579/25000,1377283/2500,23667/250;32097579/25000,1377283/2500,23667/250,67;1377283/2500, 23667/250, 67, 18/5; 23

10、667/250, 67, 18/5, 18;B=8442429/625, -767319/625, 232638/125, -14859/25;C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f 及其系数 C 如下C = 5.0911 -14.1905 6.4102 -8.2574f=716503695845759/140737488355328*x3-7988544102557579/562949953421312*x2+1804307491277693/281474976710656*x-4648521160813215/562949953421312故所求的拟合曲线为f (

11、x) 5.0911 x314.1905 x2 6.4102 x 8.2574.（4 4）编写下面的MATLAB 程序估计其误差，并作出拟合曲线和数据的图形.输入程序 xi=-2.5 -1.7-1.1-0.800.11.52.73.6;y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.126.5068.04;n=length(xi);f=5.0911.*xi.3-14.1905.*xi.2+6.4102.*xi -8.2574;x=-2.5:0.01: 3.6;F=5.0911.*x.3-14.1905.*x.2+6.4102.*x -8.2574;fy

12、=abs(f-y); fy2=fy.2; Ew=max(fy),E1=sum(fy)/n, E2=sqrt(sum(fy2)/n)plot(xi,y,r*), hold on, plot(x,F,b-), hold offlegend(数据点(xi,yi),拟合曲线y=f(x),xlabel(x), ylabel(y),title(例7.2.1的数据点(xi,yi)和拟合曲线y=f(x)的图形)运行后屏幕显示数据与(xi, yi)拟合函数 f的最大误差Ew， 平均误差 E1 和均方根误差 E2 及其数据点和(xi, yi)拟合曲线 y=f(x)的图形（略）.Ew = E1 = E2 =3.10

13、5 4 0.903 4 1.240 97.37.3函数的选取函数的选取rk(x)及其及其 MATLABMATLAB 程序程序例例7.3.17.3.1给出一组实验数据点的(xi, yi)横坐标向量为x x=(-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2.6,-2.5, -2.1,-1.5, -2.7,-3.6)，纵横坐标向量为y y=(459.26,52.81,198.27,165.60,59.17,41.66,25.92, 22.37,13.47,12.87, 11.87,6.69,14.87,24.22)，试用线性最小二乘法求拟合曲线，并用（ 7.2

14、） ， （7.3）和（7.4）式估计其误差，作出拟合曲线.解解（1 1）在 MATLAB 工作窗口输入程序x=-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2.6,-2.5,-2.1,-1.5, -2.7,-3.6;y=459.26,52.81,198.27,165.60,59.17,41.66,25.92,22.37,13.47, 12.87, 11.87,6.69,14.87,24.22;plot(x,y,r*),legend(实验数据(xi,yi)xlabel(x), ylabel(y),title(例7.3.1的数据点(xi,yi)的散点图)运行后

15、屏幕显示数据的散点图（略）.（3 3）编写下列 MATLAB程序计算在f (x)(xi, yi)处的函数值，即输入程序 syms a bx=-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2.6,-2.5,-2.1,-1.5,-2.7,-3.6; fi=a.*exp(-b.*x)运行后屏幕显示关于a和b的线性方程组fi = a*exp(17/2*b), a*exp(87/10*b), a*exp(71/10*b),a*exp(34/5*b), a*exp(51/10*b), a*exp(9/2*b), a*exp(18/5*b),a*exp(17/5*b),

16、 a*exp(13/5*b), a*exp(5/2*b), a*exp(21/10*b),a*exp(3/2*b), a*exp(27/10*b), a*exp(18/5*b)编写构造误差平方和的MATLAB 程序如下y=459.26,52.81,198.27,165.60,59.17,41.66,25.92,22.37,13.47,12.87, 11.87, 6.69,14.87,24.22;fi=a*exp(17/2*b),a*exp(87/10*b),a*exp(71/10*b),a*exp(34/5*b),a*exp(51/10*b),a*exp(9/2*b),a*exp(18/5*b

17、),a*exp(17/5*b),a*exp(13/5*b),a*exp(5/2*b),a*exp(21/10*b),a*exp(3/2*b),a*exp(27/10*b), a*exp(18/5*b);fy=fi-y;fy2=fy.2;J=sum(fy.2)运行后屏幕显示误差平方和如下J =(a*exp(17/2*b)-22963/50)2+(a*exp(87/10*b)-5281/100)2+(a*exp(71/10*b)-19827/100)2+(a*exp(34/5*b)-828/5)2+(a*exp(51/10*b)-5917/100)2+(a*exp(9/2*b)-2083/50)2

18、+(a*exp(18/5*b)-648/25)2+(a*exp(17/5*b)-2237/100)2+(a*exp(13/5*b)-1347/100)2+(a*exp(5/2*b)-1287/100)2+(a*exp(21/10*b)-1187/100)2+(a*exp(3/2*b)-669/100)2+(a*exp(27/10*b)-1487/100)2+(a*exp(18/5*b)-1211/50)2为求使达到a,bJ最小，只需利用极值的必要条件，得到关于的a,b线性方程组，这可以由下面的 MATLAB 程序完成，即输入程序 syms a bJ=(a*exp(17/2*b)-22963/5

19、0)2+(a*exp(87/10*b)-5281/100)2+(a*exp(71/10*b)-19827/100)2+(a*exp(34/5*b)-828/5)2+(a*exp(51/10*b)-5917/100)2+(a*exp(9/2*b)-2083/50)2+(a*exp(18/5*b)-648/25)2+(a*exp(17/5*b)-2237/100)2+(a*exp(13/5*b)-1347/100)2+(a*exp(5/2*b)-1287/100)2+(a*exp(21/10*b)-1187/100)2+(a*exp(3/2*b)-669/100)2+(a*exp(27/10*b)

20、-1487/100)2+(a*exp(18/5*b)-1211/50)2;Ja=diff(J,a); Jb=diff(J,b);Ja1=simple(Ja), Jb1=simple(Jb),运行后屏幕显示 J 分别对的偏导数a,b如下Ja1 =2*a*exp(3*b)+2*a*exp(17*b)+2*a*exp(87/5*b)+2*exp(68/5*b)*a+2*exp(9*b)*a+2*a*exp(34/5*b)-669/50*exp(3/2*b)-1487/50*exp(27/10*b)-2507/25*exp(18/5*b)-22963/25*exp(17/2*b)-5281/50*ex

21、p(87/10*b)-19827/50*exp(71/10*b)-2237/50*exp(17/5*b)-1656/5*exp(34/5*b)-1347/50*exp(13/5*b)-5917/50*exp(51/10*b)-1287/50*exp(5/2*b)-2083/25*exp(9/2*b)-1187/50*exp(21/10*b)+4*a*exp(36/5*b)+2*a*exp(26/5*b)+2*a*exp(71/5*b)+2*a*exp(51/5*b)+2*a*exp(5*b)+2*a*exp(21/5*b)+2*a*exp(27/5*b)Jb1 =1/500*a*(2100*a

22、*exp(21/10*b)2+8500*a*exp(17/2*b)2+6800*a*exp(34/5*b)2-10035*exp(3/2*b)-40149*exp(27/10*b)-180504*exp(18/5*b)-3903710*exp(17/2*b)-459447*exp(87/10*b)-1407717*exp(71/10*b)-76058*exp(17/5*b)-1126080*exp(34/5*b)-35022*exp(13/5*b)-301767*exp(51/10*b)-32175*exp(5/2*b)-187470*exp(9/2*b)-24927*exp(21/10*b)

23、+7100*a*exp(71/10*b)2+5100*a*exp(51/10*b)2+4500*a*exp(9/2*b)2+7200*a*exp(18/5*b)2+3400*a*exp(17/5*b)2+2600*a*exp(13/5*b)2+2500*a*exp(5/2*b)2+1500*a*exp(3/2*b)2+2700*a*exp(27/10*b)2+8700*a*exp(87/10*b)2)用解二元非线性方程组的牛顿法的MATLAB 程序求解线性方程组Ja1=0，Jb1 =0，得a = b=2.811 0 0.581 6故所求的拟合曲线（7.13）为f (x) 2.811 0e0.5

24、816 x.（7.14）（4 4）根据（7.2） ， （7.3） ， （7.4）和 （7.14）式编写下面的MATLAB 程序估计其误差，并做出拟合曲线和数据的图形.输入程序 xi=-8.5 -8.7 -7.1 -6.8 -5.10 -4.5 -3.6 -3.4 -2.6 -2.5-2.1 -1.5 -2.7 -3.6;y=459.26 52.81 198.27 165.60 59.17 41.66 25.92 22.3713.47 12.87 11.87 6.69 14.87 24.22;n=length(xi); f=2.8110.*exp(-0.5816.*xi); x=-9:0.01:

25、 -1;F=2.8110.*exp(-0.5816.*x); fy=abs(f-y); fy2=fy.2;Ew=max(fy),E1=sum(fy)/n, E2=sqrt(sum(fy2)/n), plot(xi,y,r*), hold onplot(x,F,b-), hold off,legend(数据点(xi,yi),拟合曲线y=f(x)xlabel(x), ylabel(y), title(例7.3.1的数据点(xi,yi)和拟合曲线y=f(x)的图形)运行后屏幕显示数据与(xi, yi)拟合函数 f的最大误差Ew = 390.141 5，平均误差 E1=36.942 2 和均方根误差

26、E2=106.031 7 及其数据点和拟合曲(xi, yi)线 y=f(x)的图形（略）.7.47.4多项式拟合及其多项式拟合及其 MATLABMATLAB 程序程序例例 7.4.17.4.1给出一组数据点列入表(xi, yi)7 3 中，试用线性最小二乘法求拟合曲线，并用（7.2） ， （7.3）和（7.4）式估计其误差，作出拟合曲线.表 73例 7.4.1 的一组数据(xi, yi)xiyi-2.9 -1.9 -1.1 -0.8 0 0.1 1.52.7 3.653.94 33.68 20.88 16.92 8.79 8.98 4.17 9.12 19.88解解（1 1）首先根据表）首先根

27、据表7 7 3 3 给出的数据点给出的数据点(xi, yi)，用下列，用下列 MATLABMATLAB 程序画出散点图程序画出散点图.在 MATLAB 工作窗口输入程序 x=-2.9 -1.9 -1.1 -0.8 0 0.1 1.5 2.7 3.6;y=53.94 33.68 20.88 16.92 8.79 8.98 4.17 9.1219.88;plot(x,y,r*), legend(数据点(xi,yi)xlabel(x), ylabel(y),title(例7.4.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图（略）.（3 3）用用作作线线性性最最 小小二乘二乘 拟拟合合

28、 的的多多项项式式拟拟 合合的的 MATMAT LABLAB 程程序序 求求待待定系定系 数数ak(k 1, 2,3).输入程序 a=polyfit(x,y,2)运行后输出（7.16）式的系数a =2.8302 -7.3721 9.1382故拟合多项式为f (x) 2.830 2x2 7.372 1x 9.138 2.（4 4）编写下面的MATLAB 程序估计其误差，并做出拟合曲线和数据的图形.输入程序 xi=-2.9 -1.9 -1.1 -0.8 0 0.1 1.5 2.7 3.6;y=53.94 33.68 20.88 16.92 8.79 8.98 4.17 9.12 19.88;n=l

29、ength(xi); f=2.8302.*xi.2-7.3721.*xi+9.1382x=-2.9:0.001:3.6;F=2.8302.*x.2-7.3721.*x+8.79;fy=abs(f-y); fy2=fy.2; Ew=max(fy), E1=sum(fy)/n,E2=sqrt(sum(fy2)/n), plot(xi,y,r*, x,F,b-),legend(数据点(xi,yi),拟合曲线y=f(x)xlabel(x), ylabel(y), title(例7.4.1 的数据点(xi,yi)和拟合曲线y=f(x)的图形)运行后屏幕显示数据与(xi, yi)拟合函数 f的最大误差Ew

30、，平均误差 E1 和均方根误差 E2及其数据点(xi,yi)和拟合曲线y=f(x)的图形（略）.Ew = E1 = E2 =0.745 7, 0.389 2, 0.436 37.57.5拟合曲线的线性变换及其拟合曲线的线性变换及其 MATLABMATLAB 程序程序例例 7.5.17.5.1给出一组实验数据点的(xi, yi)横坐标向量为 x x=（7.5 6.8 5.10 4.53.6 3.4 2.6 2.5 2.1 1.5 2.7 3.6） ， 纵横坐标向量为 y y= （359.26 165.60 59.1741.66 25.92 22.37 13.47 12.87 11.87 6.69

31、 14.87 24.22），试用线性变换和线性最小二乘法求拟合曲线，并用（7.2） ， （7.3）和（7.4）式估计其误差，作出拟合曲线.解解（1 1）首先根据给出的数据点(xi, yi)，用下列 MATLAB 程序画出散点图.在 MATLAB 工作窗口输入程序x=7.5 6.8 5.10 4.5 3.6 3.4 2.6 2.5 2.1 1.5 2.73.6;y=359.26 165.60 59.17 41.66 25.92 22.37 13.47 12.8711.87 6.69 14.87 24.22;plot(x,y,r*), legend(数据点(xi,yi)xlabel(x), yla

32、bel(y),title(例7.5.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图（略）.（2 2）根据数据散点图，取拟合曲线为y aebx(a 0,b 0),（7.19）其中是待定a,b系数.令Y ln y, A lna, B b，则（7.19）化为Y A Bx.在 MATLAB 工作窗口输入程序x=7.5 6.8 5.10 4.5 3.6 3.4 2.6 2.5 2.1 1.5 2.73.6;y=359.26165.6059.1741.6625.9222.3713.4712.87 11.87 6.69 14.87 24.22;Y=log(y); a=polyfit(x,Y,1

33、); B=a(1);A=a(2); b=B,a=exp(A)n=length(x); X=8:-0.01:1; Y=a*exp(b.*X); f=a*exp(b.*x);plot(x,y,r*,X,Y,b-), xlabel(x),ylabel(y)legend(数据点(xi,yi),拟合曲线y=f(x)title(例 7.5.1 的数据点(xi,yi)和拟合曲线y=f(x)的图形)fy=abs(f-y); fy2=fy.2; Ew=max(fy), E1=sum(fy)/n,E2=sqrt(sum(fy2)/n)运行后屏幕显示 e 的系y abx数b =0.624 1，a=2.703 9,数

34、据与拟合(xi, yi)函数0.6241xf 的最大误差 Ew=67.641 9，平均误差 E1=8.677 6 和均方根误差E2=20.711 3 及其数据点和拟合曲(xi, yi)线 e 的图形f (x) 2.703 9（略）.7.67.6函数逼近及其函数逼近及其 MATLABMATLAB 程序程序最佳均方逼近的最佳均方逼近的 MATLABMATLAB 主程序主程序function yy1,a,WE=zjjfbj(f,X,Y,xx)m=size(f);n=length(X);m=m(1);b=zeros(m,m); c=zeros(m,1);if n=length(Y) error( X和

35、Y的维数应该相同)endfor j=1:mfor k=1:m b(j,k)=0;for i=1:nb(j,k)=b(j,k)+feval(f(j,:),X(i)*feval(f(k,:),X(i);end endc(j)=0;for i=1:n c(j)=c(j)+feval(f(j,:),X(i)*Y(i);endenda=bc;WE=0;for i=1:n ff=0;for j=1:mff=ff+a(j)*feval(f(j,:),X(i);endWE=WE+(Y(i)-ff)*(Y(i)-ff);endif nargin=3return;endyy=;for i=1:m l=;for j

36、=1:length(xx) l=l,feval(f(i,:),xx(j);endyy=yy l;endyy=yy*a; yy1=yy; a=a;WE;2例例 7.6.17.6.1对数据 X 和Y Y, 用函数进行y 1, y x, y x逼近，用所得到的逼近函数计算在处的函x 6.5数值，并估计误差.其中X X=(13456789); Y Y=(-11-13-11-7-171729).解解在MATLAB工作窗口输入程序 X= 1 3 4 5 6 7 8 9;Y=-11 -13 -11 -7 -1 7 1729;f=fun0;fun1;fun2; yy,a,WE=zjjfbj(f,X,Y,6.5

37、)运行后屏幕显示如下yy = 2.75000000000003a = -7.00000000000010 -4.99999999999995 1.00000000000000WE = 7.172323350269439e-027例例7.6.27.6.2对数据X和Y Y， 用函数y 1, y x, y x，y sin xy cosx，y e ，进行逼近，其中 X X=（0 0.50 1.00 1.50 2.00 2.50 3.00 ） ，Y Y=（0 0.4794 0.8415 0.98150.9126 0.5985 0.1645）.解解 在MATLAB工作窗口输入程序 X= 0 0.50 1.

38、00 1.50 2.00 2.50 3.00;Y=0 0.4794 0.8415 0.9815 0.9126 0.5985 0.1645;f=fun0;fun1;fun2;fun3;fun4;fun5;xx=0:0.2:3;yy,a,WE=zjjfbj(f,X,Y, xx), plot(X,Y,ro,xx,yy,b-)运行后屏幕显示如下（图略）yy = Columns 1 through 7-0.0005 0.2037 0.3939 0.5656 0.7141 0.83480.9236 Columns 8 through 140.9771 0.9926 0.9691 0.9069 0.8080

39、 0.67660.5191 Columns 15 through 16 0.3444 0.1642a = 0.3828 0.4070 -0.3901 0.0765 -0.4598 0.5653WE = 1.5769e-004即，最佳逼近函数为y=0.3828+0.4070*x-0.3901*x2+0.0765*exp(x) -0.4598*cos(x) +0.5653*sin(x).2x7.77.7三角多项式逼近及其三角多项式逼近及其 MATLABMATLAB程序程序计算三角多项式的计算三角多项式的 MATLABMATLAB 主程序主程序function A,B,Y1,Rm=sanjiao(X

40、,Y,X1,m)n= length(X)-1;max1=fix(n-1)/2);if m max1 m=max1;endA=zeros(1,m+1);B=zeros(1,m+1);Ym=(Y(1)+Y(n+1)/2; Y(1)=Ym; Y(n+1)=Ym; A(1)=2*sum(Y)/n;for i=1:mB(i+1)=sin(i*X)*Y; A(i+1)=cos(i*X)*Y;endA=2*A/n; B=2*B/n; A(1)=A(1)/2;Y1=A(1);for k=1:mY1=Y1+A(k+1)*cos(k*X1)+ B(k+1)*sin(k*X1);Tm=A(1)+A(k+1).*co

41、s(k*X)+ B(k+1).*sin(k*X); k=k+1;endY;Tm; Rm=(sum(Y-Tm).2)/n;例例 7.7.17.7.1根据上的个, n 13, 60, 350等距横坐标点xi 2 i n(i 0,1,2,n)和函数f(x)2sinx.3（1）求的 6 阶三f (x)角多项式逼近，计算均方误差;（2）将这三个三角多项式分别与的傅里f (x)叶级数18 3nn1f(x)(1)sinnx2n19n 1的前 6 项进行比较;（3）利用三角多项式分别计算Xi= -2, 2.5 的值;（4）在同一坐标系中，画出函数f (x)，n 13, 60, 350的三角多项式和数据点的图形

42、.解解 （1 1）输入程序）输入程序X1=-pi:2*pi/13:pi;Y1=2*sin(X1/3);X1i=-2,2.5;A1,B1,Y11,Rm1=sanjiao(X1,Y1,X1i,6), X2=-pi:2*pi/60:pi;Y2=2*sin(X2/3);A2,B2,Y12,Rm2=sanjiao(X2,Y2,X1i,6)X3=-pi:2*pi/350:pi;Y3=2*sin(X3/3);A3,B3,Y13,Rm3=sanjiao(X3,Y3,X1i,6)X1i=-2,2.5;Y1=2*sin(X1i/3)for n=1:6bi=(-1)(n+1)*18*sqrt(3)*n/(pi*(9

43、*n2-1)end（2 2）画图，输入程序）画图，输入程序X1=-pi:2*pi/13:pi;Y1=2*sin(X1/3);Xi=-pi:0.001:pi; f=2*sin(Xi/3);A1,B1,Y1i,R1m=sanjiao(X1,Y1,Xi,6);X2=-pi:2*pi/60:pi;Y2=2*sin(X2/3); X3=-pi:2*pi/350:pi;Y3=2*sin(X3/3);A2,B2,Y2i,R2m=sanjiao(X2,Y2,Xi,6);A3,B3,Y3i,R3m=sanjiao(X3,Y3,Xi,6);plot(X1,Y1,r*, Xi, Y1i,b-,Xi, Y2i,g-,

44、 Xi, Y3i, m:,Xi, f, k-.)xlabel(x),ylabel(y)legend(数据点(xi,yi),n=13的三角多项式,n=60的三角多项式,n=350的三角多项式,函数f(x)title(例 7.7.1 的数据点(xi,yi)、n=13，60，350 的三角多项式 T3和函数 f(x)的图形)运行后图形（略）.7.87.8随机数据点上的二元拟合及其随机数据点上的二元拟合及其MATLABMATLAB程序程序例例7.8.17.8.1 设节点(X,Y,Z)中的X和Y分别是在区间和上的53, 3 2.5, 3.50个随机数，Z是函数Z=7-3x3e 在(X,Y)的值，拟合点(

45、XI,YI)中的 XI=-3:0.2:3,YI=-2.5:0.2:3.5.分别用二元拟合方法中最近邻内插法、三角基线性内插法、三角基三次内插法和 MATLAB 4 网格化坐标方法计算在(XI,YI)处的值，作出它们的图形，并与被拟和曲面进行比较.解解（1 1）最近邻内插法）最近邻内插法.输入程序 x=rand(50,1);y=rand(50,1); %生成50个一元均匀分布随机数x和y， x，y .X=-3+(3-(-3)*x;%利用x生成的随机变量.Y=-2.5+(3.5-(-2.5)*y; %利用y生成的随机变量.Z=7-3* X.3 .* exp(-X.2 - Y.2);%在每个随机点（

46、X,Y）处计算Z-x2-y2的值.X1=-3:0.2:3;Y1=-2.5:0.2:3.5;XI,YI = meshgrid(X1,Y1);%将坐标（XI,YI）网格化.ZI=griddata(X,Y,Z,XI,YI, nearest) %计算在每个插值点（XI,YI）处的插值ZI.mesh(XI,YI, ZI)%作二元拟合图形.xlabel(x), ylabel(y), zlabel(z),title(用最近邻内插法拟合函数z =7-3 x3 exp(-x2 - y2)的曲面和节点的图形)%legend(拟合曲面,节点(xi,yi,zi)hold on%在当前图形上添加新图形.plot3(X,

47、Y,Z, bo)%用兰色小圆圈画出每个节点(X,Y,Z).hold of%结束在当前图形上添加新图形.运行后屏幕显示用最近邻内插法拟合函数Z=7-3x3e在两组不同节点处的曲面及其插值ZI（略）.（2 2）三角基线性内插法）三角基线性内插法.输入程序 x=rand(50,1);y=rand(50,1); %生成50个一元均匀分布随机数x和y， x，y .X=-3+(3-(-3)*x;%利用x生成 上的随机变量.Y=-2.5+(3.5-(-2.5)*y; %利用y生成 上的随机变量.Z=7-3* X.3 .* exp(-X.2 - Y.2);%在每个随机点（X,Y）处计算Z的值.X1=-3:0.

48、2:3;Y1=-2.5:0.2:3.5;XI,YI = meshgrid(X1,Y1);%将坐标（XI,YI）网格化.ZI=griddata(X,Y,Z,XI,YI, linear) %计算在每个插值点（XI,YI）处的插值ZI.mesh(XI,YI, ZI)%作二元拟合图形.xlabel(x), ylabel(y), zlabel(z),title(用三角基线性内插法拟合函数z =7-3 x3 exp(-x2 - y2)的曲面和节点的图形)%legend(拟合曲面,节点(xi,yi,zi)hold on%在当前图形上添加新图形.plot3(X,Y,Z, bo)%用兰色小圆圈画出每个节点(X,

49、Y,Z).hold of%结束在当前图形上添加新图形.运行后屏幕显示用三角基线性内插法拟合函数Z=7-3x3e在两组不同节点处的曲面和节点的图形及其插值ZI（略）.（3 3）三角基三次内插法）三角基三次内插法.输入程序 x=rand(50,1);y=rand(50,1); %生成50个一元均匀分布随机数x和y， x，y .X=-3+(3-(-3)*x;%利用x生成 上的随机变量.Y=-2.5+(3.5-(-2.5)*y;%利用y生成 上的随机变量.Z=7-3* X.3 .* exp(-X.2 - Y.2);%在每个随机点（X,Y）处计算Z的值.X1=-3:0.2:3;Y1=-2.5:0.2:3

50、.5;XI,YI = meshgrid(X1,Y1);%将坐标（XI,YI）网格化.ZI=griddata(X,Y,Z,XI,YI, cubic) %计算在每个插值点（XI,YI）处的插值ZI.-x2-y2-x2-y2mesh(XI,YI, ZI)%作二元拟合图形.xlabel(x), ylabel(y), zlabel(z),title(用三角基三次内插法拟合函数z =7-3 x3 exp(-x2 - y2)的曲面和节点的图形)%legend(拟合曲面,节点(xi,yi,zi)hold on%在当前图形上添加新图形.plot3(X,Y,Z, bo)%用兰色小圆圈画出每个节点(X,Y,Z).h