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

《曲线拟合的线性最小二乘法及其MATLAB程序》由会员分享，可在线阅读，更多相关《曲线拟合的线性最小二乘法及其MATLAB程序（11页珍藏版）》请在装配图网上搜索。

1、1曲线拟合的线性最小二乘法及其MATLAB程序例7.2.1给出一组数据点（,七）列入表7 2中，试用线性最小二乘法求拟合曲线， 并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.表7-2例7.2. 1的一组数据（x , y ）xi-2.5-1.7-1.1-0.800.11顶.73.6yi-192.9-85.50-36.15-26.52-9.10-8.43 -13.126.5068.04解 (1)在MATLAB工作窗口输入程序 x=-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 -1

2、3.126.5068.04;plot(x,y,r*),legend(实验数据(xi,yi)xlabel(x), ylabel(y),title(例7.2.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图(略).(3)编写下列MATLAB程序计算f (x)在(x,七)处的函数值，即输入程序 syms a1 a2 a3 a4x=-2.5 -1.7-1.1-0.800.11.52.73.6;fi=a1.*x.八3+ a2.*x.八2+ a3.*x+ a4运行后屏幕显示关于a1,a2,。3和。4的线性方程组fi = -125/8*a1+25/4*a2-5/2*a3+a4, -4913/10

3、00*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*a3+a4编写构造误差平方和的MATLAB程序 y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.126.506

4、8.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*a2+27/10*a3+a4,5832/125*a1+324/25*a2+18/5*a3+a4;fy=fi-y; fy2 = fy2; J=sum(fy2)运行

5、后屏幕显示误差平方和如下J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)八2+(-4913/1000*a1+2 89/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*a 2+3/2*a3+a4+328/25)八2+(19683/1000*a1+

6、729/100*a2+27/10*a3+a4-13/ 2)八2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)八2为求a ,a ,a ,a使J达到最小，只需利用极值的必要条件当 =0 (k = 1,2,3,4)， 12348ak得到关于a ,a ,a ,a的线性方程组，这可以由下面的MATLAB程序完成，即输入程序 1234 syms al a2 a3 a4J=(-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/100

7、0*a1+121/100*a2-11 /10*a3+a4+723/20)八2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)八2+(a 4+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/25)八2;Ja1=diff(J,a1); Ja2=diff(J,a2);Ja3=

8、diff(J,a3);Ja4=diff(J,a4);Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3),Ja41=simple(Ja4),运行后屏幕显示J分别对aa2 ,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+67*a4 +767319/625Ja31 =1377283/2500*a1+23

9、667/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; 23667/250, 67, 18/5, 18;B=8442429/625,

10、 -767319/625, 232638/125, -14859/25;C=B/A, f=poly2sym(C)运行后屏幕显示拟合函数f及其系数C如下C =5.0911-14.19056.4102-8.2574f=716503 6958 4575 9/14 0737 4 8835532 8*xP-7 98 854 410255757 9/5 62 94 9953421312*xW+1804307491277693/281474976710656*x -4648521160813215/562949953421312 故所求的拟合曲线为f (X) = 5.0911 定-14.1905 X2 +

11、6.4102 x - 8.2574 .(4)编写下面的MATLAB程序估计其误差，并作出拟合曲线和数据的图形输入程 序 xi=-2.5 -1.7 -1.1 -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.50 68.04;n=length(xi);f=5.0911.*xi.A3-14.1905.*xi.A2+6.4102.*xi -8.2574;x=-2.5:0.01: 3.6;F=5.0911.*x.A3-14.1905.*x.A2+6.4102.*x -8.2574;fy=abs(f-y);

12、 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)的图形)运行后屏幕显示数据(x,七)与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2 及其数据点(X,七)和拟合曲线y=f(x)的图形(略).Ew =E1 =E2 =3.105 40.903 41.240 97.3函数气

13、（工）的选取及其MATLAB程序例7.3.1给出一组实验数据点(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=(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)，(7.3) 和(7.4)式估计其误差，作出拟合曲线.解 (1)在MATLAB工作窗口输入程序x=-8.5,-8.

14、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)的散点图)运行后屏幕显示数据的散点图(略).(3)编写下列MATLAB程序计算f (x)在(x,七)处的函数值，即输入程序 syms a

15、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)运行后屏幕显示关于。和屏勺线性方程组 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), a*exp(13/5*b),a*exp(5/2*b), a*exp(21/10*b),a*exp(3/2*b), a*exp(27/1

16、0*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;a*exp(71/10*b),a*exp(9/2*b), a*exp(13/5*b),a*exp(3/2*b),fi = a*exp(34/5*b), a*exp(18/5*b), a*exp(5/2*b), a*exp(27/10*b), fy=fi-y;fy2=fy.八2;J=sum(fy2)运行后屏幕显示误差平方和如下a*exp(17/

17、2*b),a*exp(87/10*b),a*exp(51/10*b),a*exp(17/5*b),a*exp(21/10*b),a*exp(18/5*b);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/1 0*b)-5917/100)八2+(a*exp(9/2*b)-2083/50)八2+(a*exp(18/5*b)-648/25 )+(a*exp(17/5*b)-2237/100)八2+(a*ex

18、p(13/5*b)-1347/100)八2+(a*ex p(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,b使J达到最小，只需利用极值的必要条件，得到关于a,b的线性方程组， 这可以由下面的MATLAB程序完成，即输入程序 syms a bJ=(a*exp(17/2*b)-22963/50)八2+(a*exp(87/10*b)-5281/100)八2 +(a*exp(71/10*b

19、)-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)-1487/100)八2+(a*exp(18/5*b)-121

20、1/50)八2;Ja=diff(J,a); Jb=diff(J,b);Ja1=simple(Ja), Jb1=simple(Jb),运行后屏幕显示J分别对a, b的偏导数如下Jal =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(2 7/10*b)-2507/25*exp(18/5*b)-22963/25*exp(17/2*b)-5281/50*exp(87 /10*b)-19827/50*exp(71/

21、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*e xp(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*exp(21/10*b)八2+8500*a*

22、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*ex p(21/10*b)+7100*a*exp(71/

23、10*b)八2+5100*a*exp(51/10*b)八2+4500*a*e xp(9/2*b)八2+7200*a*exp(18/5*b)八2+3400*a*exp(17/5*b)八2+2600*a*ex p(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程序求解线性方程组J1 =0，Jb1 =0，得a =b=2.811 00.581 6故所求的拟合曲线(7.13)为f (x) = 2.811 0e-0.581

24、6- .(7.14)(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.2652.81198.27165.6059.1741.6625.92 22.3713.4712.8711.876.6914.8724.22;n=length(xi); f=2.8110.*exp(-0.5816.*xi); x=-9:0.01: -1; F=2.8110.*exp(-0.5816.

25、*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 on plot(x,F,b-), hold off,legend(数据点(xi,yi),拟合曲线y=f(x)xlabel(x), ylabel(y),title(例7.3.1的数据点(xi,yi)和拟合曲线y=f(x)的图形)运行后屏幕显示数据(七，七)与拟合函数f的最大误差Ew = 390.141 5，平均误差 鸟=36.942 2和均方根误差E2=106.031 7及其数据点(, y)和拟合曲线y=f(x

26、)的图 形(略).Z7.4多项式拟合及其MATLAB程序例7.4.1给出一组数据点(X, y.)列入表7 3中，试用线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.表7-3例7.4. 1的一组数据（x , y.）X. 1-2.9-1.9-1.1-0.800.11.E.73.6-y53.9433.6820.8816.928.798.984.179.1219.88解 （1）首先根据表7 3给出的数据点（x., y.），用下列MATLAB程序画出散点图.在MATLAB工作窗口输入程序 x=-2.9 -1.9 -1.1 -0.8 0 0.1 1.5 2.7

27、3.6; y=53.94 33.68 20.88 16.92 8.79 8.98 4.17 9.12 19.88;plot(x,y, r*), legend(数据点(xi,yi) xlabel(x), ylabel(y),title(例7.4.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图(略).(3 )用作线性最小二乘拟合的多项式拟合的MATLAB程序求待定系数a k(k = 1, 2,3).输入程序 a=polyfit(x,y,2)运行后输出(7.16)式的系数 a =2.8302-7.37219.1382故拟合多项式为f (x) = 2.830 2x2 - 7.372 1x

28、 + 9.138 2 .(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=length(xi); f=2.8302.*xi.八2-7.3721.*xi+9.1382 x=-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

29、(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)的图形)运行后屏幕显示数据(x,七)与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2 及其数据点(x,y)和拟合曲线yfx)的图形(略). wEw =E1 =E2 =0.745 7,0.389 2,0.436 37.5拟合曲线的线性变换及其MATLAB程序例7.5.1给出一组实验数据点(x, y,)的横坐标向量为尤=(7.5 6.8 5.10 4

30、.5 3.6 3.4 2.6 2.5 2.1 1.5 2.7 3.6)，纵横坐标向量为此(359.26 165.60 59.17 41.66 25.92 22.37 13.47 12.87 11.87 6.69 14.87 24.22)，试用线性变换和线性 最小二乘法求拟合曲线，并用(7.2)，(7.3)和(7.4)式估计其误差，作出拟合曲线.解 (1)首先根据给出的数据点(x,y,)，用下列MATLAB程序画出散点图.在MATLAB工作窗口输入程序 x=7.5 6.8 5.10 4.5 3.6 3.4 2.6 2.5 2.1 1.5 2.7 3.6;y=359.26 165.60 59.17

31、 41.66 25.92 22.37 13.47 12.87 11.87 6.6914.8724.22;plot(x,y, r*), legend(数据点(xi,yi)xlabel(x), ylabel(y),title(例7.5.1的数据点(xi,yi)的散点图)运行后屏幕显示数据的散点图(略).(2) 根据数据散点图，取拟合曲线为y = a ebx (a 0,b。0),(7.19)其中a, b是待定系数.令Y = ln y, A = In a, B = b，则(7.19)化为Y = A + Bx ,在MATLAB工作窗口输入程序 x=7.56.85.104.53.63.42.62.52.1

32、1.52.73.6;y=359.26165.6059.1741.6625.9222.3713.4712.8711.87 6.6914.8724.22;Y=log(y); a=polyfit(x,Y,1); 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=ab

33、s(f-y); fy2=fy.八2; Ew=max(fy), E1=sum(fy)/n, E2=sqrt(sum(fy2)/n)运行后屏幕显示y = a ebx的系数b =0.624 1, a =2.703 9,数据(x., J.)与拟合函数f 的最大误差Ew =67.641 9，平均误差E1=8.677 6和均方根误差E2=20.711 3及其数 据点(x.，y.)和拟合曲线f (x) = 2.703 9 e0.624ix的图形(略).7.6函数逼近及其MATLAB程序最佳均方逼近的MATLAB主程序function yy1,a,WE=zjjfbj(f,X,Y,xx) m=size(f);n

34、=length(X);m=m(1);b=zeros(m,m); c=zeros(m,1); if n=length(Y)error X和Y的维数应该相同)endfor j=1:mfor k=1:mb(j,k)=0;for i=1:nb(j,k)=b(j,k)+feval(f(j,:),X(i)*feval(f(k,:),X(i);endendc(j)=0;for i=1:nc(j)=c(j)+feval(f(j,:),X(i)*Y(i);endenda=bc;WE=0;for i=1:nff=0;for j=1:mff=ff+a(j)*feval(f(j,:),X(i);endWE=WE+(Y

35、(i)-ff)*(Y(i)-ff);endif nargin=3return;endyy=;for i=1:ml=;for j=1:length(xx)l=l,feval(f(i,:),xx(j);endyy=yy l;endyy=yy*a; yy1=yy; a=a;WE;例7.6.1对数据X和Y,用函数J = 1, J = x, j = x2进行逼近，用所得到的逼近函 数计算在x = 6.5处的函数值，并估计误差.其中X=(1 3 4 567 8 9); Y=(-11 -13 -11 -7 -1 717 29).解在MATLAB工作窗口输入程序 X= 13456 789Y=-11-13-11

36、-7-1 71729;f=fun0;fun1;fun2; yy,a,WE=zjjfbj(f,X,Y,6.5) 运行后屏幕显示如下 yy = 2.75000000000003 a = -7.00000000000010-4.999999999999951.00000000000000WE = 7.172323350269439e-027例 7.6.2 对数据X和 Y，用函数 J = 1, J = x, j = x2，j = cos x，j = ex，J = sin x 进行逼近，其中 X=(0 0.50 1.00 1.50 2.00 2.50 3.00)，Y=(0 0.4794 0.8415 0

37、.9815 0.9126 0.5985 0.1645).解在MATLAB工作窗口输入程序 X= 00.501.001.502.002.503.00;Y=00.47940.84150.98150.91260.59850.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.00050.20370.39390.56560.71410.83480.9236Columns8 thr

38、ough140.97710.99260.96910.90690.80800.67660.5191Columns 15 through 16 0.34440.1642a = 0.38280.4070-0.39010.0765-0.45980.5653WE = 1.5769e-004 即，最佳逼近函数为y=0.3828+0.4070*x-0.3901*xA2+0.0765*exp(x) -0.4598*cos(x) +0.5653*sin(x)7.7三角多项式逼近及其MATLAB程序计算三角多项式的MATLAB主程序 function A,B,Y1,Rm=sanjiao(X,Y,X1,m) n=

39、length(X)-1;max1=fix(n-1)/2); if m max1 m=max1; end A=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:m B(i+1)=sin(i*X)*Y; A(i+1)=cos(i*X)*Y; end A=2*A/n; B=2*B/n; A(1)=A(1)/2;Y1=A(1); for k=1:m Y1=Y1+A(k+1)*cos(k*X1)+ B(k+1)*sin(k*X1); Tm=A(1)+A(k+1).*cos

40、(k*X)+ B(k+1).*sin(k*X); k=k+1; end Y;Tm; Rm=(sum(Y-Tm)2)/n;2 i兀 例7.7.1根据兀，兀上的n = 13, 60, 350个等距横坐标点x, =-n + -(i = 0,1,2, , n)和函数 f (x) = 2sin|.(1) 求f (x)的6阶三角多项式逼近，计算均方误差；(2) 将这三个三角多项式分别与f (x)的傅里叶级数I18 3 0n.f (x)=乙(-1) n+isin nx兀9n2 -1的前6项进行比较；回(3) 利用三角多项式分别计算X.= -2, 2.5的值；(4) 在同一坐标系中，画出函数f (x)，n =

41、 13,60,350的三角多项式和数据点的 图形.解(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

42、*sqrt(3)*n/(pi*(9*Z2-1)end(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,

43、Y2i,g-, 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-8随机数据点上的二元拟合及其MATLAB程序例7.8.1设节点(x,y,z)中的X和Y分别是在区间-3,3和-2.5, 3.5上的50个随机数，Z是函数Z=7-3x3e -x 2 -y2在(X,Y)的值，拟合点(X”YJ中的Xt=-3:0.2:3,

44、 丫工=-2.5:0.2:3.5.分别用二元拟合方法中最近邻内插法、三角基线性内插法、三角基 三次内插法和MATLAB 4网格化坐标方法计算在(X ,Y )处的值，作出它们的图形，并与 被拟和曲面进行比较.解 (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 - 丫.八2);%在每个随机点(X,Y)处计算Z的值.X1=-3:0.2:3;Y1=

45、-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 x八3 exp(-x八2 - y八2)的曲面 和节点的图形)%legend(拟合曲面，节点(xi,yi,zi)hold on%在当前图形上添加新图形.plot3(X,Y,Z, bo)%用兰色小圆圈画出每个节点(X,

46、Y,Z).hold of%结束在当前图形上添加新图形.运行后屏幕显示用最近邻内插法拟合函数Z=7-3x3e-x2-y2在两组不同节点处的曲面及其 插值4 (略).(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 - 丫.八2);%在每个随机点(X,Y)处计算Z 的值.X1=-3:0.2:3;Y1=-2.5:0.2:3.5;XI,Y

47、I = 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 x八3 exp(-x八2 - y八2 )的 曲面和节点的图形)%legend(拟合曲面，节点(xi,yi,zi)hold on%在当前图形上添加新图形.plot3(X,Y,Z, bo)%用兰色小圆圈画出每个节点(X,Y,Z).hold of%结束在

48、当前图形上添加新图形.运行后屏幕显示用三角基线性内插法拟合函数Z=7-3x3e-x2-y2在两组不同节点处的曲面和节点的图形及其插值Z/ (略).(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 - 丫.八2);%在每个随机点(X,Y)处计算Z 的值.X1=-3:0.2:3;Y1=-2.5:0.2:3.5;XI,YI = me

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

50、图形.运行后屏幕显示用三角基三次内插法拟合函数Z=7-3x3e-x2-y2在两组不同节点处的曲面 和节点的图形及其插值Z,(略).(4) MATLAB 4网格化坐标方法.输入程序 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 - 丫.八2);%在每个随机点(X,Y)处计算Z 的值.X1=-3:0.2:3; Y1=-2.5:0.2:3.5;XI,YI = mes

51、hgrid(X1,Y1);% 将坐标(XI,YI)网格化.ZI=griddata(X,Y,Z,XI,YI, v4) %计算在每个插值点(XI,YI)处的 插值ZI.mesh(XI,YI, ZI)%作二元拟合图形.xlabel(x), ylabel(y), zlabel(z),title( 用MATLAB 4网格化坐标方法拟合函数z =7-3 x八3 exp(-x八2 - y八2) 的曲面和节点的图形)%legend(拟合曲面，节点(xi,yi,zi)hold on%在当前图形上添加新图形.plot3(X,Y,Z, bo)%用兰色小圆圈画出每个节点(X,Y,Z).hold of%结束在当前图形上

52、添加新图形.运行后屏幕显示用MATLAB 4网格化坐标方法拟合函数Z=7-3x3e*-y2在两组不同节点处 的曲面和节点的图形及其插值Z/ (略).(5) 作被拟合曲面Z=7-3x佬-x2-y2和节点的图形.输入程序 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 - 丫.八2);%在每个随机点(X,Y)处计算Z的值.X1=-3.:0.1:3.;Y1=-2.5:0.1

53、:3.5;XI,YI = meshgrid(X1,Y1); % 将坐标(XI,YI)网格化.ZI=7-3* XI3 .* exp(-XI.八2 - YI2);mesh(XI,YI, ZI)%作二元拟合图形.xlabel(x), ylabel(y), zlabel(z), title(被拟合函数z =7-3 x八3 exp(-x八2 - y八2)的曲面和节点的图形 %legend(被拟合函数曲面，节点(xi,yi,zi)hold on%在当前图形上添加新图形.plot3(X,Y,Z, bo)%用兰色小圆圈画出每个节点(X,Y,Z).hold of%结束在当前图形上添加新图形.运行后屏幕显示被拟合

54、函数Z=7-3x3e-x2-y2的曲面和节点的图形及其函数值Z/ (略).7-9随机数据点上的11元拟合及其MATLAB程序例7.9.1首先利用MATLAB函数rand产生随机数据X,Y,Z,然后用线性变换 u = at + b (其中a = 5, b = 5)将随机数据X,Y ,Z1变换为节点坐标(X,Y,Z)，再用函数 W = 7 -3尤 3y( z +1) e - x 2 - y2 - z2 生成数据W,用三元最近邻内插法方法计算函数w在插值点 x.=-3:0.5:10, y,= -2:0.5:13, z广y,处拟合数据的值，并作其图形.解输入程序1 X1=-5+5*rand(10,1)

55、; Y1=-5+5*rand(10,1); Z1=Y1; X,Y,Z = meshgrid(X1,Y1,Z1); W=7-3* X.八3 .* Y.*(Z + 1).* exp(-X2 - Y.八2- Z2); xi=-3:0.5:10; yi=-2:0.5:13; zi=yi; XI,YI,ZI = meshgrid(xi,yi,zi); W1=griddata3(X, Y, Z, W, XI, YI, ZI, nearest); slice(XI,YI,ZI,W1,-2 4 9.5,9,-2 2 9), %shading flat %lighting flat xlabel(x), yla

56、bel(y), zlabel(z), title(被拟合函 数W=7-3X3Y(Z + 1)exp(-X2 - Y八2- Z八2); hold on colorbar( horiz) view(-3045)运行后屏幕显示三元线性拟合值及其图形(略).例7.9.2设节点(x,y,z,w)中的X,Y和Z分别是在区间-3,3和-2.5, 3.5, y=z 上的15个随机数，困是函数w = 2 + x e - x2 - y 2 - z 2 在(X,Y,Z)的值，拟合点(xy.z)中 的首=一3:0.2:3，y广-2.5:0.2:3.5，zy.,用linear方法计算拟合数据的值，并作 解 输入程序 x

57、=rand(15,1); y=rand(15,1);X1=-3+(3-(-3)*x;Y1=-2.5+(3.5-(-2.5)*y;Z1=Y1;X,Y,Z = meshgrid(X1,Y1,Z1);W=2+X.* exp(-X.八2 - Y.八2- Z.八2); xi=-3:0.2:3; yi=-2.5:0.2:3.5; zi=yi; X2,Y2,Z2=meshgrid(xi,yi,zi);W1=griddata3(X, Y, Z, W, X2,Y2,Z2,linear); slice(X2,Y2,Z2,W1,-1 0 1.5,2,-2 3), shading flat,lighting flat, xlabel(x), ylabel(y), zlabel(z), title(被拟合函数 W=2+X exp(-X八2 - Y八2- Z八2); hold on,colorbar(horiz), view(-3 5) 运行后屏幕显示三元线性拟合值及其图形(略).