函数迫近与曲线(面)拟合的MATLAB程序

本文格式为Word版，下载可任意编辑——函数迫近与曲线(面)拟合的MATLAB程序

函数迫近与曲线(面)拟合的MATLAB程序

7.1曲线拟合、误差及其MATLAB程序

例7.1.1已知函数yf(x)5x314x7sin2(2x)和一组数据(xi,yi)列入表7–1中，比较最大误差，平均误差，均方根误差和误差平方和.

表7–1例7.1.1的一组数据(

x,y)

解由给定的函数和数据，在MATLAB工作窗口输入

x=[-2.5,-1.7,-1.1,-0.8,0,0.1,0.5,3.6];n=length(x);

y=[-43.505.6911.3414.1601.02-6.37185.84];

f=5.*x.^3-14.*x+7.*(sin(2*pi*x)).^2;fy=abs(f-y);

fy2=fy.^2;[x',y',f',fy',fy2'],Ew=max(fy),

E1=sum(fy)/n,E2=sqrt((sum(fy2))/n),E=sum(fy2)

运行后屏幕显示如下

xyffyfy2

-2.5000-43.5000-43.12500.37500.1406

-1.70005.69005.56660.12340.0152

-1.100011.340011.16340.17660.0312

-0.800014.160014.97160.81160.6586

00000

0.10001.02001.02340.00340.0000

0.5000-6.3700-6.37500.00500.0000

3.6000185.8400185.29840.54160.2933

Ew=E1=E2=E=

0.81160.25460.37731.1390

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

例7.2.1给出一组数据点(xi,yi)列入表7–2中，试用线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.

表7–2例7.2.1的一组数据(

x,y)

解（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-13.126.5068.04];

plot(x,y,'r*'),

legend('试验数据(xi,yi)')

xlabel('x'),ylabel('y'),

title('例7.2.1的数据点(xi,yi)的散点图')

运行后屏幕显示数据的散点图（略）.

（3）编写以下MATLAB程序计算f(x)在(xi,yi)处的函数值，即输入程序

symsa1a2a3a4

x=[-2.5-1.7-1.1-0.800.11.52.73.6];

fi=a1.*x.^3+a2.*x.^2+a3.*x+a4

运行后屏幕显示关于a1,a2,a3和a4的线性方程组

函数迫近与曲线(面)拟合的MATLAB程序

-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]

编写构造误差平方和的MATLAB程序

y=[-192.9-85.50-36.15-26.52-9.10-8.43-13.126.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*a2+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+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*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,a4使J达到最小，只需利用极值的必要条件J0(k1,2,3,4)，ak

得到关于a1,a2,a3,a4的线性方程组，这可以由下面的MATLAB程序完成，即输入程序

symsa1a2a3a4

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+(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=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/625

Ja21=

32097579/25000*a1+1377283/2500*a2+23667/250*a3+67*a4

+767319/625

Ja31=

1377283/2500*a1+23667/250*a2+67*a3+18/5*a4-232638/125

Ja41=

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,-767319/625,232638/125,-14859/25];

函数迫近与曲线(面)拟合的MATLAB程序

运行后屏幕显示拟合函数f及其系数C如下

C=5.0911-14.19056.4102-8.2574

f=716503695845759/140737488355328*x^3

-7988544102557579/562949953421312*x^2

+1804307491277693/281474976710656*x

-4648521160813215/562949953421312

故所求的拟合曲线为

f(x)5.0911x314.1905x26.4102x8.2574.

（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=abs(f-y);fy2=fy.^2;Ew=max(fy),

E1=sum(fy)/n,E2=sqrt((sum(fy2))/n)

plot(xi,y,'r*'),holdon,plot(x,F,'b-'),holdoff

legend('数据点(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.10540.90341.2409

7.3函数rk(x)的选取及其MATLAB程序

例7.3.1给出一组试验数据点(xi,yi)的横坐标向量为

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.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)在(xi,yi)处的函数值，即输入程序

symsab

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];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),a*exp(13/5*b),a*exp(5/2*b),a*exp(21/10*b),

函数迫近与曲线(面)拟合的MATLAB程序

编写构造误差平方和的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),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+(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,b使J达到最小，只需利用极值的必要条件，得到关于a,b的线性方程组，这可以由下面的MATLAB程序完成，即输入程序

symsab

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/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分别对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*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/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*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)+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.81100.5816

故所求的拟合曲线（7.13）为

函数迫近与曲线(面)拟合的MATLAB程序

f(x)2.8110e0.5816x.（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.9222.37

13.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.*x);fy=abs(f-y);fy2=fy.^2;

Ew=max(fy),

E1=sum(fy)/n,E2=sqrt((sum(fy2))/n),plot(xi,y,'r*'),holdonplot(x,F,'b-'),holdoff,

legend('数据点(xi,yi)','拟合曲线y=f(x)')

xlabel('x'),ylabel('y'),

title('例7.3.1的数据点(xi,yi)和拟合曲线y=f(x)的图形')

运行后屏幕显示数据(xi,yi)与拟合函数f的最大误差Ew=390.1415，平均误差E1=36.9422和均方根误差E2=106.0317及其数据点(xi,yi)和拟合曲线y=f(x)的图形（略）.

7.4多项式拟合及其MATLAB程序

例7.4.1给出一组数据点(xi,yi)列入表7–3中，试用线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.

表7–3例7.4.1的一组数据(

x,y)

解（1）首先根据表7–3给出的数据点ii，用以下MATLAB程序画出散点图.在MATLAB工作窗口输入程序

x=[-2.9-1.9-1.1-0.800.11.52.73.6];

y=[53.9433.6820.8816.928.798.984.179.12

19.88];

plot(x,y,'r*'),legend('数据点(xi,yi)')

xlabel('x'),ylabel('y'),

title('例7.4.1的数据点(xi,yi)的散点图')

运行后屏幕显示数据的散点图（略）.

（3）用作线性最小二乘拟合的多项式拟合的MATLAB程序求待定系数ak(k1,2,3).输入程序

a=polyfit(x,y,2)

运行后输出（7.16）式的系数

a=

2.8302-7.37219.1382

故拟合多项式为

f(x)2.8302x27.3721x9.1382.

（4）编写下面的MATLAB程序估计其误差，并做出拟合曲线和数据的图形.输入程序

xi=[-2.9-1.9-1.1-0.800.11.52.73.6];

y=[53.9433.6820.8816.928.798.984.179.1219.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((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)的图形')

函数迫近与曲线(面)拟合的MATLAB程序

运行后屏幕显示数据(xi,yi)与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点(xi,yi)和拟合曲线y=f(x)的图形（略）.

Ew=E1=E2=

0.7457,0.3892,0.4363

7.5拟合曲线的线性变换及其MATLAB程序

例7.5.1给出一组试验数据点(xi,yi)的横坐标向量为x=（7.56.85.104.5

3.63.42.62.52.11.52.73.6），纵横坐标向量为y=（359.26165.6059.1741.6625.9222.3713.4712.8711.876.6914.8724.22），试用线性变换和线性最小二乘法求拟合曲线，并用（7.2），（7.3）和（7.4）式估计其误差，作出拟合曲线.

解（1）首先根据给出的数据点(xi,yi)，用以下MATLAB程序画出散点图.

在MATLAB工作窗口输入程序

x=[7.56.85.104.53.63.42.62.52.11.52.7

3.6];

y=[359.26165.6059.1741.6625.9222.3713.4712.8711.876.6914.8724.22];

plot(x,y,'r*'),legend('数据点(xi,yi)')

xlabel('x'),ylabel('y'),

title('例7.5.1的数据点(xi,yi)的散点图')

运行后屏幕显示数据的散点图（略）.

（2）根据数据散点图，取拟合曲线为

bxyae(a0,b0),（7.19）

其中a,b是待定系数.令Ylny,Alna,Bb，则（7.19）化为YABx.在MATLAB工作窗口输入程序

x=[7.56.85.104.53.63.42.62.52.11.52.7

3.6];

y=[359.26165.6059.1741.6625.9222.3713.4712.8711.876.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=abs(f-y);fy2=fy.^2;Ew=max(fy),E1=sum(fy)/n,

E2=sqrt((sum(fy2))/n)

bx运行后屏幕显示yae的系数b=0.6241，a=2.7039,数据(xi,yi)与拟合函数f

的最大误差Ew=67.6419，平均误差E1=8.6776和均方根误差E2=20.7113及其数据点(xi,yi)和拟合曲线f(x)2.7039e0.6241x的图形（略）.

7.6函数迫近及其MATLAB程序

最正确均方迫近的MATLAB主程序

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);

ifn~=length(Y)

error('X和Y的维数应当一致')

end

forj=1:m

fork=1:m

b(j,k)=0;

fori=1:n

函数迫近与曲线(面)拟合的MATLAB程序

end

end

c(j)=0;

fori=1:n

c(j)=c(j)+feval(f(j,:),X(i))*Y(i);

end

end

a=b\c;

WE=0;

fori=1:n

ff=0;

forj=1:m

ff=ff+a(j)*feval(f(j,:),X(i));

end

WE=WE+(Y(i)-ff)*(Y(i)-ff);

end

ifnargin==3

return;

end

yy=[];

fori=1:m

l=[];

forj=1:length(xx)

l=[l,feval(f(i,:),xx(j))];

end

yy=[yyl'];

end

yy=yy*a;yy1=yy';a=a';WE;

2例7.6.1对数据X和Y,用函数y1,yx,yx进行迫近，用所得到的迫近函

数计算在x6.5处的函数值，并估计误差.其中

X=(13456789);Y=(-11-13-11-7-171729).

解在MATLAB工作窗口输入程序

X=[13456789];Y=[-11-13-11-7-171729];

f=['fun0';'fun1';'fun2'];[yy,a,WE]=zjjfbj(f,X,Y,6.5)

运行后屏幕显示如下

yy=

2.75000000000003

a=

-7.00000000000010-4.999999999999951.00000000000000

WE=

7.172323350269439e-027

2ycosx，yex，ysinx例7.6.2对数据X和Y，用函数y1,yx,yx，

进行迫近，其中X=（00.501.001.502.002.503.00），Y=（00.47940.84150.98150.91260.59850.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=Columns1through7

-0.00050.20370.39390.56560.71410.83480.9236

Columns8through14

函数迫近与曲线(面)拟合的MATLAB程序

0.5191

Columns15through16

0.34440.1642

a=0.38280.4070-0.39010.0765-0.45980.5653WE=1.5769e-004

即，最正确迫近函数为

y=0.3828+0.4070*x-0.3901*x^2+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);

ifmmax1

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;fori=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);

fork=1:m

Y1=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根据[,]上的n13,60,350个等距横坐标点xi2in(i0,1,2,,n)和函数f(x)2sinx.3

（1）求f(x)的6阶三角多项式迫近，计算均方误差;

（2）将这三个三角多项式分别与f(x)的傅里叶级数

3nn1f(x)(1)sinnx2n19n1

的前6项进行比较;

（3）利用三角多项式分别计算Xi=-2,2.5的值;

（4）在同一坐标系中，画出函数f(x)，n13,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)

forn=1:6

bi=(-1)^(n+1)*18*sqrt(3)*n/(pi*(9*n^2-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;

函数迫近与曲线(面)拟合的MATLAB程序

[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--',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-3xe在(X,Y)的值，拟合点(XI,YI)中的XI=-3:0.2:3,YI=-2.5:0.2:3.5.分别用二元拟合方法中最近邻内插法、三角基线性内插法、三角基三次内插法和MATLAB4网格化坐标方法计算在(XI,YI)处的值，作出它们的图形，并与被拟和曲面进行比较.

解（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);%在每个随机点（X,Y）处计算Z的值.

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-3x^3exp(-x^2-y^2)的曲面和节点的图形')

%legend('拟合曲面','节点(xi,yi,zi)')

holdon%在当前图形上添加新图形.

plot3(X,Y,Z,'bo')%用兰色小圆圈画出每个节点(X,Y,Z).holdof%终止在当前图形上添加新图形.

运行后屏幕显示用最近邻内插法拟合函数Z=7-3x3e在两组不同节点处的曲面及其插值ZI（略）.

（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.2:3;

Y1=-2.5:0.2:3.5;

[XI,YI]=meshgrid(X1,Y1);%将坐标（XI,YI）网格化.

-x2-y23-x2-y2

函数迫近与曲线(面)拟合的MATLAB程序

ZI=griddata(X,Y,Z,XI,YI,'linear')%计算在每个插值点（XI,YI）处的插值ZI.

mesh(XI,YI,ZI)%作二元拟合图形.

xlabel('x'),ylabel('y'),zlabel('z'),

title('用三角基线性内插法拟合函数z=7-3x^3exp(-x^2-y^2)的曲面和节点的图形')

%legend('拟合曲面','节点(xi,yi,zi)')

holdon%在当前图形上添加新图形.

plot3(X,Y,Z,'bo')%用兰色小圆圈画出每个节点(X,Y,Z).holdof%终止在当前图形上添加新图形.

运行后屏幕显示用三角基线性内插法拟合函数Z=7-3x3e在两组不同节点处的曲面和节点的图形及其插值ZI（略）.

（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.5;

[XI,YI]=meshgrid(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-3x^3exp(-x^2-y^2)的曲面和节点的图形')

%legend('拟合曲面','节点(xi,yi,zi)')

holdon%在当前图形上添加新图形.

plot3(X,Y,Z,'bo')%用兰色小圆圈画出每个节点(X,Y,Z).holdof%终止在当前图形上添加新图形.

运行后屏幕显示用三角基三次内插法拟合函数Z=7-3x3e在两组不同节点处的曲面和节点的图形及其插值ZI（略）.

（4）MATLAB4网格化坐标方法.

输入程序

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.5;

[XI,YI]=meshgrid(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('用MATLAB4网格化坐标方法拟合函数z=7-3x^3exp(-x^2-y^2)的曲面和节点的图形')

%legend('拟合曲面','节点(xi,yi,zi)')

-x2-y2-x2-y2

函数迫近与曲线(面)拟合的MATLAB程序

holdon%在当前图形上添加新图形.

plot3(X,Y,Z,'bo')%用兰色小圆圈画出每个节点(X,Y,Z).holdof%终止在当前图形上添加新图形.

运行后屏幕显示用MATLAB4网格化坐标方法拟合函数Z=7-3x3e-x

的曲面和节点的图形及其插值ZI（略）.

3222-y2在两组不同节点处（5）作被拟合曲面Z=7-3xe-x-y和节点的图形.

输入程序

x=rand(50,1);

y=rand(50,1);%生成50个一元均匀分布随机数x和y，x，y.

X=-3+(3-(-3))*x;%利用