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

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

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

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

表7–2 例7.2.1的一组数据),(y x

解 (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.10 -8.43 -13.12 6.50 68.04];

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

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

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

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

运行后屏幕显示数据的散点图(略).

(3)编写下列MATLAB 程序计算)(x f 在),(i i y x 处的函数值,即输入程序

>> syms a1 a2 a3 a4

x=[-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

运行后屏幕显示关于a 1,a 2, a 3和a 4的线性方程组

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]

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

>> y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.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; 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*a 2+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

为求4321,,,a a a a 使J 达到最小,只需利用极值的必要条件0=??k

a J )4,3,2,1(=k ,

得到关于4321,,,a a a a 的线性方程组,这可以由下面的MA TLAB 程序完成,即输入程序

>> syms a1 a2 a3 a4

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+(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=diff(J,a3); Ja4=diff(J,a4);

Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3), Ja41=simple(Ja4),

运行后屏幕显示J 分别对a 1, a 2 ,a 3 ,a 4的偏导数如下

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

解线性方程组Ja 11 =0,Ja 21 =0,Ja 31 =0,Ja 41 =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];

C=B/A, f=poly2sym(C)

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

C = 5.0911 -14.1905 6.4102 -8.2574

f=716503695845759/140737488355328*x^3

-7988544102557579/562949953421312*x^2

+1804307491277693/281474976710656*x

-4648521160813215/562949953421312

故所求的拟合曲线为

8.25746.410214.19055.0911)(23-+-=x x x x f .

(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.^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*'), hold on, plot(x,F,'b-'), hold off

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

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

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

运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差E w ,平均误差E 1和均方根误差E 2及其数据点),(i i y x 和拟合曲线y =f (x )的图形(略).

Ew = E1 = E2 =

3.105 4 0.903 4 1.240 9

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

例7.3.1 给出一组实验数据点),(i i y 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.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 程序计算)(x f 在),(i i y x 处的函数值,即输入程序

>> syms a b

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), 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), 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*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

为求b a ,使J 达到最小,只需利用极值的必要条件,得到关于b a ,的线性方程组,这可以由下面的MA TLAB 程序完成,即输入程序

>> syms a 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/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

Ja=diff(J,a); Jb=diff(J,b);

Ja1=simple(Ja), Jb1=simple(Jb),

运行后屏幕显示J 分别对b a ,的偏导数如下

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*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*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/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 程序求解线性方程组J a1 =0,J b1 =0,得

a = b=

2.811 0 0.581 6

故所求的拟合曲线(7.13)为

0811.2)(=x f e x 5816.0-. (7.14)

(4)根据(7.2),(7.3),(7.4)和(7.14)式编写下面的MA TLAB 程序估计其误差,并做出拟合曲线和数据的图形.输入程序

>> 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.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 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)的图形')

运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差E w = 390.141 5,平均误差E 1=36.942 2和均方根误差E 2=106.031 7及其数据点),(i i y x 和拟合曲线y =f (x )的图形(略).

7.4 多项式拟合及其MATLAB 程序

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

表7–3 例7.4.1的一组数据),(y x

解 (1)首先根据表7–3给出的数据点i i ,用下列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('数据点(xi,yi)')

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

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

运行后屏幕显示数据的散点图(略).

(3)用作线性最小二乘拟合的多项式拟合的MATLAB 程序求待定系数k a )3,2,1(=k .输入程序

>> a=polyfit(x,y,2)

运行后输出(7.16)式的系数

a =

2.8302 -7.3721 9.1382

故拟合多项式为

2138.91372.72830.2)(2+-=x x x f .

(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((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)的图形')

运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差E w ,平均误差E1和均方根误差E 2及其数据点(x i ,y i )和拟合曲线y =f (x )的图形(略).

Ew = E1 = E2 =

0.745 7, 0.389 2, 0.436 3

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

例7.5.1 给出一组实验数据点),(i i y x 的横坐标向量为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 1

4.87 24.22),试用线性变换和线性最小二乘法求拟合曲线,并用(7.2),(7.3)和(7.4)式估计其误差,作出拟合曲线.

解 (1)首先根据给出的数据点),(i i y x ,用下列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 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.5.1的数据点(xi,yi)的散点图')

运行后屏幕显示数据的散点图(略).

(2)根据数据散点图,取拟合曲线为

a y =e bx )0,0(≠>

b a , (7.19)

其中b a ,是待定系数.令b B a A y Y ===,ln ,ln ,则(7.19)化为Bx A Y +=.在MATLAB 工作窗口输入程序

>> x=[7.5 6.8 5.10 4.5 3.6 3.4 2.6 2.5 2.1 1.5 2.7

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('数据点(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)

运行后屏幕显示a y =e bx 的系数b =0.624 1,a =2.703 9,数据),(i i y x 与拟合函数f

的最大误差Ew =67.641 9,平均误差E 1=8.677 6和均方根误差E 2=20.711 3及其数据点),(i i y x 和拟合曲线9703.2)(=x f e x 1624.0的图形(略).

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

if n~=length(Y)

error('X 和Y 的维数应该相同')

end

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

end

end

a=b\c;

WE=0;

for i=1:n

ff=0;

for j=1:m

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

end

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

end

if nargin==3

return ;

end

yy=[];

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 对数据X 和Y , 用函数2,,1x y x y y ===进行逼近,用所得到的逼近函数计算在 6.5=x 处的函数值,并估计误差.其中

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';'fun1';'fun2']; [yy,a,WE]=zjjfbj(f,X,Y,6.5)

运行后屏幕显示如下

yy =

2.75000000000003

a =

-7.00000000000010 -4.99999999999995 1.00000000000000

WE =

7.172323350269439e-027

例7.6.2 对数据X 和Y ,用函数2,,1x y x y y ===,

x y cos =,=y e x

,x y sin =进行逼近,其中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.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.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.1642

a = 0.3828 0.4070 -0.3901 0.0765 -0.4598 0.5653 WE = 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);

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(k*X)+ B(k+1).*sin(k*X); k=k+1;

end Y;Tm; Rm=(sum(Y-Tm).^2)/n;

例7.7.1 根据],[ππ-上的350,60,13=n 个等距横坐标点n

i x i π+π-=2

),,2,1,0(n i =和函数3

sin 2)(x x f =. (1)求)(x f 的6阶三角多项式逼近,计算均方误差;

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

nx n n x f n n sin 1

9)1(318)(121∑∞=+--π= 的前6项进行比较;

(3)利用三角多项式分别计算X i = -2, 2.5的值;

(4)在同一坐标系中,画出函数)(x f ,350,60,13=n 的三角多项式和数据点的图形.

解 (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: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; 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--', 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[-和]5.3,5.2[-上的50个随机数,Z 是函数Z =7-3x 3e 2

2y - -x 在(X,Y )的值,拟合点(X I ,Y I )中的X I =-3:0.2:3, Y I =-2.5:0.2:3.5.分别用二元拟合方法中最近邻内插法、三角基线性内插法、三角基三次内插法和MATLAB 4网格化坐标方法计算在(X I ,Y I )处的值,作出它们的图形,并与被拟和曲面进行比较.

解 (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;

[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,Y,Z). hold of %结束在当前图形上添加新图形.

运行后屏幕显示用最近邻内插法拟合函数Z =7-3x 3e 2

2y - -x 在两组不同节点处的曲面及其插值Z I (略).

(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 )网格化.

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 %结束在当前图形上添加新图形.

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

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

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

hold on %在当前图形上添加新图形.

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

运行后屏幕显示用三角基三次内插法拟合函数Z =7-3x 3e 22y - -x 在两组不同节点处的曲面和节点的图形及其插值Z I (略).

(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 - 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('用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 %结束在当前图形上添加新图形.

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

22y - -x 在两组不同节点处的曲面和节点的图形及其插值ZI (略).

(5)作被拟合曲面Z =7-3x 3e 22y - -x 和节点的图形.

输入程序

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

Y1=-2.5:0.1:3.5;

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

ZI=7-3* XI.^3 .* exp(-XI.^2 - YI.^2);

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 %结束在当前图形上添加新图形.

运行后屏幕显示被拟合函数Z =7-3x 3e 22y - -x 的曲面和节点的图形及其函数值ZI (略).

7.9 随机数据点上的n 元拟合及其MATLAB 程序

例7.9.1 首先利用MATLAB 函数rand 产生随机数据X 1,Y 1,Z 1,然后用线性变换b at u += (其中5,5-==b a )将随机数据X 1,Y 1 ,Z 1 变换为节点坐标(X,Y,Z ),再用函数

)1(373+-=z y x w e

2

22z y x ---生成数据W , 用三元最近邻内插法方法计算函数w 在插值点x i ,10:5.0:3-= y i ,13:5.0:2-=z i =y i 处拟合数据的值,并作其图形. 解 输入程序

>> X1=-5+5*rand(10,1);

Y1=-5+5*rand(10,1);

Z1=Y1;

[X,Y,Z] = meshgrid(X1,Y1,Z1);

W=7-3* X.^3 .* Y.*(Z+1).* exp(-X.^2 - Y.^2- Z.^2);

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'), ylabel('y'), zlabel('z'),

title('被拟合函数W=7-3X^3Y(Z+1)exp(-X^2 - Y^2- Z^2) '); hold on

colorbar('horiz')

view([-30 45])

运行后屏幕显示三元线性拟合值及其图形(略).

例7.9.2 设节点(X,Y,Z,W )中的X,Y 和Z 分别是在区间]3,3[-和]5.3,5.2[-,Y =Z

上的15个随机数,W 是函数x w +=2e 2

22z y x ---在(X,Y,Z )的值,拟合点(x i ,y i ,z i )中

的x i =-3:0.2:3, y i =-2.5:0.2:3.5,z i =y i , 用'linear '方法计算拟合数据的值,并作其图形.

解 输入程序

>> x=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])

运行后屏幕显示三元线性拟合值及其图形(略).

MATLAB中如何直接曲线拟合

MATLAB中如何直接曲线拟合,而不使用cftool的GUI 界面 我们知道在MATLAB中有个很方便的曲线拟合工具:cftool 最基本的使用方法如下,假设我们需要拟合的点集存放在两个向量X和Y中,分别储存着各离散点的横坐标和纵坐标,则在MATLAB中直接键入命令 cftool(X,Y) 就会弹出Curve Fitting Tool的GUI界面,点击界面上的fitting即可开始曲线拟合。 MATLAB提供了各种曲线拟合方法,例如:Exponential, Fourier, Gaussing, Interpolant, Polynomial, Power, Rational, Smoothing Spline, Sum of Functions, Weibull等,当然,也可以使用 Custom Equations. cftool不仅可以绘制拟合后的曲线、给出拟合参数,还能给出拟合好坏的评价 参数(Goodness of fit)如SSE, R-square, RMSE等数据,非常好用。但是如果我们已经确定了拟合的方法,只需要对数据进行计算,那么这种GUI的操作方式就不太适合了,比如在m文件中就不方便直接调用cftool。 MATLAB已经给出了解决办法,可以在cftool中根据情况生成特定的m文件,让我们直接进行特定的曲线拟合并给出参数。具体方法在帮助文件的如下文档中" \ Curve Fitting Toolbox \ Generating M-files From Curve Fitting Tool " ,以下简单举例说明: 以双色球从第125期到第145期蓝球为Y值: Y=[12 15 4 1 7 11 5 7 1 6 16 1 1 14 2 12 9 13 10 12 11]; X=1:1:21; cftool(X,Y); 点击Fitting选择最常用的多项式拟合(Polynomial),选择3次多项式拟合(cubic),然后就会出现如下拟合图形: 然后在Curve Fitting Tool窗口中点击 " \ File \ Generate M-file " 即可生成能直接曲线拟合的m函数文件,其中使用的拟合方法就是刚才使用的三次多项式拟合,文件中这条语句证明了这一点: ft_ = fittype('poly3'); 保存该m文件(默认叫做createFit.m),调用方法和通常的m文件一样,使用不同的X和Y值就能拟合出不同的曲线。但是,这种调用方法只能看到一个拟合出的图形窗口,拟合参数以及Goodness of fit参数都看不到了,因此需要在刚才的m文件中稍作修改。 找到这句话: cf_ = fit(X(ok_),Y(ok_),ft_); 修改为: [cf_,gof] = fit(X(ok_),Y(ok_),ft_); 然后将函数声明 function createFit(X,Y) 修改为 function [cf_,gof] = createFit(X,Y) ,这样我们再调用试试看: Y=[12 15 4 1 7 11 5 7 1 6 16 1 1 14 2 12 9 13 10 12 11]; X=1:1:21;

最小二乘法曲线拟合 原理及matlab实现

曲线拟合(curve-fitting ):工程实践中,用测量到的一些离散的数据},...2,1,0),,{(m i y x i i =求一个近似的函数)(x ?来拟合这组数据,要求所得的拟合曲线能最好的反映数据的基本趋势(即使)(x ?最好地逼近()x f ,而不必满足插值原则。因此没必要取)(i x ?=i y ,只要使i i i y x -=)(?δ尽可能地小)。 原理: 给定数据点},...2,1,0),,{(m i y x i i =。求近似曲线)(x ?。并且使得近似曲线与()x f 的偏差最小。近似曲线在该点处的偏差i i i y x -=)(?δ,i=1,2,...,m 。 常见的曲线拟合方法: 1.使偏差绝对值之和最小 2.使偏差绝对值最大的最小 3.使偏差平方和最小 最小二乘法: 按偏差平方和最小的原则选取拟合曲线,并且采取二项式方程为拟合曲线的方法,称为最小二乘法。 推导过程: 1. 设拟合多项式为: 2. 各点到这条曲线的距离之和,即偏差平方和如下: 3. 问题转化为求待定系数0a ...k a 对等式右边求i a 偏导数,因而我们得到 了: ....... 4、 把这些等式化简并表示成矩阵的形式,就可以得到下面的矩阵: 5. 将这个范德蒙得矩阵化简后可得到:

6. 也就是说X*A=Y,那么A = (X'*X)-1*X'*Y,便得到了系数矩阵A,同时,我们也就得到了拟合曲线。 MATLAB实现: MATLAB提供了polyfit()函数命令进行最小二乘曲线拟合。 调用格式:p=polyfit(x,y,n) [p,s]= polyfit(x,y,n) [p,s,mu]=polyfit(x,y,n) x,y为数据点,n为多项式阶数,返回p为幂次从高到低的多项式系数向量p。x必须是单调的。矩阵s包括R(对x进行QR分解的三角元素)、df(自由度)、normr(残差)用于生成预测值的误差估计。 [p,s,mu]=polyfit(x,y,n)在拟合过程中,首先对x进行数据标准化处理,以在拟合中消除量纲等影响,mu包含标准化处理过程中使用的x的均值和标准差。 polyval( )为多项式曲线求值函数,调用格式: y=polyval(p,x) [y,DELTA]=polyval(p,x,s) y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。 [y,DELTA]=polyval(p,x,s) 使用polyfit函数的选项输出s得出误差估计Y DELTA。它假设polyfit函数数据输入的误差是独立正态的,并且方差为常数。则Y DELTA将至少包含50%的预测值。 如下给定数据的拟合曲线: x=[0.5,1.0,1.5,2.0,2.5,3.0], y=[1.75,2.45,3.81,4.80,7.00,8.60]。 解:MATLAB程序如下: x=[0.5,1.0,1.5,2.0,2.5,3.0]; y=[1.75,2.45,3.81,4.80,7.00,8.60]; p=polyfit(x,y,2) x1=0.5:0.05:3.0; y1=polyval(p,x1); plot(x,y,'*r',x1,y1,'-b') 运行结果如图1 计算结果为: p =0.5614 0.8287 1.1560 即所得多项式为y=0.5614x^2+0.08287x+1.15560 图1 最小二乘法曲线拟合示例 对比检验拟合的有效性: 例:在[0,π]区间上对正弦函数进行拟合,然后在[0,2π]区间画出图形,比较拟合区间和非拟合区间的图形,考察拟合的有效性。 在MATLAB中输入如下代码: clear x=0:0.1:pi; y=sin(x); [p,mu]=polyfit(x,y,9)

最小二乘法曲线拟合原理及matlab实现

最小二乘法曲线拟合原理及m a t l a b实现 Modified by JEEP on December 26th, 2020.

曲线拟合(curve-fitting ):工程实践中,用测量到的一些离散的数据},...2,1,0),,{(m i y x i i =求一个近似的函数)(x ?来拟合这组数据,要求所得的拟合曲线能最好的反映数据的基本趋势(即使)(x ?最好地逼近()x f ,而不必满足插值原则。因此没必要取)(i x ?=i y ,只要使i i i y x -=)(?δ尽可能地小)。 原理: 给定数据点},...2,1,0),,{(m i y x i i =。求近似曲线)(x ?。并且使得近似曲线与()x f 的偏差最小。近似曲线在该点处的偏差i i i y x -=)(?δ,i=1,2,...,m 。 常见的曲线拟合方法: 1.使偏差绝对值之和最小 2.使偏差绝对值最大的最小 3.使偏差平方和最小 最小二乘法: 按偏差平方和最小的原则选取拟合曲线,并且采取二项式方程为拟合曲线的方法,称为最小二乘法。 推导过程: 1. 设拟合多项式为: 2. 各点到这条曲线的距离之和,即偏差平方和如下: 3. 问题转化为求待定系数0a ...k a 对等式右边求i a 偏导数,因而我们得到了: ....... 4、 把这些等式化简并表示成矩阵的形式,就可以得到下面的矩阵: 5. 将这个范德蒙得矩阵化简后可得到: 6. 也就是说X*A=Y ,那么A = (X'*X)-1*X'*Y ,便得到了系数矩阵A ,同时,我们也就得到了拟合曲线。

MATLAB实现: MATLAB提供了polyfit()函数命令进行最小二乘曲线拟合。 调用格式:p=polyfit(x,y,n) [p,s]= polyfit(x,y,n) [p,s,mu]=polyfit(x,y,n) x,y为数据点,n为多项式阶数,返回p为幂次从高到低的多项式系数向量p。x必须是单调的。矩阵s包括R(对x进行QR分解的三角元素)、df(自由度)、normr(残差)用于生成预测值的误差估计。 [p,s,mu]=polyfit(x,y,n)在拟合过程中,首先对x进行数据标准化处理,以在拟合中消除量纲等影响,mu包含标准化处理过程中使用的x的均值和标准差。 polyval( )为多项式曲线求值函数,调用格式: y=polyval(p,x) [y,DELTA]=polyval(p,x,s) y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。 [y,DELTA]=polyval(p,x,s) 使用polyfit函数的选项输出s得出误差估计Y DELTA。它假设polyfit函数数据输入的误差是独立正态的,并且方差为常数。则Y DELTA将至少包含50%的预测值。 如下给定数据的拟合曲线: x=[,,,,,], y=[,,,,,]。 解:MATLAB程序如下: x=[,,,,,]; y=[,,,,,]; p=polyfit(x,y,2) x1=::; y1=polyval(p,x1); plot(x,y,'*r',x1,y1,'-b') 运行结果如图1

php九九乘法表

"; for($j = 1;$j <= 9;$j++){ echo ""; for($z = 0;$z <9 - $j;$z++){ echo " "; } for($i = $j;$i >= 1;$i--){ echo "$i*$j=".$i*$j.""; } echo ""; } echo "";*/ //第三种 /*echo "

"; for($j = 9;$j >= 1;$j--){ echo ""; for($z = 1;$z <= 9 - $j;$z++ ){ echo ""; //echo "$i*$j=".$i*$j."";

} for($i = 1;$i <= $j;$i++){ echo "

"; } echo ""; } echo "
$i*$j=".$i*$j."
";*/ /*第二种 for($j = 9;$j >= 1;$j--){ for($i = 1; $i <= $j; $i++){ //echo "$i x $j = ".$i * $j; //echo " "; echo "$i x $j =".$i*$j."  "; } echo "
"; } */ /*第一种 for($j = 1;$j <= 9;$j++){

用matlab实现最小二乘递推算法辨识系统参数

用matlab实现最小二乘递推算法辨识系统参 数 自动化系统仿真实验室指导教师: 学生姓名班级计082-2 班学号撰写时间: 全文结束》》-3-1 成绩评定: 一.设计目的 1、学会用Matlab实现最小二乘法辨识系统参数。 2、进一步熟悉Matlab的界面及基本操作; 3、了解并掌握Matlab中一些函数的作用与使用;二.设计要求最小二乘递推算法辨识系统参数,利用matlab编程实现,设初始参数为零。z(k)-1、5*z(k-1)+0、7*z(k-2)=1*u(k-1)+0、5*u(k-2)+v(k); 选择如下形式的辨识模型:z(k)+a1*z(k- 1)+a2*z(k-2)=b1*u(k-1)+b2*u(k-2)+v(k);三.实验程序 m=3;N=100;uk=rand(1,N);for i=1:Nuk(i)=uk(i)*(-1)^(i-1);endyk=zeros(1,N); for k=3:N yk(k)=1、5*yk(k-1)-0、 7*yk(k-2)+uk(k-1)+0、5*uk(k-2); end%j=100;kn=0;%y=yk(m:j);%psi=[yk(m-1:j-1);yk(m-2:j-2);uk(m-1:j-1);uk(m-2:j- 2)];%pn=inv(psi*psi);%theta=(inv(psi*psi)*psi*y);theta=[0 ;0;0;0];pn=10^6*eye(4);for t=3:Nps=([yk(t-1);yk(t-

2);uk(t-1);uk(t-2)]);pn=pn- pn*ps*ps*pn*(inv(1+ps*pn*ps));theta=theta+pn*ps*(yk(t)-ps*theta);thet=theta;a1=thet(1);a2=thet(2);b1=thet(3);b2= thet(4); a1t(t)=a1;a2t(t)=a2;b1t(t)=b1;b2t(t)=b2;endt=1:N;plot(t,a 1t(t),t,a2t(t),t,b1t(t),t,b2t(t));text(20,1、 47,a1);text(20,-0、67,a2);text(20,0、97,b1);text(20,0、47,b2);四.设计实验结果及分析实验结果图:仿真结果表明,大约递推到第步时,参数辨识的结果基本到稳态状态,即a1=1、5999,b1=1,c1=0、5,d1=-0、7。五、设计感受这周的课程设计告一段落了,时间短暂,意义重大。通过这次次练习的机会,重新把matlab课本看了一遍,另外学习了系统辨识的有关内容,收获颇丰。对matlab的使用更加纯熟,也锻炼了自己在课本中搜索信息和知识的能力。在设计过程中虽然遇到了一些问题,但经过一次又一次的思考,一遍又一遍的检查终于找出了原因所在,也暴露出了前期我在这方面的知识欠缺和经验不足。同时我也进一步认识了matlab软件强大的功能。在以后的学习和工作中必定有很大的用处。

Matlab最小二乘法曲线拟合的应用实例

MATLAB机械工程 最小二乘法曲线拟合的应用实例 班级: 姓名: 学号: 指导教师:

一,实验目的 通过Matlab上机编程,掌握利用Matlab软件进行数据拟合分析及数据可视化方法 二,实验内容 1.有一组风机叶片的耐磨实验数据,如下表所示,其中X为使用时间,单位为小时h,Y为磨失质量,单位为克g。要求: 对该数据进行合理的最小二乘法数据拟合得下列数据。 x=[10000 11000 12000 13000 14000 15000 16000 17000 18000 19000 2 0000 21000 22000 23000]; y=[24.0 26.5 29.8 32.4 34.7 37.7 41.1 42.8 44.6 47.3 65.8 87.5 137.8 174. 2] 三,程序如下 X=10000:1000:23000; Y=[24.0,26.5,29.8,32.4,34.7,37.7,41.1,42.8,44.6,47.3,65.8,87.5,137.8,17 4.2] dy=1.5; %拟合数据y的步长for n=1:6 [a,S]=polyfit(x,y,n); A{n}=a;

da=dy*sqrt(diag(inv(S.R′*S.R))); Da{n}=da′; freedom(n)=S.df; [ye,delta]=polyval(a,x,S); YE{n}=ye; D{n}=delta; chi2(n)=sum((y-ye).^2)/dy/dy; end Q=1-chi2cdf(chi2,freedom); %判断拟合良好度 clf,shg subplot(1,2,1),plot(1:6,abs(chi2-freedom),‘b’) xlabel(‘阶次’),title(‘chi2与自由度’) subplot(1,2,2),plot(1:6,Q,‘r’,1:6,ones(1,6)*0.5) xlabel(‘阶次’),title(‘Q与0.5线’) nod=input(‘根据图形选择适当的阶次(请输入数值)’); elf,shg, plot(x,y,‘kx’);xlabel(‘x’),ylabel(‘y’); axis([8000,23000,20.0,174.2]);hold on errorbar(x,YE{nod},D{nod},‘r’);hold off title(‘较适当阶次的拟合’) text(10000,150.0,[‘chi2=’num2str(chi2(nod))‘~’int2str(freedom(nod))])

MATLAB实现非线性曲线拟合最小二乘法

非线性曲线拟合最小二乘法 一、问题提出 设数据(i i y x ,),(i=0,1,2,3,4).由表3-1给出,表中第四行为i i y y =ln ,可以看出数学模型为bx ae y =,用最小二乘法确定a 及b 。 i 0 1 2 3 4 i x 1.00 1.25 1.50 1.75 2.00 i y 5.10 5.79 6.53 7.45 8.46 i y 1.629 1.756 1.876 2.008 2.135 二、理论基础 根据最小二乘拟合的定义:在函数的最佳平方逼近中],[)(b a C x f ∈,如果f(x)只在一组离散点集{i x ,i=0,1,…,m},上给定,这就是科学实验中经常见到的实验数据{(i i y x ,), i=0,1,…,m}的曲线拟合,这里)(i i x f y =,i=0,1,…,m,要求一个函数)(*x S y =与所给数据{(i i y x ,),i=0,1,…,m}拟合,若记误差 i i i y x S -=)(*δ,i=0,1,…,m,T m ),,(10δδδδ, =,设)(,),(),(10x x x n ??? 是] ,[b a C 上线性无关函数族,在)}(,),(),({10x x x span n ???? =中找一函数)(*x S ,使误差平方和 ∑∑∑===∈ -=-==m i m i m i i i x S i i i y x S y x S 0 2 )(2 * 2 22 ])([])([min ? δδ , 这里 )()()()(1100x a x a x a x S n n ???+++= (n

matlab循环语句

matlab 基本语句 1.循环语句for for i=s1:s3:s2 循环语句组 end 解释:首先给i赋值s1;然后,判断i是否介于s1与s2之间;如果是,则执行循环语句组,i=i+s3(否则,退出循环.);执行完毕后,继续下一次循环。 例:求1到100的和,可以编程如下: sum=0 for i=1:1:100 sum=sum+i end 这个程序也可以用while语句编程。 注:for循环可以通过break语句结束整个for循环. 2.循环语句while 例:sum=0;i=1; while(i<=100) sum=sum+i;i=i+1; end 3.if语句 if(条件) 语句 end if(条件) 语句 else 语句 end if(条件) 语句 elseif 语句 end 4.关系表达式:

=,>,<,>=,<=,==(精确等于) 5.逻辑表达式:|(或),&(且) 6.[n,m]=size(A)(A为矩阵) 这样可以得到矩阵A的行和列数 n=length(A),可以得到向量A的分量个数;如果是矩阵,则得到矩阵A的行与列数这两个数字中的最大值。 7.!后面接Dos命令可以调用运行一个dos程序。 8.常见函数: poly():为求矩阵的特征多项式的函数,得到的为特征多项式的各个系数。如 a=[1,0,0;0,2,0;0,0,3],则poly(a)=1 -6 11 -6。相当于poly(a)=1入^3+(-6)入^2+11入+(-6)。 compan():可以求矩阵的伴随矩阵. sin()等三角函数。 MATLAB在数学建模中的应用(3) 一、程序设计概述 MATLAB所提供的程序设计语言是一种被称为第四代编程语言的高级程序设计语言,其程序简洁,可读性很强,容易调试。同时,MATLAB的编程效率比C/C++语言要高得多。 MATLAB编程环境有很多。常用的有: 1.命令窗口 2.word窗口 3.M-文件编辑器,这是最好的编程环境。 M-文件的扩展名为“.m”。M-文件的格式分为两种: ①λ M-脚本文件,也可称为“命令文件”。 ② M-函数文件。这是matlab程序设计的主流。λ 保存后的文件可以随时调用。 二、MATLAB程序结构 按照现代程序设计的观点,任何算法功能都可以通过三种基本程序结构来实现,这三种结构是:顺序结构、选择结构和循环结构。其中顺序结构是最基本的结构,它依照语句的自然顺序逐条地执行程序的各条语句。如果要根据输入数据的实际情况进行逻辑判断,对不同的结果进行不同的处理,可以使用选择结构。如果需要反复执行某些程序段落,可以使用循环结构。 1 顺序结构 顺序结构是由两个程序模块串接构成。一个程序模块是完成一项独立功能的逻辑单元,它可以是一段程序、一个函数,或者是一条语句。 看图可知,在顺序结构中,这两个程序模块是顺序执行的,即先执行<程序

MATLAB程序(线性拟合)

1、一元线性拟合 求HNO 3的正常沸点温度T b 及摩尔汽化热。 程序如下: >> t=[0 20 40 50 70 80 90 100]; >> t=t+273.15; >> p=[1919.52 6385.07 17728.9 27726.4 62251.1 89311 124902.1 170890.6] p = 1.0e+005 * 0.0192 0.0639 0.1773 0.2773 0.6225 0.8931 1.2490 1.7089 >> subplot 121 >> plot(t,p,'o',t,p) >> t1=1./t;p2=log(p); >> pp=polyfit(t1,p2,1) pp = 1.0e+003 * -4.5691 0.0243 >> subplot 122 >> plot(t1,p2,'o',t1,p2) >> gtext('p/pa'),gtext('T/K'),GTEXT('lnP/Pa'),gtext('T^-^1/K') 由克拉贝龙-克劳修斯方程式,~ ln v H P C RT ?=-+ 作1 ln ~P T -得一直线:3 1 ln 4.5691024.30P T -=-?+ 斜率为:~ 3 4.56910v H R ?-?=-

所以摩尔汽化热为:~ 314.569108.31437.99()v H kJ mol -?=??=? 并根据拟合方程,求得一大气压时 1 32.8010T --=? 则正常沸点为:357b T K = 2、多元线性拟合: 某气体混合物由四种气体组成,在常压或低压下其粘度η与各组分摩尔分数x 1,x 2,x 3,x 4之间有如下线性关系:011223344b b x b x b x b x η=++++ 试根据下表所列实验数据用最小二乘法确定上式中的各个系数,并计算其复相关系数。 Matlab 程序如下: >> a=[1.0 0.402 0.153 0.058 0.387;1.0 0.503 0.301 0.183 0.013; 1.0 0.306 0.109 0.224 0.361; 1.0 0.296 0.365 0.009 0.330; 1.0 0.309 0.405 0.109 0.177; 1.0 0.055 0.153 0.506 0.289] a = 1.0000 0.4020 0.1530 0.0580 0.3870 1.0000 0.5030 0.3010 0.1830 0.0130 1.0000 0.3060 0.1090 0.2240 0.3610 1.0000 0.2960 0.3650 0.0090 0.3300 1.0000 0.3090 0.4050 0.1090 0.1770 1.0000 0.0550 0.1530 0.5060 0.2890 >> y=[0.00625 0.00826 0.01182 0.01944 0.02372 0.03243]' y = 0.0063 0.0083 0.0118 0.0194 0.0237 0.0324 >> b=a.'*a

最小二乘法的多项式拟合matlab实现

最小二乘法的多项式拟 合m a t l a b实现 Document serial number【NL89WT-NY98YT-NC8CB-NNUUT-NUT108】

用最小二乘法进行多项式拟合(matlab 实现) 西安交通大学 徐彬华 算法分析: 对给定数据 (i=0 ,1,2,3,..,m),一共m+1个数据点,取多项式P(x),使 函数P(x)称为拟合函数或最小二乘解,令似的 使得 其中,a0,a1,a2,…,an 为待求未知数,n 为多项式的最高次幂,由此,该问题化为求 的极值问题。由多元函数求极值的必要条件: j=0,1,…,n 得到: j=0,1,…,n 这是一个关于a0,a1,a2,…,an 的线性方程组,用矩阵表示如下:

因此,只要给出数据点 及其个数m ,再给出所要拟合的参数n ,则即可求出未知数矩阵(a0,a1,a2,…,an ) 试验题1 编制以函数 为基的多项式最小二乘拟合程序,并用于对下列数据作三次多项式最小二乘拟合(取权函数wi ≡1) x i y i 总共有7个数据点,令m=6 第一步:画出已知数据的的散点图,确定拟合参数n; x=::;y=[,,,,,,]; plot(x,y,'*') xlabel 'x 轴' ylabel 'y 轴' title '散点图' hold on {} n k k x 0=

因此将拟合参数n设为3. 第二步:计算矩阵 A= 注意到该矩阵为(n+1)*(n+1)矩阵, 多项式的幂跟行、列坐标(i,j)的关系为i+j-2,由此可建立循环来求矩阵的各个元素,程序如下: m=6;n=3; A=zeros(n+1); for j=1:n+1 for i=1:n+1 for k=1:m+1 A(j,i)=A(j,i)+x(k)^(j+i-2) end end

九九乘法表的C语言代码

九九乘法表的C语言代码,黄路平编写与2012.3.6 代码一:#include int main() { int i=1,j; for (i=1,j=1;j<=9;j++) { if( j==1) printf("%d*%d=%d\n",i,j,i*j); else {for (i=1;i<=j;i++) printf("%d*%d=%d\t",i,j,i*j); printf("\n"); } } } 代码二:switch语句 #include int main() { int i=1,j; for (i=1,j=1;j<=9;j++) { switch(j) { case 1:printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 2: for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 3:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 4:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 5:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break;

case 6:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 7:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 8:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 9:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; } } }

MATLAB软件基本的曲线拟合函数命令

MATLAB软件提供了基本的曲线拟合函数的命令。 曲线拟合就是计算出两组数据之间的一种函数关系,由此可描绘其变化曲线及估计非采集数据对应的变量信息。 1.线性拟合函数:regress() 调用格式: b = regress(y,X) [b,bint,r,rint,stats] = regress(y,X) [b,bint,r,rint,stats] = regress(y,X,alpha) 说明:b=[ε; β],regress(y,X)返回X与y的最小二乘拟合的参数值β、ε,y=ε+βX。β是p′1的参数向量;ε是服从标准正态分布的随机干扰的n′1的向量;y为n′1的向量;X为n′p矩阵。 bint返回β的95%的置信区间。 r中为形状残差,rint中返回每一个残差的95%置信区间。Stats向量包含R2统计量、回归的F值和p值。 例: x=[ones(10,1) (1:10)']; y=x*[10;1]+normrnd(0,0.1,10,1); [b,bint]=regress(y,x,0.05) 结果得回归方程为:y=9.9213+1.0143x 2.多项式曲线拟合函数:polyfit() 调用格式: p = polyfit(x,y,n) [p,s] = polyfit(x,y,n) 说明:n:多项式的最高阶数; x,y:将要拟合的数据,用数组的方式输入; p:为输出参数,即拟合多项式的系数; 多项式在x处的值y可用下面程序计算: y=polyval(p,x) 例: x=1:20; y=x+3*sin(x); p=polyfit(x,y,6) xi=linspace(1,20,100); z=polyval(p,xi); % 多项式求值函数

matlab实现插值法和曲线拟合电子教案

m a t l a b实现插值法和 曲线拟合

插值法和曲线拟合 电子科技大学 摘要:理解拉格朗日多项式插值、分段线性插值、牛顿前插,曲线拟合,用matlab编程求解函数,用插值法和分段线性插值求解同一函数,比较插值余项;用牛顿前插公式计算函数,计算函数值;对于曲线拟 合,用不同曲线拟合数据。 关键字:拉格朗日插值多项式;分段线性插值;牛顿前插;曲线拟合 引言: 在数学物理方程中,当给定数据是不同散点时,无法确定函数表达式,求解函数就需要很大的计算量,我们有多种方法对给定的表格函数进行求解,我们这里,利用插值法和曲线拟合对函数进行求解,进一步了解函数性质,两种方法各有利弊,适合我们进行不同的散点函数求解。 正文: 一、插值法和分段线性插值 1拉格朗日多项式原理 对某个多项式函数,已知有给定的k + 1个取值点: 其中对应着自变量的位置,而对应着函数在这个位置的取值。 假设任意两个不同的x j都互不相同,那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为: 其中每个为拉格朗日基本多项式(或称插值基函数),其表达式为: [3] 拉格朗日基本多项式的特点是在上取值为1,在其它的点 上取值为0。 2分段线性插值原理 给定区间[a,b], 将其分割成a=x 0

曲线拟合的最小二乘法matlab举例

曲线拟合的最小二乘法 学院:光电信息学院 姓名:赵海峰 学号: 200820501001 一、曲线拟合的最小二乘法原理: 由已知的离散数据点选择与实验点误差最小的曲线 S( x) a 0 0 ( x) a 1 1(x) ... a n n ( x) 称为曲线拟合的最小二乘法。 若记 m ( j , k ) i (x i ) j (x i ) k (x i ), 0 m (f , k ) i0 (x i )f (x i ) k (x i ) d k n 上式可改写为 ( k , jo j )a j d k ; (k 0,1,..., n) 这个方程成为法方程,可写成距阵 形式 Ga d 其中 a (a 0,a 1,...,a n )T ,d (d 0,d 1,...,d n )T , 、 数值实例: 下面给定的是乌鲁木齐最近 1个月早晨 7:00左右(新疆时间 )的天气预报所得 到的温度数据表,按照数据找出任意次曲线拟合方程和它的图像。 它的平方误差为: || 2 | 2 ] x ( f

(2008 年 10 月 26~11 月 26) F 面应用Matlab 编程对上述数据进行最小二乘拟合 三、Matlab 程序代码: x=[1:1:30]; y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1]; %三次多项式拟合% %九次多项式拟合% %十五次多项式拟合% %三次多项式误差平方和 % %九次次多项式误差平方和 % %十五次多项式误差平方和 % %用*画出x,y 图像% %用红色线画出x,b1图像% %用绿色线画出x,b2图像% %用蓝色o 线画出x,b3图像% 四、数值结果: 不同次数多项式拟和误差平方和为: r1 = 67.6659 r2 = 20.1060 r3 = 3.7952 r1、r2、r3分别表示三次、九次、十五次多项式误差平方和 拟和曲线如下图: a 仁polyfit(x,y,3) a2= polyfit(x,y,9) a3= polyfit(x,y,15) b1= polyval(a1,x) b2= polyval(a2,x) b3= polyval(a3,x) r1= sum((y-b1).A 2) r2= sum((y-b2).A2) r3= sum((y-b3).A2) plot(x,y,'*') hold on plot(x,b1, 'r') hold on plot(x,b2, 'g') hold on plot(x,b3, 'b:o')

九九乘法表源代码(vb)

Private Sub Command1_Click() For i = 1 To 4 For j = 1 To 6 Print "*"; Next j Print Next i End Sub Private Sub Command2_Click() For i = 1 To 4 Print Tab(i); For j = 1 To 6 Print "*"; Next j Print Next i End Sub Private Sub Command3_Click() For i = 1 To 4 Print Tab(5 - i); For j = 1 To 6 Print "*"; Next j Print Next i End Sub Private Sub Command4_Click() For i = 1 To 9 For j = 1 To 9 Print i; "*"; j; "="; i * j; Next j Print Next i End Sub Private Sub Command5_Click() Dim se As String Print Tab(35); "乘法表" For i = 1 To 9 For j = 1 To 9 se = i & "*" & j & "=" & i * j

Print Tab((j - 1) * 9); se; Next j Picture1.Print Next i End Sub Private Sub Command6_Click() End End Sub Private Sub Command7_Click() Print Tab(35); "乘法表" For i = 1 To 9 For j = i To 9 se = i & "*" & j & "=" & i * j Print Tab((j - 1) * 9); se; Next j Print Next i End Sub Private Sub Command8_Click() Print Tab(35); "乘法表" For i = 1 To 9 For j = 1 To i se = i & "*" & j & "=" & i * j Print Tab((j - 1) * 9); se; Next j Print Next i End Sub Private Sub Command9_Click() Cls End Sub Private Sub Picture1_Click() Dim se As String Picture1.Print Tab(35); "乘法表" For i = 1 To 9 For j = 1 To 9 se = i & "*" & j & "=" & i * j Picture1.Print Tab((j - 1) * 9); se; Next j

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

3.1 曲线拟合的线性最小二乘法及其MATLAB 程序 例3.1.1 给出一组数据点),(i i y x 列入表3-1中,试用线性最小二乘法求拟合曲线,并估计其误差,作出拟合曲线. 表3-1 例3.1.1的一组数据),(y x 解 (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.10 -8.43 -13.12 6.50 68.04]; plot(x,y,'r*'), legend('实验数据(xi,yi)') xlabel('x'), ylabel('y'), title('例3.1.1的数据点(xi,yi)的散点图') 运行后屏幕显示数据的散点图(略). (3)编写下列MA TLAB 程序计算)(x f 在),(i i y x 处的函数值,即输入程序 >> syms a1 a2 a3 a4 x=[-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 运行后屏幕显示关于a 1,a 2, a 3和a 4的线性方程组 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] 编写构造误差平方和的MATLAB 程序 >> y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.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; 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*a 2+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 为求4321,,,a a a a 使J 达到最小,只需利用极值的必要条件0=??k a J )4,3,2,1(=k ,

相关文档
最新文档