1、1 曲線擬合的線性最小二乘法及其MATLAB程序例7.2.1 給出一組數(shù)據(jù)點列入表72中，試用線性最小二乘法求擬合曲線，并用（7.2），（7.3）和（7.4）式估計其誤差，作出擬合曲線.表72 例7.2.1的一組數(shù)據(jù)xi-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6yi-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04解 （1）在MATLAB工作窗口輸入程序 x=-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

2、.10 -8.43 -13.12 6.50 68.04;plot(x,y,r*),legend(實驗數(shù)據(jù)(xi,yi)xlabel(x), ylabel(y),title(例7.2.1的數(shù)據(jù)點(xi,yi)的散點圖)運行后屏幕顯示數(shù)據(jù)的散點圖（略）. （3）編寫下列MATLAB程序計算在處的函數(shù)值，即輸入程序 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.2+ a3.*x+ a4運行后屏幕顯示關(guān)于a1,a2, a3和a4的線性方程組fi = -125/8*a1+25/4*a2-5/2*a3+

3、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編寫構(gòu)造誤差平方和的MATLAB程序 y=-192.9 -85.50 -36.15 -26.52 -9.10 -8.

4、43 -13.12 6.50 68.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;

5、 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/10*a3+a4+843/100)2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)2+(

6、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為求使達到最小，只需利用極值的必要條件 ，得到關(guān)于的線性方程組，這可以由下面的MATLAB程序完成，即輸入程序 syms a1 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/1000*a1+121/100*a2-11/10*a3+a4+723/20)2+(-64/1

7、25*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/25)2; Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3); Ja4=diff(J,a4);Ja11=simple(Ja1), Ja21=si

8、mple(Ja2), Ja31=simple(Ja3), Ja41=simple(Ja4),運行后屏幕顯示J分別對a1, a2 ,a3 ,a4的偏導數(shù)如下Ja11=/10000*a1+/25000*a2+/2500*a3+23667/250*a4-/625Ja21 = /25000*a1+/2500*a2+23667/250*a3+67*a4+/625Ja31 =/2500*a1+23667/250*a2+67*a3+18/5*a4-/125Ja41 =23667/250*a1+67*a2+18/5*a3+18*a4+14859/25解線性方程組Ja11 =0，Ja21 =0，Ja31 =0，

9、Ja41 =0，輸入下列程序A=/10000, /25000, /2500, 23667/250; /25000, /2500, 23667/250, 67; /2500, 23667/250, 67, 18/5; 23667/250, 67, 18/5, 18;B=/625, -/625, /125, -14859/25; C=B/A, f=poly2sym(C)運行后屏幕顯示擬合函數(shù)f及其系數(shù)C如下C = 5.0911 -14.1905 6.4102 -8.2574f=5759/5328*x3-57579/1312*x2+77693/0656*x-13215/1312 故所求的擬合曲線為.

10、（4）編寫下面的MATLAB程序估計其誤差，并作出擬合曲線和數(shù)據(jù)的圖形.輸入程序 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.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=abs(f-y); fy2=fy.2; Ew=max(fy

11、), E1=sum(fy)/n, E2=sqrt(sum(fy2)/n)plot(xi,y,r*), hold on, plot(x,F,b-), hold offlegend(數(shù)據(jù)點(xi,yi),擬合曲線y=f(x), xlabel(x), ylabel(y),title(例7.2.1的數(shù)據(jù)點(xi,yi)和擬合曲線y=f(x)的圖形)運行后屏幕顯示數(shù)據(jù)與擬合函數(shù)f的最大誤差Ew，平均誤差E1和均方根誤差E2及其數(shù)據(jù)點和擬合曲線y=f(x)的圖形（略）.Ew = E1 = E2 =3.105 4 0.903 4 1.240 97.3 函數(shù)的選取及其MATLAB程序例7.3.1 給出一組實驗

12、數(shù)據(jù)點的橫坐標向量為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.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-2

13、.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(實驗數(shù)據(jù)(xi,yi)xlabel(x), ylabel(y),title(例7.3.1的數(shù)據(jù)點(xi,yi)的散點圖)運行后屏幕顯示數(shù)據(jù)的散點圖（略）.（3）編寫下列MATLAB程序計算在處的函數(shù)值，即輸入程序 syms a bx=-8.5,-8.7,-7.1,-6.8,-5.10,-4.5,-3.6,-3.4,-

14、2.6,-2.5,-2.1,-1.5,-2.7,-3.6; fi=a.*exp(-b.*x)運行后屏幕顯示關(guān)于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), 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)編寫構(gòu)造誤差平方和的MATLAB程序如下y

15、=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), 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/

16、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+(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

17、/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為求使達到最小，只需利用極值的必要條件，得到關(guān)于的線性方程組，這可以由下面的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)-19827/100)2+(a*exp(34/5*b)-828/5)2+(a*exp(51/10*b)-5917/1

18、00)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)-1211/50)2;Ja=diff(J,a); Jb=diff(J,b); Ja1=simple(Ja), Jb1=simple(Jb),運行后屏幕顯示J

19、分別對的偏導數(shù)如下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*exp(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/1

20、0*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*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)-*exp(18/5*b)-*exp

21、(17/2*b)-*exp(87/10*b)-*exp(71/10*b)-76058*exp(17/5*b)-*exp(34/5*b)-35022*exp(13/5*b)-*exp(51/10*b)-32175*exp(5/2*b)-*exp(9/2*b)-24927*exp(21/10*b)+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

22、/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）為e. （7.14）（4）根據(jù)（7.2），（7.3），（7.4）和（7.14）式編寫下面的MATLAB程序估計其誤差，并做出擬合曲線和數(shù)據(jù)的圖形.輸入程序 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.

23、27 165.60 59.17 41.66 25.92 22.37 13.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: -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(數(shù)據(jù)點(xi,yi),擬合曲線y=f(x)xlabel

24、(x), ylabel(y), title(例7.3.1的數(shù)據(jù)點(xi,yi)和擬合曲線y=f(x)的圖形)運行后屏幕顯示數(shù)據(jù)與擬合函數(shù)f的最大誤差Ew = 390.141 5，平均誤差E1=36.942 2和均方根誤差E2=106.031 7及其數(shù)據(jù)點和擬合曲線y=f(x)的圖形（略）.7.4 多項式擬合及其MATLAB程序例7.4.1 給出一組數(shù)據(jù)點列入表73中，試用線性最小二乘法求擬合曲線，并用（7.2），（7.3）和（7.4）式估計其誤差，作出擬合曲線.表73 例7.4.1的一組數(shù)據(jù)xi-2.9 -1.9 -1.1 -0.8 0 0.1 1.5 2.7 3.6yi53.94 33.68

25、 20.88 16.92 8.79 8.98 4.17 9.12 19.88解 （1）首先根據(jù)表73給出的數(shù)據(jù)點，用下列MATLAB程序畫出散點圖.在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.12 19.88;plot(x,y,r*), legend(數(shù)據(jù)點(xi,yi)xlabel(x), ylabel(y), title(例7.4.1的數(shù)據(jù)點(xi,yi)的散點圖)運行后屏幕顯示數(shù)據(jù)的散點圖（略）.（3）用作線性最小二乘擬合的多項式擬合

26、的MATLAB程序求待定系數(shù) .輸入程序 a=polyfit(x,y,2)運行后輸出（7.16）式的系數(shù)a =2.8302 -7.3721 9.1382故擬合多項式為 .（4）編寫下面的MATLAB程序估計其誤差，并做出擬合曲線和數(shù)據(jù)的圖形.輸入程序 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.1382x=-2.9:0.001:3.6;F=2.8302.*x.2-7.3

27、721.*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(數(shù)據(jù)點(xi,yi),擬合曲線y=f(x)xlabel(x), ylabel(y), title(例7.4.1 的數(shù)據(jù)點(xi,yi)和擬合曲線y=f(x)的圖形)運行后屏幕顯示數(shù)據(jù)與擬合函數(shù)f的最大誤差Ew，平均誤差E1和均方根誤差E2及其數(shù)據(jù)點(xi,yi)和擬合曲線y=f(x)的圖形（略）.Ew = E1 = E2 =0.745 7, 0.389 2, 0.436 37.

28、5 擬合曲線的線性變換及其MATLAB程序 例7.5.1 給出一組實驗數(shù)據(jù)點的橫坐標向量為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 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）首先根據(jù)給出的數(shù)據(jù)點，用下列MATLAB程序畫出散點圖.在MATLAB工作窗口輸入程序 x=7.5 6.8 5.10 4

29、.5 3.6 3.4 2.6 2.5 2.1 1.5 2.7 3.6;y=359.26 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(數(shù)據(jù)點(xi,yi)xlabel(x), ylabel(y),title(例7.5.1的數(shù)據(jù)點(xi,yi)的散點圖)運行后屏幕顯示數(shù)據(jù)的散點圖（略）.（2）根據(jù)數(shù)據(jù)散點圖，取擬合曲線為 e , （7.19） 其中是待定系數(shù).令，則（7.19）化為.在MATLAB工作窗口輸入程序 x=7.5 6.8 5.10 4.5 3.6 3.4

30、2.6 2.5 2.1 1.5 2.7 3.6;y=359.26 165.60 59.17 41.66 25.92 22.37 13.47 12.87 11.87 6.69 14.87 24.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(數(shù)據(jù)點(xi,yi),擬合曲線y=f(x)title(例7.5.1 的數(shù)據(jù)點(xi,y

31、i)和擬合曲線y=f(x)的圖形)fy=abs(f-y); fy2=fy.2; Ew=max(fy), E1=sum(fy)/n, E2=sqrt(sum(fy2)/n)運行后屏幕顯示e的系數(shù)b =0.624 1，a =2.703 9,數(shù)據(jù)與擬合函數(shù)f的最大誤差Ew =67.641 9，平均誤差E1=8.677 6和均方根誤差E2=20.711 3及其數(shù)據(jù)點和擬合曲線e的圖形（略）.7.6 函數(shù)逼近及其MATLAB程序最佳均方逼近的MATLAB主程序function yy1,a,WE=zjjfbj(f,X,Y,xx)m=size(f);n=length(X);m=m(1);b=zeros(m,

32、m); c=zeros(m,1);if n=length(Y) error(X和Y的維數(shù)應該相同)endfor j=1:m for k=1:m b(j,k)=0; for i=1:n b(j,k)=b(j,k)+feval(f(j,:),X(i)*feval(f(k,:),X(i); end end c(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)

33、;endif nargin=3 return;endyy=;for i=1:m l=; for j=1:length(xx) l=l,feval(f(i,:),xx(j); end yy=yy l;end yy=yy*a; yy1=yy; a=a;WE;例7.6.1 對數(shù)據(jù)X和Y, 用函數(shù)進行逼近，用所得到的逼近函數(shù)計算在處的函數(shù)值，并估計誤差.其中X=(1 3 4 5 6 7 8 9); Y=(-11 -13 -11 -7 -1 7 17 29).解 在MATLAB工作窗口輸入程序 X= 1 3 4 5 6 7 8 9; Y=-11 -13 -11 -7 -1 7 17 29;f=fun0;

34、fun1;fun2; yy,a,WE=zjjfbj(f,X,Y,6.5)運行后屏幕顯示如下yy = 2.003a = -7.010 -4.995 1.000WE = 7.9439e-027例7.6.2 對數(shù)據(jù)X和Y，用函數(shù)，e，進行逼近，其中X=（0 0.50 1.00 1.50 2.00 2.50 3.00），Y=（0 0.4794 0.8415 0.9815 0.9126 0.5985 0.1645）.解 在MATLAB工作窗口輸入程序 X= 0 0.50 1.00 1.50 2.00 2.50 3.00;Y=0 0.4794 0.8415 0.9815 0.9126 0.5985 0.1

35、645;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.8348 0.9236 Columns 8 through 14 0.9771 0.9926 0.9691 0.9069 0.8080 0.6766 0.5191 Columns 15 through 16 0.3444 0.1642a = 0.3828 0

36、.4070 -0.3901 0.0765 -0.4598 0.5653WE = 1.5769e-004即，最佳逼近函數(shù)為y=0.3828+0.4070*x-0.3901*x2+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= 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

37、(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;end A=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).*cos(k*X)+ B(k+1).*sin(k*X); k=k+1;endY;Tm; Rm=(sum(Y-Tm).2)/n;例7.7.1 根據(jù)上的個等距橫坐標點 和函數(shù).（1）求的6階三角多項式逼近，

38、計算均方誤差;（2）將這三個三角多項式分別與的傅里葉級數(shù)的前6項進行比較;（3）利用三角多項式分別計算Xi= -2, 2.5的值;（4）在同一坐標系中，畫出函數(shù)，的三角多項式和數(shù)據(jù)點的圖形.解 （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,R

39、m3=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*n2-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(X

40、2,Y2,Xi,6);A3,B3,Y3i,R3m=sanjiao(X3,Y3,Xi,6);plot(X1,Y1,r*, Xi, Y1i,b-,Xi, Y2i,g-, Xi, Y3i, m:, Xi, f, k-.)xlabel(x),ylabel(y)legend(數(shù)據(jù)點(xi,yi),n=13的三角多項式,n=60的三角多項式,n=350的三角多項式,函數(shù)f(x)title(例7.7.1 的數(shù)據(jù)點(xi,yi)、n=13，60，350的三角多項式T3和函數(shù)f(x)的圖形)運行后圖形（略）.7.8 隨機數(shù)據(jù)點上的二元擬合及其MATLAB程序例7.8.1 設(shè)節(jié)點(X,Y,Z)中的X和Y分別是在區(qū)

41、間和上的50個隨機數(shù)，Z是函數(shù)Z=7-3x3e在(X,Y)的值，擬合點(XI,YI)中的XI=-3:0.2:3, YI=-2.5:0.2:3.5.分別用二元擬合方法中最近鄰內(nèi)插法、三角基線性內(nèi)插法、三角基三次內(nèi)插法和MATLAB 4網(wǎng)格化坐標方法計算在(XI,YI)處的值，作出它們的圖形，并與被擬和曲面進行比較.解 （1）最近鄰內(nèi)插法.輸入程序 x=rand(50,1);y=rand(50,1); %生成50個一元均勻分布隨機數(shù)x和y， x，y . X=-3+(3-(-3)*x;%利用x生成的隨機變量. Y=-2.5+(3.5-(-2.5)*y; %利用y生成的隨機變量.Z=7-3* X.3

42、.* exp(-X.2 - Y.2); %在每個隨機點（X,Y）處計算Z的值.X1=-3:0.2:3;Y1=-2.5:0.2:3.5;XI,YI = meshgrid(X1,Y1); %將坐標（XI,YI）網(wǎng)格化.ZI=griddata(X,Y,Z,XI,YI, nearest) %計算在每個插值點（XI,YI）處的插值ZI.mesh(XI,YI, ZI) %作二元擬合圖形.xlabel(x), ylabel(y), zlabel(z),title(用最近鄰內(nèi)插法擬合函數(shù)z =7-3 x3 exp(-x2 - y2) 的曲面和節(jié)點的圖形)%legend(擬合曲面,節(jié)點(xi,yi,zi)hol

43、d on %在當前圖形上添加新圖形.plot3(X,Y,Z, bo) %用蘭色小圓圈畫出每個節(jié)點(X,Y,Z).hold of %結(jié)束在當前圖形上添加新圖形.運行后屏幕顯示用最近鄰內(nèi)插法擬合函數(shù)Z=7-3x3e在兩組不同節(jié)點處的曲面及其插值ZI（略）. （2）三角基線性內(nèi)插法.輸入程序 x=rand(50,1);y=rand(50,1); %生成50個一元均勻分布隨機數(shù)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); %在每個隨機

44、點（X,Y）處計算Z的值.X1=-3:0.2:3;Y1=-2.5:0.2:3.5;XI,YI = meshgrid(X1,Y1); %將坐標（XI,YI）網(wǎng)格化.ZI=griddata(X,Y,Z,XI,YI, linear) %計算在每個插值點（XI,YI）處的插值ZI.mesh(XI,YI, ZI) %作二元擬合圖形.xlabel(x), ylabel(y), zlabel(z),title(用三角基線性內(nèi)插法擬合函數(shù)z =7-3 x3 exp(-x2 - y2) 的曲面和節(jié)點的圖形)%legend(擬合曲面,節(jié)點(xi,yi,zi)hold on %在當前圖形上添加新圖形.plot3(X

45、,Y,Z, bo) %用蘭色小圓圈畫出每個節(jié)點(X,Y,Z).hold of %結(jié)束在當前圖形上添加新圖形.運行后屏幕顯示用三角基線性內(nèi)插法擬合函數(shù)Z=7-3x3e在兩組不同節(jié)點處的曲面和節(jié)點的圖形及其插值ZI（略）. （3）三角基三次內(nèi)插法.輸入程序 x=rand(50,1);y=rand(50,1); %生成50個一元均勻分布隨機數(shù)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=-

46、3:0.2:3;Y1=-2.5:0.2:3.5;XI,YI = meshgrid(X1,Y1); %將坐標（XI,YI）網(wǎng)格化.ZI=griddata(X,Y,Z,XI,YI, cubic) %計算在每個插值點（XI,YI）處的插值ZI.mesh(XI,YI, ZI) %作二元擬合圖形.xlabel(x), ylabel(y), zlabel(z),title(用三角基三次內(nèi)插法擬合函數(shù)z =7-3 x3 exp(-x2 - y2) 的曲面和節(jié)點的圖形)%legend(擬合曲面,節(jié)點(xi,yi,zi)hold on %在當前圖形上添加新圖形.plot3(X,Y,Z, bo) %用蘭色小圓圈畫

47、出每個節(jié)點(X,Y,Z).hold of %結(jié)束在當前圖形上添加新圖形.運行后屏幕顯示用三角基三次內(nèi)插法擬合函數(shù)Z=7-3x3e在兩組不同節(jié)點處的曲面和節(jié)點的圖形及其插值ZI（略）. （4）MATLAB 4網(wǎng)格化坐標方法.輸入程序 x=rand(50,1);y=rand(50,1); %生成50個一元均勻分布隨機數(shù)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=

48、-2.5:0.2:3.5;XI,YI = meshgrid(X1,Y1); %將坐標（XI,YI）網(wǎng)格化.ZI=griddata(X,Y,Z,XI,YI, v4) %計算在每個插值點（XI,YI）處的插值ZI.mesh(XI,YI, ZI) %作二元擬合圖形.xlabel(x), ylabel(y), zlabel(z),title(用MATLAB 4網(wǎng)格化坐標方法擬合函數(shù)z =7-3 x3 exp(-x2 - y2) 的曲面和節(jié)點的圖形)%legend(擬合曲面,節(jié)點(xi,yi,zi)hold on %在當前圖形上添加新圖形.plot3(X,Y,Z, bo) %用蘭色小圓圈畫出每個節(jié)點(X,Y,Z).hold of %結(jié)束在當前圖形上添加新圖形.運行后屏幕顯示用MATLAB 4網(wǎng)格化坐標方法擬合函數(shù)Z=7-3x3e在兩組不同節(jié)點處的曲面和節(jié)點的圖形及其插值ZI（略）. （5）作被擬合曲面Z=7-3x3e和節(jié)點的圖形.輸入程序 x=rand(50,1);y=rand(50,1); %生成50個一元均勻分布隨機數(shù)x和y， x，y . X=-3+(3